PyWavelets / pywt

PyWavelets - Wavelet Transforms in Python
http://pywavelets.readthedocs.org
MIT License
2.04k stars 470 forks source link

About waveletpacket reconstruct #561

Open SeventeenChen opened 4 years ago

SeventeenChen commented 4 years ago

Hello, my friends,

Why did I get the wrong data length from waveletpacket reconstruct?

For example:

x = [1, 2, 3, 4, 5, 6, 7, 8]
wp = pywt.WaveletPacket(data=x, wavelet='db2', mode='symmetric', maxlevel=3)
new_wp = pywt.WaveletPacket(data=None, wavelet='db2', mode='symmetric')
new_wp['aaa'] = wp['aaa'].data
print(new_wp.reconstruct(update=True))

the length of new_wp is not equal to x

grlee77 commented 4 years ago

This is a bit tricky, but let me try to explain what is going on. I think this is a subtlety that we have not documented well.

The first point is that for the DWT or wavelet packet transforms with all boundary modes other than 'periodization' some additional boundary coefficients must be retained to enable perfect reconstruction. During the inverse transform it is necessary to know what the original signal size was in order to properly trim any excess coefficients upon reconstruction. You might think we could guess, but it is actually not unique as the following short example shows:

import numpy as np
import pywt

c1 = pywt.dwt(np.random.randn(9), wavelet='sym4', mode='symmetric')
print(c1[0].size)
c2 = pywt.dwt(np.random.randn(10), wavelet='sym4', mode='symmetric')
print(c2[0].size)

Both length 9 or length 10 input to DWT gives 8 approximation coefficients. So, just from the size of the coefficients we cannot guess which was the original signal length.

In your wp wavelet packet object, when you created from an existing set of data of length 8, it actually stored that length in a private attribute ._data_shape. You can inspect it via:

print(wp._data_shape)`

and you will see (8,). This allows that transform to automatically trim any excess coefficients back to the original length on reconstruction. So if you call wp.reconstruct(update=True) you will get a length 8 signal.

In the second case, you created a WaveletPacket with data=None, so this shape information was unavailable and the transform kept some extra coefficients on reconstruction that it thought might be needed. This is due to the ambiguous original signal length issue discussed above. There are currently two ways to resolve this in your example:

1.) initialize new_wp using data=np.zeros(len(x)). Then your last call to reconstruct will give you a length 8 signal. 2.) Alternatively you can keep only the first 8 entries in the array returned by new_wp.reconstruct

SeventeenChen commented 4 years ago

Dear friends, 

Thanks for your detailed and patient reply! I know how to solve my problem and why I get more coefficients than the original signal length. 

Thank you very much!!! 

Best wishes,

seventeen 

seventeen <675284835@qq.com>; <cherry_c.c@outlook.com>

 

------------------ 原始邮件 ------------------ 发件人: "PyWavelets/pywt" <notifications@github.com>; 发送时间: 2020年8月14日(星期五) 凌晨2:56 收件人: "PyWavelets/pywt"<pywt@noreply.github.com>; 抄送: "陈曦"<675284835@qq.com>;"Author"<author@noreply.github.com>; 主题: Re: [PyWavelets/pywt] About waveletpacket reconstruct (#561)

This is a bit tricky, but let me try to explain what is going on. I think this is a subtlety that we have not documented well.

The first point is that for the DWT or wavelet packet transforms with all boundary modes other than 'periodization' some additional boundary coefficients must be retained to enable perfect reconstruction. During the inverse transform it is necessary to know what the original signal size was in order to properly trim any excess coefficients upon reconstruction. You might think we could guess, but it is actually not unique as the following short example shows: import numpy as np import pywt c1 = pywt.dwt(np.random.randn(9), wavelet='sym4', mode='symmetric') print(c1[0].size) c2 = pywt.dwt(np.random.randn(10), wavelet='sym4', mode='symmetric') print(c2[0].size)

Both length 9 or length 10 input to DWT gives 8 approximation coefficients. So, just from the size of the coefficients we cannot guess which was the original signal length.

In your wp wavelet packet object, when you created from an existing set of data of length 8, it actually stored that length in a private attribute ._data_shape. You can inspect it via: print(wp._data_shape)`
and you will see (8,). This allows that transform to automatically trim any excess coefficients back to the original length on reconstruction. So if you call wp.reconstruct(update=True) you will get a length 8 signal.

In the second case, you created a WaveletPacket with data=None, so this shape information was unavailable and the transform kept some extra coefficients on reconstruction that it thought might be needed. This is due to the ambiguous original signal length issue discussed above. There are currently two ways to resolve this in your example:

1.) initialize new_wp using data=np.zeros(len(x)). Then your last call to reconstruct will give you a length 8 signal. 2.) Alternatively you can keep only the first 8 entries in the array returned by new_wp.reconstruct

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

Alejandro-1996 commented 3 years ago

Hello, and thank you for the great explanation. I just have another related question. I have an image of size 512512, and many of the implemented wavelet families give larger coefficients (e.g. 258258 or 280*280. I am using pywt.dwt2()).

If I wanted to represent the wavelet transform in a traditional way (stacking the coefficients to form an image of the same size of the original, with the LL coeff in the upper left corner and the HH in the bottom right), what would be consider best practice? To form my image with a larger size? Or to trim/reshape the coefficients? If trimming, which coefficients to discard?

Thank you very much, Alejandro Noguerón

grlee77 commented 3 years ago

On Sat, Oct 10, 2020 at 9:52 AM Alejandro-1996 notifications@github.com wrote:

Hello, and thank you for the great explanation. I just have another related question. I have an image of size 512512, and many of the implemented wavelet families give larger coefficients (e.g. 258258 or 280*280. I am using pywt.dwt2()).

If I wanted to represent the wavelet transform in a traditional way (stacking the coefficients to form an image of the same size of the original, with the LL coeff in the upper left corner and the HH in the bottom right), what would be consider best practice? To form my image with a larger size? Or to trim/reshape the coefficients? If trimming, which coefficients to discard?

Hi Alejandro,

If you run the transform with wavedec, wavedec2 or wavedecn then there is already a utility called coeffs_to_array that will stack the coefficients. For example:

image = pywt.data.camera() coeffs = pywt.wavedec2(image, wavelet='sym4', mode='symmetric', level=4) stacked, coeff_slices = pywt.coeffs_to_array(coeffs)

In this case, mode='symmetric' will result in some extra coefficients, so in order to stack into a contiguous array, padding with zeros as needed is performed. This array will be slightly larger than the original image, though. For the example above, the original is (512, 512) but the stacked version is (538, 538). To get identical size to the original you would have to use mode='periodization' (with a signal that is a multiple of 2**level in length along each axis).

There is not currently an option to discard any boundary coefficients when stacking, but if you wanted to implement that, you should discard the samples from from the end, not the start. The issue with discarding coefficients is that you would no longer be able to do the reverse operation (array_to_coeffs followed by waverecn) to get back the original signal.

Thank you very much, Alejandro Noguerón

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/PyWavelets/pywt/issues/561#issuecomment-706552383, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABRZ7PKJFXPBRNHLQK73FQDSKBRJVANCNFSM4P4I3W4Q .

kelleyjbrady commented 3 years ago

Hello,

Thank you for the for the explanation. I have been searching all day for a discussion like this. How would one resolve this issue when the WaveletPacket class is not being used?

I am facing this issue when I decompose a 1D signal using wavedec, then reconstruct with wavrec. At least in the docs wacerec does not have a data argument that can be used to set the desired output size. Can I simply throw away the extra entry in my array generated by waverec as you suggest? Does one always pick the final entry in the reconstructed array as the one to throw out? I would really love to see a reference for that, please forgive my ignorance.

grlee77 commented 3 years ago

Yes, the coefficients to discard will all be at the end. You can verify by tacking the round trip transform with random data to see that it matches (to within some floating point precision)

kelleyjbrady commented 3 years ago

To clarify for others who may read this, you mean to reconstruct the signal, then discard the final entry in the reconstructed signal, right? You referred to the items to drop as coefficients in your reply, so I want to be sure you are not saying to drop the final wavelet coefficient, then reconstruct the signal with (n-1) coefficients, where n is the number of coefficients defined in the wavedec used to decompose the signal.

grlee77 commented 3 years ago

Yes, I meant drop the last samples of the reconstructed signal

rgommers commented 2 years ago

This is a bit tricky, but let me try to explain what is going on. I think this is a subtlety that we have not documented well.

The question was answered, but let's leave this open as a documentation task.