1
$\begingroup$

Background:

I am learning an interesting Fast CWT algorithm(PPCWT) by reading this paper published in 2019. The algorithm is summarized as below. The continuous wavelet transform of a signal can be calculated efficiently according to the last row of Eq.(1.4). To test this algorithm, I write a Cython module ppCWT_cymod.pyx, which is shown below the image. enter image description here

#ppCWT_cymod.pyx
import sys
import math
import numpy as np
import cython
cimport numpy as np
from libc.math cimport ceil, floor, sqrt, log, round, abs, exp

# Second derivative of Gaussian-derived wavelet
def d3GDWdx3( scale, pos ):
    x  = scale*pos
    s3 = scale**3
    p1 = s3*np.exp(-(x/2.0)**2)*x*(60.0-20.0*x**2+x**4)
    p2 = 32.0*np.sqrt(np.pi)
    return p1/p2

# Compute the Bj
def BjCoeff( freqs, Ls, N, wavelet_name="gdw", half_support=13.0, Nums=46 ):
    Dx = Ls/float(N)
    # Choose a special wavelet function
    GDW_set = {"gdw", "GDW", "Gaussian-derived wavelet", "Gaussian derived wavelet"}
    if wavelet_name in GDW_set:
        scales = Dx*freqs*2.0/np.sqrt(5.0)
        d3psi  = d3GDWdx3
        alpha  = half_support
        J      = Nums
    else:
        sys.exit("There is no wavelet function called " + wavelet_name)

    gamma  = np.linspace(-alpha,alpha,num=J+1,endpoint=True, dtype=np.float64)

    pos    = (gamma[1:]+gamma[:-1])/2.0
    ss, pp = np.meshgrid(scales,pos,indexing="ij")
    Bcoeffs= np.zeros( (len(scales),len(gamma)), dtype=np.float64)

    Bcoeffs[:,0]    = d3psi( scales, pos[0]/scales )
    Bcoeffs[:,-1]   = -d3psi( scales, pos[-1]/scales )
    Bcoeffs[:,1:-1] = d3psi( ss[:,1:], pp[:,1:]/ss[:,1:] )-d3psi( ss[:,1:], pp[:,0:-1]/ss[:,1:] )
    return gamma, Bcoeffs

# 4th-order B-spline function
cdef bspline4th( double var ):
    cdef double y
    cdef double x = abs(var)
    
    if (x<=0.5):
        y = 115.0/192.0-5.0*x**2/8.0+x**4/4.0
    elif ((x>0.5)&(x<=1.5)):
        y = (55.0+20.0*x-120.0*x**2+80.0*x**3-16.0*x**4)/96.0
    elif ((x>1.5)&(x<=2.5)):
        y = (5.0-2.0*x)**4/384.0
    elif (x>2.5):
        y = 0.0
    return y

#5th-order B-spline function
cdef bspline5th( double var ):
    cdef double y
    cdef double x = abs(var)
    if (x<=1.0):
        y = (33.0-30.0*x**2+15.0*x**4-5.0*x**5)/60.0
    elif ((x>1.0)&(x<=2.0)):
        y = (51.0+75.0*x-210.0*x**2+150.0*x**3-45.0*x**4+5.0*x**5)/120.0
    elif ((x>2.0)&(x<=3.0)):
        y = -(x-3.0)**5/120.0
    elif (x>3.0):
        y = 0.0
    return y

#6th-order B-spline function
cdef bspline6th( double var ):
    cdef double y
    cdef double x = abs(var)
    if (x<=0.5):
        y = (5887.0-4620.0*x**2+1680.0*x**4-320.0*x**6)/11520.0
    elif ((x>0.5)&(x<=1.5)):
        y = (23583.0-420.0*x-16380.0*x**2-5600.0*x**3+15120.0*x**4-\
            6720.0*x**5+960.0*x**6)/46080.0
    elif ((x>1.5)&(x<=2.5)):
        y = (4137.0+30408.0*x-59220.0*x**2+42560.0*x**3-\
            15120.0*x**4+2688.0*x**5-192.0*x**6)/23040.0
    elif ((x>2.5)&(x<=3.5)):
        y = (7.0-2.0*x)**6/46080.0
    elif (x>3.5):
        y = 0.0
    return y

#7th-order B-spline function
cdef bspline7th( double var ):
    cdef double y
    cdef double x = abs(var)
    if (x<=1.0):
        y = 151.0/315.0 - x**2/3.0+x**4/9.0-x**6/36.0+x**7/144.0
    elif ((x>1.0)&(x<=2.0)):
        y = (2472.0-392.0*x-504.0*x**2-1960.0*x**3+2520.0*x**4-\
             1176.0*x**5+252.0*x**6-21*x**7)/5040.0
    elif ((x>2.0)&(x<=3.0)):
        y = (-1112.0+12152.0*x-19320.0*x**2+13720.0*x**3-5320.0*x**4+\
            1176.0*x**5-140.0*x**6+7.0*x**7)/5040.0
    elif ((x>3.0)&(x<=4.0)):
        y = -(x-4.0)**7/5040.0
    elif (x>4.0):
        y = 0.0
    return y

# Interpolation coefficients of the signal based on 3rd-order B-spline
# This signal satisfy the zero boundary condition
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef InterpCoeffs( np.ndarray[double,ndim=1] signal ):

    cc_plus  = np.zeros_like(signal, dtype=np.double) 
    cc_minus = np.zeros_like(signal, dtype=np.double)

    cdef int    k
    cdef int    N       = len(signal)
    cdef double z1      = -2.0 + sqrt(3.0)
    cdef double z1_2    = z1**2
    cdef int    N_z1    = -<int>( round(60*log(10) / log(abs(z1))) )
    cdef double[:] y_input = signal
    cdef double[:] c_plus  = cc_plus
    cdef double[:] c_minus = cc_minus

    # Compute c_plus
    c_plus[0] = y_input[0]
    for k in range(1,N):
        c_plus[k] = y_input[k] + z1*c_plus[k-1]
    # Compute c_minus
    z1_powers    = z1_2**( np.arange(1, N_z1+1) )
    c_minus[N-1] = -6.0*( y_input[N-1]/z1 + c_plus[N-2] )*np.sum( z1_powers )
    for k in range(N-2,-1,-1):
        c_minus[k] = z1*( c_minus[k+1] - 6.0*c_plus[k] )
    c_minus_left  = c_minus[0]*( z1**np.arange(6,0,-1) )
    c_minus_right = -6.0*c_plus[N-1]*np.sum( z1_powers )*np.array([1.0, z1])
    ck            = np.concatenate((c_minus_left, c_minus, c_minus_right), axis=0)
    gk1           = np.cumsum(ck)
    gk2           = np.cumsum(gk1)
    gk3           = np.cumsum(gk2)
    gk4           = np.cumsum(gk3)
    return gk1, gk2, gk3, gk4

# 1st-th antiderivative of signal
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef AntiDerivative1st( int N, np.ndarray[double,ndim=1] InterpolationCoeffs, double new_indx ):

    cdef int k, Ak, Bk
    cdef double BSval
    cdef double out = 0.0
    cdef double[:] gk = InterpolationCoeffs[3:]

    if ( (new_indx>=0)&(new_indx<=N-1) ):
        Ak = <int>(ceil(new_indx-3))
        Bk = <int>(floor(new_indx+2))
        for k in range(Ak, Bk+1, 1):
            BSval = bspline4th(new_indx-<double>(k)-0.5)
            out   = out+gk[k+3]*BSval
        return out
    else:
        sys.exit( "Check the value of new_indx." )

# 2nd-th antiderivative of signal
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef AntiDerivative2nd( int N, np.ndarray[double,ndim=1] InterpolationCoeffs, double new_indx ):

    cdef int k, Ak, Bk
    cdef double BSval
    cdef double out = 0.0
    cdef double[:] gk = InterpolationCoeffs[2:]

    if ( (new_indx>=0)&(new_indx<=N-1) ):
        Ak = <int>(ceil(new_indx-4))
        Bk = <int>(floor(new_indx+2))
        for k in range(Ak, Bk+1, 1):
            BSval = bspline5th(new_indx-<double>(k)-1.0)
            out   = out+gk[k+4]*BSval
        return out
    else:
        sys.exit( "Check the value of new_indx." )

# 3rd-th antiderivative of signal
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef AntiDerivative3rd( int N, np.ndarray[double,ndim=1] InterpolationCoeffs, double new_indx ):

    cdef int k, Ak, Bk
    cdef double BSval
    cdef double out = 0.0
    cdef double[:] gk = InterpolationCoeffs[1:]

    if ( (new_indx>=0)&(new_indx<=N-1) ):
        Ak = <int>(ceil(new_indx-5))
        Bk = <int>(floor(new_indx+2))
        for k in range(Ak, Bk+1, 1):
            BSval = bspline6th(new_indx-<double>(k)-1.5)
            out   = out+gk[k+5]*BSval
        return out
    else:
        sys.exit( "Check the value of new_indx." )

# 4-th antiderivative of signal
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef AntiDerivative4th( int N, np.ndarray[double,ndim=1] InterpolationCoeffs, double new_indx ):

    cdef int k, Ak, Bk
    cdef double BSval
    cdef double out = 0.0
    cdef double[:] gk = InterpolationCoeffs[:]

    if ( (new_indx>=0)&(new_indx<=N-1) ):
        Ak = <int>(ceil(new_indx-6))
        Bk = <int>(floor(new_indx+2))
        for k in range(Ak, Bk+1, 1):
            BSval = bspline7th(new_indx-<double>(k)-2.0)
            out   = out+gk[k+6]*BSval
        return out
    else:
        sys.exit( "Check the value of new_indx." )

# 4th Antiderivative of the signal with zero boundaries in the whole space
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
cdef OverallAntiDerivative4th( int N, np.ndarray[double,ndim=1] InterpolationCoeffs, 
                              double C1, double C2, double C3, double C4, double new_indx ):
    cdef double x  = new_indx
    cdef double x0 = <double>(N)
    cdef double y
    if (x<0.0):
        y = 0.0
    elif ((x>=0.0)&(x<x0-1.0)):
        y = AntiDerivative4th(N, InterpolationCoeffs, x)
    elif (x>=x0-1.0):
        y = C4+C3*(x-x0+1.0)+0.5*C2*(x-x0+1.0)**2+C1*((x-x0+1.0)**3)/6.0
    return y

# CWT of the one-dimensional signal
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def ppCWT_cy( np.ndarray[double,ndim=1] freqs, np.ndarray[double,ndim=1] signal, 
              double Ls, np.ndarray[double,ndim=1] gamma, np.ndarray[double,ndim=2] Bj ):
    
    cwt_out     = np.zeros( (len(freqs),len(signal)), dtype=np.double)
    
    cdef int s, b, j
    cdef int Ns = len(freqs)
    cdef int N  = len(signal)
    cdef int Nj = len(gamma)
    cdef double C1, C2, C3, C4, F4, indx
    cdef double Dx         = Ls/<double>(N)
    cdef double[:] scales  = 2.0*freqs/sqrt(5.0)
    cdef double[:] sscales = Dx*2.0*freqs/sqrt(5.0)
    cdef double[:,:] cwt     = cwt_out

    if ( (Bj.shape[1]!=Nj)|(Bj.shape[0]!=Ns) ):
        sys.exit( "The length of Bj should be equal to that of gamma" )
    
    # Antiderivatives of signal
    gk1, gk2, gk3, gk4 = InterpCoeffs( signal )
    C1 = AntiDerivative1st( N, gk1, <double>(N)-1.0 )
    C2 = AntiDerivative2nd( N, gk2, <double>(N)-1.0 )
    C3 = AntiDerivative3rd( N, gk3, <double>(N)-1.0 )
    C4 = AntiDerivative4th( N, gk4, <double>(N)-1.0 )
 
    # Perform CWT
    for s in range(Ns):    # per scale
        for b in range(N): # per sampling point
            for j in range(Nj):
                indx     = <double>(b) - gamma[j]/sscales[s]
                F4       = OverallAntiDerivative4th( N, gk4, C1, C2, C3, C4, indx )
                cwt[s,b] = cwt[s,b] + Bj[s,j]*F4
            cwt[s,b] = Dx*sqrt(scales[s])*cwt[s,b]
    return cwt_out

Question: To check the precision of the result of my code, I choose y=exp(-x^2) as the test signal since the CWT of it has an analytic formula. The code snippet is

import numpy as np
from ppCWT_cymod import ppCWT_cy,BjCoeff

# signal
N=2048
L=10.0
Dx = L/N
freqs = np.geomspace(2*np.pi/L, N*np.pi/L, num=256)
x = Dx*np.linspace(-N/2.0,N/2.0-1,N)
y = np.exp(-x**2)   

# perform CWT
gamma, Bj = BjCoeff( freqs, L, N, wavelet_name="gdw", half_support=13, Nums=34 )
cwt = ppCWT_cy( freqs, y, L, gamma, Bj )

The result is enter image description here

Apparently, as the scale gets smaller and smaller (freqs larger), we see that wired noise arises on the right-hand side of the origin. After my examination, such a phenomenon has nothing to do with half_support, Nums. At fixed length L, the noise is more severe as sampling number N becomes larger. I really can't figure out the source of this noise. Can you help me?

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Just one small comment, always nice to have units on the axes! $\endgroup$
    – Jdip
    Commented Sep 14, 2022 at 13:09
  • $\begingroup$ @Jdip Thank you for your comment. The horizontal axis is the position, and the vertical axis is the wavelet frequency. $\endgroup$
    – Wang Yun
    Commented Sep 14, 2022 at 13:47
  • $\begingroup$ I figured! just a very general comment ;) $\endgroup$
    – Jdip
    Commented Sep 14, 2022 at 14:07
  • $\begingroup$ @Jdip Do you have any comments on the problem I've encountered? $\endgroup$
    – Wang Yun
    Commented Sep 15, 2022 at 1:23
  • $\begingroup$ FWIW the bounty wasn't automatically awarded. I get I've not fully answered your question but such detailed debugging is beyond this site in the first place. If a post is at least somewhat helpful and there aren't others, consider awarding manually. $\endgroup$
    – OverLordGoldDragon
    Commented Sep 25, 2022 at 8:25

1 Answer 1

Reset to default
2
$\begingroup$
  1. Plot the frequency response of a bad wavelet. CWT is convolution, nothing can change that, hence there exists a pf_scale such that out_scale = ifft(fft(x) * pf_scale). Also mind the sampling grid.
  2. Prefer the unit impulse for debugging: x = zeros(N); x[N//2] = 1. With zero-padding and no unpadding (or just no padding), it will yield ifft(pf_scale) for all scale.
  3. Test linearity: cwt(x1) + cwt(x2) == cwt(x1 + x2). If this fails, that's a critical flaw, and the problem is beyond wavelets: the algo isn't doing convolution. Pick x1=x2=zeros(N); x1[N//4]=1; x2[-N//4]=1.
  4. Code in Python first so it's easier to debug interactively.

A plot of pf_scale will inevitably reveal problematic behavior, which then must be traced through the steps involved in computing it. Isolate the code for wavelet generation and plot the wavelet at every step:

freqs = # ...
plot(freqs)
wav0 = func0(freqs, param0, param1)
plot(wav0)
wav1 = func1(wav0, param2)
plot(wav1)
# ...
plot(freqs)
# ...
pf_scale = func9(wav8)

then compare against a well-behaved pf_scale. Pay attention to how scale affects each step, especially non-smooth and precision-sensitive functions. I see BjCoeff divides by scales.

It's possible (and appears to be the case) that the algo doesn't use convolution in a direct form. Again, your task is to rewrite the algorithm until it does. So if the wavelet used by the algo is pf_scale_algo, then rewrite and track effects on frequency, pf_scale_algo -> pf_scale_revised0 -> pf_scale_revised1 -> ... until -> pf_scale. Worst case, visual comparison should aid considerably.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.