pywt icon indicating copy to clipboard operation
pywt copied to clipboard

About waveletpacket reconstruct

Open SeventeenChen opened this issue 4 years ago • 9 comments

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

SeventeenChen avatar Aug 12 '20 07:08 SeventeenChen

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

grlee77 avatar Aug 13 '20 18:08 grlee77

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 <[email protected]>; <[email protected]>

 

------------------ 原始邮件 ------------------ 发件人: "PyWavelets/pywt" <[email protected]>; 发送时间: 2020年8月14日(星期五) 凌晨2:56 收件人: "PyWavelets/pywt"<[email protected]>; 抄送: "陈曦"<[email protected]>;"Author"<[email protected]>; 主题: 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.

SeventeenChen avatar Aug 14 '20 02:08 SeventeenChen

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

Alejandro-1996 avatar Oct 10 '20 13:10 Alejandro-1996

On Sat, Oct 10, 2020 at 9:52 AM Alejandro-1996 [email protected] 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 .

grlee77 avatar Oct 10 '20 15:10 grlee77

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.

kelleyjbrady avatar Jul 29 '21 19:07 kelleyjbrady

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)

grlee77 avatar Jul 29 '21 22:07 grlee77

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.

kelleyjbrady avatar Jul 29 '21 22:07 kelleyjbrady

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

grlee77 avatar Jul 30 '21 01:07 grlee77

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.

rgommers avatar Nov 06 '21 21:11 rgommers