..:: 科学软件专区-emd,请教-请教：关于emd.m的问题，{Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛-(手机验证注册)}

babybear713 发表于 2009-9-26 09:46:22

请教：关于emd.m的问题

本帖最后由 babybear713 于 2009-9-28 17:00 编辑

想利用emd对一组时程数据进行处理，需要信号的瞬时频率图和相位图，但是程序进行到红色字体之前就无法继续下去了，就是出不了图。好急啊，有哪位高手能帮我看下么

% 本程序主要通过EMD和hilbert求IMF，作出HHT归一化能量谱图（三维图），边际谱图和瞬时能量图，并做完备性验证
clear all
clc
figure(10)
load('dn1.txt');
x1=dn1(:,2); t=1:1000;
save x1; plot(x1);
x=x1;
fs=50;
figure(20)
subplot(211);plot(0.02*t,x);
xlabel('t/s');
ylabel('m/s^2');
% 作出信号的傅里叶频谱图和功率谱图fft
figure(30)
subplot(211);
Y=fft(x,1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));
xlabel('frequency');
Pyy=Y.*conj(Y)/1024;figure(35);
plot(ff,Pyy(1:512));
xlabel('frequency');
%对信号作EMD分解，求出各IMF分量图和重构图及各IMF的时频图
=emd2(x,t,);
L=size(imf,1);
X=0;
for i=1:L;
X=imf(i,:)+X;
end
figure(40);subplot(211);
Y=x-X';
plot (0.02*t,Y);subplot(212);
plot(0.02*t,X);
% 画出IMF分量图
figure(50)
subplot(411);plot(0.02*t,imf(1,:));
ylabel('c1');
subplot(412);plot(0.02*t,imf(2,:));
ylabel('c2');
subplot(413);plot(0.02*t,imf(3,:));
ylabel('c3');
subplot(414);plot(0.02*t,imf(4,:));
ylabel('c4');
figure(60)
subplot(411);plot(0.02*t,imf(5,:));
ylabel('c5');
subplot(412);plot(0.02*t,imf(6,:));
ylabel('c6');
subplot(413);plot(0.02*t,imf(7,:));
ylabel('c7');
subplot(414);plot(0.02*t,imf(8,:));
ylabel('c8');

%完备性检验
figure(70);
subplot(211);
plot(0.02*t,imf(8,:));
subplot(212);
plot(0.02*t,imf(8,:)+imf(7,:));
figure(72);
subplot(211);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:));
subplot(212);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:)+imf(5,:));
figure(75);
subplot(211);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:)+imf(5,:)+imf(4,:));
subplot(212);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:)+imf(5,:)+imf(4,:)+imf(3,:));
figure(78);
subplot(211);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:)+imf(5,:)+imf(4,:)+imf(3,:)+imf(2,:));
subplot(212);
plot(0.02*t,imf(8,:)+imf(7,:)+imf(6,:)+imf(5,:)+imf(4,:)+imf(3,:)+imf(2,:)+imf(1,:));
%画出IMF分量图的频谱
figure(80);subplot(211)
Y=fft(imf(1,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(2,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(3,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(4,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));
hold off
figure(85);
subplot(212)
Y=fft(imf(5,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(6,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(7,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold on
Y=fft(imf(8,:),1024);
ff=50*(0:511)/1024;
plot(ff,abs((Y(1:512))));hold off

%作出HHT归一化能量谱图（三维图），边际谱图和瞬时能量图
N=1000;
L=size(imf,1);% Number of components in the decomposition
% loop for on each component
S=[];
clear x z m p freq
x=imf'; % now each column is a compenent
z=hilbert(x); % analytic signal
m=abs(z); % module of z
p=angle(z); % phase of z
for i=1:L
freq(:,i)=instfreq(z(:,i)); % instantaneous frequency
ceilfreq(:,i)=ceil(freq(:,i)*N);% to have integervalues-also do a smoothing given the number
   for j=1:length(x)-2
   S(ceilfreq(j,i),j+1)=m(j+1,i);
   end
end
eps=0.00001; % to avoid zero values before the log
E=20*log(S.^2+eps);% Hilbert energ spectrum
% plot S
figure(90);
t=1:length(x);% time sample
f=t/length(x)*25; % frequency %%%%%%%%%5
imagesc(t,f,E);
contour(t(1:3087),1:497,E);%画出等高线图
surf1(1:3087,1:500,S.^2); shading interp; colormap(gray);axis()
% mesh(1:3087,1:500,S.^2);
% mesh(1:2999,1:500,S); %画出三维图
set(gca,'YDir','normal');
xlabel('t/s');
ylabel('Frequency/Hz'); %%%%%%%%%%%%%%%%%%%%
% ylabel('Normalized Frequency');
% axis(); %%%%%%%%%%%%%%%%%%%%%%
S1=sum(S.^2,1);figure(100);plot(1:3087,S1); %瞬时能量
figure(110);
S2=sum(S.^2,2);plot(1:497,S2,'r') %边际能量谱
figure(120);
S2=sum(S,2);plot(1:497,S2,'g');%边际谱

dn1.txt已经加入附件

babybear713 发表于 2009-9-26 09:48:08

怎么出现这么多qq表情啊~~~:(
还好附件里面有
希望高手指教

taohe 发表于 2009-9-26 11:16:59

数据文件dn1.txt没有贴上来，应该是1000行x2列的文本文件吧。这个数据文件可能不是和该演示程序配套的，因为里面对某些矩阵的引用时所期望的矩阵维数和数据文件不符。你运行时遇到的问题也主要是因为矩阵维数不对。

问题主要出在contour和surf的调用上。下面是大概修改后的样子，把时间轴减短为一半，方便在我的电脑上测试。
...
halfTLen = length(t)/2;
numF = 497
contour(t(1: halfTLen),1: numF,E(1:numF, 1: halfTLen));%画出等高线图
numF = 500
surf(1: halfTLen,1:numF,S(1:numF, 1: halfTLen).^2); shading interp; colormap(gray);
...
另外，楼主能不能说说你准备在什么应用背景下使用EMD？

babybear713 发表于 2009-9-28 16:51:35

3# taohe
我想结合emd和hilbert变换对框架或者板的某些点的时程曲线进行分析，以此来识别损伤~~ 我贴出来的源程序是取自《爆破震动信号分析理论与技术》附录b的源程序，它的时间是t=1：3088，我当时就把这里改了一下，后面作图的就忘记改了~

babybear713 发表于 2009-9-28 16:55:28

3# taohe
在调用contour和surf和函数时我没注意，还是用的书上的~~由于matlab学得很浅，出了问题也不知道从什么地方下手检查~~真的很感谢你的帮助{:3_58:}

taohe 发表于 2009-9-28 19:01:16

原来这样，你附件中的做emd分解的程序似乎还是那个法国小组的源程序。在某些场合下，可能需要注意版权问题。如果对版权有疑问的话，可以写信问作者本人。源程序的注释中有联络信息。

我没有看过你说的那篇关于爆破的文章，不过我记得通过HHT或者emd获得的IMF并不具备正交性，这也正是这个方法美中不足的地方。但在你附件中的程序的注释中看到了有关“正交性”的信息，你知道是什么原因吗？

pasuka 发表于 2009-9-28 22:54:16

原来这样，你附件中的做emd分解的程序似乎还是那个法国小组的源程序。在某些场合下，可能需要注意版权问题。如果对版权有疑问的话，可以写信问作者本人。源程序的注释中有联络信息。

我没有看过你说的那篇关于爆破 ...
taohe 发表于 2009-9-28 19:01 http://forum.simwe.com/images/common/back.gif
法国人的程序用得人比较多，不过效果就是见仁见智了，现在还没有见过公开代码的EMD程序在处理端点上面可以达到邓拥军等在《高技术通讯》上那篇文章的效果
目前为止还没有数学上的严格证明imf的正交性，不过黄锷院士的文章中，提到了数值上验证过正交性，没啥大问题

taohe 发表于 2009-9-29 13:15:23

谢谢楼上的解释。是不是说在楼主给出的代码中关于“正交性”的部分其实就是一种关于数值验证的方式？

能不能给出关于你提到的邓拥军等在《高技术通讯》上那篇文章的详细信息，我想有机会的话找来看看。

另外，能否在推荐一些最近的比较值得看的关于HHT或者EMD的讨论的文献？

谢谢！

babybear713 发表于 2009-10-29 16:43:28

%作出HHT归一化能量谱图（三维图），边际谱图和瞬时能量图
N=1000;
L=size(imf,1);% Number of components in the decomposition
% loop for on each component
S=[];
clear x z m p freq
x=imf'; % now each column is a compenent
z=hilbert(x); % analytic signal
m=abs(z); % module of z
p=angle(z); % phase of z
for i=1：L
freq(:,i)=instfreq(z(:,i)); % instantaneous frequency
ceilfreq(:,i)=ceil(freq(:,i)*N);% to have integervalues-also do a smoothing given the number
   for j=1:length(x)-2
   S(ceilfreq(j,i),j+1)=m(j+1,i);
   end
end
eps=0.00001; % to avoid zero values before the log
E=20*log(S.^2+eps);% Hilbert energ spectrum
% plot S
figure(90);
t=1:length(x);% time sample
f=t/length(x)*25; % frequency %%%%%%%%%
imagesc(t,f,E);
contour(t(1:3087),1:497,E);%画出等高线图
surf1(1:3087,1:500,S.^2); shading interp; colormap(gray);axis()
% mesh(1:3087,1:500,S.^2);
% mesh(1:2999,1:500,S); %画出三维图
set(gca,'YDir','normal');
xlabel('t/s');
ylabel('Frequency/Hz'); %%%%%%%%%%%%%%%%%%%%
% ylabel('Normalized Frequency');
% axis(); %%%%%%%%%%%%%%%%%%%%%%
S1=sum(S.^2,1);figure(100);plot(1:3087,S1); %瞬时能量
figure(110);
S2=sum(S.^2,2);plot(1:497,S2,'r') %边际能量谱
figure(120);
S2=sum(S,2);plot(1:497,S2,'g');%边际谱

达人能否帮我解释一下上面粉色标注语句的意思啊？
我不知道这里的N代表什么意思？

页: [1]

Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛---(手机验证注册)'s Archiver

请教：关于emd.m的问题