myemd.m

新的信号处理方法－希尔伯特黄变换的EMD分解MATLAB程序 把程序看懂 对希尔伯特黄变换的原理就理解了

% 文件名为：myemd.m 
% 本算法用来对局域波进行经验模式分解
% 输入有:  - x : 该变量定义为一个行向量，即待处理信号
%          -t : 该变量为可选项，定义为一个行向量，即采样时间，默认为1到待处理信号长度值的所有整数值，
%          - stop: 该变量为可选项，定义为一个3元素的行向量，三个元素用来控制结束分解，分别为：threshold, threshold2 and tolerance 
%                              默认值分别为：0.05,0.5,0.05
%          - tst: 该变量为可选项，如果取值为1，则每迭代分解一次暂停，等待敲入回车键；如果取值为2，则不暂停
%
% 输出为: - imf : 为一个矩阵，内蕴模式函数，最后一行为残余分量
%          - ort : 正交性的索引值
%          - nbits : 为一个向量，记录每分解得到一个内蕴模式函数需要的迭代次数
%
% 调用函数:   - extr ()：统计输入序列中极值点和过0点数目
%          - io : 计算正交性的索引值

function [imf,ort,nbits] = myemd(x,t,stop,tst);

    % 设定默认的分解结束条件
    defstop = [0.05,0.5,0.05];
    
    %如果本函数被调用时，只有一个实际参数，其它缺省，则执行如下赋值运算
    if(nargin==1)
        t = 1:length(x);
        stop = defstop;
        tst = 0;
    end

    %如果本函数被调用时，只有两个实际参数，stop和tst缺省，则执行如下赋值运算
    if(nargin==2)
        stop = defstop;
        tst = 0;
    end

    %如果本函数被调用时，只有三个实际参数，tst缺省，则执行如下赋值运算
    if (nargin==3)
        tst=0;
    end

    % returns the sizes of each dimension of array x in a vector d with ndims(X) elements.
    S = size(x);
    %如果输入信号序列行数和列数都大于1，或者维数大于2，则输出出错信息，因为x必须是行向量或者列向量
    if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
        error('x must have only one row or one column')
    end

    %如果输入信号序列是列向量，则转置为行向量
    if S(1) > 1
        x = x';
    end

    % returns the sizes of each dimension of array t in a vector d with ndims(t) elements.
    S = size(t);
    %如果输入信号序列的时间间隔的行和列都大于1，或者维数大于2，则输出出错信息，因为t必须是行向量或者列向量
    if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
        error('t must have only one row or one column')
    end

    %如果输入信号序列的时间间隔向量是列向量，则转置为行向量
    if S(1) > 1
        t = t';
    end

    %如果输入信号序列的时间间隔和输入信号序列长度不一致，则输出出错信息，因为这两者必须长度一致
    if (length(t)~=length(x))
        error('x and t must have the same length')
    end

    %returns the sizes of each dimension of array stop in a vector d with ndims(stop)elements.
    S = size(stop);
    %stop中的三个元素必须是行向量或者列向量，并且其元素个数不能多于3，否则出错。(S(1) > 1) & (S(2) > 1)) |(length(S)>2)是行与列向量的判定条件
    if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 3) | (S(2) > 3) | (length(S) > 2)
        error('stop must have only one row or one column of max three elements')
    end

    %如果stop是列向量，则转置为行向量
    if S(1) > 1
        stop = stop';
        S = size(stop);
    end

    %如果行向量元素小于3，则stop向量第三个元素采用默认值
    if S(2) < 3
        stop(3)=defstop(3);
    end

    %如果行向量元素也小于2，则stop向量第二个元素也采用默认值
    if S(2) < 2
        stop(2)=defstop(2);
    end

    %把控制分解终止条件的值，分别赋值给相应变量，用来控制分解终止
    sd = stop(1);
    sd2 = stop(2);
    tol = stop(3);

    %如果tst为非0，则调用figure函数创建一个图形对象，如果tst取值为1，则每迭代分解一次暂停，等待敲入回车键；如果取值为2，则不暂停
    if tst
        figure   
    end

    % 初始化最大迭代次数
    MAXITERATIONS=2000;

    % 设定用来插值的最大均衡点，当取值为2时，为三次样条插值
    NBSYM = 2;

    %调用函数length，求得x向量的长度
    lx = length(x);

    %定义一个向量std，其元素初始化为全0
    sdt(lx) = 0;
    %sdt的值全为控制条件sd,即对每个输入序列每一点，其第一个控制分解终止条件都为值sd
    sdt = sdt+sd;
    %定义一个向量st2d，其元素初始化为全0
    sd2t(lx) = 0;
    %sd2t的值全为控制条件sd2,即对输入序列每一点，其第一个控制分解终止条件都为值sd2
    sd2t = sd2t+sd2;

    % 当前分解时，信号序列的极值点数，和过0点数，初始化为lx
    ner = lx;
    nzr = lx;

    %把输入信号赋值给r，原输入信号x保存起来，r为初始残余分量;把内蕴模态函数imf初始化为空；k用来记录分解的次数，初始化为1
    r = x;
    imf = [];
    k = 1;

    % nbit用来统计本次分解所用的迭代次数，初始化为0
    nbit=0;

    % NbIt用来统计总的分解所用的迭代次数，初始化为0
    NbIt=0;

    %循环条件：如果极值点数目大于2，则执行循环体
    while ner > 2

        % 把残余分量赋值给m,m为当前要分解的模式，r为残余分量
        m = r;

        % 保存分解前的模式到变量mp中
        mp = m;

        %第一个控制条件+1，赋值给变量sx
        sx = sd+1;

        % test用来记录极值点数是否足够，以便判断循环分解是否继续下去
        test = 0;

        %调用extr函数，三个变量分别记录极小值点序号，和极大值点序号，过0点序号
        [indmin,indmax,indzer] = extr(m);
        %极小值点的数目
        lm=length(indmin);
        %极大值点的数目
        lM=length(indmax);
        %极值点的数目
        nem=lm + lM;
        %过0点的数目
        nzm=length(indzer);

        %j初始化为1，j在后面没有使用，可以删除这一行代码
        j=1;

        %mean()函数：求向量的均值
        %循环分解条件：test等于 0，且nbit<MAXITERATIONS,且分解控制条件满足......
        while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS

            if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)
                disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])
                disp(['stop parameter mean value : ',num2str(S)])
            end

            % 处理边界点，对边界点进行延拓,边界点的处理相当不麻烦
            %处理输入序列头部边界点
            %如果第一个极大值出现在第一个极小值点前面号，则：
            if indmax(1) < indmin(1)
                %如果第一个极小值点的值小于输入信号序列的第一个值
                if m(1) > m(indmin(1))
                    %保存：indmax向量中的，从第2个极大值开始的NBSYM个极大值序号(并且反转过来)，到lmax向量中，用来与第一个极大值点进行一次插值
                    lmax = fliplr(indmax(2:min(end,NBSYM+1)));
                    %保存：indmin向量中的，从第1个极小值开始的NBSYM个极小值序号(并且反转过来)，到lmin向量中，用来与第一个极小值点进行一次插值
                    lmin = fliplr(indmin(1:min(end,NBSYM)));
                    %保存：第一个极大值点的序号到lsym中
                    lsym = indmax(1);
                %如果第一个极小值点的值不小于输入信号序列的第一个值
                else
                    %保存：indmax向量中的，从第1个极大值开始的NBSYM个极大值序号(并且反转过来)，到lmax向量中，用来与第一个极大值点进行一次插值
                    lmax = fliplr(indmax(1:min(end,NBSYM)));
                    %保存：indmin向量中的，从第2个极小值开始的NBSYM-1个极小值序号(并且反转过来)和序号1，到lmin向量中，用来与第一个极小值点进行一次插值
                    lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
                    %保存：信号序列的第一个序号到lsym中
                    lsym = 1;
                end
            %如果第一个极小值出现在第一个极大值点前面，则：
            else
                %如果第一个极大值点的值大于输入信号序列的第一个值
                if m(1) < m(indmax(1))
                    %保存：indmax向量中的，从第1个极大值开始的NBSYM个极大值序号(并且反转过来)，到lmax向量中，用来与第一个极大值点进行一次插值
                    lmax = fliplr(indmax(1:min(end,NBSYM)));
                    %保存：indmin向量中的，从第2个极小值开始的NBSYM个极小值序号(并且反转过来)，到lmin向量中，用来与第一个极小值点进行一次插值
                    lmin = fliplr(indmin(2:min(end,NBSYM+1)));
                    %保存：第一个极小值点的序号到lsym中
                    lsym = indmin(1);
                else
                    %保存：indmax向量中的，从第1个极大值开始的NBSYM个极大值序号(并且反转过来)和序号1，到lmax向量中，用来与第一个极大值点进行一次插值
                    lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
                    %保存：indmin向量中的，从第1个极小值开始的NBSYM个极小值序号(并且反转过来)，到lmin向量中，用来与第一个极小值点进行一次插值
                    lmin = fliplr(indmin(1:min(end,NBSYM)));
                    %保存：信号序列的第一个序号到lsym中
                    lsym = 1;
                end
            end
            %处理输入序列尾部边界点
            %下面的程序对尾部边界点的处理类似上面对头部边界点的处理类似：