| Title: | Kolmogorov-Zurbenko Adaptive Filters |
|---|---|
| Description: | Time Series Analysis including break detection, spectral analysis, KZ Fourier Transforms. |
| Authors: | Brian Close <[email protected]>, Igor Zurbenko <[email protected]> and Mingzeng Sun <[email protected]> |
| Maintainer: | Brian Close <[email protected]> |
| License: | GPL-3 |
| Version: | 4.1.0.1 |
| Built: | 2025-03-06 03:17:09 UTC |
| Source: | https://github.com/cran/kza |
Kolmogorov-Zurbenko low-pass linear filter.
kz(x, m, k = 3)kz(x, m, k = 3)
x |
The raw data that will be smoothed. The data can have as many as 3 dimensions. KZ will also handle a time series. |
m |
Window size for the filter. This can be up to 3 dimensions, but not more than the dimensionality of the input data x. |
k |
Number of iterations. |
KZ is an iterated moving average. The filter can be used with missing values. One iteration is equivalent to a simple moving average. Three iterations is an approximately Gaussian shaped filter.
Zurbenko, I. G., 1986: The spectral Analysis of Time Series. North-Holland, 248 pp.
## 2 dimensions set.seed(2) a <- matrix(rep(0,100*100),nrow=100) a[35:70,35:70]<-1 a <- a + matrix(rnorm(100*100,0,1),nrow=100) z<-kz(a,m=c(20,5),k=3) x <- seq(1,100) y <- x op <- par(bg = "white") c="lightblue" m="Unsmoothed" persp(x, y, a, zlab="a", ticktype="detailed", theta = 60, phi = 45, col = c, main=m) m="KZ(a,m=c(20,5),k=3)" persp(x, y, z, zlab="z", ticktype="detailed", theta = 60, phi = 45, col = c, main=m) #example t <- seq(0,20,length=20*365) set.seed(6); e <- rnorm(n = length(t), sd = 2.0) y <- sin(3*pi*t) + e z <- kz(y,30) par(mfrow=c(2,1)) plot(y,ylim=c(-5,5),type="l",main="y = sin(3*pi*t) + noise") plot(z,ylim=c(-5,5), type="l",main="KZ filter") lines(sin(3*pi*t), col="blue") par(mfrow=c(1,1))## 2 dimensions set.seed(2) a <- matrix(rep(0,100*100),nrow=100) a[35:70,35:70]<-1 a <- a + matrix(rnorm(100*100,0,1),nrow=100) z<-kz(a,m=c(20,5),k=3) x <- seq(1,100) y <- x op <- par(bg = "white") c="lightblue" m="Unsmoothed" persp(x, y, a, zlab="a", ticktype="detailed", theta = 60, phi = 45, col = c, main=m) m="KZ(a,m=c(20,5),k=3)" persp(x, y, z, zlab="z", ticktype="detailed", theta = 60, phi = 45, col = c, main=m) #example t <- seq(0,20,length=20*365) set.seed(6); e <- rnorm(n = length(t), sd = 2.0) y <- sin(3*pi*t) + e z <- kz(y,30) par(mfrow=c(2,1)) plot(y,ylim=c(-5,5),type="l",main="y = sin(3*pi*t) + noise") plot(z,ylim=c(-5,5), type="l",main="KZ filter") lines(sin(3*pi*t), col="blue") par(mfrow=c(1,1))
KZA will recover 2-dimensional or 3-dimensional image or signal buried in noise.
kza(x, m, y = NULL, k = 3, min_size = round(0.05*m), tol = 1.0e-5, impute_tails = FALSE) ## S3 method for class 'kza' plot(x, ...)kza(x, m, y = NULL, k = 3, min_size = round(0.05*m), tol = 1.0e-5, impute_tails = FALSE) ## S3 method for class 'kza' plot(x, ...)
x |
A vector of the time series or a matrix (2d) or an array (3d) of an image. |
m |
The window for the filter. |
y |
The filtered output from kz. |
k |
The number of iterations. |
min_size |
Minimum size of window q. |
tol |
The smallest value to accept as nonzero. |
impute_tails |
The default is to drop the tails. |
... |
Other parameters. |
The selection of parameters of KZA depend on the nature of the data. This function may take a long time to run, depending on the number of dimensions and the size of the dimensions.
Brian Close <[email protected]> and Igor Zurbenko <[email protected]>
I. Zurbenko, P.S. Porter, S.T. Rao, J.Y. Ku, R. Gui, R.E. Eskridge Detecting Discontinuities in Time Series of Upper-air Data: Development and Demonstration of an Adaptive Filter Technique. Journal of Climate: (1996) Vol. 9, No. 12, pp. 3548 3560. http://journals.ametsoc.org/action/doSearch?AllField=zurbenko&filter=AllField
Kevin L. Civerolo, Elvira Brankov, S. T. Rao, Igor Zurbenko Assessing the impact of the acid deposition control program. Atmospheric Environment 35 (2001) 4135-4148 http://www.elsevier.com/locate/atmosenv
J.Chen, I.Zurbenko, Nonparametric Boundary detection, Communications in Statistics, Theory and Methods, Vol.26, 12, 2999-3014, 1997.
####### # this is an example of detection of a break point in a time series ####### yrs <- 20 t <- seq(0,yrs,length=yrs*365) m <- 365 #noise e <- rnorm(n = length(t),0,1) trend <- seq(0,-1,length=length(t)) #signal bkpt <- 3452 brk <- c(rep(0,bkpt),rep(0.5,length(t)-bkpt)) signal <- trend + brk # y = seasonal + trend + break point + noise y <- sin(2*pi*t) + signal + e k.kz <- kz(y,m) # kza reconstruction of the signal k.kza <- kza(y,m,y=k.kz,min_size=10) par(mfrow=c(2,1)) plot(y,type="l", ylim=c(-3,3)) plot(signal,type="l",ylim=c(-3,3), main="Signal and KZA Reconstruction") lines(k.kza$kza, col=4) ###################### # image detection (2d) ###################### set.seed(2) a <- matrix(rep(0,100*100),nrow=100) a[35:70,35:70]<-1 a <- a + matrix(rnorm(100*100,0,1),nrow=100) y<-kz(a,m=15,k=3) v <- kza(a,m=15,y=y,k=3,impute_tails=TRUE) x <- seq(1,100) y <- x op <- par(bg = "white") ### #noise ### c="lightblue" persp(x, y, a, zlab="z", zlim=c(-5,5), ticktype="detailed", theta=30, phi=30, col=c) ### #kza filtered ### persp(x,y,v$kza,zlab="z",zlim=c(-5,5),ticktype="detailed",theta=30,phi=30,col=c) ### # another view ### par(mfrow=c(1,2)) image(a,col=gray(seq(0,1,1/255))) image(v$kza,col=gray(seq(0,1,1/255))) par(mfrow=c(1,1))####### # this is an example of detection of a break point in a time series ####### yrs <- 20 t <- seq(0,yrs,length=yrs*365) m <- 365 #noise e <- rnorm(n = length(t),0,1) trend <- seq(0,-1,length=length(t)) #signal bkpt <- 3452 brk <- c(rep(0,bkpt),rep(0.5,length(t)-bkpt)) signal <- trend + brk # y = seasonal + trend + break point + noise y <- sin(2*pi*t) + signal + e k.kz <- kz(y,m) # kza reconstruction of the signal k.kza <- kza(y,m,y=k.kz,min_size=10) par(mfrow=c(2,1)) plot(y,type="l", ylim=c(-3,3)) plot(signal,type="l",ylim=c(-3,3), main="Signal and KZA Reconstruction") lines(k.kza$kza, col=4) ###################### # image detection (2d) ###################### set.seed(2) a <- matrix(rep(0,100*100),nrow=100) a[35:70,35:70]<-1 a <- a + matrix(rnorm(100*100,0,1),nrow=100) y<-kz(a,m=15,k=3) v <- kza(a,m=15,y=y,k=3,impute_tails=TRUE) x <- seq(1,100) y <- x op <- par(bg = "white") ### #noise ### c="lightblue" persp(x, y, a, zlab="z", zlim=c(-5,5), ticktype="detailed", theta=30, phi=30, col=c) ### #kza filtered ### persp(x,y,v$kza,zlab="z",zlim=c(-5,5),ticktype="detailed",theta=30,phi=30,col=c) ### # another view ### par(mfrow=c(1,2)) image(a,col=gray(seq(0,1,1/255))) image(v$kza,col=gray(seq(0,1,1/255))) par(mfrow=c(1,1))
Kolmogorov-Zurbenko Fourier Transform is an iterated Fourier transform.
kzft(x, f=0, m=1, k=1) coeff(m, k) transfer_function(m, k, lamda = seq(-0.5,0.5,by=0.01), omega = 0 )kzft(x, f=0, m=1, k=1) coeff(m, k) transfer_function(m, k, lamda = seq(-0.5,0.5,by=0.01), omega = 0 )
x |
The raw data |
f |
The frequency that KZFT is applied at. |
m |
The window size for transform |
k |
The number of iterations for applying the KZFT |
lamda |
The frequencies used for the calculating the transfer function. |
omega |
The frequency that KZFT is applied at. |
Kolmogorov-Zurbenko Fourier Transform (KZFT) is the Fourier transform applied over every segment of length m iterated k times. The argument alg="F" will use Fast Fourier Transforms written in C (fftw library). The alg="C" is a slow Fourier Transform but has the advantage of being able to handle missing values. It currently works in one dimension. The alg="R" is an R version of KZFT for experimental purposes.The coeff function generates the coefficients for the KZFT function.
You will introduce a phase shift and decrease the fidelity of the signal if the product of f*m is not an integer.
I. G. Zurbenko, The spectral Analysis of Time Series. North-Holland, 1986.
I. G. Zurbenko, P. S. Porter, Construction of high-resolution wavelets, Signal Processing 65: 315-327, 1998.
R. Neagu, I. G. Zurbenko, Tracking and separating non-stationary multi-component chirp signals, J. Franklin Inst., 499-520, 2002.
R. H. Shumway, D. S. Stoffer, Time Series Analysis and Its Applications: With R Examples, Springer, 2006.
Wei Yang and Igor Zurbenko, kzft: Kolmogorov-Zurbenko Fourier Transform and Applications, R-Project 2007.
Igor G. Zurbenko, Amy L. Potrzeba, Tidal Waves in Atmosphere and Their Effects, Acta Geophysica Volume 58, Number 2, 356-373
# example taken from Wei Yang's KZFT package # coefficients of kzft(201,5) # function to calculate polnomial coefficients for kzft ## Not run: a<-coeff(201,5); t<-seq(1:1001)-501; z<-cos(2*pi*0.025*t); plot(z*a,type="l",xlab="Time", ylab="Coefficient", main="Coefficients of the kzft"); lines(a); lines(-1*a); ## End(Not run) # example taken from Wei Yang's KZFT package # transfer function of the kzft(201,5) at frequency 0.025 lamda<-seq(-0.1,0.1,by=0.001) tf1<-transfer_function(201,1,lamda,0.025) tf2<-transfer_function(201,5,lamda,0.025) matplot(lamda,cbind(log(tf1),log(tf2)),type="l",ylim=c(-15,0), ylab="Natural log transformation of the coefficients", xlab="Frequency (cycles/time unit)", main="Transfer function of kzft(201,5) at frequency 0.025") # example with missing values set.seed(2) period=101 f<-1/period t<-1:2000 s<-1*sin(2*pi*f*t) x<-s noise<-3*rnorm(length(t)) x<-s+noise m=101 rand_idx <- sample(t,100,replace=FALSE) x[rand_idx]<-NA x<-as.vector(na.omit(x)) system.time(z1<-kzft(x, m=m, k=1, f=f)) system.time(z2<-kzft(x, m=m, k=2, f=f)) system.time(z3<-kzft(x, m=m, k=3, f=f)) par(mfrow=c(2,2)) plot(x,type="l",main="Original time series",xlab="t", ylab="y") lines(s,col="blue") plot(2*Re(z1),type="l",main="kzft(101,1)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") plot(2*Re(z2),type="l",main="kzft(101,2)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") plot(2*Re(z3),type="l",main="kzft(101,3)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") par(mfrow=c(1,1))# example taken from Wei Yang's KZFT package # coefficients of kzft(201,5) # function to calculate polnomial coefficients for kzft ## Not run: a<-coeff(201,5); t<-seq(1:1001)-501; z<-cos(2*pi*0.025*t); plot(z*a,type="l",xlab="Time", ylab="Coefficient", main="Coefficients of the kzft"); lines(a); lines(-1*a); ## End(Not run) # example taken from Wei Yang's KZFT package # transfer function of the kzft(201,5) at frequency 0.025 lamda<-seq(-0.1,0.1,by=0.001) tf1<-transfer_function(201,1,lamda,0.025) tf2<-transfer_function(201,5,lamda,0.025) matplot(lamda,cbind(log(tf1),log(tf2)),type="l",ylim=c(-15,0), ylab="Natural log transformation of the coefficients", xlab="Frequency (cycles/time unit)", main="Transfer function of kzft(201,5) at frequency 0.025") # example with missing values set.seed(2) period=101 f<-1/period t<-1:2000 s<-1*sin(2*pi*f*t) x<-s noise<-3*rnorm(length(t)) x<-s+noise m=101 rand_idx <- sample(t,100,replace=FALSE) x[rand_idx]<-NA x<-as.vector(na.omit(x)) system.time(z1<-kzft(x, m=m, k=1, f=f)) system.time(z2<-kzft(x, m=m, k=2, f=f)) system.time(z3<-kzft(x, m=m, k=3, f=f)) par(mfrow=c(2,2)) plot(x,type="l",main="Original time series",xlab="t", ylab="y") lines(s,col="blue") plot(2*Re(z1),type="l",main="kzft(101,1)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") plot(2*Re(z2),type="l",main="kzft(101,2)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") plot(2*Re(z3),type="l",main="kzft(101,3)",xlab="t", ylab="y", ylim=c(-6,6)) lines(s,col="blue") par(mfrow=c(1,1))
Kolmogorov-Zurbenko periodogram and smoothing using DiRienzo-Zurbenko (DZ).
kzp(y, m=length(y), k=1) ## S3 method for class 'kzp' smooth(object, log=TRUE, smooth_level=0.05, method = "DZ") ## S3 method for class 'kzp' nonlinearity(x) ## S3 method for class 'kzp' variation(x) ## S3 method for class 'kzp' summary(object, digits=getOption("digits"), top=1, ...) ## S3 method for class 'kzp' plot(x, ...)kzp(y, m=length(y), k=1) ## S3 method for class 'kzp' smooth(object, log=TRUE, smooth_level=0.05, method = "DZ") ## S3 method for class 'kzp' nonlinearity(x) ## S3 method for class 'kzp' variation(x) ## S3 method for class 'kzp' summary(object, digits=getOption("digits"), top=1, ...) ## S3 method for class 'kzp' plot(x, ...)
y |
The raw data. |
m |
The width of filtering window |
k |
The number of iterations for the KZFT |
object |
Output from kzp function. |
log |
Use logarithm values for smoothing. |
smooth_level |
Percentage of smoothness to apply. |
method |
Method used for smoothing; choices are "DZ" or "NZ". |
digits |
precision of output. |
top |
list top values |
... |
Other parameters. |
x |
periodogram |
The Kolmogorov-Zurbenko Periodogram is an estimate of the spectral density using the Kolmogorov-Zurbenko Fourier Transform (KZFT).
I. G. Zurbenko, 1986: The spectral Analysis of Time Series. North-Holland, 248 pp.
I. G. Zurbenko, P. S. Porter, Construction of high-resolution wavelets, Signal Processing 65: 315-327, 1998.
A. G. DiRienzo, I. G. Zurbenko, Semi-adaptive nonparametric spectral estimation, Journal of Computational and Graphical Statistics 8(1): 41-59, 1998.
R. Neagu, I. G. Zurbenko, Algorithm for adaptively smoothing the log-periodgram, Journal of the Franklin Institute 340: 103-123, 2003.
Wei Yang and Igor Zurbenko, kzft: Kolmogorov-Zurbenko Fourier Transform and Applications, R-Project 2007.
## Not run: t<-1:6000 f1<-0.03 f2<-0.04 noise<-15*rnorm(length(t)) amp=1.5 s<-amp*sin(2*pi*f1*t)+amp*sin(2*pi*f2*t) system.time(a<-kzp(s+noise,m=500,k=3)) b<-smooth.kzp(a, smooth_level=0.01) par(mfrow=c(3,1)) plot(periodogram(s+noise),type='l') plot(a) plot(b) par(mfrow=c(1,1)) # signal/noise signal<-kzft(s+noise,m=500,k=3) print(paste("signal-to-noise ratio = ", round(sqrt(var(2*Re(signal))/var(s+noise)),4) )) summary(a, digits=2, top=2) ## End(Not run)## Not run: t<-1:6000 f1<-0.03 f2<-0.04 noise<-15*rnorm(length(t)) amp=1.5 s<-amp*sin(2*pi*f1*t)+amp*sin(2*pi*f2*t) system.time(a<-kzp(s+noise,m=500,k=3)) b<-smooth.kzp(a, smooth_level=0.01) par(mfrow=c(3,1)) plot(periodogram(s+noise),type='l') plot(a) plot(b) par(mfrow=c(1,1)) # signal/noise signal<-kzft(s+noise,m=500,k=3) print(paste("signal-to-noise ratio = ", round(sqrt(var(2*Re(signal))/var(s+noise)),4) )) summary(a, digits=2, top=2) ## End(Not run)
Kolmogorov-Zurbenko Spline
kzs(y,m=NULL,k=3,t=NULL)kzs(y,m=NULL,k=3,t=NULL)
y |
data |
m |
smooth |
k |
The number of iterations for applying the KZFT |
t |
An indexing set |
Kolmogorov-Zurbenko Spline is essentially the Kolmogorov-Zurbenko Fourier Transform at the zero frequency.
I. G. Zurbenko, The spectral Analysis of Time Series. North-Holland, 1986.
I. G. Zurbenko, P. S. Porter, Construction of high-resolution wavelets, Signal Processing 65: 315-327, 1998.
R. H. Shumway, D. S. Stoffer, Time Series Analysis and Its Applications: With R Examples, Springer, 2006.
Derek Cyr and Igor Zurbenko, kzs: Kolmogorov-Zurbenko Spatial Smoothing and Applications, R-Project 2008.
Derek Cyr and Igor Zurbenko, A Spatial Spline Algorithm and an Application to Climate Waves Over the United States, 2008 Joint Statistical Meetings.
kzft,
n <- 1000 x <- (1:n)/n true<-((exp(2.5*x)+sin(25*x))-1)/3 noise <- rnorm(n) y <- true + noise a<-kzs(y,m=60) par(mfrow=c(2,1)) plot(y,type='l') lines(true,col="red") plot(a,type='l', ylim=c(-2,4)) lines(true,col="red") par(mfrow=c(1,1)) ################ # second example ################ t <- seq(from = -round(400*pi), to = round(400*pi), by = .25) ts <- 0.5*sin(sqrt((2*pi*abs(t))/200)) signal <- ifelse(t < 0, -ts, ts) et <- rnorm(length(t), mean = 0, sd = 1) yt <- et + signal b<-kzs(yt,m=400) par(mfrow=c(2,1)) plot(yt,type='l') lines(signal,col="red") plot(b,type='l', ylim=c(-0.5,1)) lines(signal,col="red") par(mfrow=c(1,1))n <- 1000 x <- (1:n)/n true<-((exp(2.5*x)+sin(25*x))-1)/3 noise <- rnorm(n) y <- true + noise a<-kzs(y,m=60) par(mfrow=c(2,1)) plot(y,type='l') lines(true,col="red") plot(a,type='l', ylim=c(-2,4)) lines(true,col="red") par(mfrow=c(1,1)) ################ # second example ################ t <- seq(from = -round(400*pi), to = round(400*pi), by = .25) ts <- 0.5*sin(sqrt((2*pi*abs(t))/200)) signal <- ifelse(t < 0, -ts, ts) et <- rnorm(length(t), mean = 0, sd = 1) yt <- et + signal b<-kzs(yt,m=400) par(mfrow=c(2,1)) plot(yt,type='l') lines(signal,col="red") plot(b,type='l', ylim=c(-0.5,1)) lines(signal,col="red") par(mfrow=c(1,1))
Sample variance of a Kolmogorov-Zurbenko adaptive filter. You want a sigma of at least 3.
kzsv(object) ## S3 method for class 'kzsv' summary(object, digits=getOption("digits"), ...) ## S3 method for class 'kzsv' plot(x, ...)kzsv(object) ## S3 method for class 'kzsv' summary(object, digits=getOption("digits"), ...) ## S3 method for class 'kzsv' plot(x, ...)
object |
The resultant object from kza function. |
x |
The results of the kza function. |
digits |
Precision of output. |
... |
Other parameters. |
x <- c(rep(0,4000),rep(0.5,2000),rep(0,4000)) noise <- rnorm(n = 10000, sd = 1.0) # normally-distributed random variates v <- x + noise a<-kza(v, m=1000, k=3) sv<-kzsv(a)x <- c(rep(0,4000),rep(0.5,2000),rep(0,4000)) noise <- rnorm(n = 10000, sd = 1.0) # normally-distributed random variates v <- x + noise a<-kza(v, m=1000, k=3) sv<-kzsv(a)
Kolmogorov-Zurbenko Third-Order Periodogram for estimating spectrums
kztp(x,m,k,box=c(0,0.5,0,0.5))kztp(x,m,k,box=c(0,0.5,0,0.5))
x |
The signal. |
m |
The window size for the kzft filter. |
k |
The number of iterations. |
box |
The window for the application of third-order periodgram. |
The Kolmogorov-Zurbenko Third-Order Periodogram is used to estimate spectral density of a signal. The smoothing methods are adaptive allowing the bandwidth of the spectral window to vary according to the smoothness of the underlying spectral density. For details, please see to DiRienzo and Zurbenko (1998) and Neagu and Zurbenko (2003).
I. G. Zurbenko, 1986: The spectral Analysis of Time Series. North-Holland, 248 pp. I. G. Zurbenko, P. S. Porter, Construction of high-resolution wavelets, Signal Processing 65: 315-327, 1998. W. Yang, I. G. Zurbenko, A semi-adaptive smoothing algorithm in bispectrum estimation, Proceedings of the American Statistical Association, Seattle, 2006. Wei Yang and Igor Zurbenko, kzft: Kolmogorov-Zurbenko Fourier Transform and Applications, R-Project 2007.
t<-1:10000 y<-2*sin(2*pi*0.1*t)+3*sin(2*pi*0.2*t) + 10*rnorm(length(t)) a<-kztp(y,50,1) z<-log(Mod(a)) #z<-smooth.kzp(z) omega<-seq(0,1,length=51)[2:26] #filled.contour(omega,omega,s,xlab="freq",ylab="freq",main="Smoothed 3rd Order Periodogram")t<-1:10000 y<-2*sin(2*pi*0.1*t)+3*sin(2*pi*0.2*t) + 10*rnorm(length(t)) a<-kztp(y,50,1) z<-log(Mod(a)) #z<-smooth.kzp(z) omega<-seq(0,1,length=51)[2:26] #filled.contour(omega,omega,s,xlab="freq",ylab="freq",main="Smoothed 3rd Order Periodogram")
Raw periodogram.
periodogram(y)periodogram(y)
y |
The raw data. |
Periodogram is an estimate of the spectral density using FFT.
kzp,
t<-1:1000 f1<-0.3 f2<-0.4 noise<-15*rnorm(length(t)) s<-3*sin(2*pi*f1*t)+3*sin(2*pi*f2*t) plot(periodogram(s+noise),type='l')t<-1:1000 f1<-0.3 f2<-0.4 noise<-15*rnorm(length(t)) s<-3*sin(2*pi*f1*t)+3*sin(2*pi*f2*t) plot(periodogram(s+noise),type='l')
Rolling variance function, (rlv) will detect data variations of 1-dimensional vector, 2-dimensional matrix and/or 3-dimensional array, which in turn display their local discontinuities (boundaries).
rlv(inpt, krnl)rlv(inpt, krnl)
inpt |
input data. |
krnl |
length of the rolling window for detecting. |
The rlv is calculated based on a dynamic rolling window. Within a given window, directly estimating its variance allows to capture the local variations of the signal underlying analysis. Since this window is moving, when its size is optimal, the significant local variance atlas will be captured. It is neither necessary to make any assumption, nor to model anything, so it is a nonparametric statistic. However, direct application of rlv to data embedded in heavy background noises could fail to identify the significant local variations. In situations when the estimated data are embedded in background noise, necessary smoothness may be expected.
Igor G Zurbenko, Mingzeng Sun. Estimation of spatial boundaries with rolling variance and 2D KZA algorithm. Biometrics & Biostatistics International Journal. 7(4):263.270, 2018
## A. Signal/objects # to create a cylinder circleFun <- function(center = c(0,0),diameter = 1, npoints = 200){ r = diameter / 2 tt <- seq(0,2*pi,length.out = npoints) xx <- center[1] + r * cos(tt) yy <- center[2] + r * sin(tt) return(data.frame(x = xx, y = yy)) } rxo1 <- circleFun(center = c(38, 38), diameter = 38, npoints = 200) rxi1 <- circleFun(center = c(38, 38), diameter = 18, npoints = 200) rxo2 <- circleFun(center = c(58, 58), diameter = 38, npoints = 200) rxi2 <- circleFun(center = c(58, 58), diameter = 18, npoints = 200) rxo3 <- circleFun(center = c(28, 78), diameter = 18, npoints = 200) rxi3 <- circleFun(center = c(28, 78), diameter = 8, npoints = 200) set.seed(3) s2 <- matrix(rep(0,100*100),nrow=100) for(i in 1:nrow(rxo1)) { s2[rxo1[i,2],38:rxo1[i,1]]=0.5 } for(i in 1:nrow(rxi1)) { s2[rxi1[i,2],38:rxi1[i,1]]=1 } # second cylinder for(i in 1:nrow(rxo2)) { s2[rxo2[i,2],58:rxo2[i,1]]=0.5 } for(i in 1:nrow(rxi2)) { s2[rxi2[i,2],58:rxi2[i,1]]=1 } #third separated cylinder for(i in 1:nrow(rxo3)) { s2[rxo3[i,2],28:rxo3[i,1]]=0.3 } for(i in 1:nrow(rxi3)) { s2[rxi3[i,2],28:rxi3[i,1]]=0.6 } signal=s2 ## Graphing Library plotly ## Not run: plot_ly(z=signal, type="surface") ## End(Not run) ## B. Rolling variance atlas rolv=rlv(s2,3) ## Not run: plot_ly(z=rolv, type="surface") ## End(Not run)## A. Signal/objects # to create a cylinder circleFun <- function(center = c(0,0),diameter = 1, npoints = 200){ r = diameter / 2 tt <- seq(0,2*pi,length.out = npoints) xx <- center[1] + r * cos(tt) yy <- center[2] + r * sin(tt) return(data.frame(x = xx, y = yy)) } rxo1 <- circleFun(center = c(38, 38), diameter = 38, npoints = 200) rxi1 <- circleFun(center = c(38, 38), diameter = 18, npoints = 200) rxo2 <- circleFun(center = c(58, 58), diameter = 38, npoints = 200) rxi2 <- circleFun(center = c(58, 58), diameter = 18, npoints = 200) rxo3 <- circleFun(center = c(28, 78), diameter = 18, npoints = 200) rxi3 <- circleFun(center = c(28, 78), diameter = 8, npoints = 200) set.seed(3) s2 <- matrix(rep(0,100*100),nrow=100) for(i in 1:nrow(rxo1)) { s2[rxo1[i,2],38:rxo1[i,1]]=0.5 } for(i in 1:nrow(rxi1)) { s2[rxi1[i,2],38:rxi1[i,1]]=1 } # second cylinder for(i in 1:nrow(rxo2)) { s2[rxo2[i,2],58:rxo2[i,1]]=0.5 } for(i in 1:nrow(rxi2)) { s2[rxi2[i,2],58:rxi2[i,1]]=1 } #third separated cylinder for(i in 1:nrow(rxo3)) { s2[rxo3[i,2],28:rxo3[i,1]]=0.3 } for(i in 1:nrow(rxi3)) { s2[rxi3[i,2],28:rxi3[i,1]]=0.6 } signal=s2 ## Graphing Library plotly ## Not run: plot_ly(z=signal, type="surface") ## End(Not run) ## B. Rolling variance atlas rolv=rlv(s2,3) ## Not run: plot_ly(z=rolv, type="surface") ## End(Not run)