1

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  = scalepos
    s3 = scale3
    p1 = s3np.exp(-(x/2.0)2)x(60.0-20.0*x2+x*4)
    p2 = 32.0np.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 = Dxfreqs2.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.0x2+15.0x4-5.0*x5)/60.0
    elif ((x>1.0)&(x<=2.0)):
        y = (51.0+75.0x-210.0x2+150.0*x3-45.0x4+5.0x5)/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.0x2+1680.0x4-320.0*x6)/11520.0
    elif ((x>0.5)&(x<=1.5)):
        y = (23583.0-420.0x-16380.0x2-5600.0*x3+15120.0x4-

            6720.0x5+960.0*x6)/46080.0
    elif ((x>1.5)&(x<=2.5)):
        y = (4137.0+30408.0x-59220.0x2+42560.0*x3-

            15120.0x4+2688.0x5-192.0*x6)/23040.0
    elif ((x>2.5)&(x<=3.5)):
        y = (7.0-2.0x)*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 - x2/3.0+x4/9.0-x6/36.0+x7/144.0
    elif ((x>1.0)&(x<=2.0)):
        y = (2472.0-392.0x-504.0x2-1960.0*x3+2520.0x4-

             1176.0x5+252.0*x6-21x7)/5040.0
    elif ((x>2.0)&(x<=3.0)):
        y = (-1112.0+12152.0x-19320.0x2+13720.0x3-5320.0*x4+

            1176.0x5-140.0x6+7.0*x7)/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    = -&lt;int&gt;( 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&gt;=0)&amp;(new_indx&lt;=N-1) ):
    Ak = &lt;int&gt;(ceil(new_indx-3))
    Bk = &lt;int&gt;(floor(new_indx+2))
    for k in range(Ak, Bk+1, 1):
        BSval = bspline4th(new_indx-&lt;double&gt;(k)-0.5)
        out   = out+gk[k+3]*BSval
    return out
else:
    sys.exit( &quot;Check the value of new_indx.&quot; )




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&gt;=0)&amp;(new_indx&lt;=N-1) ):
    Ak = &lt;int&gt;(ceil(new_indx-4))
    Bk = &lt;int&gt;(floor(new_indx+2))
    for k in range(Ak, Bk+1, 1):
        BSval = bspline5th(new_indx-&lt;double&gt;(k)-1.0)
        out   = out+gk[k+4]*BSval
    return out
else:
    sys.exit( &quot;Check the value of new_indx.&quot; )




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&gt;=0)&amp;(new_indx&lt;=N-1) ):
    Ak = &lt;int&gt;(ceil(new_indx-5))
    Bk = &lt;int&gt;(floor(new_indx+2))
    for k in range(Ak, Bk+1, 1):
        BSval = bspline6th(new_indx-&lt;double&gt;(k)-1.5)
        out   = out+gk[k+5]*BSval
    return out
else:
    sys.exit( &quot;Check the value of new_indx.&quot; )




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&gt;=0)&amp;(new_indx&lt;=N-1) ):
    Ak = &lt;int&gt;(ceil(new_indx-6))
    Bk = &lt;int&gt;(floor(new_indx+2))
    for k in range(Ak, Bk+1, 1):
        BSval = bspline7th(new_indx-&lt;double&gt;(k)-2.0)
        out   = out+gk[k+6]*BSval
    return out
else:
    sys.exit( &quot;Check the value of new_indx.&quot; )




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.5C2(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/&lt;double&gt;(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( &quot;The length of Bj should be equal to that of gamma&quot; )

# Antiderivatives of signal
gk1, gk2, gk3, gk4 = InterpCoeffs( signal )
C1 = AntiDerivative1st( N, gk1, &lt;double&gt;(N)-1.0 )
C2 = AntiDerivative2nd( N, gk2, &lt;double&gt;(N)-1.0 )
C3 = AntiDerivative3rd( N, gk3, &lt;double&gt;(N)-1.0 )
C4 = AntiDerivative4th( N, gk4, &lt;double&gt;(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     = &lt;double&gt;(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(2np.pi/L, Nnp.pi/L, num=256)
x = Dxnp.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?

Wang Yun
  • 124
  • 1
  • 13
  • 1
    Just one small comment, always nice to have units on the axes! – Jdip Sep 14 '22 at 13:09
  • @Jdip Thank you for your comment. The horizontal axis is the position, and the vertical axis is the wavelet frequency. – Wang Yun Sep 14 '22 at 13:47
  • I figured! just a very general comment ;) – Jdip Sep 14 '22 at 14:07
  • @Jdip Do you have any comments on the problem I've encountered? – Wang Yun Sep 15 '22 at 01:23
  • 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. – OverLordGoldDragon Sep 25 '22 at 08:25

1 Answers1

2
  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.

OverLordGoldDragon
  • 8,912
  • 5
  • 23
  • 74