EMDMatlab程序

function plot_hht(x,imf,Ts) % Plot the HHT. % :: Syntax % The array x is the input signal and Ts is the sampling period. % Example on use: [x,Fs] = wavread('Hum.wav'); % plot_hht(x(1:6000),1/Fs); % Func : emd % imf = emd(x); for k = 1:length(imf) b(k) = sum(imf{k}.*imf{k}); th = unwrap(angle(hilbert(imf{k}))); % 相位 d{k} = diff(th)/Ts/(2*pi); % 瞬时频率 end [u,v] = sort(-b); b = 1-b/max(b); % 后面绘图的亮度控制 % Hilbert瞬时频率图 N = length(x); c = linspace(0,(N-2)*Ts,N-1); % 0:Ts:Ts*(N-2) for k = v(1:2) % 显示能量最大的两个IMF的瞬时频率 figure plot(c,d{k}); xlim([0 c(end)]); ylim([0 1/2/Ts]); xlabel('Time/s') ylabel('Frequency/Hz'); title(sprintf('IMF%d', k)) end % 显示各IMF M = length(imf); N = length(x); c = linspace(0,(N-1)*Ts,N); % 0:Ts:Ts*(N-1) for k1 = 0:4:M-1 figure for k2 = 1:min(4,M-k1) subplot(4,2,2*k2-1) plot(c,imf{k1+k2}) set(gca,'FontSize',8,'XLim',[0 c(end)]); title(sprintf('第%d个IMF', k1+k2)) xlabel('Time/s') ylabel(sprintf('IMF%d', k1+k2)); subplot(4,2,2*k2) [yf, f] = FFTAnalysis(imf{k1+k2}, Ts); plot(f, yf) set(gca,'FontSize',8,'XLim',[0 5000]); title(sprintf('第%d个IMF的频谱', k1+k2)) xlabel('f/Hz') ylabel('|IMF(f)|'); end end figure subplot(211) plot(c,x) set(gca,'FontSize',8,'XLim',[0 c(end)]); title('降噪信号') xlabel('Time/s') ylabel('Origin'); subplot(212) [Yf, f] = FFTAnalysis(x, Ts); plot(f, Yf) title('降噪信号的频谱') xlabel('f/Hz') ylabel('|Y(f)|');