wavelet.m

几个关于多小波的程序

function [H,Ht] = wavelet(name,varargin)

% WAVELET -- return symbols of standard wavelets
%
%        [H,Htilde] = wavelet(name,p1,p2,...,'symbolic')
%
% The allowed types and their properties are listed below. Some of
% them take extra parameters P1, P2, ... as indicated. H and Htilde
% are the symbol and dual symbol.
%
% Normally, the coefficients are returned as floating point numbers,
% but if the last argument is the string 'symbolic', they come back as
% symbolic numbers whenever possible.
%
% Orthogonal Scalar Wavelets
%
%    'daub',p,type           Daubechies with p vanishing moments
%                            type 1 (default) are the original ones,
%                            type 2 are the least asymmetric ones
%    'coif6',support,type    coiflets of length 6
%                            support is [-4,1], [-3,2], [-2,3] or [-1,4] (default is [-2,3])
%                            there are two of each, type = 1 or 2 (default is 1)
%
% Biorthogonal Scalar Wavelets
%
%    'cohen',p,ptilde        Cohen-Daubechies-Feauveau
%    'daub79'                Daubechies 7/9 pair
%
% Orthogonal Multiwavelets
%
%    'bat'                   balanced Lebrun-Vetterli BAT O1
%    'cl2',t                 Chui-Lian 2; default t = -sqrt(7)/4
%    'cl3'                   Chui-Lian 3
%    'daubbal',p,type        balanced multiwavelet version of Daubechies wavelets
%    'dghm'                  Donovan-Geronimo-Hardin-Massopust
%    'stt'                   Shen, Tan, and Tam
%
% Biorthogonal Multiwavelets
%
%     'bcl'                  Bacchelli-Cotronei-Lazzaro
%     'hc',type              Hermite cubics, 
%                            type = 'dahmen':   Dahmen's completion
%                                   'shortest:  shortest completion
%                                   'smoothest' smoothest completion of length 5
%     'hm',s                 Hardin-Marasovich
%     'jrzb',s,t,lambda,mu   Jia-Riemenschneider-Zhou biorthogonal

% The properties are: m = dilation factor, r = multiplicity,
% l = length, p = approximation order, s = Sobolev smoothness
% alpha = Holder exponent

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

% check if last input argument is the string 'symbolic'

symbolic = 0;
nargin = length(varargin) + 1;
if (nargin >= 2)
    if strcmp(lower(varargin{end}),'symbolic')
	symbolic = 1;
	varargin = varargin(1:end-1);
	nargin = nargin - 1;
    end
end

switch lower(name)
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Orthogonal Scalar Wavelets %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 case 'daub'
%  Full form: 'daub',p,type
%  Daubechies with p vanishing moments
%  Type 1 (default) are the original ones,
%  Type 2 are the least asymmetric ones
%  For n <= 3, the two are identical
%  m=2, r=1, l=n, p=n, s=0.5000 for n=1
%                        1.0000       2
%                        1.4150       3
%                        1.7756       4
%                        2.0968       5
%
%  I. Daubechies, Orthonormal Bases of Compactly Supported Wavelets,
%        Comm. Pure Appl. Math., 41 (1988), pp. 909 - 996.
%
%  I. Daubechies, Ten Lectures on Wavelets, vol. 61 of CMNS-NSF
%        Regional Conference Series in Applied Mathematics,
%        SIAM, Philadelphia, 1992.

  H = daubechies(varargin{:});
  Ht = H;

 case 'coif6'
%  Full form: 'coif6', support, type
%  Coiflets of length 6
%  support is [-4,1], [-3,2], [-2,3] or [-1,4] (default is [-2,3])
%  type = 1 or 2 (default is 1)
%
%  I. Daubechies, Ten Lectures on Wavelets, vol. 61 of CMNS-NSF
%        Regional Conference Series in Applied Mathematics,
%        SIAM, Philadelphia, 1992.
%
%  I. Daubechies, Orthonormal Bases of Compactly Supported Wavelets II:
%        Variations on a Theme, SIAM J. Math. Anal. 24 (1993), pp. 499 - 519.
%
%  C. S. Burrus, R. A. Gopinath and H. Guo, Introduction to Wavelets
%        and Wavelet Transforms: A Primer, Prentice Hall, New Your, 1998.

  if (nargin < 2)
      support = [-2,3];
  else
      support = varargin{1};
  end
  if (nargin < 3)
      type = 1;
  else
      type = varargin{2};
  end
  
  switch type
   case 1
    switch support(1)
     case {-3, -2}
      hm2 = sym('(1+sqrt(7))/32');
      hm1 = sym('(5-sqrt(7))/32');
      h0  = sym('(14-2*sqrt(7))/32');
      h1  = sym('(14+2*sqrt(7))/32');
      h2  = sym('(1+sqrt(7))/32');
      h3  = sym('(-3-sqrt(7))/32');
     case {-4,-1}
      hm1 = sym('(9-sqrt(15))/32');
      h0  = sym('(13+sqrt(15))/32');
      h1  = sym('(6+2*sqrt(15))/32');
      h2  = sym('(6-2*sqrt(15))/32');
      h3  = sym('(1-sqrt(15))/32');
      h4  = sym('(-3+sqrt(15))/32');
     otherwise
      error('unrecognized support');
    end
   case 2
    switch support(1)
     case {-3, -2}
      hm2 = sym('(1-sqrt(7))/32');
      hm1 = sym('(5+sqrt(7))/32');
      h0  = sym('(14+2*sqrt(7))/32');
      h1  = sym('(14-2*sqrt(7))/32');
      h2  = sym('(1-sqrt(7))/32');
      h3  = sym('(-3+sqrt(7))/32');
     case {-4,-1}
      hm1 = sym('(9+sqrt(15))/32');
      h0  = sym('(13-sqrt(15))/32');
      h1  = sym('(6-2*sqrt(15))/32');
      h2  = sym('(6+2*sqrt(15))/32');
      h3  = sym('(1+sqrt(15))/32');
      h4  = sym('(-3-sqrt(15))/32');
     otherwise
      error('unrecognized support');
    end
   otherwise
    error('unrecognized type');
  end
  
  switch support(1)
   case -4
    H = mpoly({[h4;-hm1],[h3;h0],[h2;-h1],[h1;h2],[h0;-h3],[hm1;h4]},-4,'symbol',2,1); 
   case -3
    H = mpoly({[h3;-hm2],[h2;hm1],[h1;-h0],[h0;h1],[hm1;-h2],[hm2;h3]},-3,'symbol',2,1); 
   case -2
    H = mpoly({[hm2;h3],[hm1;-h2],[h0;h1],[h1;-h0],[h2;hm1],[h3;-hm2]},-2,'symbol',2,1); 
   case -1
    H = mpoly({[hm1;h4],[h0;-h3],[h1;h2],[h2;-h1],[h3;h0],[h4;-hm1]},-1,'symbol',2,1); 
  end
  
  Ht = H;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Biorthogonal Scalar Wavelets %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 'cohen'
%  Full form: 'cohen', p, ptilde
%  Cohen-Daubechies biorthogonal wavelets
%  m=2, r=1, l=??
%
%  A. Cohen, I. Daubechies and J.-C. Feauveau, Biorthogonal Bases of
%        Compactly Supported Wavelets, Comm. Pure Appl. Math., 45 (1992),
%        pp. 485 - 560.
%
%  I. Daubechies, Ten Lectures on Wavelets, vol. 61 of CMNS-NSF
%        Regional Conference Series in Applied Mathematics,
%        SIAM, Philadelphia, 1992.
  [H,Ht] = cohen(varargin{:});
 
 case 'daub79'
%  Full form: 'd79'
%  Daubechies 7/9 pair
%  m=2, r=1, l=7/9, p=4/4, s=2.1226/1.4100
%
%  I. Daubechies, Ten Lectures on Wavelets, vol. 61 of CMNS-NSF
%        Regional Conference Series in Applied Mathematics,
%        SIAM, Philadelphia, 1992.

  [H,Ht] = bidaubechies(7,9,[],varargin{:});
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Orthogonal Multiwavelets %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 'bat'
%  Full form: 'bat'
%  Lebrun-Vetterli balanced multiwavelet BAT O1
%
%  J. Lebrun and M. Vetterli, Higher-Order Balanced Multiwavelets:
%        Theory, Factorization, and Design, IEEE Trans. Signal Process.
%        49 (2001), pp. 1918 - 1930.

  h0 = sym('[0,2+sqrt(7);0,2-sqrt(7)]')/8;
  h1 = sym('[3,1;1,3]')/8;
  h2 = sym('[2-sqrt(7),0;2+sqrt(7),0]')/8;
  g0 = sym('[0,-2;0,1]')/(4*sqrt(2));
  g1 = sym('[2,2;-sqrt(7),sqrt(7)]')/(4*sqrt(2));
  g2 = sym('[-2,0;-1,0]')/(4*sqrt(2));
  H = mpoly({[h0;g0],[h1;g1],[h2;g2]},0,'symbol',2,2);
  Ht = H;
  
 case 'cl2'
%  Full form: 'cl2',t
%  Chui-Lian on [0,2], -1/sqrt(2) < t < -1/2
%  symmetric, m=2, r=2, l=3, p=1, alpha = -log_2 |1/2 + 2t|
%  default t = -sqrt(7)/4 with p=2, s=1.0546
%
%  C. K. Chui and J.-A. Lian, A Study of Orthonormal Multi-Wavelets,
%        Appl. Numer. Math. 20 (1996), pp. 273 - 298.
%
%  R.-Q. Jia, S. D. Riemenschneider, and D.-X. Zhou, Vector
%        Subdivision Schemes and Multiple Wavelets, Math. Comp. 67 (1998),
%        pp. 1533 - 1563

  if (nargin < 2)
      t = sym('-sqrt(7)/4');
  else
      t = varargin{1};
  end
  if (double(t) < -1/sqrt(2) | double(t) >= -1/2)
      error('t must be between -1/sqrt(2) and -1/2');
  end
  
  mu = sqrt(2 - 4*t^2);
  h0 = [1/4, 1/4; t/2, t/2];
  h1 = [1/2,0;0,mu/2];
  h2 = [1/4,-1/4;-t/2,t/2];
  g0 = [-1/4,-1/4; mu/4, mu/4];
  g1 = [1/2, 0; 0, -t];
  g2 = [-1/4, 1/4; -mu/4, mu/4];
  
  H = mpoly({[h0;g0],[h1;g1],[h2;g2]},0,'symbol',2,2); 
  Ht = H;

 case 'cl3'
%  Full form: 'cl3'
%  Chui-Lian on [0,3]
%  symmetric, m=2, r=2, l=3, p=3, s = 1.4408
%
%  C. K. Chui and J.-A. Lian, A Study of Orthonormal Multi-Wavelets,
%        Appl. Numer. Math. 20 (1996), pp. 273 - 298.
  h0 = sym('[10-3*sqrt(10), 5*sqrt(6)-2*sqrt(15);5*sqrt(6)-3*sqrt(15), 5-3*sqrt(10)]');
  g0 = sym('[5*sqrt(6)-2*sqrt(15),-10+3*sqrt(10);-5+3*sqrt(10),5*sqrt(6)-3*sqrt(15)]');
  h1 = sym('[30+3*sqrt(10),5*sqrt(6)-2*sqrt(15);-5*sqrt(6)-7*sqrt(15),15-3*sqrt(10)]');
  g1 = sym('[-5*sqrt(6)+2*sqrt(15),30+3*sqrt(10);15-3*sqrt(10),5*sqrt(6)+7*sqrt(15)]');
  h2 = sym('[30+3*sqrt(10),-5*sqrt(6)+2*sqrt(15);5*sqrt(6)+7*sqrt(15),15-3*sqrt(10)]');
  g2 = sym('[-5*sqrt(6)+2*sqrt(15),-30-3*sqrt(10);-15+3*sqrt(10),5*sqrt(6)+7*sqrt(15)]');
  h3 = sym('[10-3*sqrt(10),-5*sqrt(6)+2*sqrt(15);-5*sqrt(6)+3*sqrt(15),5-3*sqrt(10)]');
  g3 = sym('[5*sqrt(6)-2*sqrt(15),10-3*sqrt(10) ;5-3*sqrt(10),5*sqrt(6)-3*sqrt(15)]');
  H = mpoly({[h0;g0],[h1;g1],[h2;g2],[h3;g3]},0,'symbol',2,2); 
  H = H / 80;
  Ht = H;
  
 case 'daubbal'
%  Full form: 'daubbal', p
%  Balanced Daubechies multiwavelets
%
%  J. Lebrun and M. Vetterli, Higher-Order Balanced Multiwavelets:
%        Theory, Factorization, and Design, IEEE Trans. Signal Process.
%        49 (2001), pp. 1918 - 1930.

  S = daubechies(varargin{:});
  H = S(1,:);
  G = S(2,:);
  z2 = mpoly(1,2);
  H = [H;H*z2;G;G*z2];
  H.coef = reshape(H.coef,4,2,prod(size(H.coef))/8);
  H.r = 2;
  Ht = H;
  
 case 'dghm'
%  Full form: 'dghm'
%  Donovan-Geronimo-Hardin-Massopust
%  symmetric, m=2, r=2, l=4, p=2, s=1.5000
%
%  G. C. Donovan, J. S. Geronimo, D. P. hardin, and P. R. Massopust,
%        Construction of Orthogonal Wavelets Using Fractal Interpolation
%        Functions, SIAM J. Math. Anal. 27 (1996), pp. 1158 - 1192.

  h0 = sym('[3/10, 4/(5*sqrt(2)); -1/(20*sqrt(2)), -3/20]');
  h1 = sym('[3/10, 0; 9/(20*sqrt(2)), 1/2]');
  h2 = sym('[0, 0; 9/(20*sqrt(2)), -3/20]');
  h3 = sym('[0, 0; -1/(20*sqrt(2)),0]');
  g0 = sym('[-1/(20*sqrt(2)), -3/20; 1/20, 3/(10*sqrt(2))]');
  g1 = sym('[9/(20*sqrt(2)), -1/2; -9/20, 0]');
  g2 = sym('[9/(20*sqrt(2)), -3/20; 9/20, -3/(10*sqrt(2))]');
  g3 = sym('[-1/(20*sqrt(2)), 0; -1/20, 0]');
  
  H = mpoly({[h0;g0],[h1;g1],[h2;g2],[h3;g3]},0,'symbol',2,2); 
  Ht = H;
  
 case 'stt'
%  Full form: 'stt'
%  Shen, Tan, and Tham
%  symmetric, m=2, r=2, l=4, p=1, s=0.9919
%
%  L. Shen, H. H. Tan, and J. Y. Tham, Symmetric-Antisymmetric
%        Orthonormal Multiwavelets and Related Scalar Wavelets,
%        Appl. Comput. Harmon. Anal. 8 (2000), pp. 258 - 279.

  h0 = sym('[1/(16*(4+sqrt(15))), 1/16; 1/(16*(4+sqrt(15))), -1/16]');
  h1 = sym('[(31+8*sqrt(15))/(16*(4+sqrt(15))), 1/16; -(31+8*sqrt(15))/(16*(4+sqrt(15))), 1/16]');
  h2 = sym('[(31+8*sqrt(15))/(16*(4+sqrt(15))), -1/16; (31+8*sqrt(15))/(16*(4+sqrt(15))), 1/16]');
  h3 = sym('[1/(16*(4+sqrt(15))), -1/16;  -1/(16*(4+sqrt(15))), -1/16]');
  g0 = sym('[-1/16, 1/(16*(4+sqrt(15))); -1/16, -1/(16*(4+sqrt(15)))]');
  g1 = sym('[1/16, -(31+8*sqrt(15))/(16*(4+sqrt(15))); -1/16, -(31+8*sqrt(15))/(16*(4+sqrt(15)))]');
  g2 = sym('[1/16, (31+8*sqrt(15))/(16*(4+sqrt(15))); 1/16, -(31+8*sqrt(15))/(16*(4+sqrt(15)))]');
  g3 = sym('[-1/16, -1/(16*(4+sqrt(15))); 1/16, -1/(16*(4+sqrt(15)))]');
  
  H = mpoly({[h0;g0],[h1;g1],[h2;g2],[h3;g3]},0,'symbol',2,2); 
  Ht = H;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Biorthogonal Multiwavelets %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 'bcl'
%  Full form:  'bcl'
%  Bacchelli-Cotronei-Lazzaro
%
%  S. Bacchelli, M. Cotronei, and D. Lazzaro, An Algebraic Construction
%        of k-balanced Multiwavelets via the Lifting Scheme,
%        Numer. Algorithms 23 (2000), pp. 329 - 356.

  hm1 = sym('[2,1;-1,2]')/8;
  h0 = sym('[4,0;0,4]')/8;
  h1 = sym('[2,-1;1,2]')/8;
  gm1 = sym('[-9,-2;2,-9]')/64;
  g0 = sym('[-16,4;-4,-16]')/64;
  g1 = sym('[50,0;0,50]')/64;
  g2 = sym('[-16,-4;4,-16]')/64;
  g3 = sym('[-9,2;-2,-9]')/64;
  htm2 = sym('[-9,-2;2,-9]')/64;
  htm1 = sym('[16,-4;4,16]')/64;
  ht0 = sym('[50,0;0,50]')/64;
  ht1 = sym('[16,4;-4,16]')/64;
  ht2 = sym('[-9,2;-2,-9]')/64;
  gt0 = sym('[-2,-1;1,-2]')/8;
  gt1 = sym('[4,0;0,4]')/8;
  gt2 = sym('[-2,1;-1,-2]')/8;
  zero = sym('[0,0;0,0]');

  H = mpoly({[hm1;gm1],[h0;g0],[h1;g1],[zero;g2],[zero;g3]},-1,'symbol',2,2); 
  Ht = mpoly({[htm2;zero],[htm1;zero],[ht0;gt0],[ht1;gt1],[ht2;gt2]},-2,'symbol',2,2); 
  
 case 'hc'
%  Full form: 'hc',type
%  Hermite cubics
%  type = 'dahmen':     Dahmen's completion
%  type = 'shortest':   shortest completion
%                       m=2, r=2, l=3/3, p=4/0, s=2.5000/-1.5000
%         'smoothest':  smoothest completion of support length 4
%                       m=2, r=2, l=3/5, p=4/2, s=2.5000/0.8279

  switch varargin{1}
   case 'dahmen'
%  Hermite cubics, Dahmen's completion  
%
%  W. Dahmen, B. Han, R.-Q. Jia, and A. Kunoth, Biorthogonal
%        Multiwavelets on the Interval: Cubic Hermite Splines,
%        Constr. Approx. 16 (2000), pp. 221 - 259.

    hm1 = sym('[1/4,3/8;-1/16,-1/16]');
    h0 = sym('[1/2,0;0,1/4]');
    h1 = sym('[1/4,-3/8;1/16,-1/16]');
    htm2 = sym('[-7/128,-5/128;87/256,31/128]');
    htm1 = sym('[1/4,3/32;-99/64,-37/64]');
    ht0 = sym('[39/64,0;0,15/16]');
    ht1 = sym('[1/4,-3/32;99/64,-37/64]');
    ht2 = sym('[-7/128,5/128;-87/256,31/128]');
    H = mpoly({hm1,h0,h1},-1,'symbol',2,2);
    Ht = mpoly({htm2,htm1,ht0,ht1,ht2},-2,'symbol',2,2);

   case 'shortest'
%  Hermite cubics, shortest completion
%
%  V. Strela, private communication

    h0 = sym('[1/4,-1/8;3/16,-1/16]');
    h1 = sym('[1/2,0;0,1/4]');
    h2 = sym('[1/4,1/8;-3/16,-1/8]');
    g0 = sym('[1/2,0;0,1/2]');
    zero = sym('[0,0;0,0]');
    
    ht1 = sym('[1,0;0,2]');
    gtm1 = sym('[-1/2,3/4;-1/4,1/2]');
    gt0 = sym('[1,0;0,1]');
    gt1  = sym('[-1/2,-3/4;1/4,1/4]');
    H = mpoly({[h0;g0],[h1;zero],[h2;zero]},0,'symbol',2,2);
    Ht = mpoly({[zero;gtm1],[zero;gt0],[ht1;gt1]},-1,'symbol',2,2);
    
   case 'smoothest'
%  Hermite cubics, smoothest completion of support length 4
%
%  C. Heil, G. Strang, and V. Strela, Approximation by Translates
%        of Refinable Functions, Numer. Math. 73 (1996), pp. 75 - 94.

    hm1 = sym('[1/4,3/8;-1/16,-1/16]');
    h0 = sym('[1/2,0;0,1/4]');
    h1 = sym('[1/4,-3/8;1/16,-1/16]');
    gm1 = sym('[67/480,7/480;-95/1944,-1/324]');
    g0 = sym('[-1/2,-187/120;89/486,91/162]');
    g1 = sym('[173/240,0;0,13/9]');
    g2 = sym('[-1/2,187/120;-89/486,91/162]');
    g3 = sym('[67/480,-7/480;95/1944,-1/324]');
    zero = sym('[0,0;0,0]');
    
    htm2 = sym('[-73/1296, -77/1944; 773/2160,   3229/12960]');
    htm1 = sym('[1/4, 89/972; -187/120, -91/162]');
    ht0 = sym('[397/648, 0; 0, 6091/6480]');
    ht1 = sym('[1/4, -89/972; 187/120, -91/162]');
    ht2 = sym('[-73/1296, 77/1944; -773/2160,  3229/12960]');
    gt0 = sym('[-1/8, -1/16; 3/16, 1/16]');
    gt1 = sym('[1/4, 0; 0, 1/4]');
    gt2 = sym('[-1/8, 1/16; -3/16, 1/16]');
    
    H = mpoly({[hm1;gm1],[h0;g0],[h1;g1],[zero;g2],[zero;g3]},-1,'symbol',2,2);
    Ht = mpoly({[htm2;zero],[htm1;zero],[ht0;gt0],[ht1;gt1],[ht2;gt2]},-2,'symbol',2,2);
   
   otherwise
    error('unrecognized type');
  end
  
 case 'hm'
%  Full form: 'hm', s
%  Hardin-Marasovich, -1 < s < 1/7
%  m=2, r=2, l=4/4, p=2/2
%
%  D. P. Hardin and J. A. Marasovich, Biorthogonal Multiwavelets on [-1,1],
%        Appl. Comput. Harmon. Anal. 7 (1999), pp. 34 - 53.

  s = varargin{1};
  if (symbolic)
      s = sym(s);
  end
  st = (1+2*s)/(-2+5*s);
  alpha = 3*(1-s)*(1-s*st)/(4-s-st-2*s*st);
  gamma = sqrt(6*(4-s-st-2*s*st)/(7 -4*s - 4*st + s*st));
  delta = sqrt(12*(-1+st)*(-1+s)*(-1+s*st)/(-4+s+st+2*s*st));
  a = alpha*gamma*(1-2*alpha-2*s)/(2*delta);
  b = 1/2 - alpha;
  c = alpha*gamma*(3-2*alpha-2*s)/(2*delta);
  d = alpha + s;
  e = delta/gamma;
  alphat = 3*(1-st)*(1-s*st)/(4-s-st-2*s*st);
  gammat = gamma;
  deltat = delta;
  at = alphat*gammat*(1-2*alphat-2*st)/(2*deltat);
  bt = 1/2 - alphat;
  ct = alphat*gammat*(3-2*alphat-2*st)/(2*deltat);
  dt = alphat + st;
  et = deltat/gammat;

  h0 = [0,a;0,0];
  h1 = [b,c;0,0];
  h2 = [1,c;0,d];
  h3 = [b,a;e,d];