关于瞬时频率估计，前面虽说暂时放一下，但心中始终还是念念不忘，因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多，方法很多，但问题也很多。在网上能找到的方法，简单归结如下： 相位差分法 相位建模法 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）
对希尔伯特－黄变换HHT、经验模式分解EMD及基本模式函数IMF概念稍微有点了解的人，都知道IMF必须满足的两个条件：(l)在整个数据长度中，极值点和过零点的数目必须相等或至多相差一个；(2)在任意数据点，局部极大值的包络和局部极小值的包络的均值必须为零。 对第一个条件，我几乎没加思考地理解成“所有的局部极大值都大于0，所有的局部极小值都小于0”，也就是说IMF中不应该有小于0的极大值，及大于0的极小值。我想很多人都会这样理解吧，这是很自然的。 可是今天在查看一列IMF时，突然发现了一个小于0的极大值。这使我感到非常诧异，因为EMD用的是当今被学界广泛接受的、法国人G. Rilling编的程序。现将过程记录如下。 为了简化和统一符号，将我以前小波消噪所得到的各信号序列改赋： TW4860_rd = TWpsy5den_med;%体温 MB2917_rd = MBpsy5den_low;%脉搏 GY2917_rd = GYpsy5den_low;%收缩压 DY2917_rd = DYpsy5den_low;%舒张压 JY2917_rd = JYpsy5den_low;%均压 CY2917_rd = CYpsy5den_low;%差压 作： imf_MB2917_rd = emd(MB2917_rd); x1=imf_MB2917_rd(1,:); = cenvelope(x1,1);%求上下包络线及其均值 env=env'; z=zeros(1,length(x1)) plot(x1,'r-') hold on plot(env) plot(envmoy,'k--') plot(z,'m-') xlim( ) 运行，得： 图29-1 axis( ) 运行，得： 图29-2 红线是分析信号，蓝线是下包络线，绿线是上包络线，虚线是包络线均值，紫线是0值线。那个小于0的极大值，看到了吧？ 进一步分析原因： = extr(x1);%求极值点与过零点 lmn = length(indmin)%极小值点数量 lmn = 7253 lmx = length(indmax)%极大值点数量 lmx = 7252 lzr = length(indzer)%过零点数量 lzr = 14504 说明|lmn+lmx-lzr| = 1是成立的。查看G. Rilling的程序也没有问题。再多思考一下，不禁哑然失笑。原来那个“在整个数据长度中，极值点和过零点的数目必须相等或至多相差一个”并不等价于“所有的局部极大值都大于0，所有的局部极小值都小于0”，它就是保留了在IMF中有“一个小于0的极大值或一个大于0的极小值”的可能性。 考查EMD发明人黄锷院士提出EMD的原论文："The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis"，从其中 4. Intrinsic mode functions ...... An intrinsic mode function (IMF) is a function that satisfies two conditions: （译文：一个 基本模式 函数（IMF） 是 满足这样 两个条件的函数） (1)in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one; and （译文： （1） 在 整个 数据 集 ， 极值点 数量 和 过零点数量 必须 相等或 至多 相差 一个 ，且） (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. （译文： （2） 在任何时间点 ， 由 局部极大值 定义 的包络线和由 局部极小值 定义的 包络线的 均值 为 零 。） 可知，我前面揭示的问题，产生的根源在于黄院士本人。他为什么要那样表述IMF？肯定不是为了更容易完成EMD而故意放宽IMF的条件(如果是为了这个目的，极值点过零点数量差应该跟数据长度挂钩才对，而且何必只限于一个呢？多几个岂不更容易完成EMD？)。 黄院士肯定是只考虑了数据两端的情况，因为当一列“所有的局部极大值都大于0，所有的局部极小值都小于0”的数据，(1)其最小的过零点序号小于最小的极值点序号、最大的过零点序号大于最大的极值点序号时，过零点数量比极值点数量多一个；(2)其最小的过零点序号大于最小的极值点序号、最大的过零点序号小于最大的极值点序号时，过零点数量比极值点数量少一个；(3)其余情况，过零点数量等于极值点数量。他也就不假思索地将“在整个数据长度中，极值点和过零点的数目必须相等或至多相差一个”视为与“所有的局部极大值都大于0，所有的局部极小值都小于0”等价了。跟我犯的错误是一样的啊！ 事实上，在上述“(1)其最小的过零点序号小于最小的极值点序号、最大的过零点序号大于最大的极值点序号”情况下，若数据中段存在一个小于0的极大值或一个大于0的极小值时，过零点数量比极值点数量少一个，数据是符合黄院士关于IMF的定义的。这就是问题产生的根本原因。 （本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de010167zd.html 首发时间： 2012-06-02 14:15:01 ）
标签: IMF
