EMD分解后如何求周期

hmjerry

大家好，看论坛里都是做信号处理的，我是跨专业，没接触过信号处理，现在想用EMD方法做周期分析，因为用的是气象资料，无法给出采样频率，现在已经得到了IMF图，想求出各个IMF中的周期，请问该如何编程？

请大家帮帮忙！
谢谢各位！
数据、代码和结果如下
11.txt (228 Bytes, 下载次数: 19)

2011-3-16 15:04 上传

点击文件名下载附件

1. clear
   
2. z=[textread('D:\11.txt')]
   
3. t=1982:1:2009;
   
4. x=[1982:2:2010]  
   
5. c=emd(z);  %  EMD 分解
   
6. [m,n]=size(c);  
   
7. for i=1:m-1  
   
8. %--------------------------------------------------------------------
   
9. %计算各IMF的方差贡献率
   
10. %定义：方差为平方的均值减去均值的平方
    
11. %均值的平方
    
12. %imfp2=mean(c(i,:),2).^2
    
13. %平方的均值
    
14. %imf2p=mean(c(i,:).^2,2)
    
15. %各个IMF的方差
    
16. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
    
17. end;
    
18. for i=1:m-1
    
19. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
    
20. %方差百分比，也就是方差贡献率
    
21. mseb(i)=mse(i)/sum(mse)*100;
    
22. %显示各个IMF的方差和贡献率
    
23. end;
    
24. 
    
25. figure(1)
    
26. for i=1:m-1
    
27. disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
    
28. end;
    
29. subplot(m+1,1,1)
    
30. plot(t,z)
    
31. set(gca,'fontname','times New Roman')
    
32. set(gca,'fontsize',12.0)
    
33. set(gca,'xtick',x);  
    
34. ylabel(['signal'])   
    
35. 
    
36. for i=1:m-1
    
37. subplot(m+1,1,i+1);
    
38. set(gcf,'color','w')
    
39. plot(t,c(i,:),'k')
    
40. set(gca,'fontname','times New Roman')
    
41. set(gca,'fontsize',12.0)
    
42. set(gca,'xtick',x);  
    
43. ylabel(['imf',int2str(i)])
    
44. end
    
45. subplot(m+1,1,m+1);
    
46. set(gcf,'color','w')
    
47. plot(t,c(m,:),'k')
    
48. set(gca,'fontname','times New Roman')
    
49. set(gca,'fontsize',12.0)
    
50. set(gca,'xtick',x);  
    
51. ylabel(['r',int2str(m-1)])
    
52. xlabel('Time')
    
53. 
    

复制代码

hmjerry

这么久了  怎么都没人回复呢...大家帮帮忙啊！

上山人

学习中，帮你顶。

ChaChing

回复 1 # hmjerry 的帖子

水平有限但好奇, 若如LZ所说无法给出采样频率, 可以求出各个对应的周期吗?
敢问LZ的数据是否为一年一笔数据, 若是, 那一年不就是采样週期(采样频率的倒数)

hmjerry

回复 4 # ChaChing 的帖子

谢谢你的问题   
我用的是气温数据，原数据为一天一个值，处理成年平均值，现在所用为一年一个值。
想要计算的周期是很多年的变化周期，若给出一天或一年为采样频率，则周期为一天或一年，但此周期并不能代表长年的变化周期。
有人求算出每个IMF所对应的中心频率，然后求中心频率的倒数为周期，但个人觉得还不是很准确，所为问题一直没有得到解决，若有兴趣，可以一起讨论。

xjd

我用的是轴承故障数据，滚动体有9个，出现9个冲击波形即为一个周期。

唐三否

困扰啊？如何解决？