dwt.m

时间序列分析中很用的源码,书的原名为时间序列分析的小波方法.

function [WJ, VJ, att, NJ] = dwt(X, wtf, nlevels, boundary, varargin)
%   dwt -- Compute the (partial) discrete wavelet transform (DWT).
%
%****f* wmtsa.dwt/dwt
%
% NAME
%   dwt -- Compute the (partial) discrete wavelet transform (DWT).
%
% USAGE
%   [W, V, att, NJ] = dwt(X, [wtf], [nlevels], [boundary], [{opts}])
%
% INPUTS
%   * X          -- set of observations 
%                   (vector of length N or matrix of size N x Nchan)
%   * wtf        -- (optional) wavelet transform filter name or struct 
%                   (string, case-insensitve or wtf struct).
%                   Default:  'la8'
%   * nlevels    -- (optional) maximum level J0 (integer) 
%                   or method of calculating J0 (character string).
%                   Valid values: integer>0 or a valid method name
%                   Default:  'conservative'
%   * boundary   -- (optional) boundary conditions to use (character string)
%                   Valid values: 'circular' or 'reflection'
%                   Default: 'reflection'
%   * opts       -- (optional) Additional function options.
%
% OUTPUTS
%   * WJ         --  DWT wavelet coefficents 
%                    (cell array of length J0 of NJ(:,1) x NChan matrices).
%   * VJ         --  DWT scaling coefficents 
%                    (cell array of length J0 of NJ(:,2) x NChan matrices).
%   * att        --  structure containing DWT transform attributes.
%   * NJ         --  Jx2 matrix containing number of DWT wavelet (WJ)
%                    and scaling (VJ) coefficients at each level j.
%
% SIDE EFFECTS
%   1.  wavelet is a WMTSA-supported DWT wavelet filter; otherwise error.
%
%
% DESCRIPTION
%
%
% EXAMPLE
%
%
% WARNINGS
%
%
% ERRORS
%
%
% NOTES
%
%
% BUGS
%
%
% TODO
%
%
% ALGORITHM
%
%
% REFERENCES
%
%
% SEE ALSO
%
%

% AUTHOR
%   Charlie Cornish
%
% CREATION DATE
%
%
% COPYRIGHT
%
%
% CREDITS
%
%
% REVISION
%   $Revision: 630 $
%
%***

%   $Id: dwt.m 630 2006-05-02 20:47:17Z ccornish $

defaults.wtf = 'la8';
defaults.boundary = 'reflection';
defaults.nlevels  = 'conservative';
  
opts_defaults.RetainVJ = 0;

usage_str = ['Usage:  [WJ, VJ, att, NJ] = ', mfilename, ...
             '(X, [wtf], [nlevels], [boundary], [opts])'];

%%  Check input arguments and set defaults.
error(nargerr(mfilename, nargin, '1:', nargout, [0:4], 1, usage_str, 'struct'));

set_defaults(defaults);



%% Parse and check the options.
opts = parse_opts(varargin{:});
error(validate_opts(opts, opts_defaults, 'struct'));
opts = set_opts_defaults(opts, opts_defaults);


%% Get a valid wavelet transform filter coefficients struct.
if (ischar(wtf))
  try
    [wtf_s] = dwt_filter(wtf);
  catch
    rethrow(lasterror);
  end
elseif (iswtf(wtf))
  wtf_s = wtf;
else
  error('WMTSA:invalidWaveletTransformFilter', ...
        encode_errmsg('WMTSA:invalidWaveletTransformFilter', wmtsa_err_table, 'wtf'));
end

wtfname = wtf_s.Name;
g = wtf_s.g;
h = wtf_s.h;

% If a vector, make X a column vector
if (wmtsa_isvector(X, 'nonsingleton'))
  X = X(:);
end

%%  N     = length of original series
%%  NChan = number of channels for multi-variant dataset.
[N, NChan] = size(X);

%%  If nlevels is an integer > 0, set J0 = nlevels.
%%  otherwise, select J0 based on choice method specified.
if (isa(nlevels, 'char'))
  if (strcmp(nlevels, 'conservative'))
    J0 = modwt_choose_nlevels(nlevels, wtfname, N);
  else
    error('WMTSA:invalidNLevelsValue', ...
          encode_errmsg('WMTSA:invalidNLevelsValue', wmtsa_err_table));
  end
elseif (isnumeric(nlevels))
  if (nlevels > 0)
    J0 = nlevels;
  else
    error('WMTSA:negativeJ0', ...
          ['nlevels must be an integer greater than 0.']);
  end
else
  error('WMTSA:invalidNLevelsValue', ...
         encode_errmsg('WMTSA:invalidNLevelsValue', wmtsa_err_table));
end

if (J0 < 0)
  error('WMTSA:negativeJ0', ...
        ['J0 must be greater than 0.']);
elseif (2.^J0 > N)
  error('WMTSA:DWT:TooLargeJ0', 'Level (J0) too large for sample size for DWT');
end

if (~opts.RetainVJ && ...
    (log2(N) ~= floor(log2(N))))
  error('WMTSA:DWT:nonPowerOfTwoSampleSize', ...
        ['Sample size is not a multiple of a power of 2; ', ...
         'Use RetainVJ option for odd length sample sizes']);
end

%% Initialize the scale (Vin) for first level by setting it equal to X
%% using specified  boundary conditions
switch boundary
  case 'reflection'
   Xin = cat(1, X, flipdim(X, 1));
  case {'circular', 'periodic'}
   Xin = X;
  otherwise
   error('WMTSA:invalidBoundary', ['Invalid boundary method.']);
end


%% NW = length of the extended series = number of coefficients
NW = size(Xin, 1);

NJ = zeros(J0, 2);;
NJ(:,1) = floor(NW ./ 2.^[1:J0]);
WJ = cell(J0,NChan);
VJ = cell(J0,NChan);


%% Do the DWT
for (i = 1:NChan)
  Vin = Xin(:,i);
  for (j = 1:J0)
    [W_j, Vout] = dwtj(Vin, h, g);
    WJ{j,i} = W_j;
    if (j == J0)
      VJ{J0,i} = Vout(:);
      NJ(j,2) = length(Vout);
    elseif (opts.RetainVJ)
        if (mod(length(Vout), 2))
          Vin = Vout(1:length(Vout)-1);
          VJ{j,i} = Vout(end);
          NJ(j,2) = 1;
        else
          Vin = Vout;
          NJ(j,2) = 0;
        end
    else
      Vin = Vout;
      NJ(j,2) = 0;
    end
  end
end


%% Update attributes
att.Transform = 'DWT';
if (isstruct(wtf))
  att.WTF = wtf;
else
  att.WTF = wtfname;
end
att.N = N;
att.NW = NW;
att.J0 = J0;
att.NChan = NChan;
att.Boundary = boundary;
att.Aligned = 0;
att.RetainVJ = opts.RetainVJ;

return