- 积分
- 0
- 注册时间
- 2008-8-20
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 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的时频图
[imf,ort,nbits]=emd2(x,t,[0.05,0.5,0.05]);
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 integer values-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([1 3090 1 500 0 0.004])
% 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([0 772 0 300]); %%%%%%%%%%%%%%%%%%%%%%
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已经加入附件 |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
|