4、利用instfreq函数求取瞬时频率时出现的问题.docx
《4、利用instfreq函数求取瞬时频率时出现的问题.docx》由会员分享，可在线阅读，更多相关《4、利用instfreq函数求取瞬时频率时出现的问题.docx（9页珍藏版）》请在冰点文库上搜索。
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（'老外定义的瞬时频率'）;
结果如下：
为什么会出现负频率现象，根据HHT理论，IMF满足窄带信号且是单分量信号，利用解析信号法求取瞬时频率应该不会有负频率，我查看了那个instfreq函数，是将相位的差分取
了绝对值，当然不会出现负值，可是我不知道这是什么含义，而且它难道不会影响原信号的频谱特征吗？
希望大家能够帮忙！
答：
看了几天的文献，终于搞明白了情况，下面和大家说一下，基本上做完EMD分解后，剩下的就是如何得到瞬时频率了，HUANG最一开始是由解析信号法即希尔伯特变换来求取瞬时频率的，由于频率的定义是由相位的导数来定义的，那必然是非常精确地，也就产生了一个问题，它对噪音是十分敏感的，我们通常利用差商来代替偏导数必然会造成精度的缺失，我一开始利用一阶向后差分来代替偏导数，这样出现的问题就是必须要求相位是严格递增的（解析信号的相位曲线是严格递增的，大家可以从相位曲线上看出来），注意是严格递增，也就是说，中间有噪音影响的情况下，那么也会出现后一个值大于前一个的情况，这种情况下就会出现负频率，所以解决的办法就是提高差分的精度，利用高阶差分，这样会利用到周伟多数的点，也就避免了这种情况。
当然瞬时频率的求法有很多种，解析相位法只是其中一种，HUANG已经发明了一种不利用希尔伯特变换来求取瞬时频率的方法，我就不再这复述了，希望大家利用某种数学手段分析问题的时候能够多考虑原因，
樓主所說的Huang得到瞬時頻率的方法
在Huang,N.E.,Z.Wu,S.R.Long,K.C.Arnold,X.ChenandK.Blank（2009）,Oninstantaneousfrequency,AdvanceinAdaptiveDataAnalysis.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）; %对输入信号进行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）;%设置频率点数%完全从理论公式出发。
网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解
fork=1:
size（E,1）bjp（k）=sum（E（k,:
））*1/fs;
endfigure（3）;
plot（Cenf（1,:
）*fs,bjp）;%作边际谱图进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率
xlabel（'频率/Hz'）;
ylabel（'幅值'）;
分解后的IMF
希尔伯特谱
边际谱
Hilbert谱主要是观察频率随时间的变化规律，以及能量的相对变化，而且这个HILBERT谱的幅值是对颜色进行归一化后的值，颜色的相对变化来看出它的相对高低变化的，所以你说的值并不是它的原始幅值，
6、调频信号，HHT和fft哪个正确？
对于信号：
clearfs=1;
N=1000;
t=0:
fs:
（N-1）*fs;x1=sin（3*t）;x2=cos（0.5*t+cos（0.05*t））;x3=sin（0.06*t）;x=x1+x2+x3;
HHT边际谱如下：
fft如下：
HHT和fft哪个正确？
为什么HHT边际谱调频段（即x2）的最大幅值分布在两端？
你的想法是从fft出发，这样fft肯定是正确的，但是我觉得真正的频率应该是相位对时间求导即d（0.5t+cos（0.05t））/dt=0.5+0.05sin（0.05t）的连续频带
答：
我觉得你的边际谱好像不对，调频信号它的频率是随着时间变化，它的能量应该集中在载波频率的周围。
而不应该是分布在二边。
你的求边际谱的程序可能有问题。
答：
HHT边际谱中的频率与傅立叶分析中的频率意义完全不同。
在傅立叶表达中，在某一频率W处能量的存在，代表一个正弦或余弦波在整个时间长度上都存在。
HHT中，在某一频率W处能量的存在，仅代表在数据的整个时间长度上，很可能有这样一个频率的振动波在局部出现过。
因为HHT和FFT中频率的定义不同，出现上述差异是正常的。
HHT边际谱是HHT时频谱对时间的积分，出现“为什么HHT边际谱调频段（即x2）的最大幅值分布在两端？
”这一问题，可以借助正弦信号概率分布密度函数出现类似情况来理解。
如有甚么问题，可以发信给我：
xyg_1975@
我很少上网，一般mail都会回复的
答：
边际谱从统计意义上表示了整组数据每个频率点的积累幅值分布，而傅立叶谱的某点幅值表示在整个信号里有一个含有此频率的三角函数组分，而且幅值越大只是说明在整个数据段上，局部存在的可能性越大。
再看得到的图形，FFT表示的是整个数据中，能量在一个频率上分布的可能性地描述，而边际谱表示在在每一个频率上幅值的积累，如果想知道具体时间那么就看HHT谱，这个时间-幅值-频率的三维谱。
说到瞬时频率，傅立叶变换不强调局部性，而是强调全局性。
咱们的HHT才提出一个唯一的瞬时频率的定义。
因此拿瞬时频率来衡量傅立叶变换也是不公平的。
答：
t=[0:
1/10000:
1];s=cos（10*pi*t）+2*cos（40*pi*t）;imf分量图和时频图如下，时频图中的频率对应的还是对的，第一幅对应的是20Hz，第二幅是5Hz。
imf分量图时频图
l=length（s）;a=size（c）;b=a
（1）;
form=1:
b
H（m,:
）=hilbert（c（m,:
））;
x（m,:
）=real（H（m,:
））;
y（m,:
）=imag（H（m,:
））;fi（m,:
）=atan（y（m,:
）./x（m,:
））;
end
g（m,:
）=szwf（fi（m,:
）,t）;
f（m,:
）=g（m,:
）./（2*pi）;subplot（b,1,m）;plot（t,f（m,:
））;
xlabel（'t（s）'）;
ylabel（'freq（hz）'）
title（'时频图'）;
这个是做完emd分解后，根据IMF分量作出的，就是根据我看得论文中时频图的做法，求出瞬时频率做出的。
流程是先进行hilbert变换（直接调用的），再求其虚部和实部之比的反正切，再求数值微分（这个程序是自己编的，经验证还算可以），最后得出时频图。
求出的相位要解卷绕后，才能进行数值微分。
解卷绕是求相位时用angle再unwrap就可以了。
你作的那个图不能叫时频图，没有反映幅值信息，应该叫瞬时频率图。
zhangnan3509做的那个才是时频图。
升级会员 