gauswavf.m
上传用户:haiyisale
上传日期:2013-01-09
资源大小:3246k
文件大小:3k
源码类别:

波变换

开发平台:

Matlab

  1. function [psi,X] = gauswavf(LB,UB,N,NumWAW)
  2. %GAUSWAVF Gaussian wavelet.
  3. %   [PSI,X] = GAUSWAVF(LB,UB,N,P) returns values of the Pth 
  4. %   derivative of the Gaussian function F = Cp*exp(-x^2) on 
  5. %   an N point regular grid for the interval [LB,UB]. Cp is 
  6. %   such that the 2-norm of the Pth derivative of F is equal 
  7. %   to 1.  P can be integer values from 1 to 8.
  8. %
  9. %   Output arguments are the wavelet function PSI
  10. %   computed on the grid X.
  11. %   [PSI,X] = GAUSWAVF(LB,UB,N) is equivalent to
  12. %   [PSI,X] = GAUSWAVF(LB,UB,N,1).
  13. %
  14. %   These wavelets have an effective support of [-5 5].
  15. %
  16. %   ----------------------------------------------------
  17. %   If you have access to the Extended Symbolic Toolbox,
  18. %   you may specify a value of P > 8. 
  19. %   ----------------------------------------------------
  20. %
  21. %   See also CGAUSWAVF, WAVEINFO.
  22. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
  23. %   Last Revision: 14-May-2003.
  24. %   Copyright 1995-2004 The MathWorks, Inc.
  25. %   $Revision: 1.12.4.2 $  $Date: 2004/03/15 22:40:36 $
  26. % Check arguments.
  27. %-----------------
  28. nbIn = nargin;
  29. switch nbIn
  30.   case 0 ,
  31.    if ~exist('maple') , nmax = 8; else , nmax = 45; end
  32.    psi = nmax;
  33.    % psi contains the number max for Gaussian wavelets.
  34.    % This number depends of Symbolic Toolbox
  35.    
  36.   case 2 , error('Not enough input arguments.');
  37.       
  38.   case 3 , NumWAW = 1;
  39.       
  40.   case 4 ,
  41.    if ischar(NumWAW)
  42.        [fam,NumWAW] = wavemngr('fam_num',NumWAW);
  43.        NumWAW = wstr2num(NumWAW);
  44.    end
  45.    
  46.   otherwise  , 
  47.       error('Too many input arguments.');
  48. end
  49. if errargt(mfilename,NumWAW,'int') , error('*'); end
  50. % Compute values of the Gauss wavelet.
  51. X = linspace(LB,UB,N);  % wavelet support.
  52. if find(NumWAW==[1:8])
  53.   X2 = X.^2;
  54.   F0 = (2/pi)^(1/4)*exp(-X2);
  55. end
  56. switch NumWAW
  57.   case 1
  58.     psi = -2*X.*F0;
  59.   case 2    
  60.     psi = 2/(3^(1/2)) * (-1+2*X2).*F0;
  61.   case 3
  62.     psi = 4/(15^(1/2)) * X.* (3-2*X2).*F0;
  63.                 
  64.   case 4
  65.     psi = 4/(105^(1/2)) * (3-12*X2+4*X2.^2).*F0;
  66.   case 5
  67.     psi = 8/(3*(105^(1/2))) * X.* (-15+20*X2-4*X2.^2).*F0;  
  68.   case 6
  69.     psi = 8/(3*(1155^(1/2))) * (-15+90*X2-60*X2.^2+8*X2.^3).*F0; 
  70.   case 7
  71.     psi = 16/(3*(15015^(1/2))) *X.*(105-210*X2+84*X2.^2-8*X2.^3).*F0; 
  72.   case 8
  73.     psi = 16/(45*(1001^(1/2))) * (105-840*X2+840*X2.^2-224*X2.^3+16*X2.^4).*F0;
  74.   otherwise
  75.     if ~exist('maple')
  76.         msg = '*** The Extended Symbolic Toolbox is required ***';
  77.         errargt(mfilename,msg,'msg')
  78.         error(msg);
  79.     end
  80.     y = sym('t');
  81.     f = exp(-y.*y);
  82.     d = diff(f,NumWAW);
  83.     n = sqrt(int(d*d,-Inf,Inf));
  84.     d = d/n;
  85.     for j = 1:length(X)
  86.         t = X(j);
  87.         psi(j) = eval(d);
  88.     end
  89. end
  90. switch rem(NumWAW,4)
  91.     case {2,3} , psi = -psi;
  92. end