关于EMD程序使用的问题
为什么我用的这个程序不能处理其它数据，只能处理它给定的数据，程序里带的另外三组调频或调幅的数据也处理不了，甚至连这个给定数据增加了相位后也无法处理？亟待高手解决！function s=emddec()
t=0:1:2000;
%s=sin(2*pi*0.012*t)+(1+0.2*sin(2*pi*0.075*t)).*cos(2*pi*0.15*t+0.5*sin2*(pi*0.015*t));
%s=(1+0.5*sin(2*pi*0.01*t)).*cos(2*pi*0.14*t+0.4*sin(2*pi*0.03*t))+sin(2*pi*0.02*t)+sin(2*pi*0.07*t);
s=sin(2*pi*0.08*t)+sin(2*pi*0.11*t)+sin(2*pi*0.15*t)+sin(2*pi*0.17*t);
%s1=sawtooth(0.005*t,0.5)+0.5*sin(2*pi*0.01*t);
s1=s;
sd=1;
sd1=1;
sd2=1;
sd3=1;
sd4=1;
for n=1:29
if(sd>0.3)
jd=find(diff(sign(diff(s1)))==-2)+1;%找h1极大值的位置
jx=find(diff(sign(diff(s1)))==2)+1; %找h1极小值的位置
ah=length(jd);
al=length(jx);
jdz(ah)=0;
for i=1:ah
bh=jd(i);
jdz(i)=s1(bh);
end %找出h1极大值对应函数值
jxz(al)=0;
for i=1:al
bl=jx(i);
jxz(i)=s1(bl);
end %找出h1极小值对应函数值
jsbl=interp1(jd,jdz,t,'cubic'); %极大值拟和的上包络
jxbl=interp1(jx,jxz,t,'cubic');%极小值拟和的下包络
m1=(jsbl+jxbl)/2;%上下包络均值
h1=s1-m1;
sd=sum(((s1-h1)./h1).^2);
s1=h1;
end
end
c1=h1;
figure;
subplot(6,1,1);
plot(t,s,'k')
axis();
xlabel('signal time/ms')
ylabel('amplitude')
subplot(6,1,2);
plot(t,c1,'k')
axis();
xlabel('IMF1 time/ms')
ylabel('amplitude')
r1=s-c1;
for n1=1:54
if (sd1>0.3)
jdh1=find(diff(sign(diff(r1)))==-2)+1;%找h1极大值的位置
jxh1=find(diff(sign(diff(r1)))==2)+1; %找h1极小值的位置
ahh1=length(jdh1);
alh1=length(jxh1);
jdzh(ahh1)=0;
for i=1:ahh1
bhh1=jdh1(i);
jdzh1(i)=r1(bhh1);
end %找出h1极大值对应函数值
jxzh1(alh1)=0;
for i=1:alh1
blh1=jxh1(i);
jxzh1(i)=r1(blh1);
end %找出h1极小值对应函数值
jsblh1=interp1(jdh1,jdzh1,t,'cubic'); %极大值拟和的上包络
jxblh1=interp1(jxh1,jxzh1,t,'cubic');%极小值拟和的下包络
m2=(jsblh1+jxblh1)/2;%上下包络均值
h2=r1-m2;
sd1=sum(((h2-h1)./h1).^2);
r1=h2;
end
end
c2=h2;
subplot(6,1,3);
plot(t,c2,'k')
axis();
xlabel('IMF2 time/ms')
ylabel('amplitude')
r2=s-c1-c2;
for n2=1:7
if (sd2>0.3)
jdh2=find(diff(sign(diff(r2)))==-2)+1;%找h1极大值的位置
jxh2=find(diff(sign(diff(r2)))==2)+1; %找h1极小值的位置
ahh2=length(jdh2);
alh2=length(jxh2);
jdzh2(ahh2)=0;
for i=1:ahh2
bhh2=jdh2(i);
jdzh2(i)=r2(bhh2);
end %找出h1极大值对应函数值
jxzh2(alh2)=0;
for i=1:alh2
blh2=jxh2(i);
jxzh2(i)=r2(blh2);
end %找出h1极小值对应函数值
jsblh2=interp1(jdh2,jdzh2,t,'cubic'); %极大值拟和的上包络
jxblh2=interp1(jxh2,jxzh2,t,'cubic');%极小值拟和的下包络
m3=(jsblh2+jxblh2)/2;%上下包络均值
h3=r2-m3;
sd2=sum(((h3-h2)./h2).^2);
r2=h3;
end
end
c3=h3;
subplot(6,1,4);
plot(t,c3,'k')
axis();
xlabel('IMF3 time/ms')
ylabel('amplitude')
r3=s-c1-c2-c3;
for n3=1:100
if (sd3>0.3)
jdh3=find(diff(sign(diff(r3)))==-2)+1;%找h1极大值的位置
jxh3=find(diff(sign(diff(r3)))==2)+1; %找h1极小值的位置
ahh3=length(jdh3);
alh3=length(jxh3);
jdzh3(ahh3)=0;
for i=1:ahh3
bhh3=jdh3(i);
jdzh3(i)=r3(bhh3);
end %找出h1极大值对应函数值
jxzh3(alh3)=0;
for i=1:alh3
blh3=jxh3(i);
jxzh3(i)=r3(blh3);
end %找出h1极小值对应函数值
jsblh3=interp1(jdh3,jdzh3,t,'cubic'); %极大值拟和的上包络
jxblh3=interp1(jxh3,jxzh3,t,'cubic');%极小值拟和的下包络
m4=(jsblh3+jxblh3)/2;%上下包络均值
h4=r3-m4;
sd3=sum(((h4-h3)./h3).^2);
r3=h4;
end
end
c4=h4;
subplot(6,1,5);
plot(t,c4,'k')
xlabel('IMF4 time/ms')
ylabel('amplitude')
axis();
r4=s-c1-c2-c3-c4;
subplot(6,1,6);
plot(t,r4,'k')
xlabel('r time/ms')
ylabel('amplitude')
axis();
s3=c1+c2+c3+c4;
s4=s-s3;
figure
subplot(2,1,1)
plot(t,s3,'k')
axis();
xlabel('IMF分量之和 时间/ms')
ylabel('幅度')
subplot(2,1,2)
plot(t,s4,'k')
axis();
xlabel('IMF和函数与原函数之差 时间/ms')
ylabel('幅度')
for n=1:1500
c11(n)=c1(n+300);
c21(n)=c2(n+300);
c31(n)=c3(n+300);
c41(n)=c4(n+300);
end
h=0.001;
t=0:h:0.6;
X=fft(c11);
f=t/h;
figure;
subplot(2,4,1);
plot(f(1:floor(length(f)/2)),abs(X(1:floor(length(f)/2))),'k');
ylabel('幅度')
xlabel('频率 /Hz')
X1=fft(c21);
f=t/h;
subplot(2,4,2);
plot(f(1:floor(length(f)/2)),abs(X1(1:floor(length(f)/2))),'k');
ylabel('幅度')
xlabel('频率 /Hz')
X2=fft(c31);
f=t/h;
subplot(2,4,3);
plot(f(1:floor(length(f)/2)),abs(X2(1:floor(length(f)/2))),'k');
ylabel('幅度')
xlabel('频率 /Hz')
X3=fft(c41);
f=t/h;
subplot(2,4,4);
plot(f(1:floor(length(f)/2)),abs(X3(1:floor(length(f)/2))),'k');
ylabel('幅度')
xlabel('频率 /Hz')
m=0:1000;
y=2*sin(2*pi*0.24*m)+0.2*sin(2*pi*0.17*m);
y1=2*sin(2*pi*0.17*m)+0.2*sin(2*pi*0.24*m)+0.1*sin(2*pi*0.115*m)+0.1*sin(2*pi*0.05*m);
y2=0.2*sin(2*pi*0.17*m)+0.1*sin(2*pi*0.24*m)+2*sin(2*pi*0.115*m)+0.3*sin(2*pi*0.05*m);
y3=2*sin(2*pi*0.05*m)+0.2*sin(2*pi*0.115*m)+0.1*sin(2*pi*0.17*m);
h=0.001;
t=0:h:0.6;
X=fft(y);
f=t/h;
subplot(2,4,5);
plot(f(1:floor(length(f)/2)),abs(X(1:floor(length(f)/2))),'k');
ylabel('IMF1 幅度')
xlabel('频率 /Hz')
X1=fft(y1);
f=t/h;
subplot(2,4,6);
plot(f(1:floor(length(f)/2)),abs(X1(1:floor(length(f)/2))),'k');
ylabel('IMF2 幅度')
xlabel('频率 /Hz')
X2=fft(y2);
f=t/h;
subplot(2,4,7);
plot(f(1:floor(length(f)/2)),abs(X2(1:floor(length(f)/2))),'k');
ylabel('IMF3 幅度')
xlabel('频率 /Hz')
X3=fft(y3);
f=t/h;
subplot(2,4,8);
plot(f(1:floor(length(f)/2)),abs(X3(1:floor(length(f)/2))),'k');
ylabel('IMF4 幅度')
xlabel('频率 /Hz')
:@( :@( :@( :@( :@( :@( :@( :@( :@( :@( 自己改一改里面的东西即可
回复 沙发 laoda 的帖子
程序方面的东西不是很在行。不知道是哪里出了问题。尤其是不能处理实验得到的数据。还请指教。谢谢！:@) 这个程序根本运行不了，很多地方得改我运行了楼主的程序了 有点不明白
我运行了楼主的程序了 有点不明白见附近1请问你这样做有什么意义么？也就是你程序的
m=0:1000;
y=2*sin(2*pi*0.24*m)+0.2*sin(2*pi*0.17*m);
y1=2*sin(2*pi*0.17*m)+0.2*sin(2*pi*0.24*m)+0.1*sin(2*pi*0.115*m)+0.1*sin(2*pi*0.05*m);
y2=0.2*sin(2*pi*0.17*m)+0.1*sin(2*pi*0.24*m)+2*sin(2*pi*0.115*m)+0.3*sin(2*pi*0.05*m);
y3=2*sin(2*pi*0.05*m)+0.2*sin(2*pi*0.115*m)+0.1*sin(2*pi*0.17*m);
h=0.001;
t=0:h:0.6;
X=fft(y);
f=t/h;
subplot(2,4,5);
plot(f(1:floor(length(f)/2)),abs(X(1:floor(length(f)/2))),'k');
ylabel('IMF1 幅度')
xlabel('频率 /Hz')
2原函数和imf和函数的差 这个图的意义是残余量么? 程序没有通用性，在写程序的时候没有固定的长度、极大值、极小值。。。。。。。
回复 5楼 qqvirile 的帖子
问题1：这好像是源程序，不是我的改动，不过我好像也不大明白；问题2：求差的目的就是考察分解的结果是否具备完备性，即看看分解的信号与原信号差多少，从而衡量分解信号能否重构原信号。
这是我的理解。 对我很有帮助，谢谢！ 不客气，大家一起探讨。
可是为什么没有哪位大侠解决一下我的问题？:@( 小弟又遇见了一些新问题。可能问题很弱，但确实给我带来了很大麻烦。希望哪位大侠能出手相助。
1.就是关于坐标的问题。为什么我分解出来的信号没有坐标？
2.做Hilbert边际谱总得不到正确的结果。每次都能出结果，但是好像对于任何信号得到的结果都一样，频率都是出现在450HZ-500HZ左右，是不是也是坐标出了问题？可是设置了采样频率了还是不对，如果采样设置为1000，那频率就出现在450HZ-500HZ，如果采样频率设置为10000，那频率就出现在4500HZ-5000HZ。让我感到很滑稽。
3.哪位大侠能帮我做一下HHT时频图？
如果真有好心人，就帮小弟以x=cos(2*pi*30*t+0.5*sin(2*pi*15*t))+sin(2*pi*120*t)信号为例，帮小弟做一份完整的EMD带坐标分解、边际谱和时频图。真的不胜感激。
我知道这很耽误大侠的宝贵时间，但看在我元旦还为这事头疼的份上就帮一下吧。:@) 是不是问题真的很烂？没有人回答啊。自己顶一下。:lol 问下楼主 这个程序的问题你解决了吗 我现在要作个EMD的毕设 需要他 问题解决没？
页:
[1]