Why is inverse CWT inexact / inaccurate?

Question

I'm all new to wavelet analysis. I'm trying to get a working understanding of the continuous wavelet transform and its inverse. By "working understanding", I really mean "getting some code to execute for 1D signals and have an intuitive understanding of what the result means".

I'm using the ssqueezepy library in python, by https://dsp.stackexchange.com/users/50076/overlordgolddragon

This library provide cwt and icwt routines for the continuous wavelet transform and its inverse ... but the inverse process does not seem to get the original signal back.

More precisely, using the following code, I expected to see the same signal several times:

import numpy as np
import itertools
import ssqueezepy
import matplotlib.pyplot as plt

wavelet = 'gmw'
nv = 16

n = 100
t = np.linspace(0.,1.,n)
x = np.sin(10*t*t)

x_mean = x.mean()

l1_norm_choices = [True,False]
one_int_choices = [True,False]

plt.plot(x,color='b')

for (l1_norm,one_int) in itertools.product(l1_norm_choices,one_int_choices):


    cwtmatr, freqs = ssqueezepy.cwt(x,wavelet=wavelet,nv=nv,l1_norm=l1_norm)
    x_inv = ssqueezepy.icwt(cwtmatr, wavelet=wavelet, nv=nv, one_int=one_int, l1_norm=l1_norm, x_mean=x_mean)

    plt.plot(x_inv,color='r')

plt.show()

Instead, this is my result (original signal in blue, tentative reconstructions in red) :

Can someone help me

• understand why I don't retrieve the original signal with my code
• tweak the code so that I do retrieve the original signal

UPDATE: Here is a more thorough piece of code that searches for combinations of arguments that makes icwt the inverse of cwt. Found nothing.

import numpy as np
import itertools
import ssqueezepy
import matplotlib.pyplot as plt

import os

output_folder = os.path.dirname(__file__)
basename = os.path.splitext(os.path.basename(__file__))[0]

Do_plot = True
# Do_plot = False

plot_ext = [
    '.png',
    # '.eps',
    # '.pdf',
]

eps = 1e-7

nv = 256   # Does this matter ?
n = 256    
t = np.linspace(0.,1.,n)
x = np.sin(10*t*t+1) # or whatever signal

x_mean = x.mean()

wavelet_choices = list(ssqueezepy.wavs())
l1_norm_choices = [True,False]
one_int_choices = [True,False]
scales_choices = ['log','linear','log-piecewise']
padtype_choices = ['reflect', 'symmetric', 'replicate', 'wrap', 'zero' ]
# rpadded_choices = [True,False] # True fails an assertion check
rpadded_choices = [False]

i_plot = -1
for (l1_norm,one_int,scales,wavelet,padtype,rpadded) in itertools.product(l1_norm_choices,one_int_choices,scales_choices,wavelet_choices,padtype_choices,rpadded_choices):

    cwtmatr, freqs = ssqueezepy.cwt(x,wavelet=wavelet,nv=nv,l1_norm=l1_norm,scales=scales,padtype=padtype,rpadded=rpadded)
    x_inv = ssqueezepy.icwt(cwtmatr, wavelet=wavelet, nv=nv, one_int=one_int, l1_norm=l1_norm, x_mean=x_mean,scales=scales,padtype=padtype,rpadded=rpadded)

    all_close = np.linalg.norm(x-x_inv) < eps * np.linalg.norm(x) 
    
    if (all_close):
        print('Yay, found one ')
        print(f'l1_norm={l1_norm}, one_int={one_int}, scales={scales}, wavelet={wavelet}, padtype={padtype}, rpadded={rpadded}')

    if Do_plot:
        plt.plot(x,color='b')
        plt.plot(x_inv,color='r')
        plt.title(f'l1_norm={l1_norm}, one_int={one_int}, scales={scales}, wavelet={wavelet}, padtype={padtype}, rpadded={rpadded}')

        i_plot += 1
        for ext in plot_ext:
            plt.savefig(os.path.join(output_folder,f'{basename}_{i_plot}{ext}'))

        plt.close()
```

Answer 1

Why iCWT is inexact

Refer to top of this article.

In this case you've reproduced an edge case by using a very short input signal with frequency extrema, that amplifies imperfections in a CWT filterbank.

Debugging iCWT in practice

Two plots are of interest - the filterbank, and its "energy of signal transfer function"; ssqueezepy wavelets are analytic and zero-phase so we can just plot the real part (don't do this in general):

We see the signal is attenuated for first 3 bins, and near Nyquist. I won't explain these phenomena here, only how to address them. Next we inspect the spectrum of padded signal (original also works but less accurate):

It's heavy in low frequencies, where this filterbank misbehaves. We can address this by either designing a better filterbank (change scales or wavelet), or manually zeroing original signal's frequencies and eliminating padding - latter is a terrible idea but for sake of this question it works:

Note that accurate icwt isn't an option with scipy or PyWavelets, and both have serious flaws, as explained here.

Is exact iCWT possible?

Yes, but it's not desired to do this very close to DC or Nyquist per sacrificing desired feature qualities. That is, with direct one-integral inversion; we can still manually compensate for some filter shortcomings at inversion by making the EoS TF curve flat, but in feature quality sense this is cheating, also not implemented in ssqueezepy.

ssqueezepy also performs worse here than it needs to; my auto-scales scheme has flaws. I'll "soon" release a library with better CWT, you can be notified by "Watch"-ing ssqueezepy.

Tips for accurate icwt with ssqueezepy

• Not a "tip" but a must: pass same parameters to both cwt and icwt
• Use longer input
• scales='log:maximal' is safest (but still fails in this case)
• Inspect and manually extend the scales returned by cwt if the filterbank under-tiles frequencies
• Use higher nv, but the default is already quite high
• Pass in x_mean = x.mean() to icwt

Other configs

This answer covers one-integral inversion; two-integral inversion is problematic for numerous reasons and should be avoided unless familiar with CWT and iCWT.

---

Code

import numpy as np
from numpy.fft import rfft, irfft
import matplotlib.pyplot as plt
from ssqueezepy import cwt, icwt, Wavelet, padsignal
from ssqueezepy.visuals import plot, plotscat


wavelet = Wavelet('gmw')
nv = 16

n = 100
t = np.linspace(0.,1.,n)
x = np.sin(10*t*t)

x_mean = x.mean()

plt.plot(x,color='b')

Wx, scales = cwt(x, wavelet=wavelet, nv=nv)
x_inv = icwt(Wx, wavelet=wavelet, nv=nv, x_mean=x_mean)

plt.plot(x_inv,color='r')
plt.show()

#%%
psis = wavelet._Psih
pkw = dict(w=.6, h=.8)
plot(psis.T, color='tab:blue', show=1, **pkw, title="Filterbank, real part")

eos_tfn = np.sum(psis, axis=0)**2
plotscat(eos_tfn, show=1, **pkw, title='Energy of Signal "transfer function"')

xp = padsignal(x, 'reflect')
plotscat(rfft(xp), abs=1, show=1, **pkw, title="|rfft(x_padded)|")

#%%
xf = rfft(x)
xf[1:3] = 0
xadj = irfft(xf)
x_mean = xadj.mean()

Wx, scales = cwt(xadj, wavelet=wavelet, nv=nv, padtype=None)
x_inv = icwt(Wx, wavelet=wavelet, scales=scales, nv=nv, x_mean=x_mean)

plot(xadj, color='b')
plot(x_inv, color='r', title="xadj, xadj_inv | MSE=%.3e" % np.mean((xadj-x_inv)**2))
plt.plot(xadj, color='b')
plt.plot(x_inv, color='r')
plt.show()