envelope.m

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

function [Upenvelop,Lowenvelop]=Envelope(Image,conn)
% Image 是输入图像，conn是用来检测图像极值点时所选取的像素邻域
% 对输入图像进行双向镜像延拓，双向各延拓100的长度

% ExImage=wextend('2D','symw',Image,100,'b');    % Symmetric-padding (whole-point): boundary value symmetric replication
% ExImage=wextend('2D','symh',Image,100,'b');  % Symmetric-padding (half-point): boundary value symmetric replication
% ExImage=wextend('2D','asymh',Image,100,'b'); % Antisymmetric-padding (half-point): boundary value antisymmetric replication
% ExImage=wextend('2D','asymw',Image,100,'b'); % Antisymmetric-padding (whole-point): boundary value antisymmetric replication
% ExImage=wextend('2D','ppd',Image,100,'b');   % Periodized extension (1)
% ExImage=wextend('2D','per',Image,100,'b');   % Periodized extension (2):If the signal length is odd, wextend adds an extra sample, equal to the last value, on the right and performs extension using the 'ppd' mode. Otherwise, 'per' reduces to 'ppd'. The same kind of rule stands for images.

%计算图像Image的极大值图像Imax和极小值图像Imin
[Imax,Imin]=Imextrem(Image,conn); 

%I_max=[];    %初始化存储极大值位置坐标和像素值的矩阵,三行矩阵，第一行是横坐标值，第二行是纵坐标值，第三行是坐标所对位置的像素值
%I_min=[];    %初始化存储极小值位置坐标和像素值的矩阵,三行矩阵，第一行是横坐标值，第二行是纵坐标值，第三行是坐标所对位置的像素值

% 计算I_max矩阵
[m0,n0]=size(Imax);  %求出极大值图像的大小，m代表行数，n代表列数
[m1,n1]=size(Imin);  %求出极小值图像的大小，m代表行数，n代表列数
%for i=1:m0
%    for j=1:n0
%        if(Imax(i,j)~=0)
%            I_max=[I_max,[i;j;Imax(i,j)]];
%        end
%    end
%end
% 计算I_min矩阵
%for i=1:m1
%    for j=1:n1
%       if(Imin(i,j)~=0)
%            I_min=[I_min,[i;j;Imin(i,j)]];
%        end
%    end
%end
% 计算极大值包络的插值过程
%x=I_max(1,:);   %图像极大值图像Imax的像素值不为零处的行坐标
%y=I_max(2,:);   %图像极大值图像Imax的像素值不为零处的纵坐标
%z=I_max(3,:);   %图像极大值图像Imax的不为零处像素值

[x,y,z]=find(Imax);
d=1:1:m0;       %步长
[xi,yi]=meshgrid(d,d);
zi_max=griddata(y,x,z,xi,yi,'cubic'); %三次样条插值形成的上包络图

%zi_max=griddata(y,x,z,xi,yi,'spline'); %样条插值形成的上包络图
%zi_max=griddata(y,x,z,xi,yi,'nearest'); % 最邻近插值形成的上包络图
%zi_max=griddata(y,x,z,xi,yi,'linear');  %线性插值形成的上包络图

Upenvelop=zi_max;

%figure(10);imshow(Upenvelop,[])

% 计算极小值包络的插值过程
%xx=I_min(1,:);   %图像极大值图像Imax的像素值不为零处的行坐标
%yy=I_min(2,:);   %图像极大值图像Imax的像素值不为零处的纵坐标
%zz=I_min(3,:);   %图像极大值图像Imax的不为零处像素值

[xx,yy,zz]=find(Imin);
t=1:1:m1;       %步长
[xxi,yyi]=meshgrid(t,t);
zi_min=griddata(yy,xx,zz,xxi,yyi,'cubic'); %三次样条插值形成的下包络图

%zi_min=griddata(y,x,z,xi,yi,'spline'); %样条插值形成的下包络图
%zi_min=griddata(y,x,z,xi,yi,'nearest'); % 最邻近插值形成的下包络图
%zi_min=griddata(y,x,z,xi,yi,'linear');  %线性插值形成的下包络图

Lowenvelop=zi_min;
%figure(20);imshow(Lowenvelop,[])