www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/+wavelet/+internal/morsewavelet.m
function [psihat,psi] = morsewavelet(omega,k,ga,be) % Find the time domain, psi, and frequency domain, psihat % k-th order Morse wavelet % N is the length of the signal % k is the order of the Morse wavelet, k=1,etc. % be is the \beta paramater, be = 8; % ga = the \gamma parameter, ga = 3; % [psi,psihat] = morsewavelet(omega,k,be,ga); % Algorithms due to JM Lilly % % Lilly, J. M. (2015), jLab: A data analysis package for Matlab, v. 1.6.1, % http://www.jmlilly.net/jmlsoft.html. fo = wavelet.internal.morsepeakfreq(ga,be); basicmorse =2*exp(-be.*log(fo)+fo.^ga+be.*log(abs(omega))-abs(omega).^ga).*(omega>0); %basicmorse = 2*exp(be/ga)*(ga/be)^(be/ga)*omega.^be.*exp(-abs(omega).^ga).*(omega>0); Akbg = morsenormconstant(be,ga,k); % Laguerre polynomials only needed for higher-order eigenfunctions % We are not supporting this in R2016b %Lkc = zeros(size(omega)); %Lkc(1:fix(N/2)) = wlaguerrepoly(2*abs(omega(1:fix(N/2))).^ga,k,(2*be+1)/ga-1); psihat = Akbg*basicmorse; % Only invert if the second output is requested if nargout == 2 psi = ifftshift(ifft(psihat)); end %------------------------------------------------------------- function Akbg = morsenormconstant(be,ga,k) % Returns the Morse wavelet normalization constant for the k-th order % Morse wavelet based on the parameters, \beta and \gamma % In R2016b, we are just supporting k=0, accordingly, this is just equal % to 1. r = (2*be+1)/ga; Akbg =sqrt(exp(gammaln(r)+gammaln(k+1)-gammaln(k+r))); %------------------------------------------------------------------- %------------------------------------------------------------------- %function Lkc = wlaguerrepoly(x,k,c) % Returns the k-th generalized Laguerre polynomial % c = r-1 % r = (2*be+1)/ga; May be better to just pass ga and be % the input is a vector of radian frequencies 2*omega^ga %Lkc = zeros(size(x)); % allocate an array of zeros the size of x %for m = 0:k % gammaexp = exp(gammaln(k+c+1)-gammaln(c+m+1)-gammaln(k-m+1)); % Lkc = Lkc+(-1)^m*gammaexp.*((x.^m)/gamma(m+1)); %end %-----------------------------------------------------------------