imf_sift.m

这个是EMD分解算法,相对比较简单,十分适合初学者,我开始也是从个程序入手的,觉得十分不错,与大家共享.

function IMF=IMF_sift(Image,conn,sd)
% IMF 是从图像Image中分解得到的固有模态函数
% Image 是输入图像，conn 是检测图像极值点时所用到的像素邻域，通常取conn=4或conn=8
% sd 是筛选过程中的结束条件
% imf_num 是筛选所要得到的固有模态函数IMF的个数
Ires=Image;   %复制图像给Iresi
SD=1;          % 设置迭代条件初始值

while(SD>sd)
    %Iresi=wextend('2D','symw',Iresi,100,'b'); 
    Orig=symmetric_extension(Ires,100);            %镜像延拓图像
    [m,n]=size(Orig);
    [Upenvelop,Lowenvelop]=Envelope(Orig,conn);    %计算极值点上下包络
    Emeani=(Upenvelop+Lowenvelop)./2;              %原始信号均值包络
    Iresi=Orig-Emeani;                             %第一次去均值结果h1
    IresiTrue=Iresi(101:m-100,101:n-100);          %去延拓h1图像
    Iresisym=symmetric_extension(IresiTrue,100);   %镜像延拓h1图像
    [Upenvelop,Lowenvelop]=Envelope(Iresisym,conn);  %计算h1的上下均值
    Emeanj=(Upenvelop+Lowenvelop)./2;              %求的h1的包络均值
    Iresj=Iresi-Emeanj;                            %第二次去均值结果h2
    IresjTrue=Iresj(101:m-100,101:n-100);          %去延拓h2图像
    [mm,nn]=size(IresiTrue);
    ERRup=0;
    ERRlow=0;
    for i=1:mm
        for j=1:nn
            ERRup=ERRup+(abs(IresiTrue(i,j)-IresjTrue(i,j)))^2;  
        end
    end
    for i=1:mm
        for j=1:nn
            ERRlow=ERRlow+IresiTrue(i,j)^2;
        end
    end
    SD=ERRup/ERRlow;
    Ires=IresjTrue;
end
IMF=IresjTrue;
%figure(1);imshow(IMF,[])