# Why is inverse CWT inexact / inaccurate?

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

i_plot = -1

all_close = np.linalg.norm(x-x_inv) < eps * np.linalg.norm(x)

if (all_close):
print('Yay, found one ')

if Do_plot:
plt.plot(x,color='b')
plt.plot(x_inv,color='r')

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

plt.close()
$$$$

• Were the docs helpful? help(icwt) under References: "One integral inverse CWT. dsp.stackexchange.com/a/71274/50076", and "Inverse CWT derivation. dsp.stackexchange.com/a/71148/50076" Commented Mar 17, 2023 at 11:03
• uh oh I should update that first link to instead be dsp.stackexchange.com/a/76239/50076 ; (1)-(4) there on top is all you need to debug this, I may look at your specific case later Commented Mar 17, 2023 at 11:07
• Thanks @OverLordGoldDragon for taking my question :-). Yes, I've come across those links. Let me give them a more thorough second read. Commented Mar 17, 2023 at 11:16
• Update : Reading the docs, I see that the defaults for the padtype arg of ssqueezepy.cwt and ssqueezepy.icwt is not the same. Setting it manually to be the same does not make icwt the inverse of cwt however :-( Commented Mar 17, 2023 at 11:39
• You're right on padtype, I'll note it, thanks. Note it doesn't affect one_int=True. Commented Mar 17, 2023 at 13:12

### Why iCWT is inexact

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"')

#%%
xf = rfft(x)
xf[1:3] = 0
`