有关EMD分解matlab程序

一般的的查询可在matlab里的帮助界面进行搜索，点击帮助。

打开帮助页面，左侧检索栏进行检索需要查询的语句，然后即可查看右侧谨老孝查询结果。

或者在主界面，输入help 空格+你要查询的内容，进行查询。下次你可以尝试一下。

一般程序都会有不懂得语句，或没用过的，可以在刚刚说过的帮助页面进行查询，看如何使用，输入参量什么意义，程序输出结果是什么。

你这个描述祥稿太过于简略了，含氏我只能给出这样的答案了，一般可以都给一点程序，也许会有对答题人帮助。

不过你这个程序的确有点复杂，不根据前后逻辑，和主程序和目的是很难解答的。

你可以看看主程序，再查查帮助。

希望对你有所帮助。谢谢。

%此版本为ALAN 版本的整合注释版

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% imf = emd(x)

% Func : findpeaks

x= transpose(x(:))%转置基颂为行矩阵

imf = []

while ~ismonotonic(x) %当x不是单调函数，分解终止条件

x1 = x

sd = Inf%均值

%直到x1满足IMF条件，得c1

while (sd >0.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 u1 >0, u = 0

else, u = 1end

function u = isimf(x)

%u=0表示x不是固有模式函数，u=1表示x是固有模式函数

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 = 1end

function s = getspline(x)

%三次样条函数拟合成元数据包络线

N = length(x)

p = findpeaks(x)

s = spline([0 p N+1],[0 x(p) 0],1:N)

-------------------------------------------------------------------------------

--------------------------------------------------------------------------------

function n = findpeaks(x)

% Find peaks.找到极值 ,n为极值点所在位置

% n = findpeaks(x)

n= find(diff(diff(x) >0) <0)

u= find(x(n+1) >x(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)

% Set time-frequency plots.

N = length(x)

c = linspace(0,(N-2)*Ts,N-1)

%

figure

for k = v(1:2)

plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3)hold on

set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50])%设置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-1

figure

for k2 = 1:min(4,M-k1),

subplot(4,1,k2),

plot(c,imf{k1+k2})

set(gca,'FontSize',8,'XLim',[0 c(end)])

title('EMD分解结果')

end

xlabel('Time')

end

function imf = emd(x,n)%%最好把函数名改为emd1之类的，以免和Grilling的emd冲突

%%n为你想得到的IMF的个数

c = x('% copy of the input signal (as a row vector)

N = length(x)-

% loop to decompose the input signal into n successive IMFs

imf = []% Matrix which will contain the successive IMF, and the residuefor t=1:n

% loop on successive IMFs

%-------------------------------------------------------------------------

% inner loop to find each imf

h = c% at the beginning of the sifting process, h is the signal

SD = 1% Standard deviation which will be used to stop the sifting process

while SD >0.3 % while the standard deviation is higher than 0.3 (typical value) %%筛选停止准毕坦慎则

% find local max/min points

d = diff(h)% approximate derivative %%求各点手敬导数

maxmin = []% to store the optima (min and max without distinction so far)

for i=1:N-2

if d(i)==0% we are on a zero %%导数信陪为0的点，即”驻点“，但驻点不一定都是极值点，如y=x^3的x=0处

if sign(d(i-1))~=sign(d(i+1)) % it is a maximum %%如果驻点两侧的导数异号（如一边正，一边负），那么该点为极值点

maxmin = [maxmin, i]%%找到极值点在信号中的坐标（不分极大值和极小值点）

end

elseif sign(d(i))~=sign(d(i+1)) % we are straddling a zero so%%如y=|x|在x=0处是极值点，但该点倒数不存在，所以不能用上面的判

断方法

maxmin = [maxmin, i+1] % define zero as at i+1 (not i) %%这里提供了另一类极值点的判断方法

end

end

if size(maxmin,2) <2 % then it is the residue %%判断信号是不是已经符合残余分量定义

break

end

% divide maxmin into maxes and mins %% 分离极大值点和极小值点

if maxmin(1)>maxmin(2) % first one is a max not a min

maxes = maxmin(1:2:length(maxmin))

mins = maxmin(2:2:length(maxmin))

else% is the other way around

maxes = maxmin(2:2:length(maxmin))

mins = maxmin(1:2:length(maxmin))

end% make endpoints both maxes and mins

maxes = [1 maxes N]

mins = [1 mins N]

%------------------------------------------------------------------------- % spline interpolate to get max and min envelopesform imf

maxenv = spline(maxes,h(maxes),1:N) %%用样条函数插值拟合所有的极大值点

minenv = spline(mins, h(mins),1:N)%%用样条函数插值拟合所有的极小值点

m = (maxenv + minenv)/2% mean of max and min enveloppes %%求上下包络的均值

prevh = h% copy of the previous value of h before modifying it %%h为分解前的信号

h = h - m% substract mean to h %% 减去包络均值

% calculate standard deviation

eps = 0.0000001% to avoid zero values

SD = sum ( ((prevh - h).^2) ./ (prevh.^2 + eps) )%% 计算停止准则

end

imf = [imfh]% store the extracted IMF in the matrix imf

% if size(maxmin,2)<2, then h is the residue

% stop criterion of the algo. if we reach the end before n

if size(maxmin,2) <2

break

end

c = c - h% substract the extracted IMF from the signal

end

return

---