EMD分解后如何求周期

hmjerry · 发表于 2011-3-16 15:03

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

本帖最后由 hmjerry 于 2011-3-16 15:04 编辑

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

请大家帮帮忙！
谢谢各位！
数据、代码和结果如下

11.txt (228 Bytes, 下载次数: 19)

clear
z=[textread('D:\11.txt')]
t=1982:1:2009;
x=[1982:2:2010]
c=emd(z); % EMD 分解
[m,n]=size(c);
for i=1:m-1
%--------------------------------------------------------------------
%计算各IMF的方差贡献率
%定义：方差为平方的均值减去均值的平方
%均值的平方
%imfp2=mean(c(i,:),2).^2
%平方的均值
%imf2p=mean(c(i,:).^2,2)
%各个IMF的方差
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
end;
for i=1:m-1
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
%方差百分比，也就是方差贡献率
mseb(i)=mse(i)/sum(mse)*100;
%显示各个IMF的方差和贡献率
end;
figure(1)
for i=1:m-1
disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
end;
subplot(m+1,1,1)
plot(t,z)
set(gca,'fontname','times New Roman')
set(gca,'fontsize',12.0)
set(gca,'xtick',x);
ylabel(['signal'])
for i=1:m-1
subplot(m+1,1,i+1);
set(gcf,'color','w')
plot(t,c(i,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',12.0)
set(gca,'xtick',x);
ylabel(['imf',int2str(i)])
end
subplot(m+1,1,m+1);
set(gcf,'color','w')
plot(t,c(m,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',12.0)
set(gca,'xtick',x);
ylabel(['r',int2str(m-1)])
xlabel('Time')

复制代码

hmjerry · 发表于 2011-3-17 14:42

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

上山人 · 发表于 2011-3-19 02:26

学习中，帮你顶。

ChaChing · 发表于 2011-4-30 21:31

回复 1 # hmjerry 的帖子

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

hmjerry · 发表于 2011-5-3 09:15

回复 4 # ChaChing 的帖子

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

xjd · 发表于 2011-5-3 10:35

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

唐三否 · 发表于 2013-5-17 10:41

困扰啊？如何解决？

账号		自动登录	找回密码
密码			我要加入

[HHT] EMD分解后如何求周期

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐: