关于EMD程序使用的问题

zhuxiaoxun

为什么我用的这个程序不能处理其它数据，只能处理它给定的数据，程序里带的另外三组调频或调幅的数据也处理不了，甚至连这个给定数据增加了相位后也无法处理？亟待高手解决！
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([0 400 -4 4]);
   xlabel('signal time/ms')
   ylabel('amplitude')
   subplot(6,1,2);
   plot(t,c1,'k')
   axis([0 400 -2 2]);
   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([0 400 -2 2]);
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([0 400 -1.5 1.5]);
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([0 400 -2 2]);

r4=s-c1-c2-c3-c4;
subplot(6,1,6);
plot(t,r4,'k')
xlabel('r time/ms')
   ylabel('amplitude')
     axis([0 400 -2 2]);

s3=c1+c2+c3+c4;
s4=s-s3;
figure
subplot(2,1,1)
plot(t,s3,'k')
axis([0 400 -3.5 3.5]);
xlabel('IMF分量之和 时间/ms')
ylabel('幅度')
subplot(2,1,2)
plot(t,s4,'k')
axis([0 400 -0.4 0.4]);
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

自己改一改里面的东西即可

zhuxiaoxun

程序方面的东西不是很在行。不知道是哪里出了问题。尤其是不能处理实验得到的数据。还请指教。谢谢！:@)

anycall

这个程序根本运行不了，很多地方得改

qqvirile

我运行了楼主的程序了 有点不明白见附近
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和函数的差 这个图的意义是残余量么?

吃书的老虎

程序没有通用性，在写程序的时候没有固定的长度、极大值、极小值。。。。。。。

zhuxiaoxun

问题1：这好像是源程序，不是我的改动，不过我好像也不大明白；
问题2：求差的目的就是考察分解的结果是否具备完备性，即看看分解的信号与原信号差多少，从而衡量分解信号能否重构原信号。
这是我的理解。

gzslbs2003

对我很有帮助，谢谢！

zhuxiaoxun

不客气，大家一起探讨。
可是为什么没有哪位大侠解决一下我的问题？:@(

zhuxiaoxun

小弟又遇见了一些新问题。可能问题很弱，但确实给我带来了很大麻烦。希望哪位大侠能出手相助。
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带坐标分解、边际谱和时频图。真的不胜感激。
我知道这很耽误大侠的宝贵时间，但看在我元旦还为这事头疼的份上就帮一下吧。:@)

zhuxiaoxun

是不是问题真的很烂？没有人回答啊。自己顶一下。:lol

shanzhongyao

问下楼主 这个程序的问题你解决了吗 我现在要作个EMD的毕设 需要他

shanzhongyao

问题解决没？