经验模态分解股票波动matlab,LMD经验模态分解matlab程序——原味的

曾经也用滑动平均写过LMD，其实滑动平均的EMD才是原汁原味的居于均值分解。

分享给有需要的人，程序写的不好，只是希望提供一种思路。如果谁写了更完美LMD程序，别忘了发我一份，快毕业了，一直没有把LMD写完美，对于我来说始终是个遗憾。来分完美的LMD让我也品尝下，我也无憾了~

此处没有提供测试代码，如需要可以点这里：点我

源代码如下：

%原始lmd算法，效果很不好，不知道程序哪里写错

function[PF,A,SI]=lmd(m)

c=m;

k=0

wucha1=0.001;

n_l=nengliang(m);

while 1

k=k+1;

a=1;

h=c;

[pf,a,si]=zhaochun(a,h,wucha1);

c=c-pf;

PF(k,:)=pf;

A(k,:)=a;

SI(k,:)=si;

c_pos=pos(c);

n_c=nengliang(c);

n_pf=nengliang(pf);

if

length(c_pos)<3 || n_c

length(pos(pf))

n_pf

PF(k+1,:)=c;

break

end

end

function pos=pos(y)

%功能：找序列极值点位置坐标

%y:输入序列

%pos:极值点的序列位置坐标

m = length(y);

d = diff(y);

n = length(d);

d1 = d(1:n-1);

d2 = d(2:n);

indmin = find(d1.*d2<0 &

d1<0)+1;

indmax = find(d1.*d2<0 &

d1>0)+1;

if any(d==0)

imax = [];

imin = [];

bad = (d==0);

dd = diff([0 bad 0]);

debs = find(dd == 1);

fins = find(dd == -1);

if debs(1) == 1

if

length(debs) > 1

debs = debs(2:end);

fins = fins(2:end);

else

debs = [];

fins = [];

end

end

if length(debs) > 0

if fins(end)

== m

if length(debs) > 1

debs = debs(1:(end-1));

fins = fins(1:(end-1));

else

debs = [];

fins = [];

end end

end

lc = length(debs);

if lc > 0

for k =

1:lc

if d(debs(k)-1) > 0

if d(fins(k)) < 0

imax = [imax round((fins(k)+debs(k))/2)];

end

else

if d(fins(k)) > 0

imin = [imin round((fins(k)+debs(k))/2)];

end

end

end

end

if length(imax) > 0

indmax =

sort([indmax imax]);

end

if length(imin) > 0

indmin =

sort([indmin imin]);

end

end

minmax=cat(2,indmin,indmax);

pos=sort(minmax);

end

%----------zhaochun.m

function [pf,a,si]=zhaochun(a,h,wucha1)

chun_num=0;

while 1

chun_num=chun_num+1

t=1:length(h);

h_pos=position(h);%极值点位置序列

tpoint=t(h_pos);%极值点时间值

hpoint=h(h_pos);%极值点幅度值

hpoint=bianjie(h,hpoint,1);%端点处理后的极值点，多出了2个极值点

[mi,ai]=find_miai(hpoint);%找出极值点之间的均值函数和包络函数

mi_point=junzhi(mi);%所有极值点处的均值序列(幅值)-纵坐标(点序列)

ai_point=junzhi(ai);%所有极值点处的包络序列(幅值)-纵坐标(点序列)

%此处端点处的均值和包络都是

端点和极值点之间的均值和包络值(如果端点视作极值点，这样处理是不合理的，端点都不是极值点，这样处理是正确的)

lmi_point=mi(1);%左端点的均值

rmi_point=mi(length(mi));%右端点的均值

lai_point=ai(1);%左端点的包络

rai_point=ai(length(ai));%右端点的包络

%

mi_point_d=link(lmi_point,rmi_point,mi_point);%连接端点均值及所有极值点处的 均值

(带端点的均值序列)(点序列)

ai_point_d=link(lai_point,rai_point,ai_point);%连接端点包络及所有极值点出的

包络值(带端点的包络序列)(点序列)

%

tpoint_d=link(t(1),t(length(t)),tpoint);

mi_fun=chadian1(tpoint_d,mi,mi_point_d);%包含端点和极值点和普通点的

均值序列-平缓前的均值序列

ai_fun=chadian1(tpoint_d,ai,ai_point_d);%包含端点和极值点和普通点的

包络序列-平滑前的包络序列

%以上是完整的未处理的均值函数mi_fun和包络函数ai_fun

%找出第一平滑的滑动跨度

kmax=max(diff(tpoint_d));%找出时间跨度最大的 相邻几点 间的

距离

kmax1=uint16(kmax/3);

kmax1=double(kmax1);

jiou=uint8(rem(kmax1,2));

if kmax1<3

move_k=3; %

b<3,滑动跨度=3，

elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1

move_k=(kmax1-1);

else move_k=kmax1;

%b>=3,b是奇数时，跨度=b

end

[mi_move,move_mi_num]=move(mi_fun,move_k);

[ai_move,move_ai_num]=move(ai_fun,move_k);

mi=mi_move;

ai=ai_move;

a=a.*ai;

si=(h-mi)./ai;

h=si;

ai_funmax=max(ai);

ai_funmin=min(ai);

if

(ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1)

break

end

end

pf=a.*si;

end

function [x,flag]=move(a,k)

l=length(a);

t=1:l;

% jingdu=t(2)-t(1);

flag=0;

x=a;

max_flag=10;%平滑重复的最大次数设置

max_error=0;%相邻两点间相等允许的最大差异(理论上=0才人为无差异，实际操作要考虑采样精度，所以可以设置一个允许范围)

while flag<=max_flag ||

min(abs(diff(x)))<=max_error;%这里用到abs，然后再min，因为幅值的差值会出现负值，我们的目的是找差值=0，而不是负数

if

flag==0

flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次(最大10次)

else

flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次(最大10次)

x_pos=position(x);%极值点位置序列

tpoint=t(x_pos);%极值点时间值

tpoint_d=link(t(1),t(l),tpoint);%极值点的时间轴上的位置

kmax=max(diff(tpoint_d));%找出时间跨度最大的 相邻几点 间的

距离

kmax1=uint16(kmax/3);

kmax1=double(kmax1);

jiou=uint8(rem(kmax1,2));

if kmax1<3

k=3; % b<3,滑动跨度=3，

elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1

k=(kmax1-1);

else k=kmax1; %b>=3,b是奇数时，跨度=b

end

end

y=zeros(1,l);% 初始化序列y, y是中间变量

%%%%%%%%%%%%%%%%%%%边界点的处理%%%%%%%%%%%%%%

for

i=1:(k-1)/2

v=0;

z=0;

for j=1:i*2-1

v=v+x(j);

end

y(i)=v/(i*2-1);

for

j=l:-1:l+2-i*2 %j=l:-1:l+1-(i*2-1)

z=z+x(j);

end

y(l+1-i)=z/(i*2-1);

end

%%%%%%%%%%%%%%中间点的处理 %%%%%%%%%%%%%%%%%%%%%%%%%%

for

i=1+(k-1)/2:l-(k-1)/2 %这个 (d+1)/2是跨度的中点

单边的边界点数=中点值-1=(d+1)/2-1=(d-1)/2