www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/private/waveft.m
function [wft,frequencies] = waveft(WAV,omega,scales) %WAVEFT Wavelet Fourier transform. % WAVEFT computes the wavelet values in the frequency plane. % For a wavelet defined by WAV, WFT = WAVEFT(WAV,OMEGA,SCALES) % returns the wavelet Fourier transform WFT. % % WAV can be a string, structure or a cell array. % If WAV is a string, it contains the name of the wavelet % used for the analysis. % If WAV is a structure, WAV.name and WAV.param are respectively % the name of the wavelet and, if necessary, an associated parameter. % If WAV is a cell array, WAV{1} and WAV{2} contain the name of % the wavelet and an optional parameter. % OMEGA is a vector which contains the frequencies at which the % transform was computed. % SCALES is a vector which contains the scales used for the % wavelet analysis. % WFT is a matrix of size (NbScales x Nbfreq). % % Admissible wavelets are: % - MORLET wavelet (A) - 'morl': % PSI_HAT(s) = pi^(-1/4) * exp(-(s>0)*(s-s0).^2/2) * (s>0) % Parameter: s0, default s0 = 6. % % - MORLET wavelet (B) - 'morlex': (without Heaviside function) % PSI_HAT(s) = pi^(-1/4) * exp(-(s-s0).^2/2); % Parameter: s0, default s0 = 6. % % - MORLET wavelet (C) - 'morl0': (with exact zero mean value) % PSI_HAT(s) = pi^(-1/4) * [exp(-(s-s0).^2/2) - exp(-s0.^2/2)] % Parameter: s0, default s0 = 6. % % - MEXICAN wavelet - 'mexh': % PSI_HAT(s) = (1/gamma(2+0.5))*s^2 .* exp((-s.^2)/2) % (DOG wavelet with m = 2) % % - DOG wavelet - 'dog': m-th order Derivative Of Gaussian % PSI_HAT(s) = -(1i^m/sqrt(gamma(m+0.5)))*((s)^m).*exp((-s.^2)/2) % Parameter: m (order of derivative), default m = 2. The order % m must be even. % % - PAUL wavelet - 'paul': % PSI_HAT(s) = K*s^m.*exp(-s) % The constant K is such that: % K = (2^m)/sqrt(m*fact(2*m-1)) % Parameter: m, default m = 4. % - Bump wavelet: 'bump': % PSI_HAT(s) = exp(1-(1/((s-mu)^2./sigma^2))).*(abs((s-mu)/sigma)<1) % Parameters: mu,sigma. % default: mu=5, sigma = 0.6. % Allowable parameter ranges: 3 <= mu <= 7 % 0.1 < sigma <1.2 % Normalized scales for the bump wavelet must be between 1.8 and 512. % % PSI_HAT denotes the Fourier transform of PSI. % % See also CWTFT, CWT, ICWTFT. % M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 04-Mar-2010. % Copyright 1995-2015 The MathWorks, Inc. param = []; if isstruct(WAV) wname = WAV.name; param = WAV.param; elseif iscell(WAV) wname = WAV{1}; if length(WAV)>1 , param = WAV{2}; end else wname = WAV; end StpFrq = omega(2); NbFrq = length(omega); SqrtNbFrq = sqrt(NbFrq); cfsNORM = sqrt(StpFrq)*SqrtNbFrq; NbSc = length(scales); wft = zeros(NbSc,NbFrq); switch wname case 'morl' % MORLET (A) if isempty(param) , gC = 6; else gC = param; end mul = (pi^(-0.25))*cfsNORM; for jj = 1:NbSc expnt = -(scales(jj).*omega - gC).^2/2.*(omega > 0); wft(jj,:) = mul*sqrt(scales(jj))*exp(expnt).*(omega > 0); end FourierFactor = gC/(2*pi); frequencies = FourierFactor./scales; case 'morlex' % MORLET (B) if isempty(param) , gC = 6; else gC = param; end mul = (pi^(-0.25))*cfsNORM; for jj = 1:NbSc expnt = -(scales(jj).*omega - gC).^2/2; wft(jj,:) = mul*sqrt(scales(jj))*exp(expnt); end FourierFactor = gC/(2*pi); frequencies = FourierFactor./scales; case 'morl0' % MORLET (C) if isempty(param) , gC = pi*sqrt(2/log(2)); else gC = param; end mul = (pi^(-0.25))*cfsNORM; for jj = 1:NbSc expnt = -(scales(jj).*omega - gC).^2/2; correct = -(gC^2 +(scales(jj).*omega).^2)/2; wft(jj,:) = mul*sqrt(scales(jj))*(exp(expnt)-exp(correct)); end FourierFactor = gC/(2*pi); frequencies = FourierFactor./scales; case 'mexh' mul = sqrt(scales/gamma(2+0.5))*cfsNORM; for jj = 1:NbSc scapowered = (scales(jj).*omega); expnt = -(scapowered.^2)/2; wft(jj,:) = mul(jj)*(scapowered.^2).*exp(expnt); end FourierFactor = sqrt(2+1/2)/(2*pi); frequencies = FourierFactor./scales; case 'dog' if isempty(param) , m = 2; else m = param; validateattributes(m,{'numeric'},{'scalar','even'}); end mul = -((1i^m)/sqrt(gamma(m+0.5)))*sqrt(scales)*cfsNORM; for jj = 1:NbSc scapowered = (scales(jj).*omega); expnt = -(scapowered.^2)/2; wft(jj,:) = mul(jj).*(scapowered.^m).*exp(expnt); end FourierFactor = sqrt(m+1/2)/(2*pi); frequencies = FourierFactor./scales; case 'paul' if isempty(param) , m = 4; else m = param; end mul = sqrt(scales)*(2^m/sqrt(m*prod(2:(2*m-1))))*cfsNORM; for jj = 1:NbSc expnt = -(scales(jj).*omega).*(omega > 0); daughter = mul(jj)*((scales(jj).*omega).^m).*exp(expnt); wft(jj,:) = daughter.*(omega > 0); end FourierFactor = (2*m+1)/(4*pi); frequencies = FourierFactor./scales; case 'bump' if isempty(param) mu = 5; sigma = 0.6; else mu = param(1); sigma = param(2); end for jj = 1:NbSc w = (scales(jj)*omega-mu)./sigma; expnt = -1./(1-w.^2); daughter = exp(1)*exp(expnt).*(abs(w)<1-eps(1)); daughter(isnan(daughter)) = 0; wft(jj,:) = daughter; end FourierFactor = mu/(2*pi); frequencies = FourierFactor./scales; otherwise error(message('Wavelet:FunctionInput:InvWavNam')); end end % {function}