Problems with periodization
in [Matlab]
|
From: Ivan Goremiko on 15 Jun 2010 12:25 "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7lp5$ksr$1(a)fred.mathworks.com>... > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7l6c$ekk$1(a)fred.mathworks.com>... > > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv6ch4$fek$1(a)fred.mathworks.com>... > > > I am trying to implement db4 DWT with periodazed wavelets. > > > However,my results is slightly different what matlab does. > > > As I noticed my code works fine for one-level decomposition,but > > > and higher level I've got details coefficients different from matlab code. > > > Approximation coefficients are calculating fine. > > > My algorithm in general: > > > > > > read signal > > > > > > apply circular convolution to low pass filter > > > apply circular convolution to high pass filter > > > shift last element to first (??)I have no idea why it is needed in Matlab > > > signal = details coefficients > > > repeat --. > > > > > > Can everyone tell me what is wrong with this approach and what matlab doing with > > > dwtmode set to ''per'? > > > > Hi Ivan, are you iterating on the lowpass filter output? Are you downsampling? Perhaps you're sampling a different phase of the wavelet coefficients than Matlab. > > > > Compare the following: > > > > % getting a reproducible signal > > reset(RandStream.getDefaultStream) > > x = randn(8,1); > > > > % Daubechies 4 extremal phase wavelet > > [Lo,Hi] = wfilters('db4'); > > > > % Level 1 > > extend = length(Lo)/2; > > y = wextend('1D','per',x,extend); > > A1 = conv2(y(:)',Lo(:)','valid'); > > A1 = A1(2:2:length(x)); > > D1 = conv2(y(:)',Hi(:)','valid'); > > D1 = D1(2:2:length(x)); > > > > % A1 are the level one approximation and D1 the level one details > > > > % Level 2 > > > > % input is A1 the output of the lowpass filter downsampled > > y = wextend('1D','per',A1,extend); > > A2 = conv2(y(:)',Lo(:)','valid'); > > A2 = A2(2:2:length(x)/2); > > D2 = conv2(y(:)',Hi(:)','valid'); > > D2 = D2(2:2:length(x)/2); > > > > % Compare > > [A2 D2 D1] > > > > with > > > > [C,L] = wavedec(x,2,'db4'); > > C' > > > > > > Finally, as far as "shift last element to first (??)I have no idea why it is needed in Matlab" > > > > I think you're just talking about ordering the C vector output of wavedec(). There is no magical reason why you have to do that. That is simply the way the developers decided to order the vector of coefficients. > > > > Hope that helps, > > Wayne Thanks for answer Wayne. One more question - At every level one has to extend signal periodically with "extend" elements (half of original) ? And can you write ,please, code for inverse transform?
From: Ivan Goremiko on 15 Jun 2010 12:32 "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7lp5$ksr$1(a)fred.mathworks.com>... > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7l6c$ekk$1(a)fred.mathworks.com>... > > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv6ch4$fek$1(a)fred.mathworks.com>... > > > I am trying to implement db4 DWT with periodazed wavelets. > > > However,my results is slightly different what matlab does. > > > As I noticed my code works fine for one-level decomposition,but > > > and higher level I've got details coefficients different from matlab code. > > > Approximation coefficients are calculating fine. > > > My algorithm in general: > > > > > > read signal > > > > > > apply circular convolution to low pass filter > > > apply circular convolution to high pass filter > > > shift last element to first (??)I have no idea why it is needed in Matlab > > > signal = details coefficients > > > repeat --. > > > > > > Can everyone tell me what is wrong with this approach and what matlab doing with > > > dwtmode set to ''per'? > > > > Hi Ivan, are you iterating on the lowpass filter output? Are you downsampling? Perhaps you're sampling a different phase of the wavelet coefficients than Matlab. > > > > Compare the following: > > > > % getting a reproducible signal > > reset(RandStream.getDefaultStream) > > x = randn(8,1); > > > > % Daubechies 4 extremal phase wavelet > > [Lo,Hi] = wfilters('db4'); > > > > % Level 1 > > extend = length(Lo)/2; > > y = wextend('1D','per',x,extend); > > A1 = conv2(y(:)',Lo(:)','valid'); > > A1 = A1(2:2:length(x)); > > D1 = conv2(y(:)',Hi(:)','valid'); > > D1 = D1(2:2:length(x)); > > > > % A1 are the level one approximation and D1 the level one details > > > > % Level 2 > > > > % input is A1 the output of the lowpass filter downsampled > > y = wextend('1D','per',A1,extend); > > A2 = conv2(y(:)',Lo(:)','valid'); > > A2 = A2(2:2:length(x)/2); > > D2 = conv2(y(:)',Hi(:)','valid'); > > D2 = D2(2:2:length(x)/2); > > > > % Compare > > [A2 D2 D1] > > > > with > > > > [C,L] = wavedec(x,2,'db4'); > > C' > > > > > > Finally, as far as "shift last element to first (??)I have no idea why it is needed in Matlab" > > > > I think you're just talking about ordering the C vector output of wavedec(). There is no magical reason why you have to do that. That is simply the way the developers decided to order the vector of coefficients. > > > > Hope that helps, > > Wayne Ashame on me - I was applying low pass and high pass filter on both - details and approximation coefficients.
From: Wayne King on 16 Jun 2010 06:41 "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv89hg$sh4$1(a)fred.mathworks.com>... > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7lp5$ksr$1(a)fred.mathworks.com>... > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7l6c$ekk$1(a)fred.mathworks.com>... > > > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv6ch4$fek$1(a)fred.mathworks.com>... > > > > I am trying to implement db4 DWT with periodazed wavelets. > > > > However,my results is slightly different what matlab does. > > > > As I noticed my code works fine for one-level decomposition,but > > > > and higher level I've got details coefficients different from matlab code. > > > > Approximation coefficients are calculating fine. > > > > My algorithm in general: > > > > > > > > read signal > > > > > > > > apply circular convolution to low pass filter > > > > apply circular convolution to high pass filter > > > > shift last element to first (??)I have no idea why it is needed in Matlab > > > > signal = details coefficients > > > > repeat --. > > > > > > > > Can everyone tell me what is wrong with this approach and what matlab doing with > > > > dwtmode set to ''per'? > > > > > > Hi Ivan, are you iterating on the lowpass filter output? Are you downsampling? Perhaps you're sampling a different phase of the wavelet coefficients than Matlab. > > > > > > Compare the following: > > > > > > % getting a reproducible signal > > > reset(RandStream.getDefaultStream) > > > x = randn(8,1); > > > > > > % Daubechies 4 extremal phase wavelet > > > [Lo,Hi] = wfilters('db4'); > > > > > > % Level 1 > > > extend = length(Lo)/2; > > > y = wextend('1D','per',x,extend); > > > A1 = conv2(y(:)',Lo(:)','valid'); > > > A1 = A1(2:2:length(x)); > > > D1 = conv2(y(:)',Hi(:)','valid'); > > > D1 = D1(2:2:length(x)); > > > > > > % A1 are the level one approximation and D1 the level one details > > > > > > % Level 2 > > > > > > % input is A1 the output of the lowpass filter downsampled > > > y = wextend('1D','per',A1,extend); > > > A2 = conv2(y(:)',Lo(:)','valid'); > > > A2 = A2(2:2:length(x)/2); > > > D2 = conv2(y(:)',Hi(:)','valid'); > > > D2 = D2(2:2:length(x)/2); > > > > > > % Compare > > > [A2 D2 D1] > > > > > > with > > > > > > [C,L] = wavedec(x,2,'db4'); > > > C' > > > > > > > > > Finally, as far as "shift last element to first (??)I have no idea why it is needed in Matlab" > > > > > > I think you're just talking about ordering the C vector output of wavedec(). There is no magical reason why you have to do that. That is simply the way the developers decided to order the vector of coefficients. > > > > > > Hope that helps, > > > Wayne > Thanks for answer Wayne. > One more question - > At every level one has to extend signal periodically with "extend" elements (half of original) ? > And can you write ,please, code for inverse transform? Hi Ivan, the extension is based on the length of the lowpass analysis filter. It just happened in that little example I gave you that those two lengths were the same. It's 1/2 the length of the lowpass analysis filter. Wayne
From: Ivan Goremiko on 22 Jun 2010 10:53 "Wayne King" <wmkingty(a)gmail.com> wrote in message <hva9o0$hud$1(a)fred.mathworks.com>... > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv89hg$sh4$1(a)fred.mathworks.com>... > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7lp5$ksr$1(a)fred.mathworks.com>... > > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7l6c$ekk$1(a)fred.mathworks.com>... > > > > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv6ch4$fek$1(a)fred.mathworks.com>... > > > > > I am trying to implement db4 DWT with periodazed wavelets. > > > > > However,my results is slightly different what matlab does. > > > > > As I noticed my code works fine for one-level decomposition,but > > > > > and higher level I've got details coefficients different from matlab code. > > > > > Approximation coefficients are calculating fine. > > > > > My algorithm in general: > > > > > > > > > > read signal > > > > > > > > > > apply circular convolution to low pass filter > > > > > apply circular convolution to high pass filter > > > > > shift last element to first (??)I have no idea why it is needed in Matlab > > > > > signal = details coefficients > > > > > repeat --. > > > > > > > > > > Can everyone tell me what is wrong with this approach and what matlab doing with > > > > > dwtmode set to ''per'? > > > > > > > > Hi Ivan, are you iterating on the lowpass filter output? Are you downsampling? Perhaps you're sampling a different phase of the wavelet coefficients than Matlab. > > > > > > > > Compare the following: > > > > > > > > % getting a reproducible signal > > > > reset(RandStream.getDefaultStream) > > > > x = randn(8,1); > > > > > > > > % Daubechies 4 extremal phase wavelet > > > > [Lo,Hi] = wfilters('db4'); > > > > > > > > % Level 1 > > > > extend = length(Lo)/2; > > > > y = wextend('1D','per',x,extend); > > > > A1 = conv2(y(:)',Lo(:)','valid'); > > > > A1 = A1(2:2:length(x)); > > > > D1 = conv2(y(:)',Hi(:)','valid'); > > > > D1 = D1(2:2:length(x)); > > > > > > > > % A1 are the level one approximation and D1 the level one details > > > > > > > > % Level 2 > > > > > > > > % input is A1 the output of the lowpass filter downsampled > > > > y = wextend('1D','per',A1,extend); > > > > A2 = conv2(y(:)',Lo(:)','valid'); > > > > A2 = A2(2:2:length(x)/2); > > > > D2 = conv2(y(:)',Hi(:)','valid'); > > > > D2 = D2(2:2:length(x)/2); > > > > > > > > % Compare > > > > [A2 D2 D1] > > > > > > > > with > > > > > > > > [C,L] = wavedec(x,2,'db4'); > > > > C' > > > > > > > > > > > > Finally, as far as "shift last element to first (??)I have no idea why it is needed in Matlab" > > > > > > > > I think you're just talking about ordering the C vector output of wavedec(). There is no magical reason why you have to do that. That is simply the way the developers decided to order the vector of coefficients. > > > > > > > > Hope that helps, > > > > Wayne > > Thanks for answer Wayne. > > One more question - > > At every level one has to extend signal periodically with "extend" elements (half of original) ? > > And can you write ,please, code for inverse transform? > > Hi Ivan, the extension is based on the length of the lowpass analysis filter. It just happened in that little example I gave you that those two lengths were the same. It's 1/2 the length of the lowpass analysis filter. > > Wayne Hi Wayne. Can you also explain how to do inverse transform up to third level?
From: Ivan Goremiko on 22 Jun 2010 11:00
"Wayne King" <wmkingty(a)gmail.com> wrote in message <hva9o0$hud$1(a)fred.mathworks.com>... > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv89hg$sh4$1(a)fred.mathworks.com>... > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7lp5$ksr$1(a)fred.mathworks.com>... > > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hv7l6c$ekk$1(a)fred.mathworks.com>... > > > > "Ivan Goremiko" <merlushatrash(a)gmail.com> wrote in message <hv6ch4$fek$1(a)fred.mathworks.com>... > > > > > I am trying to implement db4 DWT with periodazed wavelets. > > > > > However,my results is slightly different what matlab does. > > > > > As I noticed my code works fine for one-level decomposition,but > > > > > and higher level I've got details coefficients different from matlab code. > > > > > Approximation coefficients are calculating fine. > > > > > My algorithm in general: > > > > > > > > > > read signal > > > > > > > > > > apply circular convolution to low pass filter > > > > > apply circular convolution to high pass filter > > > > > shift last element to first (??)I have no idea why it is needed in Matlab > > > > > signal = details coefficients > > > > > repeat --. > > > > > > > > > > Can everyone tell me what is wrong with this approach and what matlab doing with > > > > > dwtmode set to ''per'? > > > > > > > > Hi Ivan, are you iterating on the lowpass filter output? Are you downsampling? Perhaps you're sampling a different phase of the wavelet coefficients than Matlab. > > > > > > > > Compare the following: > > > > > > > > % getting a reproducible signal > > > > reset(RandStream.getDefaultStream) > > > > x = randn(8,1); > > > > > > > > % Daubechies 4 extremal phase wavelet > > > > [Lo,Hi] = wfilters('db4'); > > > > > > > > % Level 1 > > > > extend = length(Lo)/2; > > > > y = wextend('1D','per',x,extend); > > > > A1 = conv2(y(:)',Lo(:)','valid'); > > > > A1 = A1(2:2:length(x)); > > > > D1 = conv2(y(:)',Hi(:)','valid'); > > > > D1 = D1(2:2:length(x)); > > > > > > > > % A1 are the level one approximation and D1 the level one details > > > > > > > > % Level 2 > > > > > > > > % input is A1 the output of the lowpass filter downsampled > > > > y = wextend('1D','per',A1,extend); > > > > A2 = conv2(y(:)',Lo(:)','valid'); > > > > A2 = A2(2:2:length(x)/2); > > > > D2 = conv2(y(:)',Hi(:)','valid'); > > > > D2 = D2(2:2:length(x)/2); > > > > > > > > % Compare > > > > [A2 D2 D1] > > > > > > > > with > > > > > > > > [C,L] = wavedec(x,2,'db4'); > > > > C' > > > > > > > > > > > > Finally, as far as "shift last element to first (??)I have no idea why it is needed in Matlab" > > > > > > > > I think you're just talking about ordering the C vector output of wavedec(). There is no magical reason why you have to do that. That is simply the way the developers decided to order the vector of coefficients. > > > > > > > > Hope that helps, > > > > Wayne > > Thanks for answer Wayne. > > One more question - > > At every level one has to extend signal periodically with "extend" elements (half of original) ? > > And can you write ,please, code for inverse transform? > > Hi Ivan, the extension is based on the length of the lowpass analysis filter. It just happened in that little example I gave you that those two lengths were the same. It's 1/2 the length of the lowpass analysis filter. > > Wayne I have approximation coefficients - A =0.3220 0 0.0097 0 0.1302 0 0.6210 0 details D =0.3335 0 -0.0926 0 -0.0200 0 -0.0344 0 on second level. Then I extend them by 2 elements periodically . Apply conv2(or just conv) and add each other. What I've got is not correct. Can somebody explain me how will be correct? |