www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/centfrq.m
function [freq,xval,recfreq] = centfrq(wname,iter,plotflag) %CENTFRQ Wavelet center frequency. % FREQ = CENTFRQ('wname') returns the center frequency in hertz % of the wavelet function 'wname' (see WAVEFUN). % % For FREQ = CENTFRQ('wname',ITER), ITER is the number % of iterations used by the WAVEFUN function to compute % the wavelet. ITER has a default value of 8 when not specified. % % [FREQ,XVAL,RECFREQ] = CENTFRQ('wname',ITER, 'plot') % returns in addition the associated center frequency based % approximation RECFREQ on the 2^ITER points grid XVAL % and plots the wavelet function and RECFREQ. % % See also SCAL2FRQ, WAVEFUN, WFILTERS. % M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 04-Mar-98. % Last Revision: 20-Oct-2014. % Copyright 1995-2014 The MathWorks, Inc. % Check arguments. narginchk(1,3) if nargin==1 iter = 8; end validateattributes(iter, {'numeric'}, ... {'integer', 'positive', 'scalar'}, 'centfrq', 'iter'); % Retrieve wavelet. wname = deblankl(wname); wtype = wavemngr('type',wname); switch wtype case 1 [~,psi,xval] = wavefun(wname,iter); case 2 [~,psi,~,~,xval] = wavefun(wname,iter); case 3 [~,psi,xval] = wavefun(wname,iter); case 4 [psi,xval] = wavefun(wname,iter); case 5 [psi,xval] = wavefun(wname,iter); end T = max(xval)-min(xval); % T is the size of the domain of psi. n = length(psi); psi = psi-mean(psi); % psi is numerically centered. psiFT = fft(psi); % computation of the modulus sp = (abs(psiFT)); % of the FT. % Compute arg max of the modulus of the FT (center frequency). Is_BIOR31 = isequal(wname,'bior3.1'); if ~Is_BIOR31 [vmax,indmax] = max(sp); if indmax > n/2 indmax = n-indmax+2; % indmax is always >= 2. end else [~,I,~] = localmax(sp(:)',1,false); indmax = I(1); % first local max vmax = sp(indmax); end per = T/(indmax-1); % period corresponding to the maximum. freq = 1/per; % associated frequency. if (nargin > 2) && strcmp(plotflag, 'plot') % Recontruct the signal from only the largest magnitude element of FFT. psiFT(sp<vmax) = 0; if Is_BIOR31 psiFT(sp>vmax) = 0; end recfreq = ifft(psiFT); recfreqreal = 0.75*max(abs(psi))*real(recfreq)/max(abs(recfreq)); recfreqimag = 0.75*max(abs(psi))*imag(recfreq)/max(abs(recfreq)); if wtype <= 4 % For real valued wavelets, show only the real part plot(xval,psi,'-b',... xval,recfreqreal,'-r'); title([getString(message('Wavelet:centfrq:PeriodLabel')) ': ' ... num2str(per) ' ' ... getString(message('Wavelet:centfrq:FreqLabel')) ': ' ... num2str(freq)]) legend({[getString(message('Wavelet:centfrq:WaveletLegend')) ... ': ' wname], ... getString(message('Wavelet:centfrq:ApproxLegend'))}) else % for complex wavelet without scale function (Type 5), show both % real and imaginary parts subplot(211) plot(xval,real(psi),'-b',... xval,recfreqreal,'-r'); title([getString(message('Wavelet:centfrq:PeriodLabel')) ': ' ... num2str(per) ' '... getString(message('Wavelet:centfrq:FreqLabel')) ': ' ... num2str(freq)]) legend({[getString(message('Wavelet:centfrq:WaveletLegend')) ... ': ' wname], ... getString(message('Wavelet:centfrq:ApproxLegend'))}) ylabel(getString(message('Wavelet:centfrq:RealLabel'))) subplot(212), plot(xval,imag(psi),'-b',... xval,recfreqimag ,'-r') legend({[getString(message('Wavelet:centfrq:WaveletLegend')) ... ': ' wname], ... getString(message('Wavelet:centfrq:ApproxLegend'))}) ylabel(getString(message('Wavelet:centfrq:ImagLabel'))) end else recfreq = []; end