EMD算法的matlab程序介绍解析.docVIP

%此版本为 ALAN 版本的整合注释版 function imf =emd(x %Empiricial Mode Decomposition (Hilbert-HuangTransform %imf =emd(x %Func :findpeaks x=transpose(x(:;%转置为行矩阵 imf =[]; while ~ismonotonic(x%当 x 不是单调函数,分解终止条件 x1=x; sd =Inf;%均值 %直到 x1满足 IMF 条件,得 c1 while (sd0.1 |~isimf(x1%当标准偏差系数 sd 大于 0.1或 x1不是固有模态函数时,分 量终止条件 s1=getspline(x1;%上包络线 s2=-getspline(-x1;%下包络线 x2=x1-(s1+s2/2;%此处的 x2为文章中的 h sd =sum((x1-x2.^2/sum(x1.^2; x1=x2; end imf{end+1}=x1; x =x-x1; end imf{end+1}=x; %FUNCTIONS function u =ismonotonic(x %u=0表示 x 不是单调函数, u=1表示 x 为单调的 u1=length(findpeaks(x*length(findpeaks(-x; if u10, u =0; else, u =1; end function u =isimf(x %u=0表示 x 不是固有模式函数, u=1表示 x 是固有模式函数 N =length(x; u1=sum(x(1:N-1.*x(2:N0; u2=length(findpeaks(x+length(findpeaks(-x; if abs(u1-u21, u =0; else, u =1; end function s =getspline(x %三次样条函数拟合成元数据包络线 N =length(x; p =findpeaks(x; s =spline([0p N+1],[0x(p0],1:N; ------------------------------------------------------------------------------- -------------------------------------------------------------------------------- function n =findpeaks(x %Find peaks. 找到极值 ,n 为极值点所在位置 %n =findpeaks(x n =find(diff(diff(x0 0; u =find(x(n+1x(n; n(u=n(u+1; ------------------------------------------------------------------------------------------ ---------------------------------------------------------------------------------------- function plot_hht00(x,Ts %双边带调幅信号的 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. clear all; close all; Ts=0.0005; t=0:Ts:10;%采样率 2000HZ %调幅信号 %x=sin(2*pi*t.*sin(40*pi*t; x=sin(2*pi*t; s1=getspline(x;%上包络线 s2=-getspline(-x;%上包络线 x1=(s1+s2/2;%此处的 x2为文章中的 h figure; plot(t,x;xlabel(Time,ylabel(Amplitude;title(双边带调幅信号 ;hold on; plot(t,s1,-r; plot(t,s2,-r; plot(t,x1,g; imf =emd(x; for k =1:length(imf b(k=sum(imf{k}.*imf{k}; th =angle(hilbert(imf{k}; d{k}=diff(th/Ts/(2*pi; end [u,v]=sort(-b; b =1-b/max(b;