Talk:Cohen-Daubechies-Feauveau wavelet

From Wikipedia, the free encyclopedia

[edit] Correct last digits and different normalizations

Hmm... I got her by clicking on "Discussion when I was Visitiong the Cohen-Daubechies-Feauveau wavelet page, and I get an error message from Wikipedia saying they dont have an article???

Hey men, the filter coeffs are incorrect; if inverted, it results in a 21/19-tap filter. The correct one with (1;1) normalization

analysis lowpass filter:	0.0267487574108098,	-0.0168641184428750,	-0.0782232665289879,	0.2668641184428723,	0.6029490182363579,	0.2668641184428723,	-0.0782232665289879,	-0.0168641184428750,	0.0267487574108098,
analysis highpass filter:	0.0456358815571247,	-0.0287717631142498,	-0.2956358815571235,	0.5575435262284970,	-0.2956358815571235,	-0.0287717631142498,	0.0456358815571247,
synthesis lowpass filter (with interleaved filtering):	0.0534975148216208,	-0.0912717631142514,	-0.1564465330579798,	0.5912717631142532,	1.2058980364727310,	0.5912717631142532,	-0.1564465330579798,	-0.0912717631142514,	0.0534975148216208,
synthesis highpass filter (with interleaved filtering):	0.0337282368857512,	-0.0575435262285022,	-0.5337282368857500,	1.1150870524570070,	-0.5337282368857500,	-0.0575435262285022,	0.0337282368857512,

--Frigo 10:39, 12 February 2006 (UTC)

I have included the correct version, with (1;1) normalization, and with inverses used in interleaved filtering (we interleave the lowpass and highpass coeffs, THEN we filter) but haven't removed the previous, despite that it is wrong; its inverse is a 21/19-tap filter (as I described above). Should it be removed? Frigo 18:51, 7 March 2006 (UTC)

Hi, Frigo, nice that You found Your way here. Sometimes, in the last year more often then acceptable, the wiki-servers have grave problems. Especially for people logged in and trying to edit. Without login, the read-only version is on different servers and should have better stability.

I will check Your concerns, the coeffient table was just pasted. Can You tell how You arrived at a 21/19-tap filter? Also, the importance of "interleaved version" needs explaination.--LutzL 08:39, 8 March 2006 (UTC)llehmann at wavelet org

Actually, those coefficients above are wrong, since the synthesis filters are just the analysis filters, multiplied by two. But 9/7 is not orthogonal. The coefficients in the old version of the article were right, and You repeated them, adding two wrong columns. Now they should be correct up to rounding and sufficiently explained to be above all doubts.--LutzL 15:54, 8 March 2006 (UTC)

Hi. The wrong version differed two coeff from the right version:

wrong:	0.0267487410809760,	-0.0168641184428750,	-0.0782232665289879,	0.2668641184428723,	0.6029490182363579,	0.2668641184428723,	-0.0782232665289879,	-0.0168641184428749,	0.0267487410809760,
right:	0.0267487574108098,	-0.0168641184428750,	-0.0782232665289879,	0.2668641184428723,	0.6029490182363579,	0.2668641184428723,	-0.0782232665289879,	-0.0168641184428750,	0.0267487574108098,

Also, I tried to invert it with gauss-jordan decomposition and with a extended cholesky decomposition (don't ask...), with matrix sizes of 64x64, 128x128, and if I remember correctly, 100x100. In every cases, the wrong version caused a "leakage", modifying the 0 coeffs of the synthesis filters. I became aware of this flaw, so I begin to search for the right coefficients of CDF 9/7. I found it in FBI's WSQ specification, with (sqrt(2),sqrt(2)) normalization (pdf file, page 49):

http://citeseer.ist.psu.edu/98590.html

I use this code, and it works 100% good with the filters I specified (the synthesis filters aren't the analysis ones multiplied with 2, I checked it.

// a_l_len=9,a_h_len=7,a_l_shift=0,a_h_shift=1

// s_l_len=9,s_h_len=7,s_l_shift=1,s_h_shift=0

void wavelet_forward_symmetric_symmetric(double input[],double temp[],int n,wavelet &w){ int x,y,h,i,j; double s;

for(x=0;x<(n>>1);x++){

h=(x<<1)+w.a_l_shift;
s=w.a_l[w.a_l_len>>1]*input[h];
for(y=1;y<=(w.a_l_len>>1);y++){
 i=abs(h+y);
 j=abs(h-y);
 if(i>=(n<<1)-2)i%=(n<<1)-2;
 if(j>=(n<<1)-2)j%=(n<<1)-2;
 if(i>=n)i=(n<<1)-2-i;
 if(j>=n)j=(n<<1)-2-j;
 s+=w.a_l[(w.a_l_len>>1)+y]*(input[i]+input[j]);};
temp[x]=s;};

for(x=0;x<(n>>1);x++){

h=(x<<1)+w.a_h_shift;
s=w.a_h[w.a_h_len>>1]*input[h];
for(y=1;y<=(w.a_h_len>>1);y++){
 i=abs(h+y);
 j=abs(h-y);
 if(i>=(n<<1)-2)i%=(n<<1)-2;
 if(j>=(n<<1)-2)j%=(n<<1)-2;
 if(i>=n)i=(n<<1)-2-i;
 if(j>=n)j=(n<<1)-2-j;
 s+=w.a_h[(w.a_h_len>>1)+y]*(input[i]+input[j]);};
temp[x+(n>>1)]=s;};

for(x=0;x<n;x++)input[x]=temp[x]; };

void wavelet_backward_symmetric_symmetric(double input[],double temp[],int n,wavelet &w){ int x,y,h,i,j; double s;

for(x=0;x<(n>>1);x++){

temp[(x<<1)+w.a_l_shift]=input[x];          // was 0. now a_l_shift. is this good?
temp[(x<<1)+w.a_h_shift]=input[x+(n>>1)];}; // was 1. now a_h_shift. is this good?

for(x=0;x<(n>>1);x++){

h=(x<<1)+w.s_l_shift;
s=w.s_l[w.s_l_len>>1]*temp[h];
for(y=1;y<=(w.s_l_len>>1);y++){
 i=abs(h+y);
 j=abs(h-y);
 if(i>=(n<<1)-2)i%=(n<<1)-2;
 if(j>=(n<<1)-2)j%=(n<<1)-2;
 if(i>=n)i=(n<<1)-2-i;
 if(j>=n)j=(n<<1)-2-j;
 s+=w.s_l[(w.s_l_len>>1)+y]*(temp[i]+temp[j]);};
input[(x<<1)+w.s_l_shift]=s;};

for(x=0;x<(n>>1);x++){

h=(x<<1)+w.s_h_shift;
s=w.s_h[w.s_h_len>>1]*temp[h];
for(y=1;y<=(w.s_h_len>>1);y++){
 i=abs(h+y);
 j=abs(h-y);
 if(i>=(n<<1)-2)i%=(n<<1)-2;
 if(j>=(n<<1)-2)j%=(n<<1)-2;
 if(i>=n)i=(n<<1)-2-i;
 if(j>=n)j=(n<<1)-2-j;
 s+=w.s_h[(w.s_h_len>>1)+y]*(temp[i]+temp[j]);};
input[(x<<1)+w.s_h_shift]=s;};

}; Frigo 23:57, 8 March 2006 (UTC)

Hi, I gave the algorithm to compute the coefficients to any precision already in the article. One only has to numerically compute the real root of a polynomial of degree 3. The table was directly computed this way in the CAS Magma. I will repeat it next week with a little more precision.--LutzL 20:49, 10 March 2006 (UTC)

Hello. I have a question: which form of the filters is the correct form, the synthesis lowpass being longer than the synthesis highpass, or vice-versa? Most paper I have seen assumed the synthesis lowpass filter is longer. A long time ago I tried both, and somewhy I choose the lowpass filter to be longer. I don't remember why I chose it and how I have implemented the other... (Perhaps I felt it is more elegant =) , f.e. in the case of interpolating filters (like 5/3), the analysis and synthesis filters are equal, with sqrt(2) normalization) Anyway, I thought it is called interleaved filtering, but I'm not sure.. And what is the most appropiate normalization for displaying on wiki? I think it would be (1;1) and/or (sqrt(2);sqrt(2)). Frigo 23:59, 10 March 2006 (UTC)

Hi, for Your convinience the filters in an almost unpractical precision

k	adual/2	bdual/2
-4	0.026748757410810088414008097968967319312004837006597093574251	0
-3	-0.016864118442874954426246480229347056133010595557498776661922	0.045635881557125045573753519770652943866989404442501223338078
-2	-0.078223266528990262508525352334563389514040539141385927620846	-0.028771763114250091147507039541305887733978808885002446676157
-1	0.266864118442874954426246480229347056133010595557498776661922	-0.295635881557125045573753519770652943866989404442501223338078
0	0.60294901823636034818903450873119214040407140426957766809319	0.557543526228500182295014079082611775467957617770004893352314
1	0.266864118442874954426246480229347056133010595557498776661922	-0.295635881557125045573753519770652943866989404442501223338078
2	-0.078223266528990262508525352334563389514040539141385927620846	-0.028771763114250091147507039541305887733978808885002446676157
3	-0.016864118442874954426246480229347056133010595557498776661922	0.045635881557125045573753519770652943866989404442501223338078
4	0.026748757410810088414008097968967319312004837006597093574251	0

I don't know a final answer to Your questions. In image processing it is preferable if the synthesis scaling filter belongs to the smoother scaling function. The errors while quantizing may be bigger with the "bumpier" analysis fitler, but the reconstructed image looks smooth. For 97 I don't know ad hoc the smoothness of both scaling functions, but for 53 resp bior-4,4, the shorter sequence belongs to the spline of order 3, which is as smooth as it can get. And no, You are wrong, the synthesis filters in 53 have also a differing length, as explained and computed in the article. For the normalization, well there has to be a factor of 1/2 somewhere, but the exact location is not important for bitplane coding etc.--LutzL 12:07, 11 March 2006 (UTC)

I noticed something strange: if I multiply the analysis highpass coefficients with -1, the synthesis lowpass filter becomes a highpass filter, and vice-versa. (basically, every 2nd of their coeffs gets multiplied by -1). So the lowpass filter becomes the shorter of the two, but it can't be better than the longer one, since we only multiplied the an. highpass coeffs with -1.

so, in the case of 5/3 filter:

original: (synthesis highpass first, lowpass second in interleave)

-0.2500000000000000,
 0.5000000000000000,
 1.5000000000000000,
 0.5000000000000000,
-0.2500000000000000,

and

-0.5000000000000000,
 1.0000000000000000,
-0.5000000000000000,

"new": (synthesis lowpass first, highpass second in interleave)

 0.2500000000000000,
 0.5000000000000000,
-1.5000000000000000,
 0.5000000000000000,
 0.2500000000000000,

and

 0.5000000000000000,
 1.0000000000000000,
 0.5000000000000000,

(And I obtain the synthesis filters via matrix inversion; analysis lowpass filter is first, analysis highpass filter is second in interleave; so I can't be wrong)

---

I noticed that the reconstruction of only the LL subband is smoother with an analysis 9/7 (L/H) filter (with 9/3 or 3/9 synthesis filter), than with an analysis 9/3 (L/H) filter (with 9/7 or 7/9 synthesis filter), the latter is full of "diamond" shape.

---

And a question: how do you reconstruct? Upsample and filter and add, or interleave and filter?

Frigo 00:19, 15 March 2006 (UTC)

[edit] The meaning of "Interleave"

Interleave and filter is not going to work. If I had to implement it efficiently, I would use lifting. But for the exposition here the first variant has to be applied, as it is explained in Fast wavelet transform.--LutzL 09:54, 16 March 2006 (UTC)

I see, Hilbert space.

But interleaved filtering can be deducted from upsampling + filtering + add, because of the add operation. And because we insert 0s, the synthesis filters can be trivially "transformed" between interleaved filtering and upsampled filtering by exchanging every 2nd coefficient of the synthesis lowpass and highpass filters.

It works. When I multiply my interleaved filtering's transform matrices, it results in a unitary matrix, when I decompose + recompose Lena or Baboon without modifying anything, it results in the original image. So it works. One can easily see this if he inverts a wavelet transform matrix.

Frigo 17:50, 17 March 2006 (UTC)

Hi, I added some headlines, so that one can edit in pieces. Of course You are right, I didn't understand the operation behind the word. If You check the "reconstruction" section in the "fast transform" piece, the abstract formula is exactly that, upsample, filter, add. By the way, the whole wavelet thing here is a total mess. The Discrete wavelet transform has less mathematical contents than the fast transform, lifting is still missing. The continuous wavelet transform is both mentioned in wavelet series, where it doesn't belong, and in wavelet, where it has too much detail, in fact more than in its own article. The latter should be a link to wavelet theory, which should be an overview with hardly more than short explainations of the list of wavelet-related transforms and the lists already present, commons:wavelet and commons:Category:wavelets are better organized;-). "discrete wavelet" should be an overview as well, linking together MRA, FWT, complex extension of both, the different classes ("haar", "daubechies", "coiflet", "biorthogonal") etc. Its all here, but in bad shape.--LutzL 19:05, 17 March 2006 (UTC)