网站模版 源码,做租号玩网站赚钱吗,做公司网站的理念,网站专项审批查询说明#xff1a;研究LMD局域均值分解有3个月左右#xff0c;能找到的相关文章也基本上看了一遍#xff0c;觉得是个很好的方法#xff0c;号称是EMD经验模态分解的改进版。但是网络上一直没有找到该算法的matlab程序#xff0c;只见文章说的天花乱坠。后来自己写了一个研究LMD局域均值分解有3个月左右能找到的相关文章也基本上看了一遍觉得是个很好的方法号称是EMD经验模态分解的改进版。但是网络上一直没有找到该算法的matlab程序只见文章说的天花乱坠。后来自己写了一个但是使用有一个问题没有解决就是分解的时候怎么去除骑行波的问题。先把自己写的程序贡献出来让大家分析下一起讨论下到底LMD的程序怎样写才能如文献中说的那样达到目的。欢迎大家热烈讨论程序仍存在不收敛的问题拿出来分享只是希望高手指点一二写的不好欢迎拍砖文件夹包含找出纯调频函数计算瞬时幅值瞬时频率的函数等%%%%%%%%%%%主程序%%%%%%%%%%%lmd1原始版%通过emd.m的三次样条镜像延拓求出上下包络及均值%局域均值函数(上包络下包络)/2%局域包络函数|上包络-下包络|/2%相关文章见《一种基于EMD的振动信号时频分析新方法研究》-胡劲松杨世锡function[PF,A,SI]lmd1(m)%最后一个PF是残余分量%A是瞬时赋值%SI是纯调频函数求它的瞬时频率就是需要的频率cm;k0;wucha10.001;n_lnengliang(m);while 1kk1;a1;hc;[pf,a,si]zhaochun1(a,h,wucha1);cc-pf;PF(k,:)pf;A(k,:)a;SI(k,:)si;c_pospos(c);n_cnengliang(c);n_pfnengliang(pf);%停止调节%1.emd用的是三次样条求包络要求至少3个极值点所以这里c的极值点个数也应该至少为3%2.如果上一个PF的极值点数比下一个PF的极值点数少说明结果也不正确(这个也可以作为停止条件考虑进去)%上面一句是否可以等价于当前PF的极值点个数一定要大于等于残量(c)的极值点个数(目前是用这个作为停止条件的一个参考写入程序)%3.当前PF分量的能量应该大于残量c的能量(这个有待商榷)%4.残余能量不能大于信号能量ifn_pfiflength(c_pos)3 ||n_cn_pfn_lPF(k1,:)c;breakendendend%%%%%%%%%%%%%%%%%%nengliang函数%%%%%%%%%%function pnengliang(y)% mymean(y);% p(y-my).*(y-my);% psum(p);psum(abs(y).^2);end%%%%%%%%%%%%%找纯函数%%%%%%%%%%%%%%%%%function [pf,a,si]zhaochun1(a,h,wucha1)chun_num0;while 1chun_numchun_num1;t1:length(h);[envmin,envmax,envmoy,indmin,indmax,indzer] envelope(t,h,spline);mi(envmaxenvmin)./2;aiabs(envmax-envmin)./2;aa.*ai;si(h-mi)./ai;hsi;ai_funmaxmax(ai);ai_funminmin(ai);if(ai_funmax1wucha1ai_funmin1-wucha1)breakendendpfa.*si;chun_numfunction [envmin, envmax,envmoy,indmin,indmax,indzer] envelope(t,x,INTERP)%computes envelopes and mean with various interpolationsNBSYM 2; % 边界延拓点数DEF_INTERP spline;if nargin 2x t;t 1:length(x);INTERP DEF_INTERP;endif nargin 2if ischar(x)INTERP x;x t;t 1:length(x);endendif ~ischar(INTERP)error(interp parameter must be linear,cubic or spline)endif ~any(strcmpi(INTERP,{linear,cubic,spline}))error(interp parameter must be linear,cubic or spline)endif min([size(x),size(t)]) 1error(x and t must be vectors)ends size(x);if s(1) 1x x;ends size(t);if s(1) 1t t;endif length(t) ~ length(x)error(x and t must have the same length)endlx length(x);[indmin,indmax,indzer] extr(x,t);%boundary conditions for interpolation[tmin,tmax,xmin,xmax] boundary_conditions(indmin,indmax,t,x,NBSYM);% definition of envelopes from interpolationenvmax interp1(tmax,xmax,t,INTERP); envmin interp1(tmin,xmin,t,INTERP);if nargout 2envmoy (envmax envmin)/2;endendfunction [tmin,tmax,xmin,xmax] boundary_conditions(indmin,indmax,t,x,nbsym)% computes the boundary conditions for interpolation (mainly mirrorsymmetry)lx length(x);% 判断极值点个数if (length(indmin) length(indmax) 3)error(not enoughextrema)end% 插值的边界条件if indmax(1) indmin(1)%第一个极值点是极大值if x(1) x(indmin(1))%以第一个极大值为对称中心lmax fliplr(indmax(2:min(end,nbsym1)));lmin fliplr(indmin(1:min(end,nbsym)));lsym indmax(1);else%如果第一个采样值小于第一个极小值则将认为该值是一个极小值以该点为对称中心lmax fliplr(indmax(1:min(end,nbsym)));lmin [fliplr(indmin(1:min(end,nbsym-1))),1];lsym 1;endelseif x(1) x(indmax(1))% 以第一个极小值为对称中心lmax fliplr(indmax(1:min(end,nbsym)));lmin fliplr(indmin(2:min(end,nbsym1)));lsym indmin(1);else%如果第一个采样值大于第一个极大值则将认为该值是一个极大值以该点为对称中心lmax [fliplr(indmax(1:min(end,nbsym-1))),1];lmin fliplr(indmin(1:min(end,nbsym)));lsym 1;endend%序列末尾情况与序列开头类似if indmax(end) indmin(end)if x(end) x(indmax(end))rmax fliplr(indmax(max(end-nbsym1,1):end));rmin fliplr(indmin(max(end-nbsym,1):end-1));rsym indmin(end);elsermax [lx,fliplr(indmax(max(end-nbsym2,1):end))];rmin fliplr(indmin(max(end-nbsym1,1):end));rsym lx;endelseif x(end) x(indmin(end))rmax fliplr(indmax(max(end-nbsym,1):end-1));rmin fliplr(indmin(max(end-nbsym1,1):end));rsym indmax(end);elsermax fliplr(indmax(max(end-nbsym1,1):end));rmin [lx,fliplr(indmin(max(end-nbsym2,1):end))];rsym lx;endend%将序列根据对称中心镜像到两边tlmin 2*t(lsym)-t(lmin);tlmax 2*t(lsym)-t(lmax);trmin 2*t(rsym)-t(rmin);trmax 2*t(rsym)-t(rmax);% in case symmetrized parts do not extend enough%如果对称的部分没有足够的极值点if tlmin(1) t(1) | tlmax(1) t(1)% 对折后的序列没有超出原序列的范围if lsym indmax(1)lmax fliplr(indmax(1:min(end,nbsym)));elselmin fliplr(indmin(1:min(end,nbsym)));endif lsym 1%这种情况不应该出现程序直接中止error(bug)endlsym 1; %直接关于第一采样点取镜像tlmin 2*t(lsym)-t(lmin);tlmax 2*t(lsym)-t(lmax);end %序列末尾情况与序列开头类似if trmin(end) t(lx) | trmax(end) t(lx)if rsym indmax(end)rmax fliplr(indmax(max(end-nbsym1,1):end));elsermin fliplr(indmin(max(end-nbsym1,1):end));endif rsym lxerror(bug)endrsym lx;trmin 2*t(rsym)-t(rmin);trmax 2*t(rsym)-t(rmax);end% 延拓点上的取值xlmax x(lmax);xlmin x(lmin);xrmax x(rmax);xrmin x(rmin);% 完成延拓tmin [tlmin t(indmin) trmin];tmax [tlmax t(indmax) trmax];xmin [xlmin x(indmin) xrmin];xmax [xlmax x(indmax) xrmax];end%---------------------------------------------------------------------------------------------------% 极值点和过零点位置提取function [indmin, indmax, indzer] extr(x,t);%extracts the indices corresponding to extremaif(nargin1)t1:length(x);endm length(x);if nargout 2x1x(1:m-1);x2x(2:m);indzer find(x1.*x20);if any(x 0)iz find( x0 );indz [];if any(diff(iz)1)zer x 0;dz diff([0 zer 0]);debz find(dz 1);finz find(dz -1)-1;indz round((debzfinz)/2);elseindz iz;endindzer sort([indzerindz]);endendd diff(x);n length(d);d1 d(1:n-1);d2 d(2:n);indmin find(d1.*d20 d10)1;indmax find(d1.*d20 d10)1;% when two or more consecutive points have the same value weconsider only one extremum in the middle of the constant area% 当连续多个采样值相同时把最中间的一个值作为极值点处理方式与连0类似if any(d0)imax [];imin [];bad (d0);dd diff([0 bad 0]);debs find(dd 1);fins find(dd -1);if debs(1) 1iflength(debs) 1debs debs(2:end);fins fins(2:end);elsedebs [];fins [];endendif length(debs) 0if fins(end) mif length(debs) 1debs debs(1:(end-1));fins fins(1:(end-1));elsedebs [];fins [];end endendlc length(debs);if lc 0for k 1:lcif d(debs(k)-1) 0if d(fins(k)) 0imax [imax round((fins(k)debs(k))/2)];endelseif d(fins(k)) 0imin [imin round((fins(k)debs(k))/2)];endendendendif length(imax) 0indmax sort([indmax imax]);endif length(imin) 0indmin sort([indmin imin]);endend endend%%%%%%%%%%%%%%%%%%pos%%%%%%%%%%%%function posspos(y)%一个输出%功能找序列极值点位置坐标%y:输入序列%pos:极值点的序列位置坐标[indmin,indmax]position(y);minmaxcat(2,indmin,indmax);posssort(minmax);end%%%%%%%%%%%%%%position%%%%%%%%返还极大值和极小值位置下标%两个输出function [indmin,indmax]position(y)m length(y);d diff(y);n length(d);d1 d(1:n-1);d2 d(2:n);indmin find(d1.*d20 d10)1;indmax find(d1.*d20 d10)1;if any(d0)imax [];imin [];bad (d0);dd diff([0 bad 0]);debs find(dd 1);fins find(dd -1);if debs(1) 1iflength(debs) 1debs debs(2:end);fins fins(2:end);elsedebs [];fins [];endendif length(debs) 0if fins(end) mif length(debs) 1debs debs(1:(end-1));fins fins(1:(end-1));elsedebs [];fins [];end endendlc length(debs);if lc 0for k 1:lcif d(debs(k)-1) 0if d(fins(k)) 0imax [imax round((fins(k)debs(k))/2)];endelseif d(fins(k)) 0imin [imin round((fins(k)debs(k))/2)];endendendendif length(imax) 0indmax sort([indmax imax]);endif length(imin) 0indmin sort([indmin imin]);endendend%说明每个程序要单独保存成m文件放在同一个文件夹下调用使用实例x(t)(2cos(90*t).*cos(500*t1800.*t.*t));fs5000;t0:1/fs:0.341;yx(t);subplot(5,1,1);plot(t,y);xlabel(原始信号);[pf,a,si]lmd1(y);subplot(5,1,2);plot(t,pf(1,:));xlabel(PF1);subplot(5,1,3);plot(t,pf(2,:));xlabel(PF2);subplot(5,1,4);plot(t,pf(3,:));xlabel(PF3);subplot(5,1,5);plot(t,pf(4,:));xlabel(残量信号);从上图可以看出成功将第一个分量和第二个分量分离出来但是残余分量存在很大的问题这是因为分解过程中的骑行波没有去处导致的至于怎样完善LMD局域均值分解算法目前个人没有时间研究了希望得到指点。
