Блог

Я пытаюсь восстановить детали и приближения на всех уровнях, используя обратное стационарное вейвлет-преобразование. Код выглядит следующим образом:

 def decompose_B(Btotal, wname, Lps, Hps):
    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]), 10)  # 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)
        A[:,-1]     = invers_res
        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)


            # *************************************************************************
        # 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):-int(pads/2)]

        swd            = keep_all[names[kk]]['swd']
        Swd[names[kk]] = swd[:,int(pads/2):-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)
    
    # 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 y,A,D, swa, swd, Coeff, Apr, Swd,extra

### for minimum reproducible example

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


y,A,D, swa, swd, Coeff, Apr, Swd,extra = decompose_B(Btotal, wname, Lps, Hps)
 

однако восстановленный сигнал, похоже, не соответствует оригиналу. Я думаю, что, возможно, придется что-то сделать с порядком, в котором библиотека pywavelets возвращает коэффициент, но я не смог найти ошибку.
Я добавил пример временных рядов для минимального воспроизводимого примера.