小波变换和多分辨率处理

偏执的太偏执、 2022-08-18 15:20 107阅读 0赞

**小波变换和多分辨率处理**

# **1.基本概念：** #

**小波变换（wavelet transform，WT）**是一种新的变换分析方法，它继承和发展了短时傅立叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点，能够提供一个随频率改变的“时间-频率”窗口，是进行信号时频分析和处理的理想工具。它的主要特点是通过变换能够充分突出问题某些方面的特征，能对时间（空间）频率的局部化分析，通过伸缩平移运算对信号（函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，解决了Fourier变换的困难问题，成为继Fourier变换以来在科学方法上的重大突破。

# 2.用到的函数： #

（1）用到matlab自带的函数：

wavemenu：Wavelet GUI工具；

Dwt, dwt2： 一维、二维小波变换；

idwt, idwt2：一维、二维小波反变换。

（2）非matlab自带的函数：

Wavefilter，Wave2gray，wavecut等。

# 3.matlab源码： #

(1) wave2gray.m

func function w = wave2gray(c,  s,  scale,  border)
    %WAVE2GRAY Display wavelet decomposition coefficients.
    %	W = WAVE2GRAY(C, S, SCALE, BORDER) displays and returns a
    %	wavelet coefficient image.
    %		EXAMPLES:
    %			wave2gray(c, s);						Display w/defaults.
    %			foo = wave2gray(c,   s);				Display and return.
    %			foo = wave2gray(c,   s, 4);				Magnify the details.
    %			foo = wave2gray(c,   s, -4);			Magnify absolute values.
    %			foo = wave2gray(c,   s, 1, 'append');	Keep border values.
    %
    %		INPUTS/OUTPUTS:
    %			[C, S] is a wavelet decomposition vector and bookkeeping
    %			matrix.
    %
    %		SCALE				Detail coefficient scaling
    %--------------------------------------------------------------------
    %		0 or 1      Maximum range (default)
    %		2, 3...     Magnify default by the scale factor
    %		-1, -2...   Magnify absolute values by abs(scale)
    %
    %		BORDER			Border between wavelet decompositions
    %--------------------------------------------------------------------
    %		'absorb'		Border replaces image (default)
    %		'append'		Border increases width of image
    %
    %		Image W:	------- ------ -------------- -------------------
    %					|      |      |              |
    %					| a(n) | h(n) |              |
    %					|      |      |              |
    %					------- ------     h(n-1)    |
    %					|      |      |              |
    %					| v(n) | d(n) |              |      h(n-2)
    %					|      |      |              |
    %					------- ------ --------------
    %					|             |              |
    %					|    v(n-1)   |    d(n-1)    |
    %					|             |              |
    %					-------------- -------------- -------------------
    %					|                            |
    %					|           v(n-2)           |      d(n-2)
    %					|                            |
    %
    %		Here,  n denotes the decomposition step scale and a, h, v, d are
    %		approximation,  horizontal,  vertical, and diagonal detail
    %		coefficients,  respectively.
    
    % Check input arguments for  reasonableness.
    
    error(nargchk(2,   4,   nargin));
    
    if (ndims(c) ~= 2)   |   (size(c, 1) ~= 1)
    	error('C must be a row vector.'); end
    
    if (ndims(s) ~= 2)   | ~isreal(s)   | ~isnumeric(s)   |   (size(s, 2) ~= 2)
    	error('S must be a real, numeric two-column array.'); end
    
    elements = prod(s, 2);
    
    if (length(c) < elements(end))   |   ...
    		~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
    	error(['[C S] must be a standard wavelet ' ...
    				'decomposition structure.']);
    end
    
    if (nargin > 2) & (~isreal(scale) | ~isnumeric(scale))
    	error('SCALE must be a real, numeric scalar.');
    end
    
    if (nargin > 3) & (~ischar(border))
    	error('BORDER must be character string.');
    end
    
    if  nargin == 2
    	scale =1;    % Default scale.
    end
    
    if nargin < 4
    	border = 'absorb'; % Default border.
    end
    
    % Scale coefficients and determine pad fill.
    absflag = scale < 0;
    scale = abs(scale);
    if scale == 0
    	scale = 1;
    end
    
    [cd, w] = wavecut('a', c, s);  w = mat2gray(w);
    cdx = max(abs(cd(:))) / scale;
    if absflag
    	cd = mat2gray(abs(cd), [0, cdx]); fill = 0;
    else
    	cd = mat2gray(cd, [-cdx, cdx]); fill = 0.5;
    end
    
    % Build gray image one decomposition at a time.
    for i = size(s, 1) - 2:-1:1
    	ws = size(w);
    	h = wavecopy('h', cd, s, i);
    	pad = ws - size(h);    frontporch = round(pad / 2);
    	h = padarray(h, frontporch, fill, 'pre');
    	h = padarray(h, pad - frontporch, fill, 'post');
    	v = wavecopy('v',   cd,   s,   i);
    	pad = ws - size(v);            frontporch = round(pad  /  2);
    	v = padarray(v,   frontporch,   fill,   'pre');
    	v = padarray(v,  pad - frontporch,  fill,   'post');
    	d = wavecopy('d',   cd,   s,   i);
    	pad = ws - size(d);            frontporch = round(pad  /  2);
    	d = padarray(d,   frontporch,   fill,   'pre');
    	d = padarray(d,   pad - frontporch,   fill,   'post');
    % Add  1   pixel white border.
    	switch  lower(border)
    		case   'append'
    			w =  padarray(w,   [1   1],   1,   'post');
    			h = padarray(h,   [1   0],   1,   'post');
    			v = padarray(v,   [0  1],   1,   'post');
    		case   'absorb'
    			w(:,   end)   =  1;         w(end,   :)   =  1;
    			h(end,   :)   =  1;          v(:,   end)   =  1;
    		otherwise
    			error('Unrecognized BORDER parameter.');
    		end
    	w = [w h; v d];                                          % Concatenate coefs.
    end
    
    if nargout   == 0
    	imshow(w);                                               % Display  result.
    end

(2) waveback.m

func function [varargout] = waveback(c, s, varargin)
    %WAVEBACK Performs a multi-level two-dimensional inverse FWT.
    %   [VARARGOUT] = WAVEBACK(C, S, VARARGIN) computes a 2D N-level
    %   partial or complete wavelet reconstruction of decomposition
    %   structure [C, S]. 
    %
    %   SYNTAX:
    %   Y = WAVEBACK(C, S, 'WNAME');  Output inverse FWT matrix Y 
    %   Y = WAVEBACK(C, S, LR, HR);   using lowpass and highpass 
    %                                 reconstruction filters (LR and 
    %                                 HR) or filters obtained by 
    %                                 calling WAVEFILTER with 'WNAME'.
    %
    %   [NC, NS] = WAVEBACK(C, S, 'WNAME', N);  Output new wavelet 
    %   [NC, NS] = WAVEBACK(C, S, LR, HR, N);   decomposition structure
    %                                           [NC, NS] after N step 
    %                                           reconstruction.
    %
    %   See also WAVEFAST and WAVEFILTER.
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
    %   $Revision: 1.5 $  $Date: 2003/10/13 01:29:36 $
    
    % Check the input and output arguments for reasonableness.
    error(nargchk(3, 5, nargin));
    error(nargchk(1, 2, nargout));
       
    if (ndims(c) ~= 2) | (size(c, 1) ~= 1)
      error('C must be a row vector.');   
    end
      
    if (ndims(s) ~= 2) | ~isreal(s) | ~isnumeric(s) | (size(s,2) ~= 2)
      error('S must be a real, numeric two-column array.');   
    end
      
    elements = prod(s, 2);
    if (length(c) < elements(end)) | ...
          ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
      error(['[C S] must be a standard wavelet ' ...
             'decomposition structure.']); 
    end
    
    % Maximum levels in [C, S].
    nmax = size(s, 1) - 2;        
    % Get third input parameter and init check flags.
    wname = varargin{1};  filterchk = 0;   nchk = 0;    
      
    switch nargin
    case 3
       if ischar(wname)   
          [lp, hp] = wavefilter(wname, 'r');   n = nmax;
       else 
          error('Undefined filter.');  
       end
       if nargout ~= 1 
          error('Wrong number of output arguments.');  
       end
    case 4
       if ischar(wname)
          [lp, hp] = wavefilter(wname, 'r');   
          n = varargin{2};    nchk = 1;
       else
          lp = varargin{1};   hp = varargin{2};   
          filterchk = 1;   n = nmax;
          if nargout ~= 1 
             error('Wrong number of output arguments.');  
          end
       end
    case 5
        lp = varargin{1};   hp = varargin{2};   filterchk = 1;
        n = varargin{3};    nchk = 1;
    otherwise
        error('Improper number of input arguments.');     
    end
      
    fl = length(lp);
    if filterchk                                        % Check filters.
      if (ndims(lp) ~= 2) | ~isreal(lp) | ~isnumeric(lp) ...
            | (ndims(hp) ~= 2) | ~isreal(hp) | ~isnumeric(hp) ...
            | (fl ~= length(hp)) | rem(fl, 2) ~= 0
          error(['LP and HP must be even and equal length real, ' ...
                 'numeric filter vectors.']); 
      end
    end
    
    if nchk & (~isnumeric(n) | ~isreal(n))          % Check scale N.
        error('N must be a real numeric.'); 
    end
    if (n > nmax) | (n < 1)
       error('Invalid number (N) of reconstructions requested.');    
    end
    if (n ~= nmax) & (nargout ~= 2) 
       error('Not enough output arguments.'); 
    end
      
    nc = c;    ns = s;    nnmax = nmax;             % Init decomposition.
    for i = 1:n
       % Compute a new approximation.
       a = symconvup(wavecopy('a', nc, ns), lp, lp, fl, ns(3, :)) + ...
           symconvup(wavecopy('h', nc, ns, nnmax), ...
                     hp, lp, fl, ns(3, :)) + ...
           symconvup(wavecopy('v', nc, ns, nnmax), ...
                     lp, hp, fl, ns(3, :)) + ...
           symconvup(wavecopy('d', nc, ns, nnmax), ...
                     hp, hp, fl, ns(3, :));
          
        % Update decomposition.
        nc = nc(4 * prod(ns(1, :)) + 1:end);     nc = [a(:)' nc];
        ns = ns(3:end, :);                       ns = [ns(1, :); ns];
        nnmax = size(ns, 1) - 2;
    end
    
    % For complete reconstructions, reformat output as 2-D.
    if nargout == 1
        a = nc;   nc = repmat(0, ns(1, :));     nc(:) = a;    
    end
      
    varargout{1} = nc;
    if nargout == 2 
       varargout{2} = ns;  
    end
    
    %-------------------------------------------------------------------%
    function z = symconvup(x, f1, f2, fln, keep)
    % Upsample rows and convolve columns with f1; upsample columns and
    % convolve rows with f2; then extract center assuming symmetrical
    % extension.
    
    y = zeros([2 1] .* size(x));      y(1:2:end, :) = x;
    y = conv2(y, f1');
    z = zeros([1 2] .* size(y));      z(:, 1:2:end) = y;
    z = conv2(z, f2);
    z = z(fln - 1:fln + keep(1) - 2, fln - 1:fln + keep(2) - 2);

(3) wavecopy.m

func function y = wavecopy(type, c, s, n) 
    %WAVECOPY Fetches coefficients of a wavelet decomposition structure. 
    %   Y = WAVECOPY(TYPE, C, S, N) returns a coefficient array based on 
    %   TYPE and N.   
    % 
    %   INPUTS: 
    %     TYPE      Coefficient category 
    %     ------------------------------------- 
    %     'a'       Approximation coefficients 
    %     'h'       Horizontal details 
    %     'v'       Vertical details 
    %     'd'       Diagonal details 
    % 
    %     [C, S] is a wavelet data structure. 
    %     N specifies a decomposition level (ignored if TYPE = 'a'). 
    % 
    %   See also WAVEWORK, WAVECUT, and WAVEPASTE. 
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins 
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004 
    %   $Revision: 1.4 $  $Date: 2003/10/13 01:20:41 $ 
    
    error(nargchk(3, 4, nargin)); 
    if nargin == 4    
       y = wavework('copy', type, c, s, n); 
    else   
       y = wavework('copy', type, c, s);   
    end

(4) wavecut.m

func function [nc, y] = wavecut(type, c, s, n)
    %WAVECUT Zeroes coefficients in a wavelet decomposition structure.
    %   [NC, Y] = WAVECUT(TYPE, C, S, N) returns a new decomposition
    %   vector whose detail or approximation coefficients (based on TYPE
    %   and N) have been zeroed. The coefficients that were zeroed are
    %   returned in Y.
    %
    %   INPUTS:
    %     TYPE      Coefficient category
    %     -------------------------------------
    %     'a'       Approximation coefficients
    %     'h'       Horizontal details
    %     'v'       Vertical details
    %     'd'       Diagonal details
    %
    %     [C, S] is a wavelet data structure.
    %     N specifies a decomposition level (ignored if TYPE = 'a').
    %
    %   See also WAVEWORK, WAVECOPY, and WAVEPASTE.
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
    %   $Revision: 1.4 $  $Date: 2003/10/13 01:20:09 $
    
    error(nargchk(3, 4, nargin));
    if nargin == 4   
       [nc, y] = wavework('cut', type, c, s, n);
    else  
       [nc, y] = wavework('cut', type, c, s);  
    end

(5) wavefast.m

function [c, s] = wavefast(x, n, varargin) 
    %WAVEFAST Perform multi-level 2-dimensional fast wavelet transform. 
    %   [C, L] = WAVEFAST(X, N, LP, HP) performs a 2D N-level FWT of 
    %   image (or matrix) X with respect to decomposition filters LP and 
    %   HP. 
    % 
    %   [C, L] = WAVEFAST(X, N, WNAME) performs the same operation but 
    %   fetches filters LP and HP for wavelet WNAME using WAVEFILTER. 
    % 
    %   Scale parameter N must be less than or equal to log2 of the 
    %   maximum image dimension.  Filters LP and HP must be even. To 
    %   reduce border distortion, X is symmetrically extended. That is, 
    %   if X = [c1 c2 c3 ... cn] (in 1D), then its symmetric extension 
    %   would be [... c3 c2 c1 c1 c2 c3 ... cn cn cn-1 cn-2 ...]. 
    % 
    %   OUTPUTS: 
    %     Matrix C is a coefficient decomposition vector: 
    % 
    %      C = [ a(n) h(n) v(n) d(n) h(n-1) ... v(1) d(1) ] 
    % 
    %     where a, h, v, and d are columnwise vectors containing 
    %     approximation, horizontal, vertical, and diagonal coefficient 
    %     matrices, respectively.  C has 3n + 1 sections where n is the 
    %     number of wavelet decompositions.  
    % 
    %     Matrix S is an (n+2) x 2 bookkeeping matrix: 
    % 
    %      S = [ sa(n, :); sd(n, :); sd(n-1, :); ... ; sd(1, :); sx ] 
    % 
    %     where sa and sd are approximation and detail size entries. 
    % 
    %   See also WAVEBACK and WAVEFILTER. 
     
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins 
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004 
    %   $Revision: 1.5 $  $Date: 2003/10/13 01:14:17 $ 
     
    % Check the input arguments for reasonableness. 
    error(nargchk(3, 4, nargin)); 
     
    if nargin == 3 
       if ischar(varargin{1})    
          [lp, hp] = wavefilter(varargin{1}, 'd'); 
       else  
          error('Missing wavelet name.');    
       end 
    else 
          lp = varargin{1};     hp = varargin{2};    
    end 
     
    fl = length(lp);      sx = size(x); 
     
    if (ndims(x) ~= 2) | (min(sx) < 2) | ~isreal(x) | ~isnumeric(x) 
       error('X must be a real, numeric matrix.');      
    end 
       
    if (ndims(lp) ~= 2) | ~isreal(lp) | ~isnumeric(lp) ... 
           | (ndims(hp) ~= 2) | ~isreal(hp) | ~isnumeric(hp) ... 
           | (fl ~= length(hp)) | rem(fl, 2) ~= 0 
       error(['LP and HP must be even and equal length real, ' ... 
              'numeric filter vectors.']);  
    end 
       
    if ~isreal(n) | ~isnumeric(n) | (n < 1) | (n > log2(max(sx))) 
       error(['N must be a real scalar between 1 and ' ... 
              'log2(max(size((X))).']);     
    end 
       
    % Init the starting output data structures and initial approximation. 
    c = [];       s = sx;     app = double(x); 
     
    % For each decomposition ... 
    for i = 1:n 
       % Extend the approximation symmetrically. 
       [app, keep] = symextend(app, fl); 
        
       % Convolve rows with HP and downsample. Then convolve columns 
       % with HP and LP to get the diagonal and vertical coefficients. 
       rows = symconv(app, hp, 'row', fl, keep); 
       coefs = symconv(rows, hp, 'col', fl, keep); 
       c = [coefs(:)' c];    s = [size(coefs); s]; 
       coefs = symconv(rows, lp, 'col', fl, keep); 
       c = [coefs(:)' c]; 
        
       % Convolve rows with LP and downsample. Then convolve columns 
       % with HP and LP to get the horizontal and next approximation 
       % coeffcients. 
       rows = symconv(app, lp, 'row', fl, keep); 
       coefs = symconv(rows, hp, 'col', fl, keep); 
       c = [coefs(:)' c]; 
       app = symconv(rows, lp, 'col', fl, keep); 
    end 
     
    % Append final approximation structures. 
    c = [app(:)' c];       s = [size(app); s]; 
     
    %-------------------------------------------------------------------% 
    function [y, keep] = symextend(x, fl) 
    % Compute the number of coefficients to keep after convolution 
    % and downsampling. Then extend x in both dimensions. 
     
    keep = floor((fl + size(x) - 1) / 2); 
    y = padarray(x, [(fl - 1) (fl - 1)], 'symmetric', 'both'); 
     
    %-------------------------------------------------------------------% 
    function y = symconv(x, h, type, fl, keep) 
    % Convolve the rows or columns of x with h, downsample, 
    % and extract the center section since symmetrically extended. 
     
    if strcmp(type, 'row') 
       y = conv2(x, h); 
       y = y(:, 1:2:end); 
       y = y(:, fl / 2 + 1:fl / 2 + keep(2)); 
    else 
       y = conv2(x, h'); 
       y = y(1:2:end, :); 
       y = y(fl / 2 + 1:fl / 2 + keep(1), :); 
    end

(6) wavefilter.m

function [varargout]=wavefilter(wname,type) 
    error(nargchk(1,2,nargin)); 
    if(nargin==1 & nargout ~=4) | (nargin==2 & nargout~=2) 
        error('Invalid number of output arguments.'); 
    end 
    if nargin==1 & ~ischar(wname) 
        error('WNAME must be a string.'); 
    end 
    if nargin==2 & ~ischar(type) 
        error('TYPE must be a string.'); 
    end  
     
    switch lower(wname) 
        case {'haar','db1'} 
            ld=[1 1]/sqrt(2); hd=[-1 1]/sqrt(2); 
            lr=ld;            hr=-hd; 
        case 'db4' 
            ld=[-1.059740178499728e-002 3.288301166698295e-002 3.084138183598697e-002 -1.870348117188811e-001 ... 
                -2.798376941698385e-002 6.308807679295904e-001 7.148465705525415e-001 2.303778133088552e-001]; 
            t=[0:7]; 
            hd=ld; hd(end:-1:1)=cos(pi*t).*ld; 
            lr=ld; lr(end:-1:1)=ld; 
            hr=cos(pi*t).*ld; 
             
        case 'sym4' 
            ld=[-7.576571478927333e-002 -2.963552764599851e-002 4.976186676320155e-001 8.037387518059161e-001 ... 
                2.978577956052774e-001  -9.921954357684722e-002 -1.260396726203783e-002 3.222310060404270e-002]; 
            t=[0:7]; 
            hd=ld; hd(end:-1:1)=cos(pi*t).*ld; 
            lr=ld; lr(end:-1:1)=ld; 
            hr=cos(pi*t).*ld; 
        case 'bior6.8' 
            ld=[0 1.908831736481291e-003 -1.9142861290088767e-003 -1.699063986760234e-002 1.1934565279772926e-002 ... 
                4.973290349094079e-002 -7.726317316720414e-002 -9.405920349573646e-002 4.207962846098268e-001 ... 
                8.259229974584023e-001 4.207962846098268e-001 -9.405920349573646e-002 -7.726317316720414e-002 ... 
                4.973290349094079e-002 1.193456527972926e-002 -1.699063986760234e-002 -1.914286129088767e-003 ... 
                1.908831736481291e-003]; 
            hd=[0 0 0 1.442628250562444e-002 -1.446750489679015e-002 -7.872200106262882e-002 4.036797903033992e-002 ... 
                4.178491091502746e-001 -7.589077294536542e-001 4.178491091502746e-001 4.036797903033992e-002 ... 
                -7.872200106262882e-002 -1.446750489679015e-002 1.442628250562444e-002 0 0 0 0]; 
            t=[0:17]; 
            lr=cos(pi*(t+1)).*hd; 
            hr=cos(pi*t).*ld; 
             
        case 'jpeg9.7' 
            ld=[0 0.02674875741080976 -0.01686411844287495 -0.07822326652898785 0.2668641184428723 ... 
                0.6029490182363579 0.2668641184428723 -0.07822326652898785 -0.01686411844287495  ... 
                0.02674875741080976]; 
            hd=[0 -0.09127176311424948 0.05754352622849957 0.5912717631142470 -1.115087052456994 ... 
                0.5912717631142470 0.05754352622849957 -0.09127176311424948 0 0]; 
            t=[0:9]; 
            lr=cos(pi*(t+1)).*hd; 
            hr=cos(pi*t).*ld; 
        otherwise  
            error('Unrecongizable wavelet name (WNAME).'); 
    end 
    if(nargin==1) 
        varargout(1:4)={ld,hd,lr,hr}; 
    else  
        switch lower(type(1)) 
            case 'd' 
                varargout={ld,hd}; 
            case 'r' 
                varargout={lr,hr}; 
            otherwise  
                error('Unrecongizable filter TYPE.'); 
        end 
    end

(7) wavepaste.m

function nc = wavepaste(type, c, s, n, x)
    %WAVEPASTE Puts coefficients in a wavelet decomposition structure.
    %   NC = WAVEPASTE(TYPE, C, S, N, X) returns the new decomposition
    %   structure after pasting X into it based on TYPE and N.
    %
    %   INPUTS:
    %     TYPE      Coefficient category
    %     -------------------------------------
    %     'a'       Approximation coefficients
    %     'h'       Horizontal details
    %     'v'       Vertical details
    %     'd'       Diagonal details
    %
    %     [C, S] is a wavelet data structure.
    %     N specifies a decomposition level (Ignored if TYPE = 'a').
    %     X is a two-dimensional approximation or detail coefficient
    %       matrix whose dimensions are appropriate for decomposition
    %       level N.
    %
    %   See also WAVEWORK, WAVECUT, and WAVECOPY.
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
    %   $Revision: 1.4 $  $Date: 2003/10/13 01:21:13 $
    
    error(nargchk(5, 5, nargin))
    nc = wavework('paste', type, c, s, n, x);

(8) wavework.m

func  function [varargout] = wavework(opcode, type, c, s, n, x)
    %WAVEWORK is used to edit wavelet decomposition structures.
    %   [VARARGOUT] = WAVEWORK(OPCODE, TYPE, C, S, N, X) gets the
    %   coefficients specified by TYPE and N for access or modification
    %   based on OPCODE.
    %
    %   INPUTS:
    %     OPCODE      Operation to perform
    %     --------------------------------------------------------------
    %     'copy'      [varargout] = Y = requested (via TYPE and N)
    %                 coefficient matrix 
    %     'cut'       [varargout] = [NC, Y] = New decomposition vector
    %                 (with requested coefficient matrix zeroed) AND 
    %                 requested coefficient matrix 
    %     'paste'     [varargout] = [NC] = new decomposition vector with
    %                 coefficient matrix replaced by X
    %
    %     TYPE        Coefficient category
    %     --------------------------------------------
    %     'a'         Approximation coefficients
    %     'h'         Horizontal details
    %     'v'         Vertical details
    %     'd'         Diagonal details
    %
    %     [C, S] is a wavelet toolbox decomposition structure.
    %     N is a decomposition level (Ignored if TYPE = 'a').
    %     X is a two-dimensional coefficient matrix for pasting.
    %
    %   See also WAVECUT, WAVECOPY, and WAVEPASTE.
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
    %   $Revision: 1.6 $  $Date: 2003/11/21 15:00:18 $
    
    error(nargchk(4, 6, nargin));
    
    if (ndims(c) ~= 2) | (size(c, 1) ~= 1)
       error('C must be a row vector.');   
    end
      
    if (ndims(s) ~= 2) | ~isreal(s) | ~isnumeric(s) | (size(s, 2) ~= 2)
       error('S must be a real, numeric two-column array.');   
    end
      
    elements = prod(s, 2);                % Coefficient matrix elements.
    if (length(c) < elements(end)) | ...
          ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
        error(['[C S] must form a standard wavelet decomposition ' ...
               'structure.']); 
    end
    
    if strcmp(lower(opcode(1:3)), 'pas') & nargin < 6
       error('Not enough input arguments.');   
    end
      
    if nargin < 5
       n = 1;      % Default level is 1.
    end             
    nmax = size(s, 1) - 2;                % Maximum levels in [C, S].
    
    aflag = (lower(type(1)) == 'a');
    if ~aflag & (n > nmax)
       error('N exceeds the decompositions in [C, S].');    
    end
      
    switch lower(type(1))                 % Make pointers into C.
    case 'a'
       nindex = 1;
       start = 1;    stop = elements(1);    ntst = nmax;
    case  {'h', 'v', 'd'}
       switch type
       case 'h', offset = 0;     % Offset to details.
       case 'v', offset = 1;
       case 'd', offset = 2;
       end
       nindex = size(s, 1) - n;      % Index to detail info.
       start = elements(1) + 3 * sum(elements(2:nmax - n + 1)) + ...
               offset * elements(nindex) + 1;
       stop = start + elements(nindex) - 1;
       ntst = n;
    otherwise
       error('TYPE must begin with "a", "h", "v", or "d".');
    end
    
    switch lower(opcode)                   % Do requested action.
    case {'copy', 'cut'}
       y = repmat(0, s(nindex, :));
       y(:) = c(start:stop);    nc = c;
       if strcmp(lower(opcode(1:3)), 'cut')
          nc(start: stop) = 0; varargout = {nc, y};
       else 
          varargout = {y};    
       end
    case 'paste'
       if prod(size(x)) ~= elements(end - ntst)
          error('X is not sized for the requested paste.');
       else
          nc = c;   nc(start:stop) = x(:);   varargout = {nc};  
       end
    otherwise
       error('Unrecognized OPCODE.');
    end

(9) wavezero.m

function [nc, g8] = wavezero(c, s, l, wname)
    %WAVEZERO Zeroes wavelet transform detail coefficients. 
    %   [NC, G8] = WAVEZERO(C, S, L, WNAME) zeroes the level L detail
    %   coefficients in wavelet decomposition structure [C, S] and
    %   computes the resulting inverse transform with respect to WNAME
    %   wavelets.
    
    %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
    %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
    %   $Revision: 1.4 $  $Date: 2003/10/13 01:31:35 $
    
    [nc, foo] = wavecut('h', c, s, l);
    [nc, foo] = wavecut('v', nc, s, l);
    [nc, foo] = wavecut('d', nc, s, l);
    i = waveback(nc, s, wname);
    g8 = im2uint8(mat2gray(i));
    figure; imshow(g8);

(10) wavetest.m

f=rgb2gray(imread('laoshan.jpg'));
    imshow(f);
    [c,s]=wavefast(f,1,'sym4');
    figure; wave2gray(c,s,-6);
    [nc,y]=wavecut('a',c,s);
    figure, wave2gray(nc,s,-6);
    edge=abs(waveback(nc,s,'sym4'));
    figure,imshow(mat2gray(edge));

# 4.处理结果： #

![Center][]

[Center]: /images/20220731/12d783464c3646c88b49729768d663e5.png