R语言实现Mallat金字塔小波分解更换D4滤波器结果不匹配求解

我正尝试从零实现Mallat金字塔算法，来获取小波分解系数。

快速小波变换（Fast wavelet transform）算法核心是包含高通、低通两个滤波器，需要迭代将观测数据、上一级分解系数与滤波器做卷积，再完成下采样操作实现分解，这正是我当前的实现目标。

我将自行计算得到的系数与wavethresh包输出的各层级系数做对比验证，Haar小波的实现结果已经可以完美匹配，效果符合预期。

针对Haar小波我已经完成了可用的代码实现，精度表现良好，代码如下：

#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)
}

验证运行方式如下：先生成长度为2的整数次幂的信号，再调用函数获取系数：

J <- 10    
x <- rnorm(2^J, 0,1)

coef <- get.coef.haar(x) 

验证代码如下，运行后所有校验项均为真，说明Haar版本实现正确：

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]))
}

但当我尝试更换为其他小波的滤波器（例如D4小波）时，得到的结果就出现了偏差，我仅修改了滤波器参数如下：

low.pass  <- c(0.4829629, 0.8365163, 0.2241439, -0.1294095)
high.pass <- c(-0.1294095, -0.2241439, 0.8365163, -0.4829629)

我推测错误出在卷积的实现逻辑上，曾尝试调用R内置的stats::filter函数实现卷积后做二倍下采样，但是结果无法匹配；使用stats::convolve也没有得到正确结果，暂时找不到解决方法。我预期的输出结果应该和下述调用的输出一致：

wave <- wavethresh::wd(data = x, filter.number =2, family = "DaubExPhase")

内容的提问来源于stack exchange，提问作者YetAnotherUsr