|||
关于瞬时频率估计,前面虽说暂时放一下,但心中始终还是念念不忘,因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多,方法很多,但问题也很多。在网上能找到的方法,简单归结如下:
相位差分法
相位建模法
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
(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de01016a4n.html
首发时间:2012-06-07 20:40:45)
1/0 | 闂備浇顕栭崜姘i幒妤婃晣闁跨噦鎷�:0 | 濠碘槅鍋撶徊楣冩偋閻樿违闁跨噦鎷� | 濠电偞鍨堕幐鎼佹晝閿濆洨鍗氶悗娑櫱滄禍婊堟煥閻曞倹瀚� | 闂佽崵濮撮幖顐︽偪閸ヮ灐褰掓晸閿燂拷 |
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2025-1-28 13:13
Powered by ScienceNet.cn
Copyright © 2007-2025 中国科学报社