EMD算法的matlab程序介紹

1、%此版本為 ALAN 版本的整合注釋版function imf =emd(x%Empiricial Mode Decomposition (Hilbert-HuangTransform%imf =emd(x%Func :findpeaksx=transpose(x(:;%轉(zhuǎn)置為行矩陣imf =;while ismonotonic(x%當(dāng) x 不是單調(diào)函數(shù),分解終止條件x1=x;sd =Inf;%均值%直到 x1滿足 IMF 條件,得 c1while (sd>0.1 |isimf(x1%當(dāng)標(biāo)準(zhǔn)偏差系數(shù) sd 大于 0.1或 x1不是固有模態(tài)函數(shù)時,分 量終止條件s1=getspline(x

2、1;%上包絡(luò)線s2=-getspline(-x1;%下包絡(luò)線x2=x1-(s1+s2/2;%此處的 x2為文章中的 hsd =sum(x1-x2.2/sum(x1.2;x1=x2;endimfend+1=x1;x =x-x1;endimfend+1=x;%FUNCTIONSfunction u =ismonotonic(x%u=0表示 x 不是單調(diào)函數(shù), u=1表示 x 為單調(diào)的u1=length(findpeaks(x*length(findpeaks(-x;if u1>0, u =0;else, u =1; endfunction u =isimf(x%u=0表示 x 不是固有模式函

3、數(shù), u=1表示 x 是固有模式函數(shù)N =length(x;u1=sum(x(1:N-1.*x(2:N<0;u2=length(findpeaks(x+length(findpeaks(-x;if abs(u1-u2>1, u =0;else, u =1; endfunction s =getspline(x%三次樣條函數(shù)擬合成元數(shù)據(jù)包絡(luò)線N =length(x;p =findpeaks(x;s =spline(0p N+1,0x(p0,1:N;-function n =findpeaks(x%Find peaks. 找到極值 ,n 為極值點所在位置%n =findpeaks(xn

4、 =find(diff(diff(x>0 <0;u =find(x(n+1>x(n;n(u=n(u+1;-function plot_hht00(x,Ts%雙邊帶調(diào)幅信號的 EMD 分解%Plot the HHT.%plot_hht(x,Ts%: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%Get HHT

5、.clear all;close all;Ts=0.0005;t=0:Ts:10;%采樣率 2000HZ%調(diào)幅信號%x=sin(2*pi*t.*sin(40*pi*t;x=sin(2*pi*t;s1=getspline(x;%上包絡(luò)線s2=-getspline(-x;%上包絡(luò)線x1=(s1+s2/2;%此處的 x2為文章中的 hfigure;plot(t,x;xlabel('Time',ylabel('Amplitude'title('雙邊帶調(diào)幅信號 'hold on;plot(t,s1,'-r'plot(t,s2,'-r

6、'plot(t,x1,'g'imf =emd(x;for k =1:length(imfb(k=sum(imfk.*imfk;th =angle(hilbert(imfk;dk=diff(th/Ts/(2*pi;endu,v=sort(-b;b =1-b/max(b;%Set time-frequency plots.N =length(x;c =linspace(0,(N-2*Ts,N-1;%figure;for k =v(1:2plot(c,dk,'k.','Color',b(kk k,'MarkerSize',3;hold on;set(gca,'FontSize',8,'XLim',0c(end,'YLim',050;%設(shè)置 x 、 y 軸句柄 xlabel('Time',ylabel('Frequency'title('原信號時頻圖 'end%Set IMF plots.M =length(imf;N =length(x;c =linspace(0,(N-1*Ts,N;for k1=0:4:M-1figurefor k2=1:min(4,M-k1,subplot(4,1,k2,plot(c,i