www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/modwt.m

    function w = modwt(x,varargin)
%MODWT Maximal overlap discrete wavelet transform.
%   W = MODWT(X) computes the maximal overlap discrete wavelet transform of
%   a 1-D real-valued, double-precision input signal, X. The signal can be
%   a row or column vector and must contain at least two samples. By
%   default, the maximal overlap discrete wavelet transform is computed
%   down to level floor(log2(length(X))) using the Daubechies
%   least-asymmetric wavelet with 4 vanishing moments ('sym4') and periodic
%   boundary handling. W is a LEV+1-by-N matrix containing the wavelet
%   coefficients and final-level scaling coefficients. LEV is the level of
%   the MODWT. The m-th row of the matrix contains the wavelet (detail)
%   coefficients for scale 2^m. The LEV+1-th row of the matrix contains the
%   scaling coefficients for scale 2^LEV.
%
%   W = MODWT(X,WNAME) computes the MODWT using the wavelet, WNAME. WNAME
%   is a character vector denoting the name of an orthogonal wavelet.
%   Orthogonal wavelets are designated as type 1 wavelets in the wavelet
%   manager. Valid built-in orthogonal wavelet families begin with
%   'haar','dbN', 'fkN', 'coifN', or 'symN' where N is the number of
%   vanishing moments for all families except 'fk'. For 'fk', N is the
%   number of filter coefficients. You can determine valid values for N by
%   using waveinfo. For example, waveinfo('db'). You can check if your
%   wavelet is orthogonal by using wavemngr('type',wname) to see if a 1 is
%   returned. For example, wavemngr('type','db2').
%
%   W = MODWT(X,Lo,Hi) computes the maximal overlap discrete wavelet
%   transform using the scaling filter, Lo, and the wavelet filter, Hi. Lo
%   and Hi are even-length row or column vectors. These filters must
%   satisfy the conditions for an orthogonal wavelet. You cannot specify
%   both WNAME and a filter pair, Lo and Hi. 
%
%   W = MODWT(...,LEV) computes the maximal overlap discrete wavelet
%   transform down to the level, LEV. LEV is a positive integer that cannot
%   exceed floor(log2(length(X))). If unspecified, LEV defaults to
%   floor(log2(length(X))).
%
%   W = MODWT(...,'reflection') uses reflection boundary handling by
%   extending the signal symmetrically at the right boundary to twice the
%   signal length,[x flip(x)], before computing the wavelet transform. The
%   number of wavelet and scaling coefficients returned are twice the
%   length of the input signal. By default, the signal is extended
%   periodically. You must enter the entire character vector 'reflection'.
%   If you added a wavelet named 'reflection' using the wavelet manager,
%   you must rename that wavelet prior to using this option. 'reflection'
%   may be placed in any position in the input argument list after X.
%
%
%   % Example 1:
%   %   Obtain the maximal overlap discrete wavelet transform of the Nile 
%   %   river minimum water level data. The data is 663 samples in length 
%   %   sampled yearly. Use the Haar wavelet and transform the data down 
%   %   to level 8. Plot the level-3 wavelet coefficients.
%    
%   load nileriverminima;
%   w = modwt(nileriverminima,'haar',8);
%   plot(w(3,:)); title('Level-3 Wavelet Coefficients');
%
%   % Example 2:
%   %   Check that the maximal overlap discrete wavelet transform 
%   %   partitions the variance of the signal by scale.
%
%   load noisdopp;
%   [~,~,Lo,Hi] = wfilters('sym8');
%   w = modwt(noisdopp,Lo,Hi,10);
%   varbylev = var(w,1,2);
%   sum(varbylev)
%   var(noisdopp,1)
%
%   See also IMODWT MODWTMRA MODWTCORR MODWTXCORR


% Check number of input arguments
narginchk(1,5);

% Validate that data is real, 1-D double-precision
% with no NaNs or Infs

validateattributes(x,{'double'},{'real','nonnan','finite'});

    

%Input must be 1-D
if (~isrow(x) && ~iscolumn(x))
    error(message('Wavelet:modwt:OneD_Input'));
end

%Input must contain at least two samples
if (numel(x)<2)
    error(message('Wavelet:modwt:LenTwo'));
end


% Convert data to row vector
x = x(:)';
% Record original data length
datalength = length(x);

%Parse input arguments
params = parseinputs(datalength,varargin{:});

%Check that the level of the transform does not exceed floor(log2(numel(x))
J = params.J;
Jmax = floor(log2(datalength));
if (J <= 0) || (J > Jmax) || (J ~= fix(J))
    error(message('Wavelet:modwt:MRALevel'));
end

boundary = params.boundary;
if (~isempty(boundary) && ~strcmpi(boundary,'reflection'))
    error(message('Wavelet:modwt:Invalid_Boundary'));
end

% increase signal length if 'reflection' is specified
if strcmpi(boundary,'reflection')
    x = [x flip(x)];
end

% obtain new signal length if needed
siglen = length(x);
Nrep = siglen;


% If wavelet specified as a string, ensure that wavelet is orthogonal
if (isfield(params,'wname') && ~isfield(params,'Lo'))
    [~,~,Lo,Hi] = wfilters(params.wname);
    wtype = wavemngr('type',params.wname);
    if (wtype ~= 1)
        error(message('Wavelet:modwt:Orth_Filt'));
    end
end

%If scaling and wavelet filters are specified as vectors, ensure they
%satisfy the orthogonality conditions

if (isfield(params,'Lo') && ~isfield(params,'wname'))
    filtercheck = CheckFilter(params.Lo,params.Hi);
    if ~filtercheck
        error(message('Wavelet:modwt:Orth_Filt'));
    end
    Lo = params.Lo;
    Hi = params.Hi;
end

% Scale the scaling and wavelet filters for the MODWT
Lo = Lo./sqrt(2);
Hi = Hi./sqrt(2);

% Ensure Lo and Hi are row vectors
Lo = Lo(:)';
Hi = Hi(:)';

% If the signal length is less than the filter length, need to 
% periodize the signal in order to use the DFT algorithm

if (siglen<numel(Lo))
    x = [x repmat(x,1,ceil(numel(Lo)-siglen))];
    Nrep = numel(x);
end

% Allocate coefficient array
w = zeros(J+1,Nrep);

% Obtain the DFT of the filters
G = fft(Lo,Nrep);
H = fft(Hi,Nrep);

%Obtain the DFT of the data
Vhat = fft(x);

% Main MODWT algorithm
for jj = 1:J
    [Vhat,What] = modwtdec(Vhat,G,H,jj);
    w(jj,:) = ifft(What);
end

w(J+1,:) = ifft(Vhat);

% Truncate data to length of boundary condition
w = w(:,1:siglen);


%----------------------------------------------------------------------
function [Vhat,What] = modwtdec(X,G,H,J)
% [Vhat,What] = modwtfft(X,G,H,J)

N = length(X);
upfactor = 2^(J-1);
Gup = G(1+mod(upfactor*(0:N-1),N));
Hup = H(1+mod(upfactor*(0:N-1),N));
Vhat = Gup.*X;
What = Hup.*X;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = CheckFilter(Lo,Hi)
% For a user-supplied scaling and wavelet filter, check that
% both correspond to an orthogonal wavelet
Lo = Lo(:);
Hi = Hi(:);
Lscaling = length(Lo);
Lwavelet = length(Hi);
evenlengthLo = 1-rem(Lscaling,2);
evenlengthHi = 1-rem(Lwavelet,2);
if all([evenlengthLo evenlengthHi])
    evenlength = 1;
else
    evenlength = 0;
end
 
if (Lscaling ~= Lwavelet)
    equallen = 0;
else
    equallen = 1;
end

normLo = norm(Lo,2);
sumLo = sum(Lo);
normHi = norm(Hi,2);
sumHi = sum(Hi);
tol = 1e-7;

if (abs(normLo-1) > tol && abs(normHi -1) > tol)
    unitnorm = 0;
else
    unitnorm = 1;
     
end

if (abs(sumLo-sqrt(2)) > tol && abs(sumHi) > tol)
    sumfilters = 0;
else
    sumfilters = 1;
end


zeroevenlags = 1;

if (Lscaling > 2)
    L = Lscaling;
    xcorrHi = conv(Hi,flipud(Hi));
    xcorrLo = conv(Lo,flipud(Lo));
    xcorrLo = xcorrLo(L+2:2:end);
    xcorrHi = xcorrHi(L+2:2:end);
    zeroevenlagsLo = 1-any(abs(xcorrLo>tol));
    zeroevenlagsHi = 1-any(abs(xcorrHi>tol));
    zeroevenlags = max(zeroevenlagsLo,zeroevenlagsHi);
end
out = all([evenlength equallen unitnorm sumfilters zeroevenlags]);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params = parseinputs(siglen,varargin)
 % Parse varargin and check for valid inputs
 
% Assign defaults 
params.boundary = [];
params.J = floor(log2(siglen));
params.wname = 'sym4';
 
% Check for 'reflection' boundary      
tfbound = strcmpi(varargin,'reflection');
 
% Determine if 'reflection' boundary is specified
if any(tfbound)
    params.boundary = varargin{tfbound>0};
    varargin(tfbound>0) = [];
end
 
 % If boundary is the only input in addition to the data, return with
 % defaults 
if isempty(varargin)
    return;
end
   
 % Only remaining char variable must be wavelet name
tfchar = cellfun(@ischar,varargin);
if (nnz(tfchar) == 1)
    params.wname = varargin{tfchar>0};
end
        
% Only scalar input must be the level
tfscalar = cellfun(@isscalar,varargin); 
 
% Check for numeric inputs 
tffilters = cellfun(@isnumeric,varargin);

% At most 3 numeric inputs are supported
if nnz(tffilters)>3
    error(message('Wavelet:modwt:Invalid_Numeric'));
end

% There's one numeric argument and it's not a wavelet level
if (nnz(tffilters)==1) && (nnz(tfscalar) == 0)
    error(message('Wavelet:FunctionInput:InvalidLoHiFilters'));
end

% If there are at least two numeric inputs, the first two must be the
% scaling and wavelet filters
if (nnz(tffilters)>1)
    idxFilt = find(tffilters,2,'first');
    params.Lo = varargin{idxFilt(1)};
    params.Hi = varargin{idxFilt(2)};
    params = rmfield(params,'wname');
    
    if (length(params.Lo) < 2 || length(params.Hi) < 2)
        error(message('Wavelet:modwt:Invalid_Filt_Length'));
    end
    
end
 
% Any scalar input must be the level
if any(tfscalar)
    params.J = varargin{tfscalar>0};
end
 
% If the user specifies a filter, use that instead of default wavelet
if (isfield(params,'Lo') && any(tfchar))
     error(message('Wavelet:FunctionInput:InvalidWavFilter'));
end
 
 
 
 
 
 

 
 
 
 
     
 
  
    
    
    
    
     
       

    
         
        



% [EOF] modwt.m