Matlab求瞬时频率（基于EMD）

你的EMD是% G. Rilling,的那个吗？

len = size(imf,1)

for k = 1:len

len1 = length(imf(k,:))

b(k) = sum(imf(k,:).*imf(k,:))/len1% 时域均方值，能量

amp(k,:) = abs(imf(k,:))

b(k) = sqrt(b(k))

th = angle(hilbert(imf(k,:)))%Hilbert变换的相位

d(k,:) = diff(th)/Ts/(2*pi)%求导，得到频率:f = (1/2*pi)*d(th)/dt

end

你的频率公式用得有点搜橘问题，求出来不闷伍应是归蚂漏或一化频率

我最近也在做这个，我用的是1kHz单频正弦波

unwrap的作用是解卷绕

没用unwrap函陆液数时，瞬时相位是这样的

加上unwrap函数之后是这样的

fs=80是你信号的采样频率

xhd1=fs*diff(xh1)/(2*pi)是要求信号的瞬弯早时频率，就是瞬时相位的导数

前面乘fs是为了跟实际频率对应上

不乘fs是这样的

纵轴频率对应的不是1kHz，但是乘了fs就变成下面这样，纵轴频率就是实际频率了

想换别的信号，就把xn换成你想改的信号，再根据你的采样频率改一下埋悉雀fs就行了

你这个程序贴错了吧，

Ts=0.001还差不多

如果是这样的话，就是fs是1000Hz。

atan在matlab里面取值范行做源围就是[-pi/2,pi/2]，因此肯定存在从pi/2到-pi/2的突变，如果采样频率是1000hz的话，1000对应2pi，那么突变-pi对应就是-500hz,，所以才会出现从50突变到-450的元素。

求瞬时频率为了避免这种突变，都是取一个中间胡御变量，如

c(i)=signal(i)./signal(i-1)或者c(i)=signal(i).*signal(i-1)' (signal(i-1)的共轭)

然后档态直接freq = atan(c)/(2*pi*Ts)

for i=1:length(siganal)

theta(i)=atan(imag(signal(i))/real(signal(i)))

end

直接写作theta=atan(imag(signal)/real(signal))

避免循环

---