基于PyWavelets的逆平稳小波变换信号重构失配问题咨询

我尝试使用Python的PyWavelets库提供的逆平稳小波变换重构所有层级的近似系数和细节系数，编写的代码如下：

import numpy as np
import pywt
import matplotlib.pyplot as plt

def UDWT(Btotal, wname, Lps, Hps, edge_eff):
    Br =  Btotal[0]; Bt =  Btotal[1]; Bn =  Btotal[2]

    ## Set parameters needed for UDWT
    samplelength=len(Br)

    # If length of data is odd, turn into even numbered sample by getting rid 
    # of one point
    if np.mod(samplelength,2)>0:
        Br = Br[0:-1]
        Bt = Bt[0:-1]
        Bn = Bn[0:-1]
        
    samplelength = len(Br)

    # edge extension mode set to periodic extension by default with this
    # routine in the rice toolbox.
    pads = 2**(np.ceil(np.log2(abs(samplelength))))-samplelength  # for edge extension, This function 
                                                                  # returns 2^{ the next power of 2 }for input: samplelength

    ## Do the UDWT decompositon and reconstruction
    keep_all = {}
    for m in range(3):
        # Gets the data size up to the next power of 2 due to UDWT restrictions
        # Although periodic extension is used for the wavelet edge handling we are
        # getting the data up to the next power of 2 here by extending the data
        # sample with a constant value
        if (m==0):
            y = np.pad(Br,pad_width = int(pads/2) ,constant_values=np.nan)
        elif (m==1):
            y = np.pad(Bt,pad_width = int(pads/2) ,constant_values=np.nan)
        else:
            y = np.pad(Bn,pad_width = int(pads/2) ,constant_values=np.nan)

        # Decompose the signal using the UDWT
        nlevel = min(pywt.swt_max_level(y.shape[-1]), 8)  # Level of decomposition, impose upper limit 10
        Coeff  = pywt.swt(y, wname, nlevel)                # List of approximation and details coefficients 
                                                           # pairs in order similar to wavedec function:
                                                           # [(cAn, cDn), ..., (cA2, cD2), (cA1, cD1)]
        # Assign approx: swa and details: swd to 
        swa  = np.zeros((len(y),nlevel))
        swd  = np.zeros((len(y),nlevel))
       
        for o in range(nlevel):
            swa[:,o]  = Coeff[o][0]
            swd[:,o]  = Coeff[o][1]

        # Reconstruct all the approximations and details at all levels
        mzero = np.zeros(np.shape(swd))
        A     = mzero

        coeffs_inverse = list(zip(swa.T,mzero.T))
        invers_res  = pywt.iswt(coeffs_inverse, wname)
        D           = mzero

        for pp in range(nlevel):
            swcfs = mzero
            swcfs[:,pp] = swd[:,pp]
            coeffs_inverse2 = list(zip(np.zeros((len(swa),1)).T , swcfs.T))
            D[:,pp]         = pywt.iswt(coeffs_inverse2, wname)
            
        for jjj in range(nlevel-1,-1,-1):
            if (jjj==nlevel-1):
                A[:,jjj] = invers_res
            else:
                A[:,jjj] = A[:,jjj+1] + D[:,jjj+1]

        # *************************************************************************
        # VERY IMPORTANT: LINEAR PHASE SHIFT CORRECTION
        # *************************************************************************
        # Correct for linear phase shift in wavelet coefficients at each level. No
        # need to do this for the low-pass filters approximations as they will be
        # reconstructed and the shift will automatically be reversed. The formula
        # for the shift has been taken from Walden's paper, or has been made up by
        # me (can't exactly remember) -- but it is verified and correct.
        # *************************************************************************
        for j in range(1,nlevel+1):
            shiftfac = Hps*(2**(j-1));  
            for l in range(1,j):
                shiftfac = int(shiftfac + Lps*(2**(l-2))*((l-2)>=0)) ;
            swd[:,j-1]   = np.roll(swd[:,j-1],shiftfac)
            
            flds = {"A": A.T,
                    "D": D.T,
                    "swd"  : swd.T,
                    }
            Btot = ['Br', 'Bt', 'Bn'] # Used Just to name files

        keep_all[str(Btot[m])] = flds 

    # 1) Put all the files together into a cell structure
    Apr = {}
    Swd = {}

    pads  = int(pads)
    names = ['Br', 'Bt', 'Bn']
    for kk in range(3):
        A              = keep_all[names[kk]]['A']
        Apr[names[kk]] = A[:,int(pads/2):len(A)-int(pads/2)]

        swd            = keep_all[names[kk]]['swd']
        Swd[names[kk]] = swd[:,int(pads/2):len(A)-int(pads/2)]

    # Returns filters list for the current wavelet in the following order
    wavelet       = pywt.Wavelet(wname)
    [h_0,h_1,_,_] = wavelet.inverse_filter_bank
    filterlength  = len(h_0)
    
    if edge_eff:
        # 2)  Getting rid of the edge effects; to keep edges skip this section
        for j in range(1,nlevel+1):
            extra = int((2**(j-2))*filterlength) # give some reasoning for this eq
            for m in range(3):
                # for approximations
                Apr[names[m]][j-1][0:extra]   = np.nan
                Apr[names[m]][j-1][-extra:-1] = np.nan
                # for details
                Swd[names[m]][j-1][0:extra]   = np.nan
                Swd[names[m]][j-1][-extra:-1] = np.nan
    
    return  Apr, Swd, pads, nlevel


aa = np.sin(np.linspace(0,2*np.pi,100000))+0.05*np.random.rand(100000)
bb = np.cos(np.linspace(0,2*np.pi,100000))+0.05*np.random.rand(100000)
cc = np.cos(np.linspace(0,4*np.pi,100000))+0.05*np.random.rand(100000)
Btotal = [aa,bb,cc]

wname     ='coif2'
Lps       = 7;          #   Low pass filter phase shift for level 1 Coiflet2
Hps       = 4;          #   High pass filter phase shift for level 1 Coiflet2
edge_eff = False

Apr, Swd, pads, nlevel = UDWT(Btotal, wname, Lps, Hps, edge_eff)

### Add the details at all levels with the highest level approximations 
##  to compare with the original timeseries. (The equation shown in website)
new = Swd['Br'][0]
for i in range(1,nlevel):
    new = Swd['Br'][i]+new
sig = Apr['Br'][-1]+new

### Now plot to comapre ##
## Reconstructed signal 1
plt.plot(sig)
### Second way to get reconstructed signal
### aa first level details with approximations
plt.plot(Apr['Br'][-1] +Swd['Br'][-1] )
### Original signal
plt.plot(aa)
plt.show()

我参照Matlab小波工具箱文档描述的流程进行实现，但得到的重构时间序列与原始信号无法完全匹配，请问该如何排查问题并解决？

• 填充值错误使用NaN
  PyWavelets的swt和iswt接口无法处理NaN数值，你使用np.nan作为填充值会导致小波变换计算全部异常。建议改为边界值填充：比如使用np.pad的mode='symmetric'对称填充，或者用0、信号均值作为常量填充值。
• 相位修正逻辑顺序错误
  你的代码中先执行了各层级细节系数的重构计算，后执行相位修正逻辑，修正后的swd完全没有参与到D的重构过程中，相位修正完全失效。需要调整代码顺序：先完成swd的相位偏移修正，再执行各层级细节系数的逆变换重构。
• 系数顺序与重构求和逻辑不匹配
  pywt.swt返回的系数顺序为[(最高层级近似cAn, 最高层级细节cDn), ..., (最低层级近似cA1, 最低层级细节cD1)]，你当前的求和逻辑没有匹配该顺序，导致重构时系数累加错误。建议先输出各层级swa、swd的数值，确认顺序后调整累加逻辑。
• 相位偏移参数与偏移方向不匹配
  你硬编码了Coif2小波的相位偏移参数Lps=7、Hps=4，需要确认该参数是否和PyWavelets内置的Coif2滤波器参数匹配，同时确认np.roll的偏移方向是否正确：正参数为右移，负参数为左移，方向错误会导致整体信号偏移。
• 填充裁剪不对称
  你计算的pads可能为奇数，左右填充长度不对称，后续裁剪时使用int(pads/2)会导致信号错位。建议调整填充逻辑确保左右填充长度一致，裁剪时分别取左右对应长度的偏移量。
• 先关闭边缘效应处理做基础验证
  排查时先将edge_eff设为False，避免边缘置NaN的逻辑干扰核心重构正确性的验证，确认核心重构逻辑正确后再开启边缘效应处理。

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