瞬时频率估计方法
關于瞬時頻率估計，前面雖說暫時放一下，但心中始終還是念念不忘，因為這是一道繞不過去的坎。在網(wǎng)上多次搜索、了解其現(xiàn)狀。感覺是關注這件事的人很多，方法很多，但問題也很多。在網(wǎng)上能找到的方法，簡單歸結如下：
?
相位差分法
相位建模法
Teager能量算子法
跨零點法
求根估計法
反余弦法
時頻分布法（譜峰檢測法？）
Shekel方法
Teager-Kaiser方法
解析信號法（HHT法）
?
????其中有些只適用于單分量信號估計，有些可適用于多分量信號估計。由于單分量信號是多分量信號的特例，所以后者包括前者。
?
????這些方法，有些我目前只能找到其名字，找不到原理論述；有些可以找到原理論述，但絕大多數(shù)方法都找不到現(xiàn)成的程序（收費網(wǎng)站內(nèi)情況未知）。可以說，這些方法大多都只是用來寫寫文章、做做腦保健操而已，實際使用的人很少。如網(wǎng)上找到的一則“反余弦法瞬時頻率估計”程序（見“附一”），它把分析信號的幅值都“規(guī)一化”了，因此它只能處理理想中的“平穩(wěn)信號”，不能處理現(xiàn)實中的“非平穩(wěn)信號”。
?
????比較起來，現(xiàn)在使用解析信號法(HHT)估計瞬時頻率的人是最多的，大有“眾望所歸”的趨勢。此方法就是本人在前幾篇博文中使用的方法。可是，它存在的問題，我在博文《18、關于IMF分量瞬時頻率跳變問題的研究》也說的很詳細了。?既然現(xiàn)在最流行、影響最大的瞬時頻率估計法都是這樣，因此我有理由對其它所有瞬時頻率估計方法，都感到?jīng)]有信心。
?
????解析信號法(HHT)估計瞬時頻率用的是hhspectrum函數(shù)，hhspectrum函數(shù)是通過調(diào)用時頻工具箱中的instfreq函數(shù)來進行瞬時頻率估計的。時頻工具箱中還有一個使用AR(2)模型進行瞬時頻率估計的函數(shù)ifestar2（因為除卻注釋部分，程序正文很短，所以我把它也附在后面，方便閱讀。見“附二”）。
?
????現(xiàn)在來看看用ifestar2函數(shù)估計脈搏序列各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”附近，如果都是這樣還比較好區(qū)分，但那些比“0.25Hz”稍高的頻率，你怎么知道它是否經(jīng)過跳變呢？如果不是，那多大的頻率是由負頻跳變而來？界限在哪里？從物理意義來講，正頻率與負頻率是截然不同的頻率，一個是正向旋轉(zhuǎn)所形成，一個是負向旋轉(zhuǎn)所形成。頻率表示方法，這是人類創(chuàng)造的理論，不管現(xiàn)實如何，用“0.49Hz”來表示“-0.01Hz”，這就叫“削足適履”。
?
????以前多次講過，現(xiàn)實信號的頻率，一般情況下應該是連續(xù)、光滑的，所以我們的頻率表示方法好不好，也應該以此為最高準則，是正頻率就表示為正頻率，是負頻率就表示為負頻率。除非能證明現(xiàn)實信號頻率的確發(fā)生了跳變，否則頻率表示當中出現(xiàn)了“跳變”，都是不對的。
?
????下面將前面兩張圖的時間軸放大看看：
?
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分量頻率，都有采樣點上沒有頻率值的情況發(fā)生(不是“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函數(shù)中的indices = find(den>0)造成的。與“反余弦法瞬時頻率估計”一樣，它能作頻率估算的地方，它就作了估算；它無法作估算的地方，它就擱那不管了。這是什么態(tài)度嘛！
?
?
?
?
?
?
?
附一：網(wǎng)上找到的“反余弦法瞬時頻率估計”程序代碼
?
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
?
?
附二 時頻工具箱中瞬時頻率估計函數(shù)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
總結
- 上一篇： linux里的run-level,lin
- 下一篇： python中单行注释采用的符号是什么_