4、利用instfreq函数求取瞬时频率时出现的问题
最近利用HHT做突变检测，遇到了如下问题：
首先根据定义求取的瞬时频率：
s1=hilbert(imf(1,:));
s2(:,1)=s1(1,:);
instphase=angle(s2);
unwrapinstphase=unwrap(instphase);
instanglefrequency=diff(unwrapinstphase);
realistfre=instanglefrequency/(2*pi);
利用函数直接求取的瞬时频率
int=instfreq(s2);

figure(1);
subplot(2,1,1);
plot(realistfre);
title('自己定义的瞬时频率');
subplot(2,1,2);
plot(inp);
title('老外定义的瞬时频率');


结果如下：
2010-3-27 13:30 上传
下载附件 (78.24 KB)
瞬时频率


为什么会出现负频率现象，根据HHT理论，IMF满足窄带信号且是单分量信号，利用解析信号法求取瞬时频率应该不会有负频率，我查看了那个instfreq函数，是将相位的差分取了绝对值，当然不会出现负值，可是我不知道这是什么含义，而且它难道不会影响原信号的频谱特征吗？希望大家能够帮忙！
答：看了几天的文献，终于搞明白了情况，下面和大家说一下，基本上做完EMD分解后，剩下的就是如何得到瞬时频率了，HUANG最一开始是由解析信号法即希尔伯特变换来求取瞬时频率的，由于频率的定义是由相位的导数来定义的，那必然是非常精确地，也就产生了一个问题，它对噪音是十分敏感的，我们通常利用差商来代替偏导数必然会造成精度的缺失，我一开始利用一阶向后差分来代替偏导数，这样出现的问题就是必须要求相位是严格递增的（解析信号的相位曲线是严格递增的，大家可以从相位曲线上看出来），注意是严格递增，也就是说，中间有噪音影响的情况下，那么也会出现后一个值大于前一个的情况，这种情况下就会出现负频率，所以解决的办法就是提高差分的精度，利用高阶差分，这样会利用到周
伟多数的点，也就避免了这种情况。当然瞬时频率的求法有很多种，解析相位法只是其中一种，HUANG已经发明了一种不利用希尔伯特变换来求取瞬时频率的方法，我就不再这复述了，希望大家利用某种数学手段分析问题的时候能够多考虑原因，
樓主所說的Huang得到瞬時頻率的方法
在Huang, N. E., Z. Wu, S. R. Long, K. C. Arnold,  X. Chen and K. Blank (2009), On instantaneous frequency, Advance in Adaptive Data Analysis . Vol.1, No.2. 177-229.
是使用正規化的的方式處理IMF,不用做HT
5、完整的EMD分解全过程，有Hilbert谱和边际谱
%示例程序
N=1000;
n=1:N;
fs=1000;
t=n/fs;
fx=10;
fy=50;
x=cos(2*pi*fx*t);
y=10*cos(2*pi*fy*t);
z=x+y;
data=z;
imf=emd(data);                        %frequency函数计算频数对输入信号进行EMD分解   
[A,f,t]=hhspectrum(imf);            %对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号
[E,t,Cenf]=toimage(A,f);            %将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率  这里横轴为时间，纵轴为频率       
                                                   %即时频图(用颜表示第三维值的大小)和三维图（三维坐标系：时间，中心频率，振幅）         
cemd_visu(data,1:length(data),imf);   %显示每个IMF分量及残余信号--------------------------------------------
disp_hhs(E);                          %希尔伯特谱----------------------------------------------------------
%画出边际谱
%N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure(3);
plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率
xlabel('频率 / Hz');
ylabel('幅值');

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。