nmdemo.m

matlab程序

      title('Two levels of decomposition');
   case 5  %%% Demo 5: Max-Lifting Wavelet %%%
      fprintf(['\n     Max-Lifting Wavelet\n\n',...
            'This morphological wavelet, designed using the lifting scheme, has a\n',...
            'decimation operator such that local maxima in a signal are preserved.\n',...
            'The decomposition algorithm is\n\n',...
            '   s[n] = x[2n]\n   d[n] = x[2n+1]\n   d[n] = d[n] - max {s[n], s[n+1]}\n',...
            '   s[n] = s[n] + max {0, d[n-1], d[n]}\n\n',...
            'The reconstruction algorithm is\n\n   s[n]    = s[n] - max {0, d[n-1], d[n]}\n',...
            '   d[n]    = d[n] + max {s[n], s[n+1]}\n   x[2n]   = s[n]\n   x[2n+1] = d[n]\n\n',...
            'The figure shows an example decomposition.  This transform is implemented\n',...
            'in MAXLIFT and MAXLIFT2.\n\n\n']);
      set(G_DEMO(5),'Name','Max-Lifting Wavelet');
      x = [1,1,1,1,1,1,1,1,1,1,1,6,4,4,4,4,4,4,4,0,0,5,4,6,0,0,0,0,0,0,0,0,0,0,0,0,...
         0,0,3,3,3,3,3,4,4,4,4,4,4,6,4,4,3,3,3,3,3,3,3,3,3,3,3,3];
      subplot(2,1,1);
      plot(0:63,x,'.-');
      xlim([0,63]);
      title('Original');
      subplot(2,1,2);
      y = maxlift(x,2);   % Perform two stages of max-lifting decomposition 
      plot(0:15,y(1:16),'b.-',16:63,y(17:64),'b.',...
         [15.5,15.5],[-4,6],'k',[31.5,31.5],[-4,6],'k');
      axis([0,63,-4,6]);
      title('Two levels of decomposition');
   case 6  %%% Demo 6: Denoising with Max-Lifting %%%
      fprintf(['\n     Denoising with Max-Lifting\n\n',...
            'The max-lifting wavelet (see previous demo) can be used to denoise signals\n',...
            'where the underlying function is piecewise constant.  In the figure, the\n',...
            'top plot shows a piecewise signal with additive Gaussian noise.  The noisy\n',...
            'signal is transformed using the max-lifting wavelet and soft thresholding is\n',...
            'applied to the transform coefficients.  The bottom plot shows the denoised\n',...
            'signal from the reconstruction.\n\n\n']);
      set(G_DEMO(5),'Name','Denoising with Max-Lifting');
      x = seqhaar([1,1,1,1,1,0,0,2,2,2,2,2,2,2,2,2,2,2,2,3,3,1,1,1,1,0,0,0,0,0,0,0,0,0,0,...
            0,0,0,3,3,3,3,3,3,3,3,0,0,0,2,2,2,2,2,0,0,0,1,1,1,1,1,1,1,zeros(1,192)],-2) ...
      +randn(1,256)/4;
      subplot(2,1,1);
      plot(0:255,x,'-');
      axis([0,255,-1,4]);
      title('Noisy Signal');
      y = maxlift(x,3);           % Perform 3 stages of decomposition
      % Use soft thresholding 
      y(129:256) = softthresh(y(129:256),.8);  % Threshold 3rd (finest) detail subband
      y(65:128) = softthresh(y(65:128),.4);    % Threshold 2nd detail subband
      y(33:64) = softthresh(y(33:64),.2);      % Threshold 1st detail subband
      subplot(2,1,2);
      plot(0:255,maxlift(y,-3),'-',0:255,x,'k:');
      axis([0,255,-1,4]);
      title('Denoised with Max Lifting');
   case 7  %%% Demo 7: ENO Interpolation %%%
      fprintf(['\n     ENO Interpolation\n\n',...
            'The Essentially Non-Oscillatory (ENO) scheme uses a clever data-dependent\n',...
            'interpolation strategy.  For every interval, the neighboring samples are\n',...
            'analyzed to find the smoothest local neighborhood.  This set of neighbors is\n',...
            'then used to interpolate the cell.\n\n',...
            'The figure shows a piecewise smooth signal interpolated with both cubic\n',...
            'interpolation and ENO interpolation.  Around the jump discontinuity, the\n',...
            'cubic interpolation has oscillation artifacts, however, the ENO interpolant\n',...
            'does not, hence its name.  ENO interpolation is implemented in ENOINT.\n\n\n']);
      set(G_DEMO(5),'Name','ENO Interpolation');
      x = [cos(1.5+(.2:.6:4.3))*0.7+0.3,0.95-(.16:.08:.88).^3];
      x(8:11) = -0.7;
      yi = interp1(x,1:.125:17,'cubic');
      yeno = enoint(enoint(enoint(x,2),2),2);  % Apply enoint three times for 8x interpolation
      subplot(2,1,1);
      set(plot(0:.125:16,yi,'b-',5:.125:6,yi(41:49),'r-',7:.125:8,yi(57:65),'r-',...
         9:.125:10,yi(73:81),'r-',11:.125:12,yi(89:97),'r-',...
         0:16,x,'k.'),'MarkerSize',16);
      xlim([0,16]);         
      title('Cubic Interpolation');
      subplot(2,1,2);
      set(plot(0:.125:16,yeno,'b-',0:16,x,'k.'),'MarkerSize',16);
      xlim([0,16]);
      title('ENO Interpolation');
   case 8  %%% Demo 8: ENO Multiresolution %%%
      fprintf(['\n     ENO Multiresolution\n\n',...
            'Using ENO interpolation (see previous demo), the ENO multiresolution scheme\n',...
            'can very sparsely represent piecewise polynomial signals.\n\n',...
            'The figure shows the decomposition of a piecewise polynomial signal.  ENO\n',...
            'multiresolution is implemented in ENOPV, ENOPV2, ENOCA, and ENOCA2.\n\n\n']);
      set(G_DEMO(5),'Name','ENO Multiresolution');
      x = polyval([0.6,-1,0,.2],0:2^-7:1);
      x2 = polyval([-.2,0,.2,0],0:2^-7:1);
      x([20:40,72:110]) = x2([20:40,72:110]);
      subplot(2,1,1);
      plot(0:2^-7:1,x,'-');
      title('Original');
      subplot(2,1,2);
      y = enopv(x,3,3);   % Decompose with 3 stages and cubic interpolation
      plot(0:16,y(1:17),'b.-',17:128,y(18:129),'b.',[16.5,16.5],[-.2,.2],...
         'k-',[32.5,32.5],[-.2,.2],'k-',[64.5,64.5],[-.2,.2],'k-');
      axis([0,128,-.2,.2]);
      title('Three levels of decomposition');
   case 9  %%% Demo 9: ENO Approximation %%%
      fprintf(['\n     ENO Approximation\n\n',...
            'ENO can be used to approximate a signal using a modified encoding \n',...
            'algorithm.  This algorithm decomposes the signal such that the\n',...
            'representation is as sparse as possible while maintaining the\n',...
            'maximum error between the reconstruction and the original signal\n',...
            'within a specified tolerance.  The figure shows an example of\n',...
            'approximation with point-value ENO.  Point-value ENO\n',...
            'approximation is implemented in ENOPV.\n\n\n']);
      set(G_DEMO(5),'Name','ENO Approximation');
      
      t = linspace(0,1,257).';
      x = 1.4*(1 - 100*(t-0.5).^2).*exp(-(t-0.5).^2*50) + 1.4;
      x(101:157) = x(101:157) - 1;
      x(129:257) = sin((129:257)/5)./((3:131)/5) + 0.6;
      
      subplot(2,2,1);
      plot(0:256,x); axis([0,256,0,2]); title('Original Signal (257 samples)');
      subplot(2,2,2);
      y = enopv(x,4,3,0.05);     % Perform approximate decomposition with tolerance 0.05
      Num = length(find(abs(y) > eps));      
      a = enopv(y,-4,3);         % Reconstruct from approximate decomposition
      plot(0:256,x,'k:',0:256,a); axis([0,256,0,2]); title('Approximation (Tolerance = 0.05)');
      subplot(2,2,3);
      plot(0:256,x - a); xlim([0,256]); title('Approximation Error');
      subplot(2,2,4);
      plot(0:256,y,'.',[16.5,16.5],[-1,2],'k',[32.5,32.5],[-1,2],'k',[64.5,64.5],[-1,2],'k', ...
         [128.5,128.5],[-1,2],'k');
      axis([0,256,-1,2]); title(sprintf('Decomposition (%d/257 components)',Num));
   case 10  %%% References %%%
      s = {'F. Arandiga and R. Donat.  ``Nonlinear Multiscale Decompositions: The',...
            'Approach of A. Harten.''''  Numerical Algorithms 23 (2000), pp. 175-216.','',...
            'A. Harten.  ``Discrete Multiresolution Analysis and Generalized',...
            'Wavelets.'''' J. Applied Numerical Mathematics (1993), vol. 12,',...
            'no. 1-3, pp. 152-192.','',...
            'A. Harten. ``Multiresolution Representation of Data II: General',...
            'Framework.'''' SIAM J. Numerical Analysis 33 (1996), pp. 1205-1256.','',...
            'H. Heijmans and J. Goutsias.  ``Nonlinear Multiresolution Signal',...
            'Decomposition Schemes--Part II: Morphological Wavelets.''''  IEEE',...
            'Trans. on Image Processing (2000), vol. 9, no. 11.','',...
            'M. Swanson and A. Tewfik.  ``A Binary Wavelet Decomposition of Binary',...
            'Images.''''  IEEE Trans. on Image Processing (1996), vol. 5, no. 12.','',...
            'R. Claypoole, G. Davis, W. Sweldens, and R. Baraniuk.  ``Nonlinear',...
            'Wavelet Transforms for Image Coding.'''' November 1997.','',''};
      fprintf('\n     References\n\n');
      disp(strvcat(s));    % Display references on console
      set(G_DEMO(5),'Name','References');
      % Also create a listbox in the figure window displaying the references
      uicontrol('Style','list','Units','Normalized','Position',[0.03,0.03,0.94,0.97], ...
         'HorizontalAlignment','Left','String',s,'BackgroundColor',[1,1,1]);
   end
end

return;


function x = softthresh(x,t)
% Applies soft thresholding to signal x with threshold t.
z = (abs(x)-t)/2;
x = sign(x).*(z+abs(z));
return;