Where to find Inverse Continuous Wavelet transform (icwt)

[yanpanlau]

Hi, May I ask how to compute Inverse Continuous Wavelet transform (icwt)? I checked the documentations but cannot find the function. The corresponding Matlab function are the following:

https://www.mathworks.com/help/wavelet/ref/icwt.html

Thanks,

Ben

[rgommers]

Yes, that's not implemented yet. I don't think there's a plan to do so - @holgern or @grlee77 please correct me if I'm wrong.

[bloyl]

Any update on this? I would be very interested in an icwt?

[vbelz]

Do you plan to implement the icwt?

[wallneradam]

This library has .reconstruction() method: https://github.com/aaren/wavelets

Maybe someone could merge the 2? It would be good to have one complete wavelet library for python.

[Saudreau]

First of all, thanks for the pywt() library, this is very useful !

In parallel to FFT analysis, I’m trying to use pywt (both cwt() and dwt()) to analyze and reconstruct a signal. To understand some features implemented in pywt, I would like to reconstruct a signal according to a given scale provided by pywt.cwt().

I went through the .reconstruction() method: https://github.com/aaren/wavelets, but it’s not really straightforward … at least for me !!

So any update on a future icwt() method ?

All the very best

[grlee77]

I think the implementation in @aaren's repo follows the description in the paper by Torrence and Compo. There is a link to the text here. See the section title "Reconstruction". The following two lines appear to correspond to Eq. 11 in that paper:

https://github.com/aaren/wavelets/blob/a213d7c307e0a6b5a40cc648d4e3f69c2e77c12b/wavelets/transform.py#L430-L431

[Saudreau]

Thanks grlee77 for your reply... I will have a look at the Torrence and Compo paper and see if I can reconstruct my signal.

All the very best

[oseibrefo]

is there any update on the inverse cwt implementation?

[jspaezp]

Well... This chunk seems to work for me ...

If I have a singa that i transform like this

my_singal = ...
scales = np.array([2 ** x for x in range(7)])
coef, freqs = pywt.cwt(
    data=my_singal, scales=scales, wavelet="morl"
)

I can reconstruct it using the following ...

mwf = pywt.ContinuousWavelet('morl').wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (1 / y_0)

Original:

Reconstructed:

Overlaid section:

I feel like this would be a great addition to the library but I understand if it is not in the interest of the people who actually have to/have volunteer to maintain it.

[egaznep]

Well... This chunk seems to work for me ...

If I have a singa that i transform like this

my_singal = ...
scales = np.array([2 ** x for x in range(7)])
coef, freqs = pywt.cwt(
    data=my_singal, scales=scales, wavelet="morl"
)

I can reconstruct it using the following ...

mwf = pywt.ContinuousWavelet('morl').wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (1 / y_0)

Original:

Reconstructed:

Overlaid section:

I feel like this would be a great addition to the library but I understand if it is not in the interest of the people who actually have to/have volunteer to maintain it.

Thank you for sharing this. It indeed works with morlets, however scale is off with Mexican-hat wavelets. I am trying to find a solution and I'll edit my reply in case I succeed.

I respect the decision of not adding something like this to the library, however I cannot say that I understand the reasoning behind. Neither SciPy, nor PyWavelets nor few others have inbuilt inverse continuous wavelet transform, and if I cannot reconstruct, then what is this all been about, what are we working toward...

[bingjo]

Well... This chunk seems to work for me ... If I have a singa that i transform like this

my_singal = ...
scales = np.array([2 ** x for x in range(7)])
coef, freqs = pywt.cwt(
    data=my_singal, scales=scales, wavelet="morl"
)

I can reconstruct it using the following ...

mwf = pywt.ContinuousWavelet('morl').wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (1 / y_0)

Original: Reconstructed: Overlaid section: I feel like this would be a great addition to the library but I understand if it is not in the interest of the people who actually have to/have volunteer to maintain it.

Thank you for sharing this. It indeed works with morlets, however scale is off with Mexican-hat wavelets. I am trying to find a solution and I'll edit my reply in case I succeed.

I respect the decision of not adding something like this to the library, however I cannot say that I understand the reasoning behind. Neither SciPy, nor PyWavelets nor few others have inbuilt inverse continuous wavelet transform, and if I cannot reconstruct, then what is this all been about, what are we working toward...

Did you find a solution?

[rgommers]

respect the decision of not adding something like this to the library, however I cannot say that I understand the reasoning behind. Neither SciPy, nor PyWavelets nor few others have inbuilt inverse continuous wavelet transform, and if I cannot reconstruct, then what is this all been about, what are we working toward...

I don't think anyone said that? It seems like a good thing to add. Let me give some context here. At the moment the only two somewhat-active maintainers are @grlee77 and myself. We're both quite busy though, given we're also maintainer things like scikit-image, SciPy and NumPy - so PyWavelets is not high on our agendas. PyWavelets is quite stable with not many bugs, and we have usually found the time to review and merge PRs. And @grlee77 has been adding features as he needs them. I'm no longer a user, maintaining this package is basically public service:) Adding continuous wavelets was done by other previous maintainers mostly, and icwt never made it.

tl;dr if someone submits a clean implementation, we will review and get it merged.

@grlee77 please correct me if I missed something here.

[prevayTST]

Well... This chunk seems to work for me ... If I have a singa that i transform like this

my_singal = ...
scales = np.array([2 ** x for x in range(7)])
coef, freqs = pywt.cwt(
    data=my_singal, scales=scales, wavelet="morl"
)

I can reconstruct it using the following ...

mwf = pywt.ContinuousWavelet('morl').wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (1 / y_0)

Original: Reconstructed: Overlaid section: I feel like this would be a great addition to the library but I understand if it is not in the interest of the people who actually have to/have volunteer to maintain it.

Thank you for sharing this. It indeed works with morlets, however scale is off with Mexican-hat wavelets. I am trying to find a solution and I'll edit my reply in case I succeed. I respect the decision of not adding something like this to the library, however I cannot say that I understand the reasoning behind. Neither SciPy, nor PyWavelets nor few others have inbuilt inverse continuous wavelet transform, and if I cannot reconstruct, then what is this all been about, what are we working toward...

Did you find a solution?

I believe that the approach to reconstructing the signal that @jspaezp provided is a special case that works for the chosen wavelet and scales.

The .reconstruction() method from the wavelets repo provides a more generic approach:

The reconstructed time series is found as the sum of the
real part of the wavelet transform over all scales,

x_n = (dj * dt^(1/2)) / (C_d * Y_0(0)) \
        * Sum_(j=0)^J { Re(W_n(s_j)) / s_j^(1/2) }
where the factor C_d comes from the reconstruction of a delta
function from its wavelet transform using the wavelet
function Y_0. This C_d is a constant for each wavelet
function.

So the approach provided by @jspaezp assumes $dj = dt = C_d = 1$

Here is what Torrence and Compo have to say about the quantities $dj$ and $C_d$ (or in their notation $\delta_j$ and $C_\delta$):

It is convenient to write the scales as fractional powers of two: $$s_j = s_02^{j\delta j}, \quad j = 0, 1, \dots, J$$ where $s_0$ is the smallest resolvable scale and $J$ determines the largest scale.

Therefore, we have $\delta j = \frac{\log s_j - \log s_0}{\log 2^j}$.

For $C_d$ they write the following:

To derive $C_\delta$ for a new wavelet function, first assume a time series a $\delta$ function at time $n=0$, given by $x_n = \delta_{n0}$. This time series has a Fourier transform $\hat{x}_k = N^{-1}$, constant for all $k$. Substituting $\hat{x}_k$ into [the equation for the wavelet transform express in terms of the wavelet's Fourier transform] at time $n=0$ (the peak), the wavelet transform becomes:

$$W_\delta (s) = \frac{1}{N}\sum_{k=0}^{N-1}\hat{\psi}*(s\omega_k).$$

The reconstruction [equation] then gives

$$C_\delta = \frac{\delta j \delta t^{1/2}}{\psi_0(0)}\sum_{j=0}^J\frac{\mathfrak{R} \big[ W_\delta (s_j) \big] }{s_j^{1/2}}$$

I believe this gives all the information we need to define a general reconstruction method.

Using a similar setup as before:

wt = ...
max_j = ...
my_signal = ...
scales = np.array([2 ** j for j in range(max_j)])
my_coefs, my_freqs = pywt.cwt(
    data=my_signal, scales=scales, wavelet=wt
)

we reconstruct the signal like this (I'm assuming $dt=1$ here):

mwf = pywt.ContinuousWavelet(wt).wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

def W_delta(s, n, wf, freqs):
    psi_hat_star = np.conjugate(np.fft.fft(mwf[0]))
    return np.sum([psi_hat_star[np.argmin(
        np.abs(mwf[0] - s * pywt.frequency2scale(freqs[k]))
        )] for k in range(n)]) / n

def dj(s_arr):
    return ((log(s_arr[1]) - log(s_arr[0])) / log(2))

kwargs = dict(
    n=len(scales),
    wf=mwf,
    freqs=my_freqs
)

d_j = dj(scales)
C_d = d_j / y_0 * np.sum([np.real(W_delta(s, **kwargs))/sqrt(s) for s in scales])

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (d_j / (C_d * y_0))

# finally add the mean of the data back on ('reconstructed' is an anomaly time series)
reconstructed += my_signal.mean()

For what it's worth, I haven't tested this yet, but I will try it out and if this works, I'll open up a pull request.

[bigleaffrog]

Well... This chunk seems to work for me ...

If I have a singa that i transform like this

my_singal = ...
scales = np.array([2 ** x for x in range(7)])
coef, freqs = pywt.cwt(
    data=my_singal, scales=scales, wavelet="morl"
)

I can reconstruct it using the following ...

mwf = pywt.ContinuousWavelet('morl').wavefun()
y_0 = mwf[0][np.argmin(np.abs(mwf[1]))]

r_sum = np.transpose(np.sum(np.transpose(coef)/ scales ** 0.5, axis=-1))
reconstructed = r_sum * (1 / y_0)

Original:

Reconstructed:

Overlaid section:

I feel like this would be a great addition to the library but I understand if it is not in the interest of the people who actually have to/have volunteer to maintain it.

I can't use the method to rebuild my signal, what's going on? Can you give a more complete code?

[Yecaergi]

Hello! I need a code for icwt that reconstruct exactly my signal. Any help? Im loose