《EMD算法的matlab程序介绍解析》由会员分享,可在线阅读,更多相关《EMD算法的matlab程序介绍解析(6页珍藏版)》请在人人文库网上搜索。
1、此版本为 ALAN 版本的整合注释版function imf =emd(x%Empiricial Mode Decomposition (Hilbert-HuangTransform%imf =emd(x%Func :findpeaksx=transpose(x(:;%转置为行矩阵imf =;while ismonotonic(x%当 x 不是单调函数,分解终止条件x1=x;sd =Inf;%均值%直到 x1满足 IMF 条件,得 c1while (sd0.1 |isimf(x1%当标准偏差系数 sd 大于 0.1或 x1不是固有模态函数时,分 量终止条件s1=getspline(x1;%上包。
2、络线s2=-getspline(-x1;%下包络线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 不是单调函数, u=1表示 x 为单调的u1=length(findpeaks(x*length(findpeaks(-x;if u10, u =0;else, u =1; endfunction u =isimf(x%u=0表示 x 不是固有模式函数, u=1表示 。
3、x 是固有模式函数N =length(x;u1=sum(x(1:N-1.*x(2:N1, u =0;else, u =1; endfunction 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(xn =find(diff(diff(x0 x(n;n(u=n(u+1;-function plot_hht00(x,Ts%双边带调幅信号的 EMD 。
4、分解%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。
5、(x;%上包络线s2=-getspline(-x;%上包络线x1=(s1+s2/2;%此处的 x2为文章中的 hfigure;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(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-frequenc。
6、y 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;%设置 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,imfk1+k2;set(gca,FontSize,8,XLim,0c(end;title(EMD分解结果 ;endxlabel(Time;end。