有关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

---