fields (version 5.02)

Wtransform: Quadratic W wavelet transform for 1-d vectors or rectangular or cylindrical images

Description

Finds the forward or inverse discrete wavelet transform using the quadratic W basis.

Usage

Wtransform(x, inv = FALSE, transpose = FALSE, cut.min = 8)
Wtransform.image(x, inv=FALSE, transpose=FALSE, cut.min=8)
Wtransform.cylinder.image(x, inv=FALSE, transpose=FALSE, 
byX=TRUE,cut.min=8)

Value

A matrix the same size as x.

References

Nychka,D. , C. Wikle and J.A. Royle, (2002) Multi-resolution models for nonstationary spatial covariance functions Statistical Modelling 2 315-332.

See also the original report on W matrices:

Man Kam Kwong and P. T. Peter Tang, "W-Matrices, Nonorthogonal Multiresolution Analysis, and Finite Signals of Arbitrary Length," Preprint MCS-P449-0794, September 1994

ftp://info.mcs.anl.gov/pub/tech_reports/reports/P449.ps.Z

Details

The wavelet transform can be thought as matrix multiplication W %*% vec(x) where vec(x) is the matrix x stacked by columns. The inverse transform is inv(W) %*% vec(x) and transpose is t(W) %*% vec( x). One can interpret the columns of inv(W) as basis functions and they follow the usual pattern of translations and dilations of mother and father wavelets. (See example below.)

Another way of thinking of the transformation is by recursion: apply a smoother and a "rougher" to a vector and taking every other value. Now apply the same operation to the smooth results, now half the length of the previous vector. At each step one reduced the vector by a half and cut.min specifies the size when to stop. The function WQS performs one step of the recursion with smoother weights (-1, 3, 3, -1) and roughening weights (-1, 3, -3, 1) away from the edges. WQS actually works for a matrix where each column is transformed in this manner. Boundary adjustments are made to preserve the shape of the basis functions. (See example below.) By convention the returned matrix has the smoothed results in the first half of the columns and the rough results in the second half. The discrete wavelet transform performs this operation recursively on the smoothed results until the smoothed vector is less than a set size. Because each step reduces the size of the vector by 2 it only makes sense to apply this algorithm to vectors whose length is dyadic or the product of a small integer and a dyadic (e.g. 96=3*32). The precise tests are done by dyadic.check and dyadic.2check and if n is the dimension the constraint it that n=p*2^L where L is less or equal to than cut.min. The transform will result in p scale (father) basis functions if cut.min is equal to p.

At the end of the day this recursive algorithm defines a linear transformation of the original image to something that we call the wavelet decomposition. This is the full matrix W mentioned above. It is also possible to express multiplication of inv(W) and transposes by a similar recursive scheme and related sets of filter weights. (Try WQSi( WQS(x)) as a test.) The reader is referred to WQSi WQS.T and WQSi.T for the filtering primitives. Finally it should be noted that the particular wavelets chosen here are not orthogonal by have nice smooth properties and the father wavelet resembles a Gaussian density while the mother looks like its derivative. (See example below for some plots of the implied basis functions.) Note that Wtransform is "vectorized" so that with little extra overhead one can do transforms for many separate 1-d series with one call. In particular, Wtransform(diag(1,128),inv=TRUE) will neatly generate the W matrix given above.

For two dimensions for the basic step one applies WQS to the columns of the matrix and then to the rows. e.g. t(WQS(t(WQS(x))). This primitive step is implemented in WQS2d. The final algorithm calls WQS2d or its variants repeatedly on a matrix that decreases by a factor of two in size along each dimension at each iteration.

If the cylinder variant is specified the transform uses periodic boundaries on one of the coordinates. This is useful for data on zonal section of a sphere where a constant line of latitude should be periodic. For byX=T the wavelet basis functions wrap on the x-axis when an image plot is made. (See example below.) This convention may cause some confusion because R experts will know that the image plot rotates the the matrix so that the (1,1) element is at the lower left corner. Thus enforcing periodicity along the X-coordinate of the image pertains to the columns of the matrix used to represent the image in R. (Compare matrix( 1:10,2,5) to image( matrix( 1:10,2,5)).

See Also

Wtransform.image, WQS, WQSi

Examples

# W transform of a simple function
x<- seq( 0,1,256)
y<- x*(1-x)**3 + ifelse( x>.3, x,0)

Wtransform( y)-> out

# plot the implied wavelet basis functions
ID<- diag( 1, 256)
WQS.basis( 256)-> Wbasis
set.panel(2,2)
matplot( 1:256, Wbasis[,1:8], type="l", lty=c(1,2), col=2)
title("Father")
matplot( 1:256, Wbasis[,9:16], type="l", lty=c(1,2), col=2)
title("Mother")
matplot( 1:256, Wbasis[,17:32], type="l", lty=c(1,2), col=2)
title("Mother scaled by 2")
matplot( 1:256, Wbasis[,33:64], type="l", lty=c(1,2), col=2)
title("Mother scaled by 4")

set.panel( 1,1)

# test that the transform works


# Precise definition of what the transform is doing in terms of 
# explicit matrix multiplication  all of 
# these should be machine zero
# Note that the direct matrix multiplications will be substantially slower
# for large vectors.
# y<- rnorm( 256)
# y<- y /sqrt(mean( y**2))


#sqrt(mean( c( Wtransform(y, inv=TRUE) - Wbasis%*%y)**2))
#sqrt(mean( c(Wtransform(y, inv=TRUE, transpose=TRUE) - t(Wbasis)%*%y)**2))
#sqrt(mean( c(Wtransform(y) - solve(Wbasis)%*%y)**2))
#sqrt(mean( c(Wtransform(y, transpose=TRUE) - t(solve(Wbasis))%*%y)**2))

#
## 2-d examples
#


# Wtransform of John Lennon image
data(lennon)
look<- Wtransform.image( lennon)
# 
### take a look: 
 image.plot( look)

 
# a diagonal detail basis function 
# we find this by just multipling W by a unit vector!

temp<- matrix(0, nrow=32, ncol=32) 
temp[8,5]<- 1 
look<- Wtransform.image( temp , inv=TRUE, cut.min=4)
image( look)
title("diagonal detail W-wavelet")

#just for fun: redo this example for all indices up to 8!
#
#set.panel( 8,8)
#par( mar=c(0,0,0,0))
#for (  k in (1:8)){
#for ( j in (1:8)){
#temp<- matrix( 0 , nx,ny)
#temp[k,j] <- 1
#Wtransform.image( temp, inv=T, cut.min=cut.min)-> look
#image( look, axes=FALSE, xlab="", ylab="")
#box()
#}
#}


# examine a basis function to see periodic condition enforces along 
# X axis of image. 
temp<- matrix(0, nrow=32, ncol=32)
temp[8,5]<- 1
image( Wtransform.cylinder.image( temp , inv=TRUE, cut.min=4))
# now along Y-axis
image( Wtransform.cylinder.image( temp , inv=TRUE, cut.min=4, byX=FALSE))

# reset panel 
set.panel( 1,1)