Daubechies-1 is the Haar wavelet. It is a combination of the low-pass filter h
and a high-pass filter g
>> s = sqrt(2);
>> h = [1 1] / s;
>> g = [1 -1] / s;
To see the dwt in action you can make a signal
>> x = (1:64) / 64;
>> y = humps(x) - humps(0);
and as said by @kl3755 you need to apply it 3 times:
- every iteration, the dwt returns
a lowpass filtered signal (approximation)
a highpass filtered signal (detail)
- every next iteration is applied to the previous approximation
The term dwt
is ambiguous though -- it is often used for the fast wavelet transform fwt
, which downsamples approximation and detail by a factor 2 at each iteration. As it is the most used version, let's do the FWT here:
>> c1 = filter (h, s, y); % level 1 approximation
>> d1 = filter (g, s, y); % level 1 detail
>> c1 = c1 (2:2:end); d1 = d1 (2:2:end); % downsample
>> c2 = filter (h, s, c1); % level 2 approximation
>> d2 = filter (g, s, c1); % level 2 detail
>> c2 = c2 (2:2:end); d2 = d2 (2:2:end); % downsample
>> c3 = filter (h, s, c2); % level 3 approximation
>> d3 = filter (g, s, c2); % level 3 detail
>> c3 = c3 (2:2:end); d3 = d3 (2:2:end); % downsample
It's easy to see how you would program this recursion.
The fwt
output typically only uses the final approximation (c3) and the detail signals:
>> fwt_y_3 = [c3 d3 d2 d1];
The 'magic' of the fwt
representation is that you can reconstruct the original from just that, by filtering and upsampling in the same way as above, after reversing the filters:
>> g=reverse(g); % h is symmetric
>> d3 (2,:) = 0; d3 = d3 (:)'; % upsample d3
>> c3 (2,:) = 0; c3 = c3 (:)'; % upsample c3
>> r2 = filter (h, 1/s, c3) + filter (g, 1/s, d3); % level 2 reconstruction
>> d2 (2,:) = 0; d2 = d2 (:)'; % upsample d2
>> r2 (2,:) = 0; r2 = r2 (:)'; % upsample r2
>> r1 = filter (h, 1/s, r2) + filter (g, 1/s, d2); % level 1 reconstruction
>> d1 (2,:) = 0; d1 = d1 (:)'; % upsample d1
>> r1 (2,:) = 0; r1 = r1 (:)'; % upsample r1
>> ry = filter (h, 1/s, r1) + filter (g, 1/s, d1); % reconstruction of y
check if all is as it should be:
>> subplot(2,2,1);plot([y' 80+ry']);
>> subplot(2,2,2);plot([c1' 80+r1(1:2:end)']);
>> subplot(2,2,3);plot([c2' 80+r2(1:2:end)']);
>> subplot(2,2,4);plot(fwt_y_3);hold on;plot(80+c3(1:2:end));hold off
The name dwt
may refer to different versions of the undecimated wavelet transform. The fwt
is much faster and not redundant, but its main drawback is it's not shift invariant: you cannot shift y by reconstructing a shifted fwt
of y. Undecimated wavelet transforms are shift-invariant:
1. The continuous wavelet transform cwt
is as above, but without down- and upsampling
see dl.acm.org/citation.cfm?id=603242.603246
2. The 'cycle spinning' or 'à trous' transform does the subsampling, but performs it for all possible shifts at each level.
see dl.acm.org/citation.cfm?id=1746851