Spectral energy with Wavelet Packet Analysis
Adrien Depaillat
I need to extract the spectral energy of a signal contained within given frequency subbands ([0.0033 - 0.04 Hz], [0.04 - 0.15 Hz], [0.15 - 0.40 Hz]). As the signal is only 30 seconds long, at a sampling rate of 16Hz, I chose to use Wavelet Packet Analysis.
However, to be able to isolate these subbands, I have to perform the analysis at level 12 at least. And I saw by computing the sum of all the squared coefficients of the "wpspectrum" that this total energy varies with the level of analysis.. Though i thought that the signal energy was not altered by the wavelet packet decomposition..
Anyway, I can't compute the energy in these subbands, and I have been trying to do that for a month now, so it would be just great if some one could rescue me, as the teachers of my school can't help me on this.
Thanks for your attention
Wayne King
Hi Adrien, the problem is that wpspectrum() repeats the spaced wavelet packet coefficients to create a smooth image of the wavelet packet spectrum. This is described in the Algorithms section on the reference page for wpspectrum. See
>>doc wpspectrum
If you want to decompose your signal and look at percentage of signal energy in the various terminal nodes, use wpdec() and wenergy(). For example:
x = randn(1024,1);
dwtmode('per');
T = wpdec(x,5,'haar');
E = wenergy(T);
sum(E)
Hope that helps,
Wayne
Adrien Depaillat
>
> >>doc wpspectrum
>
> If you want to decompose your signal and look at percentage of signal energy in the various terminal nodes, use wpdec() and wenergy(). For example:
>
> x = randn(1024,1);
> dwtmode('per');
> T = wpdec(x,5,'haar');
> E = wenergy(T);
> sum(E)
>
> Hope that helps,
> Wayne
Hi,
Thanks for such a quick answwer! Indeed this seems to be an appropriate solution, but as there isn't any explicit use of the sampling frequency in the use of wpdec, I don't know how to recognize the correspondance between the pseudo-frequency and the nodes of the tree..
Could you please tell me how to find out the pseudo-frequencies attributed to the Nth node, so I can choose the nodes to consider to compute the energy in the subbands mentionned precedently?
I thought it would be written in the tree, but it doesn't seem to be explicitly given (and I apologize if I am being wrong).
Thank you for your help,
Adrien Depaillat
Wayne King
Hi Adrien, you can use otnodes() to return the frequeny-ordered terminal nodes of the wavelet packet tree. Then realize at the J-th level, the wavelet packet transform divides the frequency interval from [0,Nyquist] into intervals of [n*Fs/2^{J+1}, (n+1)*Fs/2^{J+1}]
n=0,1,2,...2^J-1
Hope that helps,
Wayne
Adrien Depaillat
> n=0,1,2,...2^J-1
>
>
> Hope that helps,
> Wayne
Hi Wayne,
Thanks for letting me know about otnodes. If I get this right, Tn_Seq(i)=j means that the node corresponding to the frequency i*Fs/2^{J+1} is the terminal node number J as ordered in the WPTree..
I thus modified my script and tried it with 2 different sampling frequencies. If the total energy does not change over these tests (100 : I guess it is then a pourcentage), the energy computed for each subbands (cf first message) does significantly change! For example, in the lowest frequency band [0.0033 to 0.04 Hz], when I change Fs from 16 to 32Hz, the energy goes from about 2 to 4.17.. Whereas for the band [0.15 to 0.4 Hz], it stays at 3.6..
The signal computed is a pure sinusoid, with f = 4 Hz. The energy computed for the node containing the information (f=4Hz) is 17.8 with Fs=16Hz, and 5.7 with Fs=32 Hz.. I thought it would be more than that, for both cases given this is the real frequency of the signal, which a synthetised one (so it is a "perfect" one : no noise..)
Alose, I thought that provided Fs was high enough to respect Shannon (at least twice the maximum frequency contained in the signal), the analysis would give me the same repartition of energy in the different subbands for different reasonnable sampling frequencies. Am I wrong?
Thank for your help!
Adrien Depaillat
Wayne King
Hi Adrien, you should include code examples so people can reproduce what you're talking about. Changing the sampling frequency changes the width of the approximate octave bands.
t = linspace(0,1,64);
x = cos(2*pi*4*t);
y = cos(2*pi*20*t);
dwtmode('per','nodisplay');
wptx = wpdec(x,2,'db4');
wpty = wpdec(y,2,'db4');
[~,tn_Seq] = otnodes(wpt);
wenergy(wptx)
wenergy(wpty)
you see that the overwhelming percentage of energy in both cases is contained in one terminal node, the terminal node that contains the single frequency in the x and y waveforms. The terminal node ordering in frequency order is:
3 4 6 5
these nodes contain energy from the approximate octave bands (in Hz):
[0,8] [8,16],[16,24],[24,32]
If you look at wenergy(wptx) you get:
98.8065 1.1822 0.0067 0.0047
If you look at wenergy(wpty) you get:
1.4428 10.3655 1.4455 86.7461
Terminal node 6 (the last one) corresponds approximately to [16,24] Hz so we get exactly the result we expect in terms of where the largest percentage of energy is.
Hope that helps,
Wayne
Adrien Depaillat
>
> t = linspace(0,1,64);
> x = cos(2*pi*4*t);
> y = cos(2*pi*20*t);
> dwtmode('per','nodisplay');
> wptx = wpdec(x,2,'db4');
> wpty = wpdec(y,2,'db4');
> [~,tn_Seq] = otnodes(wpt);
> wenergy(wptx)
> wenergy(wpty)
>
> you see that the overwhelming percentage of energy in both cases is contained in one terminal node, the terminal node that contains the single frequency in the x and y waveforms. The terminal node ordering in frequency order is:
> 3 4 6 5
> these nodes contain energy from the approximate octave bands (in Hz):
> [0,8] [8,16],[16,24],[24,32]
>
> If you look at wenergy(wptx) you get:
> 98.8065 1.1822 0.0067 0.0047
>
> If you look at wenergy(wpty) you get:
> 1.4428 10.3655 1.4455 86.7461
>
> Terminal node 6 (the last one) corresponds approximately to [16,24] Hz so we get exactly the result we expect in terms of where the largest percentage of energy is.
>
> Hope that helps,
> Wayne
Hi Wayne,
I was finally able to solve my problem! Thanks for everything,
I really appreciated your help, which was more than useful.
Cordially,
Adrien Depaillat