Question

I'm trying to implement some watermarking algorithm found in a paper1. This is a line of the paper:

Do the H-level DWT for all the renumbered segments.

Then in simulation section the author explains the wavelet used for experiments.

DWT transform adopted the common wavelet "Daubechies-1" and level H = 3.

I just don't get what does H means, how do I input H=3 in matlab DWT function?

My actual code is:

[cA,cD] = dwt(audio,'db3');

Can someone help me?


1Ji, Y. & Kim, J. A Quantified Audio Watermarking Algorithm Based on DWT-DCT. Multimedia, Computer Graphics and Broadcasting 339–344 (2011)

Was it helpful?

Solution

1. Q: "What does H (level) mean?"

A: Wikipedia describes this concept nicely, but I'll attempt to summarize. For each level, the data (original data for level 1, otherwise approximation data from previous level) is decomposed into approximation and detail data. The result is coefficient which describe the data in different frequency bins.

2. Q: How do I input H=3 in matlab DWT function?

A: As you point out, they are using db1. To extract the correct coefficients for level H=3, we'll need to implement this cascading algorithm. Here is a rough sketch of the code.

nLevels = 3;

% Get the coefficients for level 1
[cA,cD{1}] = dwt(audio,'db1');

% Continue to cascade to get additional coefficients at each level
for n = 2:nLevels
   [cA,cD{n}] = dwt(cA,'db1');
end

% Final coefficients are cA from highest level and cD from each level

OTHER TIPS

Your question has been answered very nicely by kl3755, albeit the provided solution could be further improved. For a multilevel wavelet transform, instead of using the dwt command, use wavedec:

H = 3;
[cA, cD] = wavedec(audio, H, 'db1');

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

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top