MATLAB EM2D 算法

function e2md = e2md( input_args )%——————本函数输入一灰度图象，输出该图象的2维EMD结果imf1 imf2 im3%——————二维插值use Delaunay triangulation and then cubic interpolation%——————input_args为图象路径[X,map]=imread(input_args);%X为原始图象 XMAX为极大值图象，非极大值处值为0 XMIN 为极小值图象,非极小值处值为255C=double(X);X=double(X);x_widthtrue=size(X,1);%-----复制图象，延拓100行（列）%---上下各做延拓for i=1:100    X=[X(2*i+1,:);X];endfor i=1:100    X=[X;X(x_widthtrue+101-i,:)];end%---上下延拓完毕，将做左右延拓for i=1:100    X=[X(:,2*i+1),X];endfor i=1:100    X=[X,X(:,x_widthtrue+101-i)];end%---上下延拓完毕%-----复制图象I=X;      %I描述原图象的双精度形式orig=X;orig=double(orig);      %orig描述原图象的双精度形式,用来在筛分结束的时候表示残差%————————————————————————————————x_width=size(X,1);       %x_width表示延拓之后图象的大小(长=宽)%***************************************************************for imfcount=1:3%for#01             % NL=NL_num+imfcount^2;              X=orig;              count=0;              SD=100;%SD初始化while(SD>0.3)%while#01%————————————————————————搜索极值:1——————%————————建立极值图象 极大值图象XMAX 极小值图象XMIN—————XMAX=X;XMAX=X;%————————End：建立极值图象 极大值图象XMAX 极小值图象XMIN———count=count+1;xl_max=[];          %极大值3行matrix，第一二三行分别为行数列数和灰度值xl_min=[];          %极小值3行matrix，第一二三行分别为行数列数和灰度值          %这样有利于构造包络for i=2:x_width-1    for j=2:x_width-1        if(X(i,j)>X(i-1,j))&(X(i,j)>X(i+1,j))&(X(i,j)>X(i-1,j+1))&(X(i,j)>X(i+1,j+1))&(X(i,j)>X(i,j+1))&(X(i,j)>X(i-1,j-1))&(X(i,j)>X(i+1,j-1))&(X(i,j)>X(i,j-1))            XMAX(i,j)=X(i,j);        else            XMAX(i,j)=-500;        end        if(X(i,j)<X(i-1,j))&(X(i,j)<X(i+1,j))&(X(i,j)< X(i-1,j+1))&(X(i,j)<X(i+1,j+1))&(X(i,j)<X(i,j+1))&(X(i,j)<X(i-1,j-1))&(X(i,j)<X(i+1,j-1))&(X(i,j)<X(i,j-1))            XMIN(i,j)=X(i,j);        else            XMIN(i,j)=500;        end    endendfor i=2:x_width-1    for j=2:x_width-1        if(XMAX(i,j)~=-500)            xl_max=[xl_max,[i;j;X(i,j)]];        end        if(XMIN(i,j)~=500)            xl_min=[xl_min,[i;j;X(i,j)]];        end    endend%——————————————————————End:搜索极值:2——————%————————————————End:用NL的值来进行极值的筛选————%——————————————————————构造上下包络——————x=xl_max(1,:);y=xl_max(2,:);z=xl_max(3,:);ti=1:1:x_width;[xi yi]=meshgrid(ti,ti);zi_max=griddata(y,x,z,xi,yi,'cubic');            %对极大值点进行插值            %x,y为极值点的坐标，z为该处灰度的值，            %zi显示插值结果,是一个width*width 的matrixx=xl_min(1,:);y=xl_min(2,:);z=xl_min(3,:);ti=1:1:x_width;[xi yi]=meshgrid(ti,ti);zi_min=griddata(y,x,z,xi,yi,'cubic');            %对极大值点进行插值            %x,y为极值点的坐标，z为该处灰度的值，            %zi显示插值结果,是一个width*width 的matrix%————————————————————End:构造上下包络——————zi=(zi_max+zi_min)/2;            %求得均值包络SD=max(abs(zi(101:x_widthtrue+100,101:x_widthtrue+100)))/max(abs(orig(101:x_widthtrue+100,101:x_widthtrue+100)));X=X-zi;     %将图象减去均值包络end%while#01减均值的过程imf=X;if(imfcount==1)    imf1=imf(101:x_widthtrue+100,101:x_widthtrue+100);    else if (imfcount==2)    imf2=imf(101:x_widthtrue+100,101:x_widthtrue+100);    else if(imfcount==3)        imf3=imf(101:x_widthtrue+100,101:x_widthtrue+100);        else if(imfcount==4)                imf4=imf(101:x_widthtrue+100,101:x_widthtrue+100);            else                imf5=imf(101:x_widthtrue+100,101:x_widthtrue+100);            end        end    endendorig=orig-imf;end%for#01得到IMF个数%————————————————以上得到IMF1和IMF2%********************************************************************%★★★★★在这里设置breakpoint可以得到imf1、imf2、orig(残差)★★★★★★function isitn = isitn( input_args2,isitn_n )%判断input_args2中的值是否全为isitn_n,是则返回0，否则返回1isitn=0;size_input_args2=size(input_args2,1);for i_isitn=1:size_input_args2    for j_isitn=1:size_input_args2        if (input_args2(i_isitn,j_isitn)~=isitn_n)&(i_isitn~=(floor(size_input_args2/2)+1))&(j_isitn~=(floor(size_input_args2/2)+1))            isitn=1;        end    endendfunction puguji=puguji(fxxm,m)%p为阶数,初始化为原始数据列的列数%m表示向后预测的数据的个数%fxxm为原始信号p=size(fxxm,2);for j=1:m%m表示向后预测m个    a=lpc(fxxm,p+j-1);    fxxm(p+j)=0;    for i=1:(p+j-1)        fxxm(p+j)=fxxm(p+j)-fxxm(p-i+j)*a(i+1);    endendpuguji=fxxm;function puguji2=puguji2(fxxm,m)%p为阶数,初始化为原始数据列的列数%m表示向后预测的数据的个数%fxxm为原始信号p=size(fxxm,2);fxxm=fxxm(end:-1:1);for j=1:m%m表示向后预测m个    a=lpc(fxxm,p+j-1);    fxxm(p+j)=0;    for i=1:(p+j-1)        fxxm(p+j)=fxxm(p+j)-fxxm(p-i+j)*a(i+1);    endendfxxm=fxxm(end:-1:1);puguji2=fxxm;