瞬时频率估计方法

关于瞬时频率估计，前面虽说暂时放一下，但心中始终还是念念不忘，因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多，方法很多，但问题也很多。在网上能找到的方法，简单归结如下：

 

相位差分法

相位建模法

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);

[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

 

    运行 得

    图30-1

 

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

 

    运行得

    图30-2

 

    可见，IMF除了第10、第13分量以外，瞬时频率都有很严重的“跳变”，比HHT法更严重。从接近“0Hz”频率的地方跳到“0.5Hz”附近，如果都是这样还比较好区分，但那些比“0.25Hz”稍高的频率，你怎么知道它是否经过跳变呢？如果不是，那多大的频率是由负频跳变而来？界限在哪里？从物理意义来讲，正频率与负频率是截然不同的频率，一个是正向旋转所形成，一个是负向旋转所形成。频率表示方法，这是人类创造的理论，不管现实如何，用“0.49Hz”来表示“-0.01Hz”，这就叫“削足适履”。

 

    以前多次讲过，现实信号的频率，一般情况下应该是连续、光滑的，所以我们的频率表示方法好不好，也应该以此为最高准则，是正频率就表示为正频率，是负频率就表示为负频率。除非能证明现实信号频率的确发生了跳变，否则频率表示当中出现了“跳变”，都是不对的。
 

    下面将前面两张图的时间轴放大看看：

 

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

 

    运行 得

    图30-3

 

 

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 

    图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

     33143     34589     34843     34901     34945     34961     34980

  Columns 8 through 13

     34993     34993     35001     35000     34999     35001

lt =

  Columns 1 through 7

     33143     34589     34843     34901     34945     34961     34980

  Columns 8 through 13

     34993     34993     35001     35000     34999     35001

rj =

  Columns 1 through 8

   0.9469   0.9882   0.9955   0.9971   0.9984   0.9989   0.9994   0.9998

  Columns 9 through 13

   0.9998   1.0000   1.0000   0.9999   1.0000

    lf、lt是各IMF分量作了频率估算的频率点、采样点长度，rj是作了频率估算的采样点长度占整个采样点长度的比例。

 

    这是由ifestar2函数中的indices = find(den>0)造成的。与“反余弦法瞬时频率估计”一样，它能作频率估算的地方，它就作了估算；它无法作估算的地方，它就搁那不管了。这是什么态度嘛！

 

 

 

 

 

 

 

附一：网上找到的“反余弦法瞬时频率估计”程序代码

 

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-
%  FA

% Kenneth Arnold (NASA GSFC)  Summer 2003, Initial

%----- 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  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;
[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