关于瞬时频率估计,前面虽说暂时放一下,但心中始终还是念念不忘,因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多,方法很多,但问题也很多。在网上能找到的方法,简单归结如下: 相位差分法 相位建模法 Teager能量算子法 跨零点法 求根估计法 反余弦法 时频分布法(谱峰检测法?) Shekel方法 Teager-Kaiser方法 解析信号法(HHT法) 其中有些只适用于单分量信号估计,有些可适用于多分量信号估计。由于单分量信号是多分量信号的特例,所以后者包括前者。 这些方法,有些我目前只能找到其名字,找不到原理论述;有些可以找到原理论述,但绝大多数方法都找不到现成的程序(收费网站内情况未知)。可以说,这些方法大多都只是用来写写文章、做做脑保健操而已,实际使用的人很少。如网上找到的一则“反余弦法瞬时频率估计”程序(见“附一”),它把分析信号的幅值都“规一化”了,因此它只能处理理想中的“平稳信号”,不能处理现实中的“非平稳信号”。 比较起来,现在使用解析信号法(HHT)估计瞬时频率的人是最多的,大有“众望所归”的趋势。此方法就是本人在前几篇博文中使用的方法。可是,它存在的问题,我在博文《18、关于IMF分量瞬时频率跳变问题的研究》也说的很详细了。既然现在最流行、影响最大的瞬时频率估计法都是这样,因此我有理由对其它所有瞬时频率估计方法,都感到没有信心。 解析信号法(HHT)估计瞬时频率用的是hhspectrum函数,hhspectrum函数是通过调用时频工具箱中的instfreq函数来进行瞬时频率估计的。时频工具箱中还有一个使用AR(2)模型进行瞬时频率估计的函数ifestar2(因为除却注释部分,程序正文很短,所以我把它也附在后面,方便阅读。见“附二”)。 现在来看看用ifestar2函数估计脉搏序列各IMF的瞬时频率: imfMBrd=emd(MB2917_rd); imfMBrd=imfMBrd'; lx=size(imfMBrd,1); for i=1:size(imfMBrd,2)-1; x=imfMBrd(:,i); =ifestar2(x); M{i,1}=fnorm; M{i,2}=t; M{i,3}=rejratio; end fid1=figure(1); set(fid1,'position', ) li=6; lj=1; for i=1:li for j=1:lj n=(i-1)*lj+j; s(n)=subplot(li,lj,n); set(s(n),'position', ); plot(M{n,2},M{n,1}) xlim( ) title( ) end end 运行 得 图30-1 fid2=figure(2); set(fid2,'position', ) li=7;% lj=1; for i=1:li for j=1:lj n=(i-1)*lj+j; s(n)=subplot(li,lj,n); set(s(n),'position', ); plot(M{n+6,2},M{n+6,1}) xlim( ) title( ) end end 运行得 图30-2 可见,IMF除了第10、第13分量以外,瞬时频率都有很严重的“跳变”,比HHT法更严重。从接近“0Hz”频率的地方跳到“0.5Hz”附近,如果都是这样还比较好区分,但那些比“0.25Hz”稍高的频率,你怎么知道它是否经过跳变呢?如果不是,那多大的频率是由负频跳变而来?界限在哪里?从物理意义来讲,正频率与负频率是截然不同的频率,一个是正向旋转所形成,一个是负向旋转所形成。频率表示方法,这是人类创造的理论,不管现实如何,用“0.49Hz”来表示“-0.01Hz”,这就叫“削足适履”。 以前多次讲过,现实信号的频率,一般情况下应该是连续、光滑的,所以我们的频率表示方法好不好,也应该以此为最高准则,是正频率就表示为正频率,是负频率就表示为负频率。除非能证明现实信号频率的确发生了跳变,否则频率表示当中出现了“跳变”,都是不对的。 下面将前面两张图的时间轴放大看看: fid3=figure(3); set(fid3,'position', ) li=6; lj=1; for i=1:li for j=1:lj n=(i-1)*lj+j; s(n)=subplot(li,lj,n); set(s(n),'position', ); stem(M{n,2},M{n,1}) xlim( ) title( ) end end 运行 得 图30-3 fid4=figure(4); set(fid4,'position', ) li=7; lj=1; for i=1:li for j=1:lj n=(i-1)*lj+j; s(n)=subplot(li,lj,n); set(s(n),'position', ); stem(M{n+6,2},M{n+6,1}) xlim( ) title( ) end end 图30-4 可以看到,图30-3的第1、第2分量频率、图30-4的第1分量频率,都有采样点上没有频率值的情况发生(不是“0Hz”情况哦)。其实每一个分量中都有这样的“瞬时频率无定义点”,只是因为在同一个绘图窗口将它们都绘出来不方便,所以就免了。 for i=1:size(imfMBrd,2)-1; lf(i)=length(M{i,1}) lt(i)=length(M{i,2}) rj(i)=M{i,3} end lf lt rj lf = Columns 1 through 7 33143345893484334901349453496134980 Columns 8 through 13 349933499335001350003499935001 lt = Columns 1 through 7 33143345893484334901349453496134980 Columns 8 through 13 349933499335001350003499935001 rj = Columns 1 through 8 0.94690.9882 0.9955 0.9971 0.9984 0.9989 0.9994 0.9998 Columns 9 through 13 0.99981.0000 1.0000 0.9999 1.0000 lf、lt是各IMF分量作了频率估算的频率点、采样点长度,rj是作了频率估算的采样点长度占整个采样点长度的比例。 这是由ifestar2函数中的indices = find(den0)造成的。与“反余弦法瞬时频率估计”一样,它能作频率估算的地方,它就作了估算;它无法作估算的地方,它就搁那不管了。这是什么态度嘛! 附一:网上找到的“反余弦法瞬时频率估计”程序代码 function = faacos(data, dt) % The function FAACOS generates an arccosine frequency and amplitude % of previously normalized data(n,k), where n is the length of the % time series and k is the number of IMFs. % % FAACOS finds the frequency by applying % the arccosine function to the normalized data and % checking the points where slope of arccosine phase changes. % Nature of the procedure suggests not to use the function % to process the residue component. % % Calling sequence- % = faacos(data, dt) % % Input- % data - 2-D matrix of IMF components % dt - time increment per point % Output- % f - 2-D matrix f(n,k) that specifies frequency % a - 2-D matrix a(n,k) that specifies amplitude % % Used by- % FA % Kenneth Arnold (NASA GSFC) Summer 2003, Initial %----- Get dimensions = size(data); %----- Flip data if necessary flipped=0; if (npts nimf) flipped=1; data=data'; = size(data); end %----- Input is normalized, so assume that amplitude is always 1 a = ones(npts,nimf); %----- Mark points that are above 1 as invalid (Not a Number) data(find(abs(data)1)) = NaN; %----- Process each IMF for c=1:nimf %----- Compute "phase" by arccosine acphase = acos(data(:,c)); %----- Mark points where slope of arccosine phase changes as invalid for i=2:npts-1 prev = data(i-1,c); cur = data(i,c); next = data(i+1,c); if (prev cur cur next) | (prev cur cur next) acphase(i) = NaN; end end %----- Get phase differential frequency acfreq = abs(diff(acphase))/(2*pi*dt); %----- Mark points with negative frequency as invalid acfreq(find(acfreq0)) = NaN; %----- Fill in invalid points using a spline legit = find(~isnan(acfreq)); if (length(legit) npts) f(:,c) = spline(legit, acfreq(legit), 1:npts)'; else f(:,c) = acfreq; end end %----- Flip again if data was flipped at the beginning if (flipped) f=f'; a=a'; end end 附二 时频工具箱中瞬时频率估计函数ifestar2 function =ifestar2(x,t); %IFESTAR2 Instantaneous frequency estimation using AR2 modelisation. % =IFESTAR2(X,T) computes an estimate of the % instantaneous frequency of the real signal X at time %instant(s) T. The result FNORM lies between 0.0 and 0.5. This %estimate is based only on the 4 last signal points, and has %therefore an approximate delay of 2.5 points. % %X : real signal to be analyzed. %T : Time instants (must be greater than 4) %(default : 4:length(X)). %FNORM : Output (normalized) instantaneous frequency. %T2 : Time instants coresponding to FNORM. Since the %algorithm can not always give a value, T2 is %different of T. % RATIO : proportion of instants where the algorithm yields %an estimation % %Examples : % =fmlin(50,0.05,0.3,5); x=real(x); =ifestar2(x); % plot(t,if(t),t,if2); % % N=1100; =fmconst(N,0.05); deter=real(deter); % noise=randn(N,1); NbSNR=101; SNR=linspace(0,100,NbSNR); % for iSNR=1:NbSNR, % sig=sigmerge(deter,noise,SNR(iSNR)); % =ifestar2(sig); % EQM(iSNR)=norm(if(t)-if2)^2 / length(t) ; % end; % subplot(211); plot(SNR,-10.0 * log10(EQM)); grid; % xlabel('SNR'); ylabel('-10 log10(EQM)'); % subplot(212); plot(SNR,ratio); grid; % xlabel('SNR'); ylabel('ratio'); % % See also INSTFREQ, KAYTTH, SGRPDLAY. %F. Auger, April 1996. % This estimator is the causal version of the estimator called % "4 points Prony estimator" in the article "Instantaneous %frequency estimation using linear prediction with comparisons %to the dESAs", IEEE Signal Processing Letters, Vo 3, No 2, p %54-56, February 1996. %Copyright (c) 1996 by CNRS (France). % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA if (nargin == 0), error('At least one parameter required'); end; = size(x); if (xcol~=1), error('X must have only one column'); end x = real(x); if (nargin == 1), t=4:xrow; end; = size(t); if (trow~=1), error('T must have only one row'); elseif min(t)4, error('The smallest value of T must be greater than 4'); end; Kappa = x(t-1) .* x(t-2) - x(t ) .* x(t-3) ; psi1 = x(t-1) .* x(t-1) - x(t ) .* x(t-2) ; psi2 = x(t-2) .* x(t-2) - x(t-1) .* x(t-3) ; den = psi1 .* psi2 ; indices = find(den0); arg=0.5*Kappa(indices)./sqrt(den(indices)); indarg=find(abs(arg)1); arg(indarg)=sign(arg(indarg)); fnorm = acos(arg)/(2.0*pi); rejratio = length(indices)/length(t); t = t(indices); end (本文首发于: http://blog.sina.com.cn/s/blog_6ad0d3de01016a4n.html 首发时间:2012-06-07 20:40:45)