关于瞬时频率估计,前面虽说暂时放一下,但心中始终还是念念不忘,因为这是一道绕不过去的坎。在网上多次搜索、了解其现状。感觉是关注这件事的人很多,方法很多,但问题也很多。在网上能找到的方法,简单归结如下:
相位差分法
相位建模法
Teager能量算子法
跨零点法
求根估计法
反余弦法
时频分布法(谱峰检测法?)
Shekel方法
Teager-Kaiser方法
解析信号法(HHT法)
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 =
lt =
rj =
附一:网上找到的“反余弦法瞬时频率估计”程序代码
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-
%
%
% Output-
%
%
%
% Used by-
%
% Kenneth Arnold (NASA GSFC)
%----- 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.
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
if (nargin == 0),
end;
[xrow,xcol] = size(x);
if (xcol~=1),
end
x = real(x);
if (nargin == 1),
end;
[trow,tcol] = size(t);
if (trow~=1),
elseif min(t)<4,
end;
Kappa = x(t-1) .* x(t-2) - x(t
psi1
psi2
den
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