| Title: | Basic Wavelet Routines for One-, Two-, and Three-Dimensional Signal Processing |
|---|---|
| Description: | Basic wavelet routines for time series (1D), image (2D) and array (3D) analysis. The code provided here is based on wavelet methodology developed in Percival and Walden (2000); Gencay, Selcuk and Whitcher (2001); the dual-tree complex wavelet transform (DTCWT) from Kingsbury (1999, 2001) as implemented by Selesnick; and Hilbert wavelet pairs (Selesnick 2001, 2002). All figures in chapters 4-7 of GSW (2001) are reproducible using this package and R code available at the book website(s) below. |
| Authors: | Brandon Whitcher |
| Maintainer: | Brandon Whitcher <[email protected]> |
| License: | BSD_3_clause + file LICENSE |
| Version: | 1.8.5 |
| Built: | 2024-11-01 05:54:44 UTC |
| Source: | https://github.com/bjw34032/waveslim |
The autocovariance and autocorrelation sequences from the time series model in Figures 8, 9, 10, and 11 of Andel (1986). They were obtained through numeric integration of the spectral density function.
data(acvs.andel8) data(acvs.andel9) data(acvs.andel10) data(acvs.andel11)data(acvs.andel8) data(acvs.andel9) data(acvs.andel10) data(acvs.andel11)
A data frame with 4096 rows and three columns: lag, autocovariance sequence, autocorrelation sequence.
Andel, J. (1986) Long memory time series models, Kypernetika, 22, No. 2, 105-123.
Analysis and synthesis filter banks used in dual-tree wavelet algorithms.
afb(x, af) afb2D(x, af1, af2 = NULL) afb2D.A(x, af, d) sfb(lo, hi, sf) sfb2D(lo, hi, sf1, sf2 = NULL) sfb2D.A(lo, hi, sf, d)afb(x, af) afb2D(x, af1, af2 = NULL) afb2D.A(x, af, d) sfb(lo, hi, sf) sfb2D(lo, hi, sf1, sf2 = NULL) sfb2D.A(lo, hi, sf, d)
x |
vector or matrix of observations |
af |
analysis filters. First element of the list is the low-pass filter, second element is the high-pass filter. |
af1, af2
|
analysis filters for the first and second dimension of a 2D array. |
sf |
synthesis filters. First element of the list is the low-pass filter, second element is the high-pass filter. |
sf1, sf2
|
synthesis filters for the first and second dimension of a 2D array. |
d |
dimension of filtering (d = 1 or 2) |
lo |
low-frequecy coefficients |
hi |
high-frequency coefficients |
The functions afb2D.A and sfb2D.A implement the convolutions,
either for analysis or synthesis, in one dimension only. Thus, they are the
workhorses of afb2D and sfb2D. The output for the analysis
filter bank along one dimension (afb2D.A) is a list with two elements
low-pass subband
high-pass subband
where
the dimension of analysis will be half its original length. The output for
the synthesis filter bank along one dimension (sfb2D.A) will be the
output array, where the dimension of synthesis will be twice its original
length.
In one dimension the output for the analysis filter bank
(afb) is a list with two elements
lo |
Low frequecy output |
hi |
High frequency output |
and the output for the synthesis filter
bank (sfb) is the output signal.
In two dimensions the output for the analysis filter bank (afb2D) is
a list with four elements
lo |
low-pass subband |
hi[[1]] |
'lohi' subband |
hi[[2]] |
'hilo' subband |
hi[[3]] |
'hihi' subband |
and
the output for the synthesis filter bank (sfb2D) is the output array.
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
## EXAMPLE: afb, sfb af = farras()$af sf = farras()$sf x = rnorm(64) x.afb = afb(x, af) lo = x.afb$lo hi = x.afb$hi y = sfb(lo, hi, sf) err = x - y max(abs(err)) ## EXAMPLE: afb2D, sfb2D x = matrix(rnorm(32*64), 32, 64) af = farras()$af sf = farras()$sf x.afb2D = afb2D(x, af, af) lo = x.afb2D$lo hi = x.afb2D$hi y = sfb2D(lo, hi, sf, sf) err = x - y max(abs(err)) ## Example: afb2D.A, sfb2D.A x = matrix(rnorm(32*64), 32, 64) af = farras()$af sf = farras()$sf x.afb2D.A = afb2D.A(x, af, 1) lo = x.afb2D.A$lo hi = x.afb2D.A$hi y = sfb2D.A(lo, hi, sf, 1) err = x - y max(abs(err))## EXAMPLE: afb, sfb af = farras()$af sf = farras()$sf x = rnorm(64) x.afb = afb(x, af) lo = x.afb$lo hi = x.afb$hi y = sfb(lo, hi, sf) err = x - y max(abs(err)) ## EXAMPLE: afb2D, sfb2D x = matrix(rnorm(32*64), 32, 64) af = farras()$af sf = farras()$sf x.afb2D = afb2D(x, af, af) lo = x.afb2D$lo hi = x.afb2D$hi y = sfb2D(lo, hi, sf, sf) err = x - y max(abs(err)) ## Example: afb2D.A, sfb2D.A x = matrix(rnorm(32*64), 32, 64) af = farras()$af sf = farras()$sf x.afb2D.A = afb2D.A(x, af, 1) lo = x.afb2D.A$lo hi = x.afb2D.A$hi y = sfb2D.A(lo, hi, sf, 1) err = x - y max(abs(err))
Simulated AR(1) series used in Gencay, Selcuk and Whitcher (2001).
data(ar1)data(ar1)
A vector containing 200 observations.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Computes the band-pass variance for fractional difference (FD) or seasonal persistent (SP) processes using numeric integration of their spectral density function.
bandpass.fdp(a, b, d) bandpass.spp(a, b, d, fG) bandpass.spp2(a, b, d1, f1, d2, f2) bandpass.var.spp(delta, fG, J, Basis, Length)bandpass.fdp(a, b, d) bandpass.spp(a, b, d, fG) bandpass.spp2(a, b, d1, f1, d2, f2) bandpass.var.spp(delta, fG, J, Basis, Length)
fG, f1, f2
|
Gegenbauer frequency. |
J |
Depth of the wavelet transform. |
Basis |
Logical vector representing the adaptive basis. |
Length |
Number of elements in Basis. |
a |
Left-hand boundary for the definite integral. |
b |
Right-hand boundary for the definite integral. |
d, delta, d1, d2
|
Fractional difference parameter. |
See references.
Band-pass variance for the FD or SP process between and
.
B. Whitcher
McCoy, E. J., and A. T. Walden (1996) Wavelet analysis and synthesis of stationary long-memory processes, Journal for Computational and Graphical Statistics, 5, No. 1, 26-56.
Whitcher, B. (2001) Simulating Gaussian stationary processes with unbounded spectra, Journal for Computational and Graphical Statistics, 10, No. 1, 112-134.
The Barbara image comes from Allen Gersho's lab at the University of California, Santa Barbara.
data(barbara)data(barbara)
A 256 256 matrix.
Internet.
Produce a vector of zeros and ones from a vector of basis names.
basis(x, basis.names)basis(x, basis.names)
x |
Output from the discrete wavelet package transfrom (DWPT). |
basis.names |
Vector of character strings that describe leaves on the DWPT basis tree. See the examples below for appropriate syntax. |
None.
Vector of zeros and ones.
dwpt.
data(acvs.andel8) ## Not run: x <- hosking.sim(1024, acvs.andel8[,2]) x.dwpt <- dwpt(x, "la8", 7) ## Select orthonormal basis from wavelet packet tree x.basis <- basis(x.dwpt, c("w1.1","w2.1","w3.0","w4.3","w5.4","w6.10", "w7.22","w7.23")) for(i in 1:length(x.dwpt)) x.dwpt[[i]] <- x.basis[i] * x.dwpt[[i]] ## Resonstruct original series using selected orthonormal basis y <- idwpt(x.dwpt, x.basis) par(mfrow=c(2,1), mar=c(5-1,4,4-1,2)) plot.ts(x, xlab="", ylab="", main="Original Series") plot.ts(y, xlab="", ylab="", main="Reconstructed Series") ## End(Not run)data(acvs.andel8) ## Not run: x <- hosking.sim(1024, acvs.andel8[,2]) x.dwpt <- dwpt(x, "la8", 7) ## Select orthonormal basis from wavelet packet tree x.basis <- basis(x.dwpt, c("w1.1","w2.1","w3.0","w4.3","w5.4","w6.10", "w7.22","w7.23")) for(i in 1:length(x.dwpt)) x.dwpt[[i]] <- x.basis[i] * x.dwpt[[i]] ## Resonstruct original series using selected orthonormal basis y <- idwpt(x.dwpt, x.basis) par(mfrow=c(2,1), mar=c(5-1,4,4-1,2)) plot.ts(x, xlab="", ylab="", main="Original Series") plot.ts(y, xlab="", ylab="", main="Reconstructed Series") ## End(Not run)
data(blocks)data(blocks)
A vector containing 512 observations.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
Sets the first wavelet coefficients to NA.
brick.wall(x, wf, method = "modwt") dwpt.brick.wall(x, wf, n.levels, method = "modwpt") brick.wall.2d(x, method = "modwt")brick.wall(x, wf, method = "modwt") dwpt.brick.wall(x, wf, n.levels, method = "modwpt") brick.wall.2d(x, method = "modwt")
x |
DWT/MODWT/DWPT/MODWPT object |
wf |
Character string; name of wavelet filter |
method |
Either |
n.levels |
Specifies the depth of the decomposition. This must be a number less than or equal to log(length(x),2). |
The fact that observed time series are finite causes boundary issues. One
way to get around this is to simply remove any wavelet coefficient computed
involving the boundary. This is done here by replacing boundary wavelet
coefficients with NA.
Same object as x only with some missing values.
B. Whitcher
Lindsay, R. W., D. B. Percival and D. A. Rothrock (1996). The discrete wavelet transform and the scale anlaysis of the surface properties of sea ice, IEEE Transactions on Geoscience and Remote Sensing, 34, No. 3, 771-787.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Use the Fast Fourier Transform to perform convolutions between a sequence and each column of a matrix.
convolve2D(x, y, conj = TRUE, type = c("circular", "open"))convolve2D(x, y, conj = TRUE, type = c("circular", "open"))
x |
MxN matrix. |
y |
numeric sequence of length N. |
conj |
logical; if |
type |
character; one of For |
This is a corrupted version of convolve made by replacing fft
with mvfft in a few places. It would be nice to submit this to the R
Developers for inclusion.
B. Whitcher
Monthly U.S. consumer price index from 1948:1 to 1999:12.
data(cpi)data(cpi)
A vector containing 624 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Dual-tree complex 2D discrete wavelet transform (DWT).
cplxdual2D(x, J, Faf, af) icplxdual2D(w, J, Fsf, sf)cplxdual2D(x, J, Faf, af) icplxdual2D(w, J, Fsf, sf)
x |
2D array. |
J |
number of stages. |
Faf |
first stage analysis filters for tree i. |
af |
analysis filters for the remaining stages on tree i. |
w |
wavelet coefficients. |
Fsf |
last stage synthesis filters for tree i. |
sf |
synthesis filters for the preceeding stages. |
For the analysis of x, the output is
w |
wavelet
coefficients indexed by |
For the synthesis of
w, the output is
y |
output signal. |
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
FSfarras, farras, afb2D,
sfb2D.
## Not run: ## EXAMPLE: cplxdual2D x = matrix(rnorm(32*32), 32, 32) J = 5 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = cplxdual2D(x, J, Faf, af) y = icplxdual2D(w, J, Fsf, sf) err = x - y max(abs(err)) ## End(Not run)## Not run: ## EXAMPLE: cplxdual2D x = matrix(rnorm(32*32), 32, 32) J = 5 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = cplxdual2D(x, J, Faf, af) y = icplxdual2D(w, J, Fsf, sf) err = x - y max(abs(err)) ## End(Not run)
Miscellaneous functions for dual-tree wavelet software.
cshift(x, m) cshift2D(x, m) pm(a, b)cshift(x, m) cshift2D(x, m) pm(a, b)
x |
N-point vector |
m |
amount of shift |
a, b
|
input parameters |
y |
vector |
u |
(a + b) / sqrt(2) |
v |
(a - b) / sqrt(2) |
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
A wavelet packet tree, from the discrete wavelet packet transform (DWPT), is tested node-by-node for white noise. This is the first step in selecting an orthonormal basis for the DWPT.
cpgram.test(y, p = 0.05, taper = 0.1) css.test(y) entropy.test(y) portmanteau.test(y, p = 0.05, type = "Box-Pierce")cpgram.test(y, p = 0.05, taper = 0.1) css.test(y) entropy.test(y) portmanteau.test(y, p = 0.05, type = "Box-Pierce")
y |
wavelet packet tree (from the DWPT) |
p |
significance level |
taper |
weight of cosine bell taper ( |
type |
|
Top-down recursive testing of the wavelet packet tree is
Boolean vector of the same length as the number of nodes in the wavelet packet tree.
B. Whitcher
Brockwell and Davis (1991) Time Series: Theory and Methods, (2nd. edition), Springer-Verlag.
Brown, Durbin and Evans (1975) Techniques for testing the constancy of regression relationships over time, Journal of the Royal Statistical Society B, 37, 149-163.
Percival, D. B., and A. T. Walden (1993) Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press.
data(mexm) J <- 6 wf <- "la8" mexm.dwpt <- dwpt(mexm[-c(1:4)], wf, J) ## Not implemented yet ## plot.dwpt(x.dwpt, J) mexm.dwpt.bw <- dwpt.brick.wall(mexm.dwpt, wf, 6, method="dwpt") mexm.tree <- ortho.basis(portmanteau.test(mexm.dwpt.bw, p=0.025)) ## Not implemented yet ## plot.basis(mexm.tree)data(mexm) J <- 6 wf <- "la8" mexm.dwpt <- dwpt(mexm[-c(1:4)], wf, J) ## Not implemented yet ## plot.dwpt(x.dwpt, J) mexm.dwpt.bw <- dwpt.brick.wall(mexm.dwpt, wf, 6, method="dwpt") mexm.tree <- ortho.basis(portmanteau.test(mexm.dwpt.bw, p=0.025)) ## Not implemented yet ## plot.basis(mexm.tree)
A digital photograph of Ingrid Daubechies taken at the 1993 AMS winter meetings in San Antonio, Texas. The photograph was taken by David Donoho with a Canon XapShot video still frame camera.
data(dau)data(dau)
A 256 256 matrix.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
Perform simple de-noising of an image using the two-dimensional discrete wavelet transform.
denoise.dwt.2d( x, wf = "la8", J = 4, method = "universal", H = 0.5, noise.dir = 3, rule = "hard" )denoise.dwt.2d( x, wf = "la8", J = 4, method = "universal", H = 0.5, noise.dir = 3, rule = "hard" )
x |
input matrix (image) |
wf |
name of the wavelet filter to use in the decomposition |
J |
depth of the decomposition, must be a number less than or equal to log(minM,N,2) |
method |
character string describing the threshold applied, only
|
H |
self-similarity or Hurst parameter to indicate spectral scaling, white noise is 0.5 |
noise.dir |
number of directions to estimate background noise standard deviation, the default is 3 which produces a unique estimate of the background noise for each spatial direction |
rule |
either a |
See Thresholding.
Image of the same dimension as the original but with high-freqency fluctuations removed.
B. Whitcher
See Thresholding for references concerning
de-noising in one dimension.
## Xbox image data(xbox) n <- nrow(xbox) xbox.noise <- xbox + matrix(rnorm(n*n, sd=.15), n, n) par(mfrow=c(2,2), cex=.8, pty="s") image(xbox.noise, col=rainbow(128), main="Original Image") image(denoise.dwt.2d(xbox.noise, wf="haar"), col=rainbow(128), zlim=range(xbox.noise), main="Denoised image") image(xbox.noise - denoise.dwt.2d(xbox.noise, wf="haar"), col=rainbow(128), zlim=range(xbox.noise), main="Residual image") ## Daubechies image data(dau) n <- nrow(dau) dau.noise <- dau + matrix(rnorm(n*n, sd=10), n, n) par(mfrow=c(2,2), cex=.8, pty="s") image(dau.noise, col=rainbow(128), main="Original Image") dau.denoise <- denoise.modwt.2d(dau.noise, wf="d4", rule="soft") image(dau.denoise, col=rainbow(128), zlim=range(dau.noise), main="Denoised image") image(dau.noise - dau.denoise, col=rainbow(128), main="Residual image")## Xbox image data(xbox) n <- nrow(xbox) xbox.noise <- xbox + matrix(rnorm(n*n, sd=.15), n, n) par(mfrow=c(2,2), cex=.8, pty="s") image(xbox.noise, col=rainbow(128), main="Original Image") image(denoise.dwt.2d(xbox.noise, wf="haar"), col=rainbow(128), zlim=range(xbox.noise), main="Denoised image") image(xbox.noise - denoise.dwt.2d(xbox.noise, wf="haar"), col=rainbow(128), zlim=range(xbox.noise), main="Residual image") ## Daubechies image data(dau) n <- nrow(dau) dau.noise <- dau + matrix(rnorm(n*n, sd=10), n, n) par(mfrow=c(2,2), cex=.8, pty="s") image(dau.noise, col=rainbow(128), main="Original Image") dau.denoise <- denoise.modwt.2d(dau.noise, wf="d4", rule="soft") image(dau.denoise, col=rainbow(128), zlim=range(dau.noise), main="Denoised image") image(dau.noise - dau.denoise, col=rainbow(128), main="Residual image")
data(doppler)data(doppler)
A vector containing 512 observations.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
This is now a wrapper to the function multitaper::dpss().
dpss.taper(n, k, nw = 4)dpss.taper(n, k, nw = 4)
n |
length of data taper(s) |
k |
number of data tapers; 1, 2, 3, ... (do not use 0!) |
nw |
product of length and half-bandwidth parameter (w) |
v |
matrix of data tapers (cols = tapers) |
eigen |
eigenvalue associated with each data taper, discarded |
B. Whitcher
Percival, D. B. and A. T. Walden (1993) Spectral Estimation for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press.
Kingsbury's Q-filters for the dual-tree complex DWT.
dualfilt1()dualfilt1()
These cofficients are rounded to 8 decimal places.
af |
List (i=1,2) - analysis filters for tree i |
sf |
List (i=1,2) - synthesis filters for tree i |
Note: af[[2]] is the reverse
of af[[1]].
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
Kingsbury, N.G. (2000). A dual-tree complex wavelet transform with improved orthogonality and symmetry properties, Proceedings of the IEEE Int. Conf. on Image Proc. (ICIP).
One- and two-dimensional dual-tree complex discrete wavelet transforms developed by Kingsbury and Selesnick et al.
dualtree(x, J, Faf, af) idualtree(w, J, Fsf, sf) dualtree2D(x, J, Faf, af) idualtree2D(w, J, Fsf, sf)dualtree(x, J, Faf, af) idualtree(w, J, Fsf, sf) dualtree2D(x, J, Faf, af) idualtree2D(w, J, Fsf, sf)
x |
N-point vector or MxN matrix. |
J |
number of stages. |
Faf |
analysis filters for the first stage. |
af |
analysis filters for the remaining stages. |
w |
DWT coefficients. |
Fsf |
synthesis filters for the last stage. |
sf |
synthesis filters for the preceeding stages. |
In one dimension is divisible by and
.
In two dimensions, these two conditions must hold for both and
.
For the analysis of x, the output is
w |
DWT coefficients. Each wavelet scale is a list containing the real and imaginary parts. The final scale (J+1) contains the low-pass filter coefficients. |
For the synthesis of w, the output is
y |
output signal |
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
FSfarras, farras,
convolve, cshift, afb,
sfb.
## EXAMPLE: dualtree x = rnorm(512) J = 4 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = dualtree(x, J, Faf, af) y = idualtree(w, J, Fsf, sf) err = x - y max(abs(err)) ## Example: dualtree2D x = matrix(rnorm(64*64), 64, 64) J = 3 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = dualtree2D(x, J, Faf, af) y = idualtree2D(w, J, Fsf, sf) err = x - y max(abs(err)) ## Display 2D wavelets of dualtree2D.m J <- 4 L <- 3 * 2^(J+1) N <- L / 2^J Faf <- FSfarras()$af Fsf <- FSfarras()$sf af <- dualfilt1()$af sf <- dualfilt1()$sf x <- matrix(0, 2*L, 3*L) w <- dualtree2D(x, J, Faf, af) w[[J]][[1]][[1]][N/2, N/2+0*N] <- 1 w[[J]][[1]][[2]][N/2, N/2+1*N] <- 1 w[[J]][[1]][[3]][N/2, N/2+2*N] <- 1 w[[J]][[2]][[1]][N/2+N, N/2+0*N] <- 1 w[[J]][[2]][[2]][N/2+N, N/2+1*N] <- 1 w[[J]][[2]][[3]][N/2+N, N/2+2*N] <- 1 y <- idualtree2D(w, J, Fsf, sf) image(t(y), col=grey(0:64/64), axes=FALSE)## EXAMPLE: dualtree x = rnorm(512) J = 4 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = dualtree(x, J, Faf, af) y = idualtree(w, J, Fsf, sf) err = x - y max(abs(err)) ## Example: dualtree2D x = matrix(rnorm(64*64), 64, 64) J = 3 Faf = FSfarras()$af Fsf = FSfarras()$sf af = dualfilt1()$af sf = dualfilt1()$sf w = dualtree2D(x, J, Faf, af) y = idualtree2D(w, J, Fsf, sf) err = x - y max(abs(err)) ## Display 2D wavelets of dualtree2D.m J <- 4 L <- 3 * 2^(J+1) N <- L / 2^J Faf <- FSfarras()$af Fsf <- FSfarras()$sf af <- dualfilt1()$af sf <- dualfilt1()$sf x <- matrix(0, 2*L, 3*L) w <- dualtree2D(x, J, Faf, af) w[[J]][[1]][[1]][N/2, N/2+0*N] <- 1 w[[J]][[1]][[2]][N/2, N/2+1*N] <- 1 w[[J]][[1]][[3]][N/2, N/2+2*N] <- 1 w[[J]][[2]][[1]][N/2+N, N/2+0*N] <- 1 w[[J]][[2]][[2]][N/2+N, N/2+1*N] <- 1 w[[J]][[2]][[3]][N/2+N, N/2+2*N] <- 1 y <- idualtree2D(w, J, Fsf, sf) image(t(y), col=grey(0:64/64), axes=FALSE)
All possible filtering combinations (low- and high-pass) are performed to decompose a vector or time series. The resulting coefficients are associated with a binary tree structure corresponding to a partitioning of the frequency axis.
dwpt(x, wf = "la8", n.levels = 4, boundary = "periodic") idwpt(y, y.basis)dwpt(x, wf = "la8", n.levels = 4, boundary = "periodic") idwpt(y, y.basis)
x |
a vector or time series containing the data be to decomposed. This must be a dyadic length vector (power of 2). |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
n.levels |
Specifies the depth of the decomposition.This must be a
number less than or equal to
|
boundary |
Character string specifying the boundary condition. If
|
y |
Object of S3 class |
y.basis |
Vector of character strings that describe leaves on the DWPT basis tree. |
The code implements the one-dimensional DWPT using the pyramid algorithm (Mallat, 1989).
Basically, a list with the following components
w?.? |
Wavelet coefficient vectors. The first index is associated with the scale of the decomposition while the second is associated with the frequency partition within that level. |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7), 674–693.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Wickerhauser, M. V. (1994) Adapted Wavelet Analysis from Theory to Software, A K Peters.
dwt, modwpt, wave.filter.
data(mexm) J <- 4 mexm.mra <- mra(log(mexm), "mb8", J, "modwt", "reflection") mexm.nomean <- ts( apply(matrix(unlist(mexm.mra), ncol=J+1, byrow=FALSE)[,-(J+1)], 1, sum), start=1957, freq=12) mexm.dwpt <- dwpt(mexm.nomean[-c(1:4)], "mb8", 7, "reflection")data(mexm) J <- 4 mexm.mra <- mra(log(mexm), "mb8", J, "modwt", "reflection") mexm.nomean <- ts( apply(matrix(unlist(mexm.mra), ncol=J+1, byrow=FALSE)[,-(J+1)], 1, sum), start=1957, freq=12) mexm.dwpt <- dwpt(mexm.nomean[-c(1:4)], "mb8", 7, "reflection")
All possible filtering combinations (low- and high-pass) are performed to decompose a matrix or image. The resulting coefficients are associated with a quad-tree structure corresponding to a partitioning of the two-dimensional frequency plane.
dwpt.2d(x, wf = "la8", J = 4, boundary = "periodic") idwpt.2d(y, y.basis)dwpt.2d(x, wf = "la8", J = 4, boundary = "periodic") idwpt.2d(y, y.basis)
x |
a matrix or image containing the data be to decomposed. This ojbect must be dyadic (power of 2) in length in each dimension. |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
J |
Specifies the depth of the decomposition. This must be a number
less than or equal to |
boundary |
Character string specifying the boundary condition. If
|
y |
|
y.basis |
Boolean vector, the same length as |
The code implements the two-dimensional DWPT using the pyramid algorithm of Mallat (1989).
Basically, a list with the following components
w?.?-w?.? |
Wavelet coefficient matrices (images). The first index is
associated with the scale of the decomposition while the second is
associated with the frequency partition within that level. The left and
right strings, separated by the dash ‘-’, correspond to the first |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, No. 7, 674-693.
Wickerhauser, M. V. (1994) Adapted Wavelet Analysis from Theory to Software, A K Peters.
dwt.2d, modwt.2d,
wave.filter.
An adaptive orthonormal basis is selected in order to perform the naive bootstrap within nodes of the wavelet packet tree. A bootstrap realization of the time series is produce by applying the inverse DWPT.
dwpt.boot(y, wf, J = log(length(y), 2) - 1, p = 1e-04, frac = 1)dwpt.boot(y, wf, J = log(length(y), 2) - 1, p = 1e-04, frac = 1)
y |
Not necessarily dyadic length time series. |
wf |
Name of the wavelet filter to use in the decomposition. See
|
J |
Depth of the discrete wavelet packet transform. |
p |
Level of significance for the white noise testing procedure. |
frac |
Fraction of the time series that should be used in constructing the likelihood function. |
A subroutines is used to select an adaptive orthonormal basis for the piecewise-constant approximation to the underlying spectral density function (SDF). Once selected, sampling with replacement is performed within each wavelet packet coefficient vector and the new collection of wavelet packet coefficients are reconstructed into a bootstrap realization of the original time series.
Time series of length $N$, where $N$ is the length of y.
B. Whitcher
Percival, D.B., S. Sardy, A. Davision (2000) Wavestrapping Time Series: Adaptive Wavelet-Based Bootstrapping, in B.J. Fitzgerald, R.L. Smith, A.T. Walden, P.C. Young (Eds.) Nonlinear and Nonstationary Signal Processing, pp. 442-471.
Whitcher, B. (2001) Simulating Gaussian Stationary Time Series with Unbounded Spectra, Journal of Computational and Graphical Statistics, 10, No. 1, 112-134.
Whitcher, B. (2004) Wavelet-Based Estimation for Seasonal Long-Memory Processes, Technometrics, 46, No. 2, 225-238.
A seasonal persistent process may be characterized by a spectral density
function with an asymptote occuring at a particular frequency in
. It's time domain representation was first
noted in passing by Hosking (1981). Although an exact time-domain approach
to simulation is possible, this function utilizes the discrete wavelet
packet transform (DWPT).
dwpt.sim(N, wf, delta, fG, M = 2, adaptive = TRUE, epsilon = 0.05)dwpt.sim(N, wf, delta, fG, M = 2, adaptive = TRUE, epsilon = 0.05)
N |
Length of time series to be generated. |
wf |
Character string for the wavelet filter. |
delta |
Long-memory parameter for the seasonal persistent process. |
fG |
Gegenbauer frequency. |
M |
Actual length of simulated time series. |
adaptive |
Logical; if |
epsilon |
Threshold for adaptive basis selection. |
Two subroutines are used, the first selects an adaptive orthonormal basis
for the true spectral density function (SDF) while the second computes the
bandpass variances associated with the chosen orthonormal basis and SDF.
Finally, when a uniform random variable is
generated in order to select a random piece of the simulated time series.
For more details see Whitcher (2001).
Time series of length N.
B. Whitcher
Hosking, J. R. M. (1981) Fractional Differencing, Biometrika, 68, No. 1, 165-176.
Whitcher, B. (2001) Simulating Gaussian Stationary Time Series with Unbounded Spectra, Journal of Computational and Graphical Statistics, 10, No. 1, 112-134.
hosking.sim for an exact time-domain method and
wave.filter for a list of available wavelet filters.
## Generate monthly time series with annual oscillation ## library(ts) is required in order to access acf() x <- dwpt.sim(256, "mb16", .4, 1/12, M=4, epsilon=.001) par(mfrow=c(2,1)) plot(x, type="l", xlab="Time") acf(x, lag.max=128, ylim=c(-.6,1)) data(acvs.andel8) lines(acvs.andel8$lag[1:128], acvs.andel8$acf[1:128], col=2)## Generate monthly time series with annual oscillation ## library(ts) is required in order to access acf() x <- dwpt.sim(256, "mb16", .4, 1/12, M=4, epsilon=.001) par(mfrow=c(2,1)) plot(x, type="l", xlab="Time") acf(x, lag.max=128, ylim=c(-.6,1)) data(acvs.andel8) lines(acvs.andel8$lag[1:128], acvs.andel8$acf[1:128], col=2)
This function performs a level decomposition of the input vector or
time series using the pyramid algorithm (Mallat 1989).
dwt(x, wf = "la8", n.levels = 4, boundary = "periodic") dwt.nondyadic(x) idwt(y)dwt(x, wf = "la8", n.levels = 4, boundary = "periodic") dwt.nondyadic(x) idwt(y)
x |
a vector or time series containing the data be to decomposed. This must be a dyadic length vector (power of 2). |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
n.levels |
Specifies the depth of the decomposition. This must be a number less than or equal to log(length(x),2). |
boundary |
Character string specifying the boundary condition. If
|
y |
An object of S3 class |
The code implements the one-dimensional DWT using the pyramid algorithm (Mallat, 1989). The actual transform is performed in C using pseudocode from Percival and Walden (2001). That means convolutions, not inner products, are used to apply the wavelet filters.
For a non-dyadic length vector or time series, dwt.nondyadic pads
with zeros, performs the orthonormal DWT on this dyadic length series and
then truncates the wavelet coefficient vectors appropriately.
Basically, a list with the following components
d? |
Wavelet coefficient vectors. |
s? |
Scaling coefficient vector. |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Daubechies, I. (1992) Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM: Philadelphia.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7), 674–693.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
## Figures 4.17 and 4.18 in Gencay, Selcuk and Whitcher (2001). data(ibm) ibm.returns <- diff(log(ibm)) ## Haar ibmr.haar <- dwt(ibm.returns, "haar") names(ibmr.haar) <- c("w1", "w2", "w3", "w4", "v4") ## plot partial Haar DWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(a)") for(i in 1:4) plot.ts(up.sample(ibmr.haar[[i]], 2^i), type="h", axes=FALSE, ylab=names(ibmr.haar)[i]) plot.ts(up.sample(ibmr.haar$v4, 2^4), type="h", axes=FALSE, ylab=names(ibmr.haar)[5]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) ## LA(8) ibmr.la8 <- dwt(ibm.returns, "la8") names(ibmr.la8) <- c("w1", "w2", "w3", "w4", "v4") ## must shift LA(8) coefficients ibmr.la8$w1 <- c(ibmr.la8$w1[-c(1:2)], ibmr.la8$w1[1:2]) ibmr.la8$w2 <- c(ibmr.la8$w2[-c(1:2)], ibmr.la8$w2[1:2]) for(i in names(ibmr.la8)[3:4]) ibmr.la8[[i]] <- c(ibmr.la8[[i]][-c(1:3)], ibmr.la8[[i]][1:3]) ibmr.la8$v4 <- c(ibmr.la8$v4[-c(1:2)], ibmr.la8$v4[1:2]) ## plot partial LA(8) DWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(b)") for(i in 1:4) plot.ts(up.sample(ibmr.la8[[i]], 2^i), type="h", axes=FALSE, ylab=names(ibmr.la8)[i]) plot.ts(up.sample(ibmr.la8$v4, 2^4), type="h", axes=FALSE, ylab=names(ibmr.la8)[5]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))## Figures 4.17 and 4.18 in Gencay, Selcuk and Whitcher (2001). data(ibm) ibm.returns <- diff(log(ibm)) ## Haar ibmr.haar <- dwt(ibm.returns, "haar") names(ibmr.haar) <- c("w1", "w2", "w3", "w4", "v4") ## plot partial Haar DWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(a)") for(i in 1:4) plot.ts(up.sample(ibmr.haar[[i]], 2^i), type="h", axes=FALSE, ylab=names(ibmr.haar)[i]) plot.ts(up.sample(ibmr.haar$v4, 2^4), type="h", axes=FALSE, ylab=names(ibmr.haar)[5]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) ## LA(8) ibmr.la8 <- dwt(ibm.returns, "la8") names(ibmr.la8) <- c("w1", "w2", "w3", "w4", "v4") ## must shift LA(8) coefficients ibmr.la8$w1 <- c(ibmr.la8$w1[-c(1:2)], ibmr.la8$w1[1:2]) ibmr.la8$w2 <- c(ibmr.la8$w2[-c(1:2)], ibmr.la8$w2[1:2]) for(i in names(ibmr.la8)[3:4]) ibmr.la8[[i]] <- c(ibmr.la8[[i]][-c(1:3)], ibmr.la8[[i]][1:3]) ibmr.la8$v4 <- c(ibmr.la8$v4[-c(1:2)], ibmr.la8$v4[1:2]) ## plot partial LA(8) DWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(b)") for(i in 1:4) plot.ts(up.sample(ibmr.la8[[i]], 2^i), type="h", axes=FALSE, ylab=names(ibmr.la8)[i]) plot.ts(up.sample(ibmr.la8$v4, 2^4), type="h", axes=FALSE, ylab=names(ibmr.la8)[5]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))
Performs a separable two-dimensional discrete wavelet transform (DWT) on a matrix of dyadic dimensions.
dwt.2d(x, wf, J = 4, boundary = "periodic") idwt.2d(y)dwt.2d(x, wf, J = 4, boundary = "periodic") idwt.2d(y)
x |
input matrix (image) |
wf |
name of the wavelet filter to use in the decomposition |
J |
depth of the decomposition, must be a number less than or equal to log(minM,N,2) |
boundary |
only |
y |
an object of class |
See references.
List structure containing the sub-matrices from the
decomposition.
B. Whitcher
Mallat, S. (1998) A Wavelet Tour of Signal Processing, Academic Press.
Vetterli, M. and J. Kovacevic (1995) Wavelets and Subband Coding, Prentice Hall.
## Xbox image data(xbox) xbox.dwt <- dwt.2d(xbox, "haar", 3) par(mfrow=c(1,1), pty="s") plot.dwt.2d(xbox.dwt) par(mfrow=c(2,2), pty="s") image(1:dim(xbox)[1], 1:dim(xbox)[2], xbox, xlab="", ylab="", main="Original Image") image(1:dim(xbox)[1], 1:dim(xbox)[2], idwt.2d(xbox.dwt), xlab="", ylab="", main="Wavelet Reconstruction") image(1:dim(xbox)[1], 1:dim(xbox)[2], xbox - idwt.2d(xbox.dwt), xlab="", ylab="", main="Difference") ## Daubechies image data(dau) par(mfrow=c(1,1), pty="s") image(dau, col=rainbow(128)) sum(dau^2) dau.dwt <- dwt.2d(dau, "d4", 3) plot.dwt.2d(dau.dwt) sum(plot.dwt.2d(dau.dwt, plot=FALSE)^2)## Xbox image data(xbox) xbox.dwt <- dwt.2d(xbox, "haar", 3) par(mfrow=c(1,1), pty="s") plot.dwt.2d(xbox.dwt) par(mfrow=c(2,2), pty="s") image(1:dim(xbox)[1], 1:dim(xbox)[2], xbox, xlab="", ylab="", main="Original Image") image(1:dim(xbox)[1], 1:dim(xbox)[2], idwt.2d(xbox.dwt), xlab="", ylab="", main="Wavelet Reconstruction") image(1:dim(xbox)[1], 1:dim(xbox)[2], xbox - idwt.2d(xbox.dwt), xlab="", ylab="", main="Difference") ## Daubechies image data(dau) par(mfrow=c(1,1), pty="s") image(dau, col=rainbow(128)) sum(dau^2) dau.dwt <- dwt.2d(dau, "d4", 3) plot.dwt.2d(dau.dwt) sum(plot.dwt.2d(dau.dwt, plot=FALSE)^2)
Three-dimensional separable discrete wavelet transform (DWT).
dwt.3d(x, wf, J = 4, boundary = "periodic") idwt.3d(y)dwt.3d(x, wf, J = 4, boundary = "periodic") idwt.3d(y)
x |
input array |
wf |
name of the wavelet filter to use in the decomposition |
J |
depth of the decomposition, must be a number less than or equal to log(minZ,Y,Z,2) |
boundary |
only |
y |
an object of class |
B. Whitcher
The discrete Hilbert wavelet transforms (DHWTs) for seasonal and time-varying time series analysis. Transforms include the usual orthogonal (decimated), maximal-overlap (non-decimated) and maximal-overlap packet transforms.
dwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...) dwt.hilbert.nondyadic(x, ...) idwt.hilbert(y) modwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...) imodwt.hilbert(y) modwpt.hilbert(x, wf, n.levels = 4, boundary = "periodic")dwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...) dwt.hilbert.nondyadic(x, ...) idwt.hilbert(y) modwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...) imodwt.hilbert(y) modwpt.hilbert(x, wf, n.levels = 4, boundary = "periodic")
x |
Real-valued time series or vector of observations. |
wf |
Hilbert wavelet pair |
n.levels |
Number of levels (depth) of the wavelet transform. |
boundary |
Boundary treatment, currently only |
... |
Additional parametes to be passed on. |
y |
An object of S3 class |
Hilbert wavelet transform object (list).
B. Whitcher
Selesnick, I. (200X). IEEE Signal Processing Magazine
Selesnick, I. (200X). IEEE Transactions in Signal Processing
Whither, B. and P.F. Craigmile (2004). Multivariate Spectral Analysis Using Hilbert Wavelet Pairs, International Journal of Wavelets, Multiresolution and Information Processing, 2(4), 567–587.
Monthly foreign exchange rates for the Deutsche Mark - U.S. Dollar (DEM-USD) and Japanese Yen - U.S. Dollar (JPY-USD) starting in 1970.
data(exchange)data(exchange)
A bivariate time series containing 348 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Parameter estimation for a fractional difference (long-memory, self-similar) process is performed via maximum likelihood on the wavelet coefficients.
fdp.mle(y, wf, J = log(length(y), 2))fdp.mle(y, wf, J = log(length(y), 2))
y |
Dyadic length time series. |
wf |
Name of the wavelet filter to use in the decomposition. See
|
J |
Depth of the discrete wavelet transform. |
The variance-covariance matrix of the original time series is approximated
by its wavelet-based equivalent. A Whittle-type likelihood is then
constructed where the sums of squared wavelet coefficients are compared to
bandpass filtered version of the true spectrum. Minimization occurs only
for the fractional difference parameter , while variance is estimated
afterwards.
List containing the maximum likelihood estimates (MLEs) of
and , along with the value of the likelihood for those
estimates.
B. Whitcher
M. J. Jensen (2000) An alternative maximum likelihood estimator of long-memory processes using compactly supported wavelets, Journal of Economic Dynamics and Control, 24, No. 3, 361-387.
McCoy, E. J., and A. T. Walden (1996) Wavelet analysis and synthesis of stationary long-memory processes, Journal for Computational and Graphical Statistics, 5, No. 1, 26-56.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
## Figure 5.5 in Gencay, Selcuk and Whitcher (2001) fdp.sdf <- function(freq, d, sigma2=1) sigma2 / ((2*sin(pi * freq))^2)^d dB <- function(x) 10 * log10(x) per <- function(z) { n <- length(z) (Mod(fft(z))**2/(2*pi*n))[1:(n %/% 2 + 1)] } data(ibm) ibm.returns <- diff(log(ibm)) ibm.volatility <- abs(ibm.returns) ibm.vol.mle <- fdp.mle(ibm.volatility, "d4", 4) freq <- 0:184/368 ibm.vol.per <- 2 * pi * per(ibm.volatility) ibm.vol.resid <- ibm.vol.per/ fdp.sdf(freq, ibm.vol.mle$parameters[1]) par(mfrow=c(1,1), las=0, pty="m") plot(freq, dB(ibm.vol.per), type="l", xlab="Frequency", ylab="Spectrum") lines(freq, dB(fdp.sdf(freq, ibm.vol.mle$parameters[1], ibm.vol.mle$parameters[2]/2)), col=2)## Figure 5.5 in Gencay, Selcuk and Whitcher (2001) fdp.sdf <- function(freq, d, sigma2=1) sigma2 / ((2*sin(pi * freq))^2)^d dB <- function(x) 10 * log10(x) per <- function(z) { n <- length(z) (Mod(fft(z))**2/(2*pi*n))[1:(n %/% 2 + 1)] } data(ibm) ibm.returns <- diff(log(ibm)) ibm.volatility <- abs(ibm.returns) ibm.vol.mle <- fdp.mle(ibm.volatility, "d4", 4) freq <- 0:184/368 ibm.vol.per <- 2 * pi * per(ibm.volatility) ibm.vol.resid <- ibm.vol.per/ fdp.sdf(freq, ibm.vol.mle$parameters[1]) par(mfrow=c(1,1), las=0, pty="m") plot(freq, dB(ibm.vol.per), type="l", xlab="Frequency", ylab="Spectrum") lines(freq, dB(fdp.sdf(freq, ibm.vol.mle$parameters[1], ibm.vol.mle$parameters[2]/2)), col=2)
Draws the spectral density functions (SDFs) for standard long-memory processes including fractional difference (FD), seasonal persistent (SP), and seasonal fractional difference (SFD) processes.
fdp.sdf(freq, d, sigma2 = 1) spp.sdf(freq, d, fG, sigma2 = 1) spp2.sdf(freq, d1, f1, d2, f2, sigma2 = 1) sfd.sdf(freq, s, d, sigma2 = 1)fdp.sdf(freq, d, sigma2 = 1) spp.sdf(freq, d, fG, sigma2 = 1) spp2.sdf(freq, d1, f1, d2, f2, sigma2 = 1) sfd.sdf(freq, s, d, sigma2 = 1)
freq |
vector of frequencies, normally from 0 to 0.5 |
d, d1, d2
|
fractional difference parameter |
sigma2 |
innovations variance |
fG, f1, f2
|
Gegenbauer frequency |
s |
seasonal parameter |
The power spectrum from an FD, SP or SFD process.
B. Whitcher
dB <- function(x) 10 * log10(x) fdp.main <- expression(paste("FD", group("(",d==0.4,")"))) sfd.main <- expression(paste("SFD", group("(",list(s==12, d==0.4),")"))) spp.main <- expression(paste("SPP", group("(",list(delta==0.4, f[G]==1/12),")"))) freq <- 0:512/1024 par(mfrow=c(2,2), mar=c(5-1,4,4-1,2), col.main="darkred") plot(freq, dB(fdp.sdf(freq, .4)), type="l", xlab="frequency", ylab="spectrum (dB)", main=fdp.main) plot(freq, dB(spp.sdf(freq, .4, 1/12)), type="l", xlab="frequency", ylab="spectrum (dB)", font.main=1, main=spp.main) plot(freq, dB(sfd.sdf(freq, 12, .4)), type="l", xlab="frequency", ylab="spectrum (dB)", main=sfd.main)dB <- function(x) 10 * log10(x) fdp.main <- expression(paste("FD", group("(",d==0.4,")"))) sfd.main <- expression(paste("SFD", group("(",list(s==12, d==0.4),")"))) spp.main <- expression(paste("SPP", group("(",list(delta==0.4, f[G]==1/12),")"))) freq <- 0:512/1024 par(mfrow=c(2,2), mar=c(5-1,4,4-1,2), col.main="darkred") plot(freq, dB(fdp.sdf(freq, .4)), type="l", xlab="frequency", ylab="spectrum (dB)", main=fdp.main) plot(freq, dB(spp.sdf(freq, .4, 1/12)), type="l", xlab="frequency", ylab="spectrum (dB)", font.main=1, main=spp.main) plot(freq, dB(sfd.sdf(freq, 12, .4)), type="l", xlab="frequency", ylab="spectrum (dB)", main=sfd.main)
Subroutine for use in simulating seasonal persistent processes using the discrete wavelet packet transform.
find.adaptive.basis(wf, J, fG, eps)find.adaptive.basis(wf, J, fG, eps)
wf |
Character string; name of the wavelet filter. |
J |
Depth of the discrete wavelet packet transform. |
fG |
Gegenbauer frequency. |
eps |
Threshold for the squared gain function. |
The squared gain functions for a Daubechies (extremal phase or least asymmetric) wavelet family are used in a filter cascade to compute the value of the squared gain function for the wavelet packet filter at the Gengenbauer frequency. This is done for all nodes of the wavelet packet table.
The idea behind this subroutine is to approximate the relationship between the discrete wavelet transform and long-memory processes, where the squared gain function is zero at frequency zero for all levels of the DWT.
Boolean vector describing the orthonormal basis for the DWPT.
B. Whitcher
Used in dwpt.sim.
Farras nearly symmetric filters for orthogonal 2-channel perfect reconstruction filter bank and Farras filters organized for the dual-tree complex DWT.
FSfarras()FSfarras()
af |
List (i=1,2) - analysis filters for tree i |
sf |
List (i=1,2) - synthesis filters for tree i |
Matlab: S. Cai, K. Li and I. Selesnick; R port: B. Whitcher
A. F. Abdelnour and I. W. Selesnick. “Nearly symmetric orthogonal wavelet bases”, Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing (ICASSP), May 2001.
data(heavisine)data(heavisine)
A vector containing 512 observations.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
Converts name of Hilbert wavelet pair to filter coefficients.
hilbert.filter(name)hilbert.filter(name)
name |
Character string of Hilbert wavelet pair, see acceptable names
below (e.g., |
Simple switch statement selects the appropriate HWP. There are two
parameters that define a Hilbert wavelet pair using the notation of
Selesnick (2001,2002), and . Currently, the only implemented
combinations are (3,3), (3,5), (4,2) and (4,4).
List containing the following items:
L |
length of the wavelet filter |
h0, g0
|
low-pass filter coefficients |
h1, g1
|
high-pass filter coefficients |
B. Whitcher
Selesnick, I.W. (2001). Hilbert transform pairs of wavelet bases. IEEE Signal Processing Letters 8(6), 170–173.
Selesnick, I.W. (2002). The design of approximate Hilbert transform pairs of wavelet bases. IEEE Transactions on Signal Processing 50(5), 1144–1152.
hilbert.filter("k3l3") hilbert.filter("k3l5") hilbert.filter("k4l2") hilbert.filter("k4l4")hilbert.filter("k3l3") hilbert.filter("k3l5") hilbert.filter("k4l2") hilbert.filter("k4l4")
Uses exact time-domain method from Hosking (1984) to generate a simulated time series from a specified autocovariance sequence.
hosking.sim(n, acvs)hosking.sim(n, acvs)
n |
Length of series. |
acvs |
Autocovariance sequence of series with which to generate, must
be of length at least |
Length n time series from true autocovariance sequence
acvs.
B. Whitcher
Hosking, J. R. M. (1984) Modeling persistence in hydrological time series using fractional differencing, Water Resources Research, 20, No. 12, 1898-1908.
Percival, D. B. (1992) Simulating Gaussian random processes with specified spectra, Computing Science and Statistics, 22, 534-538.
dB <- function(x) 10 * log10(x) per <- function (z) { n <- length(z) (Mod(fft(z))^2/(2 * pi * n))[1:(n%/%2 + 1)] } spp.sdf <- function(freq, delta, omega) abs(2 * (cos(2*pi*freq) - cos(2*pi*omega)))^(-2*delta) data(acvs.andel8) n <- 1024 ## Not run: z <- hosking.sim(n, acvs.andel8[,2]) per.z <- 2 * pi * per(z) par(mfrow=c(2,1), las=1) plot.ts(z, ylab="", main="Realization of a Seasonal Long-Memory Process") plot(0:(n/2)/n, dB(per.z), type="l", xlab="Frequency", ylab="dB", main="Periodogram") lines(0:(n/2)/n, dB(spp.sdf(0:(n/2)/n, .4, 1/12)), col=2) ## End(Not run)dB <- function(x) 10 * log10(x) per <- function (z) { n <- length(z) (Mod(fft(z))^2/(2 * pi * n))[1:(n%/%2 + 1)] } spp.sdf <- function(freq, delta, omega) abs(2 * (cos(2*pi*freq) - cos(2*pi*omega)))^(-2*delta) data(acvs.andel8) n <- 1024 ## Not run: z <- hosking.sim(n, acvs.andel8[,2]) per.z <- 2 * pi * per(z) par(mfrow=c(2,1), las=1) plot.ts(z, ylab="", main="Realization of a Seasonal Long-Memory Process") plot(0:(n/2)/n, dB(per.z), type="l", xlab="Frequency", ylab="dB", main="Periodogram") lines(0:(n/2)/n, dB(spp.sdf(0:(n/2)/n, .4, 1/12)), col=2) ## End(Not run)
Daily IBM stock prices spanning May~17, 1961 to November~2, 1962.
data(ibm)data(ibm)
A vector containing 369 observations.
Box, G. E.P. and Jenkins, G. M. (1976) Time Series Analysis: Forecasting and Control, Holden Day, San Francisco, 2nd edition.
Quarterly Japanese gross national product from 1955:1 to 1996:4.
data(japan)data(japan)
A vector containing 169 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Hecq, A. (1998) Does seasonal adjustment induce common cycles?, Empirical Economics, 59, 289-297.
data(jumpsine)data(jumpsine)
A vector containing 512 observations.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
Seismograph (vertical acceleration, nm/sq.sec) of the Kobe earthquake, recorded at Tasmania University, HobarTRUE, Australia on 16 January 1995 beginning at 20:56:51 (GMTRUE) and continuing for 51 minutes at 1 second intervals.
data(kobe)data(kobe)
A vector containing 3048 observations.
Data management centre, Washington University.
data(linchirp)data(linchirp)
A vector containing 512 observations.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.
Perform wavelet shrinkage using data-analytic, hybrid SURE, manual, SURE, or universal thresholding.
da.thresh(wc, alpha = .05, max.level = 4, verbose = FALSE, return.thresh = FALSE) hybrid.thresh(wc, max.level = 4, verbose = FALSE, seed = 0) manual.thresh(wc, max.level = 4, value, hard = TRUE) sure.thresh(wc, max.level = 4, hard = TRUE) universal.thresh(wc, max.level = 4, hard = TRUE) universal.thresh.modwt(wc, max.level = 4, hard = TRUE)da.thresh(wc, alpha = .05, max.level = 4, verbose = FALSE, return.thresh = FALSE) hybrid.thresh(wc, max.level = 4, verbose = FALSE, seed = 0) manual.thresh(wc, max.level = 4, value, hard = TRUE) sure.thresh(wc, max.level = 4, hard = TRUE) universal.thresh(wc, max.level = 4, hard = TRUE) universal.thresh.modwt(wc, max.level = 4, hard = TRUE)
wc |
wavelet coefficients |
max.level |
maximum level of coefficients to be affected by threshold |
value |
threshold value (only utilized in |
hard |
Boolean value, if |
alpha |
level of the hypothesis tests |
verbose |
if |
seed |
sets random seed (only utilized in |
return.thresh |
if |
An extensive amount of literature has been written on wavelet shrinkage. The functions here represent the most basic approaches to the problem of nonparametric function estimation. See the references for further information.
The default output is a list structure, the same length as was input, containing only those wavelet coefficients surviving the threshold.
B. Whitcher (some code taken from R. Todd Ogden)
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Ogden, R. T. (1996) Essential Wavelets for Statistical Applications and Data Analysis, Birkhauser.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Vidakovic, B. (1999) Statistical Modeling by Wavelets, John Wiley and Sons.
Percentage changes in monthly Mexican money supply.
data(mexm)data(mexm)
A vector containing 516 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Performs time-varying or seasonal coherence and phase anlaysis between two time seris using the maximal-overlap discrete Hilbert wavelet transform (MODHWT).
modhwt.coh(x, y, f.length = 0) modhwt.phase(x, y, f.length = 0) modhwt.coh.seasonal(x, y, S = 10, season = 365) modhwt.phase.seasonal(x, y, season = 365)modhwt.coh(x, y, f.length = 0) modhwt.phase(x, y, f.length = 0) modhwt.coh.seasonal(x, y, S = 10, season = 365) modhwt.phase.seasonal(x, y, season = 365)
x |
MODHWT object. |
y |
MODHWT object. |
f.length |
Length of the rectangular filter. |
S |
Number of "seasons". |
season |
Length of the "season". |
The idea of seasonally-varying spectral analysis (SVSA, Madden 1986) is
generalized using the MODWT and Hilbert wavelet pairs. For the seasonal
case, seasons are used to produce a consistent estimate of the
coherence and phase. For the non-seasonal case, a simple rectangular
(moving-average) filter is applied to the MODHWT coefficients in order to
produce consistent estimates.
Time-varying or seasonal coherence and phase between two time
series. The coherence estimates are between zero and one, while the phase
estimates are between and .
B. Whitcher
Madden, R.A. (1986). Seasonal variation of the 40–50 day oscillation in the tropics. Journal of the Atmospheric Sciences 43(24), 3138–3158.
Whither, B. and P.F. Craigmile (2004). Multivariate Spectral Analysis Using Hilbert Wavelet Pairs, International Journal of Wavelets, Multiresolution and Information Processing, 2(4), 567–587.
This function performs a level decomposition of the input vector
using the non-decimated discrete wavelet transform. The inverse transform
performs the reconstruction of a vector or time series from its maximal
overlap discrete wavelet transform.
modwt(x, wf = "la8", n.levels = 4, boundary = "periodic") imodwt(y)modwt(x, wf = "la8", n.levels = 4, boundary = "periodic") imodwt(y)
x |
a vector or time series containing the data be to decomposed. There is no restriction on its length. |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
n.levels |
Specifies the depth of the decomposition. This must be a number less than or equal to log(length(x),2). |
boundary |
Character string specifying the boundary condition. If
|
y |
an object of class |
The code implements the one-dimensional non-decimated DWT using the pyramid algorithm. The actual transform is performed in C using pseudocode from Percival and Walden (2001). That means convolutions, not inner products, are used to apply the wavelet filters.
The MODWT goes by several names in the statistical and engineering literature, such as, the “stationary DWT”, “translation-invariant DWT”, and “time-invariant DWT”.
The inverse MODWT implements the one-dimensional inverse transform using the pyramid algorithm (Mallat, 1989).
Basically, a list with the following components
d? |
Wavelet coefficient vectors. |
s? |
Scaling coefficient vector. |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Percival, D. B. and P. Guttorp (1994) Long-memory processes, the Allan variance and wavelets, In Wavelets and Geophysics, pages 325-344, Academic Press.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
## Figure 4.23 in Gencay, Selcuk and Whitcher (2001) data(ibm) ibm.returns <- diff(log(ibm)) # Haar ibmr.haar <- modwt(ibm.returns, "haar") names(ibmr.haar) <- c("w1", "w2", "w3", "w4", "v4") # LA(8) ibmr.la8 <- modwt(ibm.returns, "la8") names(ibmr.la8) <- c("w1", "w2", "w3", "w4", "v4") # shift the MODWT vectors ibmr.la8 <- phase.shift(ibmr.la8, "la8") ## plot partial MODWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(a)") for(i in 1:5) plot.ts(ibmr.haar[[i]], axes=FALSE, ylab=names(ibmr.haar)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(b)") for(i in 1:5) plot.ts(ibmr.la8[[i]], axes=FALSE, ylab=names(ibmr.la8)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))## Figure 4.23 in Gencay, Selcuk and Whitcher (2001) data(ibm) ibm.returns <- diff(log(ibm)) # Haar ibmr.haar <- modwt(ibm.returns, "haar") names(ibmr.haar) <- c("w1", "w2", "w3", "w4", "v4") # LA(8) ibmr.la8 <- modwt(ibm.returns, "la8") names(ibmr.la8) <- c("w1", "w2", "w3", "w4", "v4") # shift the MODWT vectors ibmr.la8 <- phase.shift(ibmr.la8, "la8") ## plot partial MODWT for IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(a)") for(i in 1:5) plot.ts(ibmr.haar[[i]], axes=FALSE, ylab=names(ibmr.haar)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.returns, axes=FALSE, ylab="", main="(b)") for(i in 1:5) plot.ts(ibmr.la8[[i]], axes=FALSE, ylab=names(ibmr.la8)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))
Performs a separable two-dimensional maximal overlap discrete wavelet transform (MODWT) on a matrix of arbitrary dimensions.
modwt.2d(x, wf, J = 4, boundary = "periodic") imodwt.2d(y)modwt.2d(x, wf, J = 4, boundary = "periodic") imodwt.2d(y)
x |
input matrix |
wf |
name of the wavelet filter to use in the decomposition |
J |
depth of the decomposition |
boundary |
only |
y |
an object of class |
See references.
List structure containing the sub-matrices from the
decomposition.
B. Whitcher
Liang, J. and T. W. Parks (1994) A two-dimensional translation invariant wavelet representation and its applications, Proceedings ICIP-94, Vol. 1, 66-70.
Liang, J. and T. W. Parks (1994) Image coding using translation invariant wavelet transforms with symmetric extensions, IEEE Transactions on Image Processing, 7, No. 5, 762-769.
## Xbox image data(xbox) xbox.modwt <- modwt.2d(xbox, "haar", 2) ## Level 1 decomposition par(mfrow=c(2,2), pty="s") image(xbox.modwt$LH1, col=rainbow(128), axes=FALSE, main="LH1") image(xbox.modwt$HH1, col=rainbow(128), axes=FALSE, main="HH1") frame() image(xbox.modwt$HL1, col=rainbow(128), axes=FALSE, main="HL1") ## Level 2 decomposition par(mfrow=c(2,2), pty="s") image(xbox.modwt$LH2, col=rainbow(128), axes=FALSE, main="LH2") image(xbox.modwt$HH2, col=rainbow(128), axes=FALSE, main="HH2") image(xbox.modwt$LL2, col=rainbow(128), axes=FALSE, main="LL2") image(xbox.modwt$HL2, col=rainbow(128), axes=FALSE, main="HL2") sum((xbox - imodwt.2d(xbox.modwt))^2) data(dau) par(mfrow=c(1,1), pty="s") image(dau, col=rainbow(128), axes=FALSE, main="Ingrid Daubechies") sum(dau^2) dau.modwt <- modwt.2d(dau, "d4", 2) ## Level 1 decomposition par(mfrow=c(2,2), pty="s") image(dau.modwt$LH1, col=rainbow(128), axes=FALSE, main="LH1") image(dau.modwt$HH1, col=rainbow(128), axes=FALSE, main="HH1") frame() image(dau.modwt$HL1, col=rainbow(128), axes=FALSE, main="HL1") ## Level 2 decomposition par(mfrow=c(2,2), pty="s") image(dau.modwt$LH2, col=rainbow(128), axes=FALSE, main="LH2") image(dau.modwt$HH2, col=rainbow(128), axes=FALSE, main="HH2") image(dau.modwt$LL2, col=rainbow(128), axes=FALSE, main="LL2") image(dau.modwt$HL2, col=rainbow(128), axes=FALSE, main="HL2") sum((dau - imodwt.2d(dau.modwt))^2)## Xbox image data(xbox) xbox.modwt <- modwt.2d(xbox, "haar", 2) ## Level 1 decomposition par(mfrow=c(2,2), pty="s") image(xbox.modwt$LH1, col=rainbow(128), axes=FALSE, main="LH1") image(xbox.modwt$HH1, col=rainbow(128), axes=FALSE, main="HH1") frame() image(xbox.modwt$HL1, col=rainbow(128), axes=FALSE, main="HL1") ## Level 2 decomposition par(mfrow=c(2,2), pty="s") image(xbox.modwt$LH2, col=rainbow(128), axes=FALSE, main="LH2") image(xbox.modwt$HH2, col=rainbow(128), axes=FALSE, main="HH2") image(xbox.modwt$LL2, col=rainbow(128), axes=FALSE, main="LL2") image(xbox.modwt$HL2, col=rainbow(128), axes=FALSE, main="HL2") sum((xbox - imodwt.2d(xbox.modwt))^2) data(dau) par(mfrow=c(1,1), pty="s") image(dau, col=rainbow(128), axes=FALSE, main="Ingrid Daubechies") sum(dau^2) dau.modwt <- modwt.2d(dau, "d4", 2) ## Level 1 decomposition par(mfrow=c(2,2), pty="s") image(dau.modwt$LH1, col=rainbow(128), axes=FALSE, main="LH1") image(dau.modwt$HH1, col=rainbow(128), axes=FALSE, main="HH1") frame() image(dau.modwt$HL1, col=rainbow(128), axes=FALSE, main="HL1") ## Level 2 decomposition par(mfrow=c(2,2), pty="s") image(dau.modwt$LH2, col=rainbow(128), axes=FALSE, main="LH2") image(dau.modwt$HH2, col=rainbow(128), axes=FALSE, main="HH2") image(dau.modwt$LL2, col=rainbow(128), axes=FALSE, main="LL2") image(dau.modwt$HL2, col=rainbow(128), axes=FALSE, main="HL2") sum((dau - imodwt.2d(dau.modwt))^2)
Three-dimensional separable maximal overlap discrete wavelet transform (MODWT).
modwt.3d(x, wf, J = 4, boundary = "periodic") imodwt.3d(y)modwt.3d(x, wf, J = 4, boundary = "periodic") imodwt.3d(y)
x |
input array |
wf |
name of the wavelet filter to use in the decomposition |
J |
depth of the decomposition |
boundary |
only |
y |
an object of class |
B. Whitcher
This function performs a level additive decomposition of the input
vector or time series using the pyramid algorithm (Mallat 1989).
mra(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")mra(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")
x |
A vector or time series containing the data be to decomposed. This
must be a dyadic length vector (power of 2) for |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
J |
Specifies the depth of the decomposition. This must be a number less than or equal to log(length(x), 2). |
method |
Either |
boundary |
Character string specifying the boundary condition. If
|
This code implements a one-dimensional multiresolution analysis introduced by Mallat (1989). Either the DWT or MODWT may be used to compute the multiresolution analysis, which is an additive decomposition of the original time series.
Basically, a list with the following components
D? |
Wavelet detail vectors. |
S? |
Wavelet smooth vector. |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, No. 7, 674-693.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
## Easy check to see if it works... x <- rnorm(32) x.mra <- mra(x) sum(x - apply(matrix(unlist(x.mra), nrow=32), 1, sum))^2 ## Figure 4.19 in Gencay, Selcuk and Whitcher (2001) data(ibm) ibm.returns <- diff(log(ibm)) ibm.volatility <- abs(ibm.returns) ## Haar ibmv.haar <- mra(ibm.volatility, "haar", 4, "dwt") names(ibmv.haar) <- c("d1", "d2", "d3", "d4", "s4") ## LA(8) ibmv.la8 <- mra(ibm.volatility, "la8", 4, "dwt") names(ibmv.la8) <- c("d1", "d2", "d3", "d4", "s4") ## plot multiresolution analysis of IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.volatility, axes=FALSE, ylab="", main="(a)") for(i in 1:5) plot.ts(ibmv.haar[[i]], axes=FALSE, ylab=names(ibmv.haar)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.volatility, axes=FALSE, ylab="", main="(b)") for(i in 1:5) plot.ts(ibmv.la8[[i]], axes=FALSE, ylab=names(ibmv.la8)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))## Easy check to see if it works... x <- rnorm(32) x.mra <- mra(x) sum(x - apply(matrix(unlist(x.mra), nrow=32), 1, sum))^2 ## Figure 4.19 in Gencay, Selcuk and Whitcher (2001) data(ibm) ibm.returns <- diff(log(ibm)) ibm.volatility <- abs(ibm.returns) ## Haar ibmv.haar <- mra(ibm.volatility, "haar", 4, "dwt") names(ibmv.haar) <- c("d1", "d2", "d3", "d4", "s4") ## LA(8) ibmv.la8 <- mra(ibm.volatility, "la8", 4, "dwt") names(ibmv.la8) <- c("d1", "d2", "d3", "d4", "s4") ## plot multiresolution analysis of IBM data par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.volatility, axes=FALSE, ylab="", main="(a)") for(i in 1:5) plot.ts(ibmv.haar[[i]], axes=FALSE, ylab=names(ibmv.haar)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368)) par(mfcol=c(6,1), pty="m", mar=c(5-2,4,4-2,2)) plot.ts(ibm.volatility, axes=FALSE, ylab="", main="(b)") for(i in 1:5) plot.ts(ibmv.la8[[i]], axes=FALSE, ylab=names(ibmv.la8)[i]) axis(side=1, at=seq(0,368,by=23), labels=c(0,"",46,"",92,"",138,"",184,"",230,"",276,"",322,"",368))
This function performs a level additive decomposition of the input
matrix or image using the pyramid algorithm (Mallat 1989).
mra.2d(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")mra.2d(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")
x |
A matrix or image containing the data be to decomposed. This must
be have dyadic length in both dimensions (but not necessarily the same) for
|
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
J |
Specifies the depth of the decomposition. This must be a number less than or equal to log(length(x),2). |
method |
Either |
boundary |
Character string specifying the boundary condition. If
|
This code implements a two-dimensional multiresolution analysis by performing the one-dimensional pyramid algorithm (Mallat 1989) on the rows and columns of the input matrix. Either the DWT or MODWT may be used to compute the multiresolution analysis, which is an additive decomposition of the original matrix (image).
Basically, a list with the following components
LH? |
Wavelet detail image in the horizontal direction. |
HL? |
Wavelet detail image in the vertical direction. |
HH? |
Wavelet detail image in the diagonal direction. |
LLJ |
Wavelet smooth image at the coarsest resolution. |
J |
Depth of the wavelet transform. |
wavelet |
Name of the wavelet filter used. |
boundary |
How the boundaries were handled. |
B. Whitcher
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, No. 7, 674-693.
Mallat, S. G. (1998) A Wavelet Tour of Signal Processing, Academic Press.
## Easy check to see if it works... ## -------------------------------- x <- matrix(rnorm(32*32), 32, 32) # MODWT x.mra <- mra.2d(x, method="modwt") x.mra.sum <- x.mra[[1]] for(j in 2:length(x.mra)) x.mra.sum <- x.mra.sum + x.mra[[j]] sum((x - x.mra.sum)^2) # DWT x.mra <- mra.2d(x, method="dwt") x.mra.sum <- x.mra[[1]] for(j in 2:length(x.mra)) x.mra.sum <- x.mra.sum + x.mra[[j]] sum((x - x.mra.sum)^2)## Easy check to see if it works... ## -------------------------------- x <- matrix(rnorm(32*32), 32, 32) # MODWT x.mra <- mra.2d(x, method="modwt") x.mra.sum <- x.mra[[1]] for(j in 2:length(x.mra)) x.mra.sum <- x.mra.sum + x.mra[[j]] sum((x - x.mra.sum)^2) # DWT x.mra <- mra.2d(x, method="dwt") x.mra.sum <- x.mra[[1]] for(j in 2:length(x.mra)) x.mra.sum <- x.mra.sum + x.mra[[j]] sum((x - x.mra.sum)^2)
This function performs a level additive decomposition of the input
array using the pyramid algorithm (Mallat 1989).
mra.3d(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")mra.3d(x, wf = "la8", J = 4, method = "modwt", boundary = "periodic")
x |
A three-dimensional array containing the data be to decomposed.
This must be have dyadic length in all three dimensions (but not necessarily
the same) for |
wf |
Name of the wavelet filter to use in the decomposition. By
default this is set to |
J |
Specifies the depth of the decomposition. This must be a number
less than or equal to |
method |
Either |
boundary |
Character string specifying the boundary condition. If
|
This code implements a three-dimensional multiresolution analysis by performing the one-dimensional pyramid algorithm (Mallat 1989) on each dimension of the input array. Either the DWT or MODWT may be used to compute the multiresolution analysis, which is an additive decomposition of the original array.
List structure containing the filter triplets associated with the multiresolution analysis.
B. Whitcher
Mallat, S. G. (1989) A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, No. 7, 674-693.
Mallat, S. G. (1998) A Wavelet Tour of Signal Processing, Academic Press.
This is the major subroutine for testing.hov, providing the
workhorse algorithm to recursively test and locate multiple variance changes
in so-called long memory processes.
mult.loc(dwt.list, modwt.list, wf, level, min.coef, debug)mult.loc(dwt.list, modwt.list, wf, level, min.coef, debug)
dwt.list |
List of wavelet vector coefficients from the |
modwt.list |
List of wavelet vector coefficients from the |
wf |
Name of the wavelet filter to use in the decomposition. |
level |
Specifies the depth of the decomposition. |
min.coef |
Minimum number of wavelet coefficients for testing purposes. |
debug |
Boolean variable: if set to |
For details see Section 9.6 of Percival and Walden (2000) or Section 7.3 in Gencay, Selcuk and Whitcher (2001).
Matrix.
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Computes the autocovariance function (ACF) for a time series or the cross-covariance function (CCF) between two time series.
my.acf(x) my.ccf(a, b)my.acf(x) my.ccf(a, b)
x, a, b
|
time series |
The series is zero padded to twice its length before the discrete Fourier transform is applied. Only the values corresponding to nonnegative lags are provided (for the ACF).
The autocovariance function for all nonnegative lags or the cross-covariance function for all lags.
B. Whitcher
data(ibm) ibm.returns <- diff(log(ibm)) plot(1:length(ibm.returns) - 1, my.acf(ibm.returns), type="h", xlab="lag", ylab="ACVS", main="Autocovariance Sequence for IBM Returns")data(ibm) ibm.returns <- diff(log(ibm)) plot(1:length(ibm.returns) - 1, my.acf(ibm.returns), type="h", xlab="lag", ylab="ACVS", main="Autocovariance Sequence for IBM Returns")
Yearly minimal water levels of the Nile river for the years 622 to 1281, measured at the Roda gauge near Cairo (Tousson, 1925, p. 366-385). The data are listed in chronological sequence by row.
data(nile)data(nile)
A length 663 vector.
The original Nile river data supplied by Beran only contained only 500 observations (622 to 1121). However, the book claimed to have 660 observations (622 to 1281). The remaining observations from the book were added, by hand, but the series still only contained 653 observations (622 to 1264).
Note, now the data consists of 663 observations (spanning the years 622-1284) as in original source (Toussoun, 1925).
Toussoun, O. (1925) M\'emoire sur l'Histoire du Nil, Volume 18 in M\'emoires a l'Institut d'Egypte, pp. 366-404.
Beran, J. (1994) Statistics for Long-Memory Processes, Chapman Hall: Englewood, NJ.
An orthonormal basis for the discrete wavelet transform may be characterized
via a disjoint partitioning of the frequency axis that covers
. This subroutine produces an orthonormal
basis from a full wavelet packet tree.
ortho.basis(xtree)ortho.basis(xtree)
xtree |
is a vector whose entries are associated with a wavelet packet tree. |
A wavelet packet tree is a binary tree of Boolean variables. Parent nodes are removed if any of their children exist.
Boolean vector describing the orthonormal basis for the DWPT.
B. Whitcher
data(japan) J <- 4 wf <- "mb8" japan.mra <- mra(log(japan), wf, J, boundary="reflection") japan.nomean <- ts(apply(matrix(unlist(japan.mra[-(J+1)]), ncol=J, byrow=FALSE), 1, sum), start=1955, freq=4) japan.nomean2 <- ts(japan.nomean[42:169], start=1965.25, freq=4) plot(japan.nomean2, type="l") japan.dwpt <- dwpt(japan.nomean2, wf, 6) japan.basis <- ortho.basis(portmanteau.test(japan.dwpt, p=0.01, type="other")) # Not implemented yet # par(mfrow=c(1,1)) # plot.basis(japan.basis)data(japan) J <- 4 wf <- "mb8" japan.mra <- mra(log(japan), wf, J, boundary="reflection") japan.nomean <- ts(apply(matrix(unlist(japan.mra[-(J+1)]), ncol=J, byrow=FALSE), 1, sum), start=1955, freq=4) japan.nomean2 <- ts(japan.nomean[42:169], start=1965.25, freq=4) plot(japan.nomean2, type="l") japan.dwpt <- dwpt(japan.nomean2, wf, 6) japan.basis <- ortho.basis(portmanteau.test(japan.dwpt, p=0.01, type="other")) # Not implemented yet # par(mfrow=c(1,1)) # plot.basis(japan.basis)
Computation of the periodogram via the Fast Fourier Transform (FFT).
per(z)per(z)
z |
time series |
Author: Jan Beran; modified: Martin Maechler, Date: Sep 1995.
Wavelet coefficients are circularly shifted by the amount of phase shift induced by the wavelet transform.
phase.shift(z, wf, inv = FALSE) phase.shift.packet(z, wf, inv = FALSE)phase.shift(z, wf, inv = FALSE) phase.shift.packet(z, wf, inv = FALSE)
z |
DWT object |
wf |
character string; wavelet filter used in DWT |
inv |
Boolean variable; if |
The center-of-energy argument of Hess-Nielsen and Wickerhauser (1996) is used to provide a flexible way to circularly shift wavelet coefficients regardless of the wavelet filter used. The results are not identical to those used by Percival and Walden (2000), but are more flexible.
phase.shift.packet is not yet implemented fully.
DWT (DWPT) object with coefficients circularly shifted.
B. Whitcher
Hess-Nielsen, N. and M. V. Wickerhauser (1996) Wavelets and time-frequency analysis, Proceedings of the IEEE, 84, No. 4, 523-540.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Wavelet coefficients are circularly shifted by the amount of phase shift induced by the discrete Hilbert wavelet transform.
phase.shift.hilbert(x, wf)phase.shift.hilbert(x, wf)
x |
Discete Hilbert wavelet transform (DHWT) object. |
wf |
character string; Hilbert wavelet pair used in DHWT |
The "center-of-energy" argument of Hess-Nielsen and Wickerhauser (1996) is used to provide a flexible way to circularly shift wavelet coefficients regardless of the wavelet filter used.
DHWT (DHWPT) object with coefficients circularly shifted.
B. Whitcher
Hess-Nielsen, N. and M. V. Wickerhauser (1996) Wavelets and time-frequency analysis, Proceedings of the IEEE, 84, No. 4, 523-540.
Organizes the wavelet coefficients from a 2D DWT into a single matrix and plots it. The coarser resolutions are nested within the lower-lefthand corner of the image.
## S3 method for class 'dwt.2d' plot(x, cex.axis = 1, plot = TRUE, ...)## S3 method for class 'dwt.2d' plot(x, cex.axis = 1, plot = TRUE, ...)
x |
input matrix (image) |
cex.axis |
|
plot |
if |
... |
additional graphical parameters if necessary |
The wavelet coefficients from the DWT object (a list) are reorganized into a single matrix of the same dimension as the original image and the result is plotted.
Image plot.
B. Whitcher
Computes the quadrature mirror filter from a given filter.
qmf(g, low2high = TRUE)qmf(g, low2high = TRUE)
g |
Filter coefficients. |
low2high |
Logical, default is |
None.
Quadrature mirror filter.
B. Whitcher
Any basic signal processing text.
## Haar wavelet filter g <- wave.filter("haar")$lpf qmf(g)## Haar wavelet filter g <- wave.filter("haar")$lpf qmf(g)
Provides the normalized cumulative sums of squares from a sequence of coefficients with the diagonal line removed.
rotcumvar(x)rotcumvar(x)
x |
vector of coefficients to be cumulatively summed (missing values excluded) |
The rotated cumulative variance, when plotted, provides a qualitative way to study the time dependence of the variance of a series. If the variance is stationary over time, then only small deviations from zero should be present. If on the other hand the variance is non-stationary, then large departures may exist. Formal hypothesis testing may be performed based on boundary crossings of Brownian bridge processes.
Vector of coefficients that are the sumulative sum of squared input coefficients.
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Compute phase shifts for wavelet sub-matrices based on the “center of energy” argument of Hess-Nielsen and Wickerhauser (1996).
shift.2d(z, inverse = FALSE)shift.2d(z, inverse = FALSE)
z |
Two-dimensional MODWT object |
inverse |
Boolean value on whether to perform the forward or inverse operation. |
The "center of energy" technique of Wickerhauser and Hess-Nielsen (1996) is employed to find circular shifts for the wavelet sub-matrices such that the coefficients are aligned with the original series. This corresponds to applying a (near) linear-phase filtering operation.
Two-dimensional MODWT object with circularly shifted coefficients.
B. Whitcher
Hess-Nielsen, N. and M. V. Wickerhauser (1996) Wavelets and time-frequency analysis, Proceedings of the IEEE, 84, No. 4, 523-540.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
n <- 512 G1 <- G2 <- dnorm(seq(-n/4, n/4, length=n)) G <- 100 * zapsmall(outer(G1, G2)) G <- modwt.2d(G, wf="la8", J=6) k <- 50 xr <- yr <- trunc(n/2) + (-k:k) par(mfrow=c(3,3), mar=c(1,1,2,1), pty="s") for (j in names(G)[1:9]) { image(G[[j]][xr,yr], col=rainbow(64), axes=FALSE, main=j) } Gs <- shift.2d(G) for (j in names(G)[1:9]) { image(Gs[[j]][xr,yr], col=rainbow(64), axes=FALSE, main=j) }n <- 512 G1 <- G2 <- dnorm(seq(-n/4, n/4, length=n)) G <- 100 * zapsmall(outer(G1, G2)) G <- modwt.2d(G, wf="la8", J=6) k <- 50 xr <- yr <- trunc(n/2) + (-k:k) par(mfrow=c(3,3), mar=c(1,1,2,1), pty="s") for (j in names(G)[1:9]) { image(G[[j]][xr,yr], col=rainbow(64), axes=FALSE, main=j) } Gs <- shift.2d(G) for (j in names(G)[1:9]) { image(Gs[[j]][xr,yr], col=rainbow(64), axes=FALSE, main=j) }
Computes sinusoidal data tapers directly from equations.
sine.taper(n, k)sine.taper(n, k)
n |
length of data taper(s) |
k |
number of data tapers |
See reference.
A vector or matrix of data tapers (cols = tapers).
B. Whitcher
Riedel, K. S. and A. Sidorenko (1995) Minimum bias multiple taper spectral estimation, IEEE Transactions on Signal Processing, 43, 188-195.
Computes wavelet cross-covariance or cross-correlation between two time series.
spin.covariance(x, y, lag.max = NA) spin.correlation(x, y, lag.max = NA)spin.covariance(x, y, lag.max = NA) spin.correlation(x, y, lag.max = NA)
x |
first time series |
y |
second time series, same length as |
lag.max |
maximum lag to compute cross-covariance (correlation) |
See references.
List structure holding the wavelet cross-covariances (correlations) according to scale.
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Whitcher, B., P. Guttorp and D. B. Percival (2000) Wavelet analysis of covariance with application to atmospheric time series, Journal of Geophysical Research, 105, No. D11, 14,941-14,962.
wave.covariance, wave.correlation.
## Figure 7.9 from Gencay, Selcuk and Whitcher (2001) data(exchange) returns <- diff(log(exchange)) returns <- ts(returns, start=1970, freq=12) wf <- "d4" demusd.modwt <- modwt(returns[,"DEM.USD"], wf, 8) demusd.modwt.bw <- brick.wall(demusd.modwt, wf) jpyusd.modwt <- modwt(returns[,"JPY.USD"], wf, 8) jpyusd.modwt.bw <- brick.wall(jpyusd.modwt, wf) n <- dim(returns)[1] J <- 6 lmax <- 36 returns.cross.cor <- NULL for(i in 1:J) { blah <- spin.correlation(demusd.modwt.bw[[i]], jpyusd.modwt.bw[[i]], lmax) returns.cross.cor <- cbind(returns.cross.cor, blah) } returns.cross.cor <- ts(as.matrix(returns.cross.cor), start=-36, freq=1) dimnames(returns.cross.cor) <- list(NULL, paste("Level", 1:J)) lags <- length(-lmax:lmax) lower.ci <- tanh(atanh(returns.cross.cor) - qnorm(0.975) / sqrt(matrix(trunc(n/2^(1:J)), nrow=lags, ncol=J, byrow=TRUE) - 3)) upper.ci <- tanh(atanh(returns.cross.cor) + qnorm(0.975) / sqrt(matrix(trunc(n/2^(1:J)), nrow=lags, ncol=J, byrow=TRUE) - 3)) par(mfrow=c(3,2), las=1, pty="m", mar=c(5,4,4,2)+.1) for(i in J:1) { plot(returns.cross.cor[,i], ylim=c(-1,1), xaxt="n", xlab="Lag (months)", ylab="", main=dimnames(returns.cross.cor)[[2]][i]) axis(side=1, at=seq(-36, 36, by=12)) lines(lower.ci[,i], lty=1, col=2) lines(upper.ci[,i], lty=1, col=2) abline(h=0,v=0) }## Figure 7.9 from Gencay, Selcuk and Whitcher (2001) data(exchange) returns <- diff(log(exchange)) returns <- ts(returns, start=1970, freq=12) wf <- "d4" demusd.modwt <- modwt(returns[,"DEM.USD"], wf, 8) demusd.modwt.bw <- brick.wall(demusd.modwt, wf) jpyusd.modwt <- modwt(returns[,"JPY.USD"], wf, 8) jpyusd.modwt.bw <- brick.wall(jpyusd.modwt, wf) n <- dim(returns)[1] J <- 6 lmax <- 36 returns.cross.cor <- NULL for(i in 1:J) { blah <- spin.correlation(demusd.modwt.bw[[i]], jpyusd.modwt.bw[[i]], lmax) returns.cross.cor <- cbind(returns.cross.cor, blah) } returns.cross.cor <- ts(as.matrix(returns.cross.cor), start=-36, freq=1) dimnames(returns.cross.cor) <- list(NULL, paste("Level", 1:J)) lags <- length(-lmax:lmax) lower.ci <- tanh(atanh(returns.cross.cor) - qnorm(0.975) / sqrt(matrix(trunc(n/2^(1:J)), nrow=lags, ncol=J, byrow=TRUE) - 3)) upper.ci <- tanh(atanh(returns.cross.cor) + qnorm(0.975) / sqrt(matrix(trunc(n/2^(1:J)), nrow=lags, ncol=J, byrow=TRUE) - 3)) par(mfrow=c(3,2), las=1, pty="m", mar=c(5,4,4,2)+.1) for(i in J:1) { plot(returns.cross.cor[,i], ylim=c(-1,1), xaxt="n", xlab="Lag (months)", ylab="", main=dimnames(returns.cross.cor)[[2]][i]) axis(side=1, at=seq(-36, 36, by=12)) lines(lower.ci[,i], lty=1, col=2) lines(upper.ci[,i], lty=1, col=2) abline(h=0,v=0) }
Parameter estimation for a seasonal persistent (seasonal long-memory) process is performed via maximum likelihood on the wavelet coefficients.
spp.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, frac = 1) spp2.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, dyadic = TRUE, frac = 1)spp.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, frac = 1) spp2.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, dyadic = TRUE, frac = 1)
y |
Not necessarily dyadic length time series. |
wf |
Name of the wavelet filter to use in the decomposition. See
|
J |
Depth of the discrete wavelet packet transform. |
p |
Level of significance for the white noise testing procedure. |
frac |
Fraction of the time series that should be used in constructing the likelihood function. |
dyadic |
Logical parameter indicating whether or not the original time series is dyadic in length. |
The variance-covariance matrix of the original time series is approximated
by its wavelet-based equivalent. A Whittle-type likelihood is then
constructed where the sums of squared wavelet coefficients are compared to
bandpass filtered version of the true spectral density function.
Minimization occurs for the fractional difference parameter and the
Gegenbauer frequency , while the innovations variance is
subsequently estimated.
List containing the maximum likelihood estimates (MLEs) of
, and , along with the value of the
likelihood for those estimates.
B. Whitcher
Whitcher, B. (2004) Wavelet-based estimation for seasonal long-memory processes, Technometrics, 46, No. 2, 225-238.
Computes the variance of a seasonal persistent (SP) process using a hypergeometric series expansion.
spp.var(d, fG, sigma2 = 1) Hypergeometric(a, b, c, z)spp.var(d, fG, sigma2 = 1) Hypergeometric(a, b, c, z)
d |
Fractional difference parameter. |
fG |
Gegenbauer frequency. |
sigma2 |
Innovations variance. |
a, b, c, z
|
Parameters for the hypergeometric series. |
See Lapsa (1997). The subroutine to compute a hypergeometric series was taken from Numerical Recipes in C.
The variance of an SP process.
B. Whitcher
Lapsa, P.M. (1997) Determination of Gegenbauer-type random process models. Signal Processing 63, 73-90.
Press, W.H., S.A. Teukolsky, W.T. Vetterling and B.P. Flannery (1992) Numerical Recipes in C, 2nd edition, Cambridge University Press.
Produces the modulus squared of the Fourier transform for a given filtering sequence.
squared.gain(wf.name, filter.seq = "L", n = 512)squared.gain(wf.name, filter.seq = "L", n = 512)
wf.name |
Character string of wavelet filter. |
filter.seq |
Character string of filter sequence. |
n |
Length of zero-padded filter. Frequency resolution will be
|
Uses cascade subroutine to compute the squared gain function from a
given filtering sequence.
Squared gain function.
B. Whitcher
par(mfrow=c(2,2)) f.seq <- "H" plot(0:256/512, squared.gain("d4", f.seq), type="l", ylim=c(0,2), xlab="frequency", ylab="L = 4", main="Level 1") lines(0:256/512, squared.gain("fk4", f.seq), col=2) lines(0:256/512, squared.gain("mb4", f.seq), col=3) abline(v=c(1,2)/4, lty=2) legend(-.02, 2, c("Daubechies", "Fejer-Korovkin", "Minimum-Bandwidth"), lty=1, col=1:3, bty="n", cex=1) f.seq <- "HL" plot(0:256/512, squared.gain("d4", f.seq), type="l", ylim=c(0,4), xlab="frequency", ylab="", main="Level 2") lines(0:256/512, squared.gain("fk4", f.seq), col=2) lines(0:256/512, squared.gain("mb4", f.seq), col=3) abline(v=c(1,2)/8, lty=2) f.seq <- "H" plot(0:256/512, squared.gain("d8", f.seq), type="l", ylim=c(0,2), xlab="frequency", ylab="L = 8", main="") lines(0:256/512, squared.gain("fk8", f.seq), col=2) lines(0:256/512, squared.gain("mb8", f.seq), col=3) abline(v=c(1,2)/4, lty=2) f.seq <- "HL" plot(0:256/512, squared.gain("d8", f.seq), type="l", ylim=c(0,4), xlab="frequency", ylab="", main="") lines(0:256/512, squared.gain("fk8", f.seq), col=2) lines(0:256/512, squared.gain("mb8", f.seq), col=3) abline(v=c(1,2)/8, lty=2)par(mfrow=c(2,2)) f.seq <- "H" plot(0:256/512, squared.gain("d4", f.seq), type="l", ylim=c(0,2), xlab="frequency", ylab="L = 4", main="Level 1") lines(0:256/512, squared.gain("fk4", f.seq), col=2) lines(0:256/512, squared.gain("mb4", f.seq), col=3) abline(v=c(1,2)/4, lty=2) legend(-.02, 2, c("Daubechies", "Fejer-Korovkin", "Minimum-Bandwidth"), lty=1, col=1:3, bty="n", cex=1) f.seq <- "HL" plot(0:256/512, squared.gain("d4", f.seq), type="l", ylim=c(0,4), xlab="frequency", ylab="", main="Level 2") lines(0:256/512, squared.gain("fk4", f.seq), col=2) lines(0:256/512, squared.gain("mb4", f.seq), col=3) abline(v=c(1,2)/8, lty=2) f.seq <- "H" plot(0:256/512, squared.gain("d8", f.seq), type="l", ylim=c(0,2), xlab="frequency", ylab="L = 8", main="") lines(0:256/512, squared.gain("fk8", f.seq), col=2) lines(0:256/512, squared.gain("mb8", f.seq), col=3) abline(v=c(1,2)/4, lty=2) f.seq <- "HL" plot(0:256/512, squared.gain("d8", f.seq), type="l", ylim=c(0,4), xlab="frequency", ylab="", main="") lines(0:256/512, squared.gain("fk8", f.seq), col=2) lines(0:256/512, squared.gain("mb8", f.seq), col=3) abline(v=c(1,2)/8, lty=2)
Stack plot of an object. This function attempts to mimic a function called
stack.plot in S+WAVELETS.
stackPlot( x, plot.type = c("multiple", "single"), panel = lines, log = "", col = par("col"), bg = NA, pch = par("pch"), cex = par("cex"), lty = par("lty"), lwd = par("lwd"), ann = par("ann"), xlab = "Time", main = NULL, oma = c(6, 0, 5, 0), layout = NULL, same.scale = 1:dim(x)[2], ... )stackPlot( x, plot.type = c("multiple", "single"), panel = lines, log = "", col = par("col"), bg = NA, pch = par("pch"), cex = par("cex"), lty = par("lty"), lwd = par("lwd"), ann = par("ann"), xlab = "Time", main = NULL, oma = c(6, 0, 5, 0), layout = NULL, same.scale = 1:dim(x)[2], ... )
x |
|
plot.type, panel, log, col, bg, pch, cex, lty, lwd, ann, xlab, main, oma, ...
|
See
|
layout |
Doublet defining the dimension of the panel. If not specified, the dimensions are chosen automatically. |
same.scale |
Vector the same length as the number of series to be plotted. If not specified, all panels will have unique axes. |
Produces a set of plots, one for each element (column) of x.
B. Whitcher
A recursive algorithm for detecting and locating multiple variance change points in a sequence of random variables with long-range dependence.
testing.hov(x, wf, J, min.coef = 128, debug = FALSE)testing.hov(x, wf, J, min.coef = 128, debug = FALSE)
x |
Sequence of observations from a (long memory) time series. |
wf |
Name of the wavelet filter to use in the decomposition. |
J |
Specifies the depth of the decomposition. This must be a number
less than or equal to |
min.coef |
Minimum number of wavelet coefficients for testing purposes. Empirical results suggest that 128 is a reasonable number in order to apply asymptotic critical values. |
debug |
Boolean variable: if set to |
For details see Section 9.6 of Percival and Walden (2000) or Section 7.3 in Gencay, Selcuk and Whitcher (2001).
Matrix whose columns include (1) the level of the wavelet transform where the variance change occurs, (2) the value of the test statistic, (3) the DWT coefficient where the change point is located, (4) the MODWT coefficient where the change point is located. Note, there is currently no checking that the MODWT is contained within the associated support of the DWT coefficient. This could lead to incorrect estimates of the location of the variance change.
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
dwt, modwt, rotcumvar,
mult.loc.
Quarterly U.S. tourism figures from 1960:1 to 1999:4.
data(tourism)data(tourism)
A vector containing 160 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Monthly U.S. unemployment figures from 1948:1 to 1999:12.
data(unemploy)data(unemploy)
A vector containing 624 observations.
Unknown.
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Upsamples a given vector.
up.sample(x, f, y = NA)up.sample(x, f, y = NA)
x |
vector of observations |
f |
frequency of upsampling; e.g, 2, 4, etc. |
y |
value to upsample with; e.g., NA, 0, etc. |
A vector twice its length.
B. Whitcher
Any basic signal processing text.
Converts name of wavelet filter to filter coefficients.
wave.filter(name)wave.filter(name)
name |
Character string of wavelet filter. |
Simple switch statement selects the appropriate filter.
List containing the following items:
L |
Length of the wavelet filter. |
hpf |
High-pass filter coefficients. |
lpf |
Low-pass filter coefficients. |
B. Whitcher
Daubechies, I. (1992) Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM: Philadelphia.
Doroslovacki (1998) On the least asymmetric wavelets, IEEE Transactions for Signal Processing, 46, No. 4, 1125-1130.
Morris and Peravali (1999) Minimum-bandwidth discrete-time wavelets, Signal Processing, 76, No. 2, 181-193.
Nielsen, M. (2000) On the Construction and Frequency Localization of Orthogonal Quadrature Filters, Journal of Approximation Theory, 108, No. 1, 36-52.
Produces an estimate of the multiscale variance, covariance or correlation along with approximate confidence intervals.
wave.variance(x, type = "eta3", p = 0.025) wave.covariance(x, y) wave.correlation(x, y, N, p = 0.975)wave.variance(x, type = "eta3", p = 0.025) wave.covariance(x, y) wave.correlation(x, y, N, p = 0.975)
x |
first time series |
type |
character string describing confidence interval calculation;
valid methods are |
p |
(one minus the) two-sided p-value for the confidence interval |
y |
second time series |
N |
length of time series |
The time-independent wavelet variance is basically the average of the squared wavelet coefficients across each scale. As shown in Percival (1995), the wavelet variance is a scale-by-scale decomposition of the variance for a stationary process, and certain non-stationary processes.
Matrix with as many rows as levels in the wavelet transform object. The first column provides the point estimate for the wavelet variance, covariance, or correlation followed by the lower and upper bounds from the confidence interval.
B. Whitcher
Gencay, R., F. Selcuk and B. Whitcher (2001) An Introduction to Wavelets and Other Filtering Methods in Finance and Economics, Academic Press.
Percival, D. B. (1995) Biometrika, 82, No. 3, 619-631.
Percival, D. B. and A. T. Walden (2000) Wavelet Methods for Time Series Analysis, Cambridge University Press.
Whitcher, B., P. Guttorp and D. B. Percival (2000) Wavelet Analysis of Covariance with Application to Atmospheric Time Series, Journal of Geophysical Research, 105, No. D11, 14,941-14,962.
## Figure 7.3 from Gencay, Selcuk and Whitcher (2001) data(ar1) ar1.modwt <- modwt(ar1, "haar", 6) ar1.modwt.bw <- brick.wall(ar1.modwt, "haar") ar1.modwt.var2 <- wave.variance(ar1.modwt.bw, type="gaussian") ar1.modwt.var <- wave.variance(ar1.modwt.bw, type="nongaussian") par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1) matplot(2^(0:5), ar1.modwt.var2[-7,], type="b", log="xy", xaxt="n", ylim=c(.025, 6), pch="*LU", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="") matlines(2^(0:5), as.matrix(ar1.modwt.var)[-7,2:3], type="b", pch="LU", lty=1, col=3) axis(side=1, at=2^(0:5)) legend(1, 6, c("Wavelet variance", "Gaussian CI", "Non-Gaussian CI"), lty=1, col=c(1,4,3), bty="n") ## Figure 7.8 from Gencay, Selcuk and Whitcher (2001) data(exchange) returns <- diff(log(as.matrix(exchange))) returns <- ts(returns, start=1970, freq=12) wf <- "d4" J <- 6 demusd.modwt <- modwt(returns[,"DEM.USD"], wf, J) demusd.modwt.bw <- brick.wall(demusd.modwt, wf) jpyusd.modwt <- modwt(returns[,"JPY.USD"], wf, J) jpyusd.modwt.bw <- brick.wall(jpyusd.modwt, wf) returns.modwt.cov <- wave.covariance(demusd.modwt.bw, jpyusd.modwt.bw) par(mfrow=c(1,1), las=0, mar=c(5,4,4,2)+.1) matplot(2^(0:(J-1)), returns.modwt.cov[-(J+1),], type="b", log="x", pch="*LU", xaxt="n", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="Wavelet Covariance") axis(side=1, at=2^(0:7)) abline(h=0) returns.modwt.cor <- wave.correlation(demusd.modwt.bw, jpyusd.modwt.bw, N = dim(returns)[1]) par(mfrow=c(1,1), las=0, mar=c(5,4,4,2)+.1) matplot(2^(0:(J-1)), returns.modwt.cor[-(J+1),], type="b", log="x", pch="*LU", xaxt="n", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="Wavelet Correlation") axis(side=1, at=2^(0:7)) abline(h=0)## Figure 7.3 from Gencay, Selcuk and Whitcher (2001) data(ar1) ar1.modwt <- modwt(ar1, "haar", 6) ar1.modwt.bw <- brick.wall(ar1.modwt, "haar") ar1.modwt.var2 <- wave.variance(ar1.modwt.bw, type="gaussian") ar1.modwt.var <- wave.variance(ar1.modwt.bw, type="nongaussian") par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1) matplot(2^(0:5), ar1.modwt.var2[-7,], type="b", log="xy", xaxt="n", ylim=c(.025, 6), pch="*LU", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="") matlines(2^(0:5), as.matrix(ar1.modwt.var)[-7,2:3], type="b", pch="LU", lty=1, col=3) axis(side=1, at=2^(0:5)) legend(1, 6, c("Wavelet variance", "Gaussian CI", "Non-Gaussian CI"), lty=1, col=c(1,4,3), bty="n") ## Figure 7.8 from Gencay, Selcuk and Whitcher (2001) data(exchange) returns <- diff(log(as.matrix(exchange))) returns <- ts(returns, start=1970, freq=12) wf <- "d4" J <- 6 demusd.modwt <- modwt(returns[,"DEM.USD"], wf, J) demusd.modwt.bw <- brick.wall(demusd.modwt, wf) jpyusd.modwt <- modwt(returns[,"JPY.USD"], wf, J) jpyusd.modwt.bw <- brick.wall(jpyusd.modwt, wf) returns.modwt.cov <- wave.covariance(demusd.modwt.bw, jpyusd.modwt.bw) par(mfrow=c(1,1), las=0, mar=c(5,4,4,2)+.1) matplot(2^(0:(J-1)), returns.modwt.cov[-(J+1),], type="b", log="x", pch="*LU", xaxt="n", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="Wavelet Covariance") axis(side=1, at=2^(0:7)) abline(h=0) returns.modwt.cor <- wave.correlation(demusd.modwt.bw, jpyusd.modwt.bw, N = dim(returns)[1]) par(mfrow=c(1,1), las=0, mar=c(5,4,4,2)+.1) matplot(2^(0:(J-1)), returns.modwt.cor[-(J+1),], type="b", log="x", pch="*LU", xaxt="n", lty=1, col=c(1,4,4), xlab="Wavelet Scale", ylab="Wavelet Correlation") axis(side=1, at=2^(0:7)) abline(h=0)
Produces an estimate of the multiscale variance with approximate confidence intervals using the 2D MODWT.
wave.variance.2d(x, p = 0.025)wave.variance.2d(x, p = 0.025)
x |
image |
p |
(one minus the) two-sided p-value for the confidence interval |
The wavelet variance is basically the average of the squared wavelet coefficients across each scale and direction of an image. As shown in Mondal and Percival (2012), the wavelet variance is a scale-by-scale decomposition of the variance for a stationary spatial process, and certain non-stationary spatial processes.
Data frame with 3J+1 rows.
B. Whitcher
Mondal, D. and D. B. Percival (2012). Wavelet variance analysis for random fields on a regular lattice. IEEE Transactions on Image Processing 21, 537–549.
Create a wavelet filter at arbitrary scale.
wavelet.filter(wf.name, filter.seq = "L", n = 512)wavelet.filter(wf.name, filter.seq = "L", n = 512)
wf.name |
Character string of wavelet filter. |
filter.seq |
Character string of filter sequence. |
n |
Length of zero-padded filter. Frequency resolution will be
|
Uses cascade subroutine to compute higher-order wavelet coefficient
vector from a given filtering sequence.
Vector of wavelet coefficients.
B. Whitcher
Bruce, A. and H.-Y. Gao (1996). Applied Wavelet Analysis with S-PLUS, Springer: New York.
Doroslovacki, M. L. (1998) On the least asymmetric wavelets, IEEE Transactions on Signal Processing, 46, No. 4, 1125-1130.
Daubechies, I. (1992) Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM: Philadelphia.
Morris and Peravali (1999) Minimum-bandwidth discrete-time wavelets, Signal Processing, 76, No. 2, 181-193.
Nielsen, M. (2001) On the Construction and Frequency Localization of Finite Orthogonal Quadrature Filters, Journal of Approximation Theory, 108, No. 1, 36-52.
## Figure 4.14 in Gencay, Selcuk and Whitcher (2001) par(mfrow=c(3,1), mar=c(5-2,4,4-1,2)) f.seq <- "HLLLLL" plot(c(rep(0,33), wavelet.filter("mb4", f.seq), rep(0,33)), type="l", xlab="", ylab="", main="D(4) in black, MB(4) in red") lines(c(rep(0,33), wavelet.filter("d4", f.seq), rep(0,33)), col=2) plot(c(rep(0,35), -wavelet.filter("mb8", f.seq), rep(0,35)), type="l", xlab="", ylab="", main="D(8) in black, -MB(8) in red") lines(c(rep(0,35), wavelet.filter("d8", f.seq), rep(0,35)), col=2) plot(c(rep(0,39), wavelet.filter("mb16", f.seq), rep(0,39)), type="l", xlab="", ylab="", main="D(16) in black, MB(16) in red") lines(c(rep(0,39), wavelet.filter("d16", f.seq), rep(0,39)), col=2)## Figure 4.14 in Gencay, Selcuk and Whitcher (2001) par(mfrow=c(3,1), mar=c(5-2,4,4-1,2)) f.seq <- "HLLLLL" plot(c(rep(0,33), wavelet.filter("mb4", f.seq), rep(0,33)), type="l", xlab="", ylab="", main="D(4) in black, MB(4) in red") lines(c(rep(0,33), wavelet.filter("d4", f.seq), rep(0,33)), col=2) plot(c(rep(0,35), -wavelet.filter("mb8", f.seq), rep(0,35)), type="l", xlab="", ylab="", main="D(8) in black, -MB(8) in red") lines(c(rep(0,35), wavelet.filter("d8", f.seq), rep(0,35)), col=2) plot(c(rep(0,39), wavelet.filter("mb16", f.seq), rep(0,39)), type="l", xlab="", ylab="", main="D(16) in black, MB(16) in red") lines(c(rep(0,39), wavelet.filter("d16", f.seq), rep(0,39)), col=2)
data(xbox)data(xbox)
A 128 128 matrix.
S+WAVELETS.
Bruce, A., and H.-Y. Gao (1996) Applied Wavelet Analysis with S-PLUS, Springer: New York.