急求HHT运行的EMD分解出的图没有坐标显示？
本帖最后由 zt861217 于 2011-4-14 16:26 编辑以下是HHT程序，用到的txt文件在附件里。求各位高手给想想办法，我是新手，还有急用。
% EMD_FMSIN.M
%
% P. Flandrin, Mar. 13, 2003 - modified Mar. 2, 2006
%
% computes and displays EMD for the sum of 2 sinusoidal
% FM's + 1 Gaussian logon
%
% displays reassigned spectrograms of the sum signal and of
% the 3 first modes extracted
%
% produces Figure 1 in
%
% G. Rilling, P. Flandrin and P. Gon峚lv弒
% "On Empirical Mode Decomposition and its algorithms"
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
clear;
NN=660;
%fid=fopen('g2204.375.txt','rt');
%fid=fopen('2174.5.TXT','rt');%泥岩
%fid=fopen('2174.75.txt','rt');%泥岩
%fid=fopen('2174.625.txt','rt');%泥岩
%fid=fopen('2270.TXT','rt');%油水层
%fid=fopen('2270.125.txt','rt');%油水层
%fid=fopen('2270.25.txt','rt');%油水层
%fid=fopen('s2176.875.txt','rt');
fid=fopen('g2204.25.txt','rt');
%NN=360;
%fid=fopen('460.txt','rt');
%%fid=fopen('495a.txt','rt');
%fid=fopen('Test_520ab.txt','rt');
%fid=fopen('Test_600.txt','rt');
%fid=fopen('Test_605.txt','rt');
%fid=fopen('670.txt','rt');
%fid=fopen('842.txt','rt');
%fid=fopen('1000.txt','rt');
%fid=fopen('1140.txt','rt');
x=fscanf(fid,'%d',NN);
fclose(fid);
x=x-x(1);
T = 1:1:NN;
t = 1:NN;
%N = 2000;% # of data samples
%p = N/2;% period of the 2 sinusoidal FM's
% sinusoidal FM 1
%fmin1 = 1/64;% min frequency
%fmax1 = 1.5*1/8;% max frequency
%x1 = fmsin(N,fmin1,fmax1,p,N/2,fmax1);
% sinusoidal FM 1
%fmin2 = 1/32;% min frequency
%fmax2 = 1.5*1/4;% max frequency
%x2 = fmsin(N,fmin2,fmax2,p,N/2,fmax2);
% logon
%f0 = 1.5*1/16;% center frequency
%x3 = amgauss(N,N/2,N/8).*fmconst(N,f0);
%a1 = 1;
%a2 = 1;
%a3 = 1;
%x = real(a1*x1+a2*x2+a3*x3);
%x = x/max(abs(x));
= emd(x); % <-------- adapted to current emd.m
% 可以在这里进行滤波
%imf(6,:)=0;
%imf(5,:)=0;
%imf(4,:)=0;
%imf(3,:)=0;
%imf(2,:)=0;
%imf(1,:)=0;
%=size(imf);
%imf(Nimf-4,:)=0;
%imf(Nimf-3,:)=0;
%imf(Nimf-2,:)=0;
%imf(Nimf-1,:)=0;
%imf(Nimf,:)=0;
emd_visu(x,t,imf,1);
figure(1)
% = hhspectrum(imf,t,l,aff)
= hhspectrum(imf,t,2);
%SUCCESS3=XLSWRITE('HHT_IF.xls',f');
% = toimage(A,f,t,splx,sply)
splx = 178;
% number of columns of the output im (time resolution).
% If different from length(t), works only for uniform sampling.
sply = 178;
% number of rows of the output im (frequency resolution).
= toimage(A,f,tt,splx,sply);
disp_hhs(im,Ttt);
DeltaFF = 2* pi * 12 / 1000.;
BB=sum(im');
BBA = sum(BB);
BB = (BB / BBA);
=size(BB);
%h = SIZEOF(BB)
ABBA = 0;
for i = (Mbb-28):-1:(Mbb-50)
ABBA = ABBA + BB(i);
end
ABBA
figure;
plot(Ttt*DeltaFF,fliplr(BB)');
%plot(fliplr(BB)');
%axis()
xlabel('频率 / kHz')
%xlabel('时间 / us')
ylabel(['幅值'])
%ylabel(['HHT的瞬时能量'])
title('Hilbert边际谱')
%title('Hilbert瞬时能量谱')
imm = im.^2;
AAA = sum(imm);
AAB = sum(AAA);
AAA =( AAA / AAB);
%h = SIZEOF(AAA)
figure;
plot(Ttt*12,(AAA)');
%xlabel('频率 / kHz')
xlabel('时间 / μs')
%ylabel(['幅值'])
ylabel(['HHT的瞬时能量'])
%title('Hilbert边际谱')
title('Hilbert瞬时能量谱')
AAAA = sum(imm');
AAAB = sum(AAAA);
AAAA =( AAAA / AAAB);
figure;
%h = SIZEOF(AAAA)
%xxmin = 0;
%xxmax = Ttt*DeltaFF;
%yymin = 0;
%yymax = 20000000;
plot(Ttt*DeltaFF,fliplr(AAAA)');
%plot(Ttt,fliplr(AAAA)')
%axis()
xlabel('频率 / kHz')
%xlabel('时间 / us')
ylabel(['能量'])
%ylabel(['HHT的瞬时能量'])
%title('Hilbert边际谱')
title('Hilbert能量谱')
clear
我才是新手，都不知道HHT怎么用MATLAB实现 我的也没有坐标急求 这个要坐标自己加一下就行了，跟初始的振动信号坐标一样
用axis自己加到EMD文件里去
不过这没什么意思
加了之后图出来就不好看了
看个趋势就够了
看看怎么实现！ {:{10}:}{:{10}:}
页:
[1]