0

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)

Eitan T
  • 32,660
  • 14
  • 72
  • 109
Jorge Zapata
  • 2,316
  • 1
  • 30
  • 57

3 Answers3

4

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
klurie
  • 1,060
  • 8
  • 16
  • I also thought that H=3 meant the 3rd wavelet of Dauchebies but in fact they use Dauchebies-1 (in simulation chapter) which I think is db1 in matlab, isn't it? – Jorge Zapata Mar 27 '13 at 03:35
  • 1
    +1: Good answer, but MATLAB has `wavedec` for a multilevel DWT. – Eitan T Mar 27 '13 at 13:32
  • It's not obvious from the paper that an undecimated dwt is used -- what is called dwt in the paper is most likely to be (in matlab terms) an fwt. This [paper](ftp://216-58-54-46.cpe.distributel.net/AiDisk_a1/MDM/Skule/4th_Year/FALL/ECE462/Projects/References/%5B7%5DA%20Novel%20Synchronization%20Invarient%20Audio%20Watermarking%20Scheme%20based%20on%20DWT%20and%20DCT.pdf) decribes a similar algorithm which does decimate the higher-level decompositions. – alle_meije Mar 27 '13 at 13:42
  • Quite nice explanation, now it's clear for me. Thanks for your time and answer1 – Jorge Zapata Apr 01 '13 at 15:48
3

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');
Eitan T
  • 32,660
  • 14
  • 72
  • 109
1

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

alle_meije
  • 2,424
  • 1
  • 19
  • 40