I'm trying to implement "by the scratch" Mallat's pyramid algorithm to obtain the wavelet decomposition coefficients.
The Fast wavelet transform algorithm is described in MathWorks in the hyperlink. Essentially you have two filters, a high-pass filter and a low-pass filter, around which you will convolve your observed data with the filters and the previous coefficients iteratively. The algorithm is supposed to perform the convolution and downsampling for the decomposition, which is what I'm trying to do.
For means of comparison, I'm comparing my coefficients which I calculated from the scratch with the coefficients in each level of the output of "wavethresh" package. The results are a perfect match to it, which is nice.
For the Haar wavelet, I managed to write the code and it worked with good precision. Here's my code:
#Mallat Pyramid algorithm for Haar system
get.coef.haar <- function(x){
#x is my signal
#High pass filter
high.pass <- c(1/sqrt(2),-1/sqrt(2))
#Low pass filter
low.pass <- c(1/sqrt(2), 1/sqrt(2))
N <- length(x)
if(is.pow.2(N)==FALSE){stop("Signal is not of a power of two")}
J <- log(N,base = 2)
#ncol = j, nrow=k
#C and D matrices for the Mallat pyramid algorithm
c <- matrix(NA,nrow = 2^(J-1), ncol = J)
d <- matrix(NA,nrow = 2^(J-1), ncol = J)
#L for the convolution
L <- length(high.pass)
#Double for to perform the pyramid algorithm
#Basically convolves and performs dyadic downsampling
for(j in 1:J){
for(k in 1:2^(J-j)){
if(j==1){
m <- 1:L-1
c[k,j] <- x[2*k-m] %*% low.pass
d[k,j] <- x[2*k-m] %*% high.pass
}else{
m <- 1:L-1
c[k,j] <- c[2*k-m,j-1] %*% low.pass
d[k,j] <- c[2*k-m,j-1] %*% high.pass
}
}
}
colnames(c) <- c(paste0("c",J:1-1))
colnames(d) <- c(paste0("d",J:1-1))
mat <- list(C=c, D = (-1)*d)
return(mat)
}
To check how it works, basically you can create a signal with length power of two, 2^J for an integer J, then get the coefficients.
J <- 10
x <- rnorm(2^J, 0,1)
coef <- get.coef.haar(x)
To check that it works, you can do:
wave <- wavethresh::wd(data = x, filter.number =1, family = "DaubExPhase")
for(i in 0:(J-1)){
print(all(accessC(wave, level=i)==coef$C[1:2^i,ncol(coef$C)-i]))
}
for(i in 0:(J-1)){
print(all(accessD(wave, level=i)==coef$D[1:2^i,ncol(coef$D)-i]))
}
which is all true to me.
However when I try to use another filter for another wavelet (say D4), my results differ. Running the same code above, but changing the filters:
low.pass <- c(0.4829629, 0.8365163, 0.2241439, -0.1294095)
high.pass <- c(-0.1294095, -0.2241439, 0.8365163, -0.4829629)
The error I assume is in the way that I wrote the code to perform the convolution. I tried (unsuccessfully) to use the R function stats::filter and do the dyadic downsampling, but the result was not a match. Also stats::convolve didn't bring me the same result, and I don't know how to solve this.
The expected results should be:
wave <- wavethresh::wd(data = x, filter.number =2, family = "DaubExPhase")
Updated answer:
After thorough review of the algorithm, it turns out your index is incorrect. m should be 1:L-2. Also, index for C and D are different. The correct indexes with periodic adjustment are:
c_index <- ifelse(2*k+m>0,2*k+m,2^(J)+2*k+m)
c_index <- ifelse(c_index<2^J+1,c_index,c_index-2^(J))
d_index <- ifelse(2*k+m-2>0,2*k+m-2,2^(J)+2*k+m-2)
d_index <- ifelse(d_index<2^J+1,d_index,d_index-2^(J))
For j>1, the indexes become:
c_index <- ifelse(2*k+m>0,2*k+m,2^(J-j+1)+2*k+m)
c_index <- ifelse(c_index<2^(J-j+1)+1,c_index,c_index-2^(J-j+1))
d_index <- ifelse(2*k+m-2>0,2*k+m-2,2^(J-j+1)+2*k+m-2)
d_index <- ifelse(d_index<2^(J-j+1)+1,d_index,d_index-2^(J-j+1))
This takes the size of each level into account.
so the entire for-loop is:
for(j in 1:J){
for(k in 1:2^(J-j)){
if(j==1){
m <- 1:L-2
c_index <- ifelse(2*k+m>0,2*k+m,2^(J)+2*k+m)
c_index <- ifelse(c_index<2^J+1,c_index,c_index-2^(J))
d_index <- ifelse(2*k+m-2>0,2*k+m-2,2^(J)+2*k+m-2)
d_index <- ifelse(d_index<2^J+1,d_index,d_index-2^(J))
c[k,j] <- x[c_index] %*% low.pass
d[k,j] <- x[d_index] %*% high.pass
}else{
m <- 1:L-2
c_index <- ifelse(2*k+m>0,2*k+m,2^(J-j+1)+2*k+m)
c_index <- ifelse(c_index<2^(J-j+1)+1,c_index,c_index-2^(J-j+1))
d_index <- ifelse(2*k+m-2>0,2*k+m-2,2^(J-j+1)+2*k+m-2)
d_index <- ifelse(d_index<2^(J-j+1)+1,d_index,d_index-2^(J-j+1))
c[k,j] <- c[c_index,j-1] %*% low.pass
d[k,j] <- c[d_index,j-1] %*% high.pass
}
}
}
For comparison, there are some small difference due to your filter does not have enough precision, I am able to match to at least 1e-5 (Note, it is plus sign instead of minus sign for D because you have D = (-1)*d):
# mat is the output
for(i in 0:(J-1)){
print(all(abs(accessC(wave, level=i)-mat$C[1:2^i,ncol(mat$C)-i])<1e-5))
}
for(i in 0:(J-1)){
print(all(abs(accessD(wave, level=i)+mat$D[1:2^i,ncol(mat$D)-i])<1e-5))
}
This returns all TRUE for the following:
wave <- wavethresh::wd(data = x, filter.number =1, family = "DaubExPhase")
wave <- wavethresh::wd(data = x, filter.number =2, family = "DaubExPhase")
Original answer:
The code fails because sometimes the index is negative, which is not acceptable in R.
Using the filter:
low.pass <- c(0.4829629, 0.8365163, 0.2241439, -0.1294095)
high.pass <- c(-0.1294095, -0.2241439, 0.8365163, -0.4829629)
L <- length(high.pass) # This is 4
At first iteration, this is what happened (error generated)
j <- 1
k <- 1
m <- 1:L-1
x[2*k-m]
Error in x[2 * k - m] : only 0's may be mixed with negative subscripts
2*k-m
[1] 2 1 0 -1
There should not be 0 and negative in the index.
I think what you want to do it to add the following code at the beginning of for loop to ensure no subscript out of bound:
if(2^k<L){
next
}
So the for loop becomes
for(j in 1:J){
for(k in 1:2^(J-j)){
if(2^k<L){
next
}
if(j==1){
m <- 1:L-1
c[k,j] <- x[2*k-m] %*% low.pass
d[k,j] <- x[2*k-m] %*% high.pass
}else{
m <- 1:L-1
c[k,j] <- c[2*k-m,j-1] %*% low.pass
d[k,j] <- c[2*k-m,j-1] %*% high.pass
}
}
}