Different results with approximation coefficients with ndwt and swt

Should one expect that the results generated from swt and ndwt to match?

I am attempting the following but as one can see, the results are different. Is there anyway to obtain similar denoising characteristics using ndwt as one can obtain with swt?

load noissin;
x = noissin(1:512);
W = ndwt(x,4,'db4','mode','per');
a = indwt(W, 'a', 4);
[swa,swd] = swt(x,4,'db4');
mzero = zeros(size(swd));
b = iswt(swa,mzero,'db4');
subplot(211),plot(a);
subplot(212),plot(b);

Answer 1

[Wayne King]

Dave, what you can do is symmetrically extend the time series, then operate on those wavelet coefficients and then keep the "central" part.

For example:

      load Data_GlobalIdx1; 
      dat = cumsum(diff(log(Dataset.SP)));
      N = length(dat);
      L = nextpow2(N);
      ext = 2^L;
      difflen = 2^L-N;
      ext = round(difflen/2);
      SIG = wextend('ar','sym',dat,ext);

Now I'm going to do something a bit complicated here but, I'm just going to denoise dat before I invert it. This is just the level-dependent universal threshold (with soft thresholding). I'll invert at the end and keep the "central" part and plot the result.

   wname = 'db4';
    level = 5;

    sorh = 's';    % Specified soft or hard thresholding
    thrSettings =  {...
      [...
      1   4096     0.387; ...
      ]; ...
      [...
      1   4096      0.478; ...
      ]; ...
      [...
      1   4096      0.647; ...
      ]; ...
      [...
      1   4096      0.966; ...
      ]; ...
      [...
      1   4096      1.41; ...
      ]; ...
      };

    % Decompose using SWT.
    wDEC = swt(SIG,level,wname);

    % Denoise.
    %---------
    len = length(SIG);
     for k = 1:level
      thr_par = thrSettings{k};
      if ~isempty(thr_par)
          NB_int = size(thr_par,1);
          x      = [thr_par(:,1) ; thr_par(NB_int,2)];
          x      = round(x);
          x(x<1) = 1;
          x(x>len) = len;
          thr = thr_par(:,3);
          for j = 1:NB_int
              if j==1 , d_beg = 0; else d_beg = 1; end
              j_beg = x(j)+d_beg;
              j_end = x(j+1);
              j_ind = (j_beg:j_end);
              wDEC(k,j_ind) = wthresh(wDEC(k,j_ind),sorh,thr(j));
          end
      end
  end

% Reconstruct the denoise signal using ISWT.
%-------------------------------------------
sigDEN = iswt(wDEC,wname);

% now discard the extension

   sigDEN = sigDEN(ext+1:ext+length(dat));
   plot(dat,'r');
   hold on;
   plot(sigDEN,'linewidth',2);
   legend('Original Data','Denoised Signal');

[Dave]

Thank you! This is great. Can you explain at all how you came up with the values for thr. I was under the impression that ddencmp was usually used for this.

[Dave]

I also attempted using:

[sigden,coefs,thrParamsOut] = cmddenoise(dat,'db4',5);

But was unable to come up with similar threshold values.

Answer 2

[Wayne King]

Hi Dave, there is a bug in ndwt.m. This has been reported. swt() is correct. Are you not able to use swt() due to length constraints?

The bug in ndwt.m is that the subband filtering only works correctly for the first stage and not subsequently. You can see the results of this bug in the fact that your 'a' vector above has high frequency oscillations when in fact, the approximation coefficients at that level, level 4, should be very smooth (containing no high frequencies)

[Dave]

That is correct, my input is arbitrary in length.

What you are describing is exactly what I am seeing. What can one do to get around this until a fix is made available?

[Wayne King]

The thing to do here Dave is extend your data vector to a length compatible with the swt() and then truncate it upon output. You have to be a little smart about the way you extend it. Do you have a sample data file you can provide? If you do that and tell me the level of the multiresolution analysis you want and the wavelet you want to use, I can show you an example of how to do that.

[Dave]

The following data should be reasonable.

load Data_GlobalIdx1 dat = cumsum(diff(log(Dataset.SP)));

I am using a 'db4' wavelet and level depth of 4 or 5.

When using a wavelet over a sliding window (new data points arriving in real time) how much of an issue is the right most end point?