wtfilter.m

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

function [wtf] = wtfilter(wtfname, transform)
% wtfilter -- Define wavelet transform filter coefficients.
%
%****f* wmtsa.dwt/wtfilter
%
% NAME
%   wtfilter -- Define wavelet transform filter coefficients.
%
% SYNOPSIS
%  [wtf] = wtfilter(wtfname, transform)
%
% INPUTS
%   * wtfname    -- name of wavelet transform filter (string, case-insenstive).
%   * transform  -- name of wavelet transform  (string, case-insenstive).
%
% OUTPUTS
%   * wtf        -- wavelet tranform filter struct (wtf_s).
%
% SIDE EFFECTS
%
%
% DESCRIPTION
%   wtfilter returns a wtf_s struct with the wavelet (high-pass) and
%   scaling (low-pass) filter coefficients, and associated attributes. 
%
%   The wtf_s struct has fields:
%   * g         -- scaling (low-pass) filter coefficients (vector).
%   * h         -- wavelet (high-pass) filter coefficients (vector).
%   * L         -- filter length (= number of coefficients) (integer).
%   * Name      -- name of wavelet filter (character string).
%   * WTFclass  -- class of wavelet filters (character string).
%   * Transform -- name of transform (character string).
%
%   Typing wtfilter('list') displays a list of supported filters.
% 
%   Typing wtfilter('all') returns a struct array of wtf_s of all
%   supported filters.
%
%   The MODWT filter coefficients are calculated from the DWT filter 
%   coefficients:
%
%      ht = h / sqrt(2)
%      gt = g / sqrt(2)
%
%   The wavelet filter coefficients (h) are calculated from the scaling
%   filter coefficients via the QMF function (wmtsa_qmf).
% 
% USAGE
%  % Return a wtf struct.
%  wtf = wtfilter(wtfname, transform)   
%
%  % Display and return a list of available filter ames.
%  wtfnames = wtfilter('list')
%  
%  % Display and return a struct array of wtf_s for available filters.
%  wtfs = wtfilter('all', transform)
%
% ERRORS  
%   WMTSA:InvalidNumArguments
%   WMTSA:MissingRequiredArgument
%   WMTSA:InvalidArgumentValue
%
% EXAMPLE
%    wtf = wtfilter('LA8', 'modwt')
%    wtf = wtfilter('haar', 'dwt')
%
% ALGORITHM
%
%
% REFERENCES
%   Percival, D. B. and A. T. Walden (2000) Wavelet Methods for
%     Time Series Analysis. Cambridge: Cambridge University Press.
%
% SEE ALSO
%  wtf_s, wtf_qmf
%
% TOOLBOX
%   wmtsa/dwt
%
% CATEGORY
%   Filters: Filters
%

% AUTHOR
%   Charlie Cornish
%
% CREATION DATE
%   2005-02-28
%
% COPYRIGHT
%   (c) Charles R. Cornish 2005
%
% CREDITS
%   Based on the original function (myfilter.m) by Brandon Whitcher.
%
% REVISION
%   $Revision: 612 $
%
%***

% $Id: wtfilter.m 612 2005-10-28 21:42:24Z ccornish $

usage_str = ['Usage:  [wtf] = ', mfilename, ...
             ' (wtfname, transform)'];

ext_usage_str = [usage_str, '\n', ...
                 'Type ', mfilename, ...
                 '(''list'') for list of available filters.', '\n', ...
                 'Type ', mfilename, ...
                 '(''all'') for struct array of all available filters.', '\n' ...
                ];

%%   Check Input Arguments
error(nargerr(mfilename, nargin, [1:2], nargout, [0:1], ...
              1, sprintf(ext_usage_str), 'struct'));

%% Process 'list' option
if (strcmpi(wtfname, 'list'))
  wtfilters = available_wtfilters;
  wtfilter_names = {wtfilters.Name};
  disp('Wavelet transform filters are available for the following:');
  for ( i = 1:length(wtfilter_names))
    disp(wtfilter_names{i});
  end
  if (nargout > 0)
    wtf = wtfilter_names;
  end
  return
end

if (~exist('transform', 'var') || isempty(transform))
  error(['WMTSA:missingRequiredArgument'], ...
        encode_errmsg('WMTSA:missingRequiredArgument', wmtsa_err_table, 'transform'));
else
  if (isempty(strmatch(upper(transform), {'DWT', 'MODWT'}, 'exact')))
    error(['WMTSA:invalidArgumentValue'], ...
          encode_errmsg('WMTSA:invalidArgumentValue', wmtsa_err_table, 'transform', ...
                        transform));
  end
end

%% Process 'all' option
if (strcmpi(wtfname, 'all'))
  wtfilters = available_wtfilters;
  if (strcmpi(transform, 'MODWT'))
    for (i = 1:length(wtfilters))
      wtfilters(i).Transform = 'MODWT';
      wtfilters(i).g = wtfilters(i).g ./ sqrt(2);
      wtfilters(i).h = wtfilters(i).h ./ sqrt(2);
    end
  end
  wtf = wtfilters;
  return
end

hwtf = str2func(['wtf_', lower(wtfname)]);

try
  wtf = feval(hwtf);
catch
  error('WMTSA:invalidWaveletTransformFilter', ...
        ['Unrecognized or unsupported wavelet transform filter name ' , ...
         '(', wtfname, ').']);
end

if (strcmpi(transform, 'MODWT'))
  wtf.Transform = 'MODWT';
  wtf.g = wtf.g ./ sqrt(2);
  wtf.h = wtf.h ./ sqrt(2);
end

return

function wtfs = available_wtfilters
%   avaliable_filters -- Returns a struct array of the available filters.
  wtf_hlist = ...
  {
      @wtf_haar, ...
      @wtf_d4, ...
      @wtf_d6, ...
      @wtf_d8, ...
      @wtf_d10, ...
      @wtf_d12, ...
      @wtf_d14, ...
      @wtf_d16, ...
      @wtf_d18, ...
      @wtf_d20, ...
      @wtf_la8, ...
      @wtf_la10, ...
      @wtf_la12, ...
      @wtf_la14, ...
      @wtf_la16, ...
      @wtf_la18, ...
      @wtf_la20, ...
      @wtf_bl14, ...
      @wtf_bl18, ...
      @wtf_bl20, ...
      @wtf_c6, ...
      @wtf_c12, ...
      @wtf_c18, ...
      @wtf_c24
   };
  
  wtf = struct;
  for (i = 1:length(wtf_hlist))
    wtf = feval(wtf_hlist{i});
    wtfs(i) = wtf;
  end
  
return

% ******************************************************************************
% *** Extremal phase filters ***
% ******************************************************************************
function wtf = wtf_haar
  wtf.Name = 'Haar';
  wtf.g = [0.7071067811865475 0.7071067811865475];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 2;
  wtf.WTFClass = 'ExtremalPhase';
  wtf.Transform = 'DWT';
return
  
function wtf = wtf_d4
  wtf.Name = 'D4';
  wtf.g = [0.4829629131445341 0.8365163037378077 0.2241438680420134 ...
           -0.1294095225512603];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 4;
  wtf.WTFClass = 'ExtremalPhase';
  wtf.Transform = 'DWT';
  
function wtf = wtf_d6
  wtf.Name = 'D6';
  wtf.g = [0.3326705529500827 0.8068915093110928 0.4598775021184915 ...
           -0.1350110200102546 -0.0854412738820267 0.0352262918857096];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 6;
  wtf.WTFClass = 'ExtremalPhase';
  wtf.Transform = 'DWT';
return

function wtf = wtf_d8
  wtf.Name ='D8';
  wtf.g = [0.2303778133074431 0.7148465705484058 0.6308807679358788 ...
           -0.0279837694166834 -0.1870348117179132 0.0308413818353661 ...
           0.0328830116666778 -0.0105974017850021];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 8;
  wtf.WTFClass = 'ExtremalPhase';
  wtf.Transform = 'DWT';
return

function wtf = wtf_d10
  wtf.Name ='D10';
  wtf.g = [0.1601023979741930 0.6038292697971898 0.7243085284377729 ...
           0.1384281459013204 -0.2422948870663824 -0.0322448695846381 ...
           0.0775714938400459 -0.0062414902127983 -0.0125807519990820 ...
           0.0033357252854738];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 10;
  wtf.WTFClass = 'ExtremalPhase';
  wtf.Transform = 'DWT';
return

function wtf = wtf_d12
  wtf.Name ='D12';
  wtf.g = [0.1115407433501094 0.4946238903984530 0.7511339080210954 ...
           0.3152503517091980 -0.2262646939654399 -0.1297668675672624 ...
           0.0975016055873224 0.0275228655303053 -0.0315820393174862 ...
           0.0005538422011614 0.0047772575109455 -0.0010773010853085];
  wtf.h = wmtsa_qmf(wtf.g);
  wtf.L = 12;