LMD经验模态分解matlab程序要点.docx
- 文档编号:12939378
- 上传时间:2023-06-09
- 格式:DOCX
- 页数:38
- 大小:402.93KB
LMD经验模态分解matlab程序要点.docx
《LMD经验模态分解matlab程序要点.docx》由会员分享,可在线阅读,更多相关《LMD经验模态分解matlab程序要点.docx(38页珍藏版)》请在冰点文库上搜索。
LMD经验模态分解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);
while1
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);
iflength(c_pos)<3||n_c PF(k+1,: )=c; break end end functionpos=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; ifany(d==0) imax=[]; imin=[]; bad=(d==0); dd=diff([0bad0]); debs=find(dd==1); fins=find(dd==-1); ifdebs (1)==1 iflength(debs)>1 debs=debs(2: end); fins=fins(2: end); else debs=[]; fins=[]; end end iflength(debs)>0 iffins(end)==m iflength(debs)>1 debs=debs(1: (end-1)); fins=fins(1: (end-1)); else debs=[]; fins=[]; end end end lc=length(debs); iflc>0 fork=1: lc ifd(debs(k)-1)>0 ifd(fins(k))<0 imax=[imaxround((fins(k)+debs(k))/2)]; end else ifd(fins(k))>0 imin=[iminround((fins(k)+debs(k))/2)]; end end end end iflength(imax)>0 indmax=sort([indmaximax]); end iflength(imin)>0 indmin=sort([indminimin]); end end minmax=cat(2,indmin,indmax); pos=sort(minmax); end %----------zhaochun.m function[pf,a,si]=zhaochun(a,h,wucha1) chun_num=0; while1 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)); ifkmax1<3 move_k=3;%b<3,滑动跨度=3, elseifjiou==0%b>=3,当b是偶数时,跨度=b-1 move_k=(kmax1-1); elsemove_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才人为无差异,实际操作要考虑采样精度,所以可以设置一个允许范围) whileflag<=max_flag||min(abs(diff(x)))<=max_error;%这里用到abs,然后再min,因为幅值的差值会出现负值,我们的目的是找差值=0,而不是负数 ifflag==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)); ifkmax1<3 k=3;%b<3,滑动跨度=3, elseifjiou==0%b>=3,当b是偶数时,跨度=b-1 k=(kmax1-1); elsek=kmax1;%b>=3,b是奇数时,跨度=b end end y=zeros(1,l);%初始化序列y, y是中间变量 %%%%%%%%%%%%%%%%%%%边界点的处理%%%%%%%%%%%%%% fori=1: (k-1)/2 v=0; z=0; forj=1: i*2-1 v=v+x(j); end y(i)=v/(i*2-1); forj=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 %%%%%%%%%%%%%%中间点的处理%%%%%%%%%%%%%%%%%%%%%%%%%% fori=1+(k-1)/2: l-(k-1)/2%这个(d+1)/2是跨度的中点单边的边界点数=中点值-1=(d+1)/2-1=(d-1)/2 w=x(i); forj=1: (k-1)/2 w=w+x(i-j)+x(i+j); end y(i)=w/k; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x=y; end %k % 调试的时候查看 %flag% 调试的时候查看 end functionhpp=bianjie(hh,hp,style) %hh: 需要边界处理的序列 %hp: 需要边界处理序列的极值点的幅度值 %style: 边界处理类型1: 端点全部看作极值点2: 左右两端各增加1个幅值为0的极值点 %hpp: 边界处理后,极值点幅值多处了两个,因为左边右边各人为加上了一个 ifstyle==1 a=hh (1); b=hh(length(hh)); end ifstyle==2 a=0; b=0; end hpp=link(a,b,hp); end functionpos=position(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; ifany(d==0) imax=[]; imin=[]; bad=(d==0); dd=diff([0bad0]); debs=find(dd==1); fins=find(dd==-1); ifdebs (1)==1 iflength(debs)>1 debs=debs(2: end); fins=fins(2: end); else debs=[]; fins=[]; end end iflength(debs)>0 iffins(end)==m iflength(debs)>1 debs=debs(1: (end-1)); fins=fins(1: (end-1)); else debs=[]; fins=[]; end end end lc=length(debs); iflc>0 fork=1: lc ifd(debs(k)-1)>0 ifd(fins(k))<0 imax=[imaxround((fins(k)+debs(k))/2)]; end else ifd(fins(k))>0 imin=[iminround((fins(k)+debs(k))/2)]; end end end end iflength(imax)>0 indmax=sort([indmaximax]); end iflength(imin)>0 indmin=sort([indminimin]); end end minmax=cat(2,indmin,indmax); pos=sort(minmax); end function[mi,ai]=find_miai(ypoint) % Lyp=length(ypoint); al=wkeep(ypoint,Lyp-1,'l'); ar=wkeep(ypoint,Lyp-1,'r'); mi=(al+ar)/2; ai=abs(ar-al)/2; end functionc=junzhi(a) l=length(a); al=wkeep(a,l-1,'l'); ar=wkeep(a,l-1,'r'); c=(al+ar)/2; end functiond=link(a,b,c) d=[a';c';b']'; %头: 尾: 中间 end functionf_value=chadian1(a,b,c) %chadian1把端点及极值点处对应的总坐标值插入,原来的均值函数的方波序列中 %输入参数a: 点序列(行向量)(包含端点和极值点以在时间上的位置-横坐标)(n个)(点序列-横坐标) %输入参数b: 段序列(行向量)(点序列a每两点之间的纵坐标的值-纵坐标)(n-1)-(段序列-纵坐标) %输入参数c: 点序列(行向量)(包含端点和极值点在对应时间上的幅值-纵坐标)(n个)-(点序列-纵坐标) %输入参数d: 原始数据的采样精度 %输出参数f_value(行向量): 点序列a插入段序列的值之后,以c的精度的值(对应于横坐标,纵坐标的值) %精度是0.001 l=length(b); al=wkeep(a,l,'l'); ar=wkeep(a,l,'r'); d={[]};%d={0}这样是为了初始化元胞数组 fori=1: l %采样精度0.001 d{i}=ones(1,(ar(i)-al(i)-1))*b(i);%ones函数参数要为整数,unint16就是数据强制类型转换, end %这里没有使用到单独为uint16((ar(i)-al(i))*1000)-1)这个自变量赋值,所以只是个中间变量,对数据不会产生污染 y=c (1); fori=1: l y=link(y,c(i+1),d{i}); end f_value=y; end end %------ functionp=nengliang(y) %my=mean(y); %p=(y-my).*(y-my); %p=sum(p); p=sum(abs(y).^2); end %-------- 求一维序列的信息熵(香浓熵)的matlab程序实例 对于一个二维信号,比如灰度图像,灰度值的范围是0-255,因此只要根据像素灰度值(0-255)出现的概率,就可以计算出信息熵。 但是,对于一个一维信号,比如说心电信号,数据值的范围并不是确定的,不会是(0-255)这么确定,如果进行域值变换,使其转换到一个整数范围的话,就会丢失数据,请高手指点,怎么计算。 比如数字信号是x(n),n=1~N (1)先用Hist函数对x(n)的赋值范围进行分块,比如赋值范围在0~10的对应第 一块,10~20的第二块,以此类推。 这之前需要对x(n)做一些归一化处理 (2)统计每一块的数据个数,并求出相应的概率 (3)用信息熵公式求解 以上求解方法获得的虽然是近似的信息熵,但是一般认为,这么做是没有问题的 求一维序列的信息熵的matlab程序代码如下: (已写成调用的函数形式) 测试程序: fs=12000; N=12000; T=1/fs; t=(0: N-1)*T; ff=104; sig=0.5*(1+sin(2*pi*ff*t)).*sin(2*pi*3000*t)+rand(1,length(t)); Hx=yyshang(sig,10) %———————求一维离散序列信息熵matlab代码 functionHx=yyshang(y,duan) %不以原信号为参考的时间域的信号熵 %输入: maxf: 原信号的能量谱中能量最大的点 %y: 待求信息熵的序列 %duan: 待求信息熵的序列要被分块的块数 %Hx: y的信息熵 %duan=10;%将序列按duan数等分,如果duan=10,就将序列分为10等份 x_min=min(y); x_max=max(y); maxf (1)=abs(x_max-x_min); maxf (2)=x_min; duan_t=1.0/duan; jiange=maxf (1)*duan_t; %fori=1: 10 %pnum(i)=length(find((y_p>=(i-1)*jiange)&(y_p %end pnum (1)=length(find(y (2)+jiange)); fori=2: duan-1 pnum(i)=length(find((y>=maxf (2)+(i-1)*jiange)&(y (2)+i*jiange))); end pnum(duan)=length(find(y>=maxf (2)+(duan-1)*jiange)); %sum(pnum) ppnum=pnum/sum(pnum);%每段出现的概率 %sum(ppnum) Hx=0; fori=1: duan ifppnum(i)==0 Hi=0; else Hi=-ppnum(i)*log2(ppnum(i)); end Hx=Hx+Hi; end end %---------------- 扩展阅读: 实验一: 计算离散信源的熵 一、实验设备: 1、计算机 2、软件: Matlab 二、实验目的: 1、熟悉离散信源的特点; 2、学习仿真离散信源的方法 3、学习离散信源平均信息量的计算方法 4、熟悉 Matlab 编程; 三、实验内容: 1、写出计算自信息量的Matlab 程序 2、写出计算离散信源平均信息量的Matlab 程序。 3、掌握二元离散信源的最大信息量与概率的关系。 4、将程序在计算机上仿真实现,验证程序的正确性并完成习题。 四、实验报告要求 简要总结离散信源的特点及离散信源平均信息量的计算,写出习题的MATLAB实现语句。 信息论基础: 自信息的计算公式 Matlab实现: I=log2(1/p) 或I=-log2(p) 熵(平均自信息)的计算公式 Matlab实现: HX=sum(-x.*log2(x));或者h=h-x(i)*log2(x(i)); 习题: 1. 甲地天气预报构成的信源空间为: 乙地信源空间为: 求此两个信源的熵。 求各种天气的自信息量。 案: 运行程序: p1=[1/2,1/4,1/8,1/8];%p1代表甲信源对应的概率 p2=[7/8,1/8];%p2代表乙信源对应的概率 H1=0.0; H2=0.0; I=[]; J=[]; fori=1: 4 H1=H1+p1(i)*log2(1/p1(i)); I(i)=log2(1/p1(i)); end disp('自信息量分别为: '); I disp('H1信源熵为: '); H1 forj=1: 2 H2=H2+p2(j)*log2(1/p2(j)); J(j)=log2(1/p2(j)); end disp('自信息量分别为: '); J disp('H2信源熵为: '); H2 emd经验模态分解matla
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- LMD 经验 分解 matlab 程序 要点