关于瞬时频率估计，前面虽说暂时放一下，但心中始终还是念念不忘，因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多，方法很多，但问题也很多。在网上能找到的方法，简单归结如下：
相位差分法
相位建模法
Teager能量算子法
跨零点法
求根估计法
反余弦法
时频分布法（谱峰检测法？）
Shekel方法
Teager-Kaiser方法
解析信号法（HHT法）
比较起来，现在使用解析信号法(HHT)估计瞬时频率的人是最多的，大有“众望所归”的趋势。此方法就是本人在前几篇博文中使用的方法。可是，它存在的问题，我在博文《18、关于IMF分量瞬时频率跳变问题的研究》也说的很详细了。 既然现在最流行、影响最大的瞬时频率估计法都是这样，因此我有理由对其它所有瞬时频率估计方法，都感到没有信心。
imfMBrd=emd(MB2917_rd);
imfMBrd=imfMBrd';
lx=size(imfMBrd,1);
for i=1:size(imfMBrd,2)-1;
x=imfMBrd(:,i);
[fnorm,t,rejratio]=ifestar2(x);
M{i,1}=fnorm;
M{i,2}=t;
M{i,3}=rejratio;
end
fid1=figure(1);
set(fid1,'position',[50,10,700,1050])
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',[0.06+(j-1)/lj,0.005+(li-0.95*(0.999^i)*i)/li,0.92/lj,0.72/li]);
plot(M{n,2},M{n,1})
xlim([0 lx+100])
title(['fnorm',int2str(n)])
end
end
fid2=figure(2);
set(fid2,'position',[50,10,700,1050])
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',[0.06+(j-1)/lj,(li-0.94*(0.999^i)*i)/li,0.92/lj,0.65/li]);
plot(M{n+6,2},M{n+6,1})
xlim([0 lx+10])
title(['fnorm',int2str(n+6)])
end
end 运行得
以前多次讲过，现实信号的频率，一般情况下应该是连续、光滑的，所以我们的频率表示方法好不好，也应该以此为最高准则，是正频率就表示为正频率，是负频率就表示为负频率。除非能证明现实信号频率的确发生了跳变，否则频率表示当中出现了“跳变”，都是不对的。
fid3=figure(3);
set(fid3,'position',[50,10,700,1050])
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',[0.06+(j-1)/lj,0.005+(li-0.95*(0.999^i)*i)/li,0.92/lj,0.70/li]);
stem(M{n,2},M{n,1})
xlim([850 950])
title(['fnorm',int2str(n)])
end
end 
fid4=figure(4);
set(fid4,'position',[50,10,700,1050])
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',[0.06+(j-1)/lj,(li-0.94*(0.999^i)*i)/li,0.92/lj,0.60/li]);
stem(M{n+6,2},M{n+6,1})
xlim([7530 7580])
title(['fnorm',int2str(n+6)])
end
end 
附一：网上找到的“反余弦法瞬时频率估计”程序代码
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 =
33143 34589 34843 34901 34945 34961 34980
34993 34993 35001 35000 34999 35001
lt =
33143 34589 34843 34901 34945 34961 34980
34993 34993 35001 35000 34999 35001
rj =
0.9469 0.9882 0.9955 0.9971 0.9984 0.9989 0.9994
0.9998 1.0000 1.0000 0.9999 1.0000
function [f,a] = 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-
% [f,a] = 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-
%
%----- Get dimensions
[npts,nimf] = size(data); %----- Flip data if necessary
flipped=0;
if (npts < nimf)
flipped=1;
data=data';
[npts,nimf] = 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(acfreq<0)) = 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 [fnorm,t,rejratio]=ifestar2(x,t);
%IFESTAR2 Instantaneous frequency estimation using AR2 modelisation.
% [FNORM,T2,RATIO]=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 :
% [x,if]=fmlin(50,0.05,0.3,5); x=real(x); [if2,t]=ifestar2(x);
% plot(t,if(t),t,if2);
%
% N=1100; [deter,if]=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));
% [if2,t,ratio(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 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 if (nargin == 0),
error('At least one parameter required');
end;
[xrow,xcol] = size(x);
if (xcol~=1),
error('X must have only one column');
end
x = real(x); if (nargin == 1),
t=4:xrow;
end; [trow,tcol] = 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(den>0);
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本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
