怀疑instfreq.m程序的正确性

zhangnan3509

eight又发怒了！:@o  大家好好想想，自己先琢磨一下。还是靠自己，并在大家的帮助下，才能有进步。:handshake

eight

原帖由 zhangnan3509 于 2007-3-21 10:38 发表
eight又发怒了！:@o  大家好好想想，自己先琢磨一下。还是靠自己，并在大家的帮助下，才能有进步。:handshake

也没有发怒，只是觉得大家太依赖别人了。说实在的，很多问题我也想请教别人，但是最后也还是靠自己摸索解决的，而且我最怕别人催我（一个早上不到的时间几个帖子就可以看作“催”了）


粗略看了一下，EMD的最新版本主要有两个改进：
1. 提供了快速算法 —— 把c程序编译成 dll ，然后在matlab中调用。我以前做过类似的，速度快80倍左右
2. 提出了Bivariate Empirical Mode Decomposition，还没有细看，不过估计没有论文的话很难理解，有耐心和精力的可以不妨看看，跟以往的很不一样

[ 本帖最后由 eight 于 2007-3-21 11:13 编辑 ]

zhangnan3509

心情我理解，我也要注意类似的问题。:loveliness:  也觉得别人帮忙就省事很多，但是毕竟是别人嚼过的馒头不甜呀。

eight

原帖由 zhangnan3509 于 2007-3-21 11:11 发表
心情我理解，我也要注意类似的问题。:loveliness:  也觉得别人帮忙就省事很多，但是毕竟是别人嚼过的馒头不甜呀。

建议版主回帖时候点击“引用”或者下图第一行的“回复”，不然有时候不知道你在回答谁或者哪层楼，当然了，没有歧义的话也可以点击第二行的“回复”:loveliness:

zhangnan3509

收到！:@L

eight

前两天看了一下，问题的确存在，这是由于 instfreq 函数对瞬时频率的求解方式引起的。

和师弟研究了一下，发现使用如下代码即可实现瞬时频率的近似：

1. 
   
2. Himf = hilbert(imf(1,:));
   
3. f1 = diff(phase(Himf))*fs/(2*pi); % fs 是采样频率，2*pi 可以去掉，这样得到的就是原信号中 cos 或 sin 后面的那个频率值
   
4. 
   
5. Himf = hilbert(imf(2,:));
   
6. f2 = diff(phase(Himf))*fs/(2*pi);
   
7. 
   
8. t1 = t(1:end-1);
   
9. t2 = t(1:end-1);
   
10. 
    
11. figure;
    
12. plot(t1,f1,'r',t2,f2,'b');

复制代码

对 x = cos(t)+0.5*cos(t/2) 和 x = sin(t) + sin(3*t) 这两个信号测试成功。

不过仍然存在以下问题：
1. “近似”表现在 diff 函数的不稳定上，即求得的瞬时频率值可能很大
2. 没有对信号两端特殊处理（故此时得到的在两端的瞬时频率值没有意义），且得到的 f1 或者 f2 比原数据的长度少了1
3. diff 是最简单、最粗糙的一种近似，可以用二阶差分来改善，另外，无论是一阶还是二阶差分，本质上仍然是微分，因此并非数值稳定。要追求更稳定的结果请用数值积分取代，这可以通过使用卷积核来实现

以下是 x = sin(t) + sin(3*t) 的测试结果：

[ 本帖最后由 eight 于 2007-5-2 11:53 编辑 ]

huangyong87

与这里的问题一样
http://forum.vibunion.com/thread-41028-1-1.html请参见6、8、9、12楼的意见。
楼主若将[fnor,t]=instfreq(x');变为[fnor,t]=instfreq(x.');就能得到正确答案。因为x是复数，x'为复共轭转置，x.'为非共轭转置。具体修改可参见下文：
Himf1=hilbert(imf(1,:));
[fnormhat1,t1]=instfreq(Himf1.');         
Himf2=hilbert(imf(2,:));
[fnormhat2,t2]=instfreq(Himf2.');
plot(t1,fnormhat1*fs,'r',t2,fnormhat2*fs,'b');

eight

原帖由 huangyong87 于 2007-5-2 15:51 发表
与这里的问题一样http://forum.vibunion.com/thread-41028-1-1.html请参见6、8、9、12楼的意见。
楼主若将[fnor,t]=instfreq(x');变为[fnor,t]=instfreq(x.');就能得到正确答案。因为x是复数，x'为复共 ...

这样可以得到正确瞬时频率吗？如果可以的话，那就没有问题了

huangyong87

这样可以得到正确瞬时频率吗？如果可以的话，那就没有问题了

x = cos(t)+0.5*cos(t/2)的运算结果。

lsclsc27

哈，instfreq函数是归一化了的，取消其归一化就可以了，我用的是一阶差分:victory:

liujing123

问题是怎样取消归一化?instfreq.m我怎么也研究不透彻,想改进也没办法/就取消归一化大家能不能给指点一下,我也是想用这个解决一下频率随时间变化的单分量信号，也就是每个时间对应一个频率。

green

我按照5楼的程序运行了一下，得到的图并不是全部红线都在蓝线下面，而是前半段红线在下，但后半段是红线在上

fancy78066994

看了两天的关于HHT的帖子，想必你应该是最了解HHT的，呵呵！我有这么几个关于HHt的问题：
1对信号进行EMD分解成IMF后，每一个IMF应该是窄带的单分量信号，那么利用希尔伯特变换求出的瞬时频率就不应该有负的频率存在，况且解析信号的频谱的负频率应该都是0，而我在分析一段实际数据的时候，利用相位的导数定义的瞬时频率求取时，会出现负频率，我看了instfreq程序，相位的差分它取了绝对值，那么频率是一定不会出现负的了，但是我不知道这到底是什么物理含义？不知你在做实际数据分析的时候，有没有碰到过负频率？
2、EMD分解使其成为IMF的条件2，（均值为0）是利用HUANG等人的SD在0.2与0.3之间定义的，我不知道这个问题是否出在这？
3、不知如何方便和你联系！呵呵