About decimation in DWT

Hi! I'm studying DWT by this link http://slideplayer.com/slide/4998644/. And I have 2 question (slide on 18.07):

1. We should do decimation anyway: after LP-filter and after HP-filter. After LP-filter it's ok, but after HP-filter we will have an aliasing? Am I right? For ex: We have a signal with f1=100 Hz and f2=400 Hz. Fs=1000 Hz. After LP (0-250 Hz) we have subband with f1=100 Hz, Fs= 500. It's ok (f1<Fs/2). After HP (250-500 Hz) we have subband with f2=400 Hz, Fs= 500. It's not ok. Aliasing (f2>Fs/2)?

2. How we can use only 2 filters in DWT? I don't understand, why the filters change their cut-off frequency when we use decimated signal? For ex: "Therefore, if you calculated, say, a low-pass filter with a cutoff frequency of 10 MHz at a sampling frequency of 100 MHz, and then samples taken at a frequency of 10 GHz were applied to the input of this filter, the filter cutoff frequency would be 100 times larger, i.e., 1 GHz ."

Answer 1

[Wayne King]

Hi Alexander,

1. Yes, you are correct that the downsampling causes aliasing (or can cause aliasing). However, the analysis filters and synthesis filters in the wavelet transform are designed in such a way that any aliasing is canceled upon reconstruction.
2. You are only using two filters because you are using the so-called pyramid algorithm. You can think of the equivalent filters in terms of successive convolutions, but you should keep in the mind the so called Noble identities from multirate signal processing. That allows you to see the effect of interchanging downsampling with filtering (or on the synthesis, the effect of interchanging upsampling with filtering). Read up on the Noble identities, for one example, if you lowpass filter with the scaling filter, G(f), and then downsample by two followed by filtering with the wavelet filter, H(f) that is equivalent to the second-level wavelet coefficients (prior to decimating those by two). Schematically, that is X(f) -> G(f) -> downsample by two -> H(f). The noble identities will tell you that is equivalent to X(f) -> G(f)H(2f) -> downsample by two so the equivalent filter for the second-level wavelet coefficients is G(f)H(2f). Now let's look at that filter in MATLAB.

N = 64;
[G,H] = wfilters('sym4');
Gdft = fft(G,N);
Hdft = fft(H,N);
H2 = Hdft(1+mod(2*(0:N-1),N));
H2 = Gdft.*H2;

Now let's plot Hdft (1st-level wavelet filter response) and H2 (second-level wavelet filter response)

   f = 0:1/N:1/2;
     plot(f,abs(Hdft(1:N/2+1)));
     hold on;
     plot(f,abs(H2(1:N/2+1)),'r');
     legend('First Level Wavelet','Second Level Wavelet');
     xlim([0 1/2]);
     grid on;
     xlabel('Cycles/Sample');

[Alexander Voznesensky]

Thank you! I was trying to understand what you said. Here is the plot. I see, that we have a scaling effect. For ex, f=100 Hz (after deimation) corresponds to f=200 Hz (before decimation). I supose, this is the reason of using only 2 filters. We should multiply new frequences with factor M*level_number and then match them to original f-axis where we have filter's amplitude-frequency response. Formally this is understandable, but... Looks like a trick. Could you recommend some stuff about it?

Answer 2

[Wayne King]

Hi Alexander,

I don't think there is a trick here. The filter responses scale as you go down in resolution but that is because the discrete wavelet transform uses L2 normalization so that the L2 norm is preserved. That is good for lots of applications, but not for others. For example, in the CWT we (and many others) believe that L1 normalization is better and that is why the CWT in MATLAB from 16b on uses L1 normalization.

One thing to keep in mind is that the downsampling at each level turns essentially half-band signals into full band signals, so we can say that (this is an approximation) the scaling coefficients at level 1 capture aspects of the data in the frequency interval (-1/4,1/4) while the wavelet coefficients capture information in (-1/2 -1/4) union (1/4,1/2). After downsampling both of these become full-band signals so they both have content over (-1/2, 1/2) but while the scaling coefficients maintain the frequency order, the wavelet coefficients actually map the frequencies in (1/4, 1/2) in reverse order. So a frequency in (1/4,1/2) gets mapped essentially to 2*(1/2-f) where f is the original signal. For example:

% for the classic orthonormal DWT use a power of two and dwtmode 'per'
dwtmode('per');
n = 0:255;
x = cos(2*pi*7/16*n);  % frequency is 7/16 cycles/sample almost at Nyquist
f = -1/2:1/256:1/2-1/256;
plot(f,fftshift(abs(fft(x))))

If you zoom in, you see the frequency is 7/16 as expected. Now we expect that the wavelet coefficients at level one will map this to approximately 2*(1/2-7/16)

[A,D] = dwt(x,'sym4');
fnew = -1/2:1/128:1/2-1/128;
plot(fnew,fftshift(abs(fft(D))))

Now the frequency is 1/8 as expected.

I think this is cover in a lot detail in Don Percival's and Andrew Walden's textbook on wavelets in time series analysis

"Wavelet Methods for Time Series Analysis".