myemd.m

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

            if indmax(end) < indmin(end)
                if m(end) < m(indmax(end))
                    rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
                    rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
                    rsym = indmin(end);
                else
                    rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
                    rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
                    rsym = lx;
                end
            else
                if m(end) > m(indmin(end))
                    rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
                    rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
                    rsym = indmax(end);
                else
                    rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
                    rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
                    rsym = lx;
                end
            end

            %给出各延拓的边界点的采样时间
            tlmin = 2*t(lsym)-t(lmin);
            tlmax = 2*t(lsym)-t(lmax);
            trmin = 2*t(rsym)-t(rmin);
            trmax = 2*t(rsym)-t(rmax);

            % 如果头部延拓的边界点在原采样时间范围内，则作如下处理
            if tlmin(1) > t(1) | tlmax(1) > t(1)
                %如果lsym值为第一个极大值的序号，则修改lmax向量为：
                if lsym == indmax(1)
                    lmax = fliplr(indmax(1:min(end,NBSYM)));
                %如果lsym值为第一个极小值的序号，则修改lmin向量为：
                else
                    lmin = fliplr(indmin(1:min(end,NBSYM)));
                end
                %如果lsym为1，则输出出错信息，原程序存在bug
                if lsym == 1
                    error('bug')
                end
                %lsym为原信号的序号1
                lsym = 1;
                %修改tlmin和tlmax，使其处入采样时间范围之外
                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-NBSYM+1,1):end));
                else
                    rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
                end
                if rsym == lx
                    error('bug')
                end
                rsym = lx;
                trmin = 2*t(rsym)-t(rmin);
                trmax = 2*t(rsym)-t(rmax);
             end 

            %得到延拓的边界点对应的信号值
            mlmax =m(lmax); 
            mlmin =m(lmin);
            mrmax =m(rmax); 
            mrmin =m(rmin);
            % 由极值点及延拓的头部、尾部边界极值点通过三次样条插值分别得到上下包络线，
            envmax = interp1([tlmax t(indmax) trmax],[mlmax m(indmax) mrmax],t,'spline');
            envmin = interp1([tlmin t(indmin) trmin],[mlmin m(indmin) mrmin],t,'spline');
            %求包络均值
            envmoy = (envmax + envmin)/2;
            %求内蕴模式函数
            m = m - envmoy;

            %调用extr函数，记录分离出的内蕴模式函数的，极小值点序号，和极大值点序号，过0点序号到向量中
            [indmin,indmax,indzer] = extr(m);
            %极小值点的数目
            lm=length(indmin);
            %极大值点的数目
            lM=length(indmax);
            %极值点的数目
            nem = lm + lM;
            %过0点的数目
            nzm = length(indzer);
            
            % 求包络均值与包络差值的比值的2倍：均值越小，差值越大，该比值越小
            sx=2*(abs(envmoy))./(abs(envmax-envmin));
            %求各点比值的均值
            s = mean(sx);

            % display       
            %如果调用时tst为非0，则执行下述代码进行显示
            if tst
                %将显示区域分为一列，四行的小区域，且后一个subplot之前的显示将在第1区域里绘制
                subplot(4,1,1)
                %绘制分解前的原序列
                plot(t,mp);
                %The hold function determines whether new graphics objects are added to the graph or replace objects in the graph.
                hold on;
                %绘制分解前的原序列的上、下包络及均值
                %plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
                %定义本图形小区域的标题：内蕴模式函数第k个，迭代次数nbit，此为迭代前的序列
                title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
                %gca里将得到相应坐标轴句柄
                set(gca,'XTick',[])
                %The hold function determines whether new graphics objects are added to the graph or replace objects in the graph.
                hold  off

                %下一个subplot之前的显示将在第2区域里绘制
                subplot(4,1,2)
                %绘制分解后的控制条件3
                plot(t,sx)
                %The hold function determines whether new graphics objects are added to the graph or replace objects in the graph.
                hold on
                %绘制分解后的控制条件1
                plot(t,sdt,'--r')
                %绘制分解后的控制条件2
                plot(t,sd2t,':k')
                %定义本图形小区域的标题：进行分解后的参数
                title('stop parameter')
                %gca里将得到相应坐标轴句柄
                set(gca,'XTick',[])
                %The hold function determines whether new graphics objects are added to the graph or replace objects in the graph.
                hold off

                %下一个subplot之前的显示将在第3区域里绘制
                subplot(4,1,3)
                %绘制本次分解得到的内蕴模式函数
                plot(t,m)
                %定义本小区域的标题：内蕴模式函数第k个，迭代次数nbit，此为分解后的序列
                title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
                %gca里将得到相应坐标轴句柄
                set(gca,'XTick',[])

                %下一个subplot之前的显示将在第4区域里绘制
                subplot(4,1,4);
                %绘制本次分解后得到的残余分量
                plot(t,r-m)
                %定义本小区域的标题：residue，残余分量
                title('residue');
                %显示本次分解后的s值：stop parameter mean value :s
                disp(['stop parameter mean value : ',num2str(s)])
                %如果本函数调用时tst参数值为2，则暂停0.01秒继续;否则暂停，直到敲回车键继续
                if tst == 2
                    pause(0.01)
                else
                    pause
                end     
            end
            
            %如果极值点的个数小于3，则循环终止条件满足，分解将结束
            if nem < 3
                test = 1;
            end
            
            % 记录当前得到的内蕴模式函数，作为下次分解时的分解前模式（用于显示）
            mp = m;
            %本次循环的迭代次数计数，总的循环迭代次数计数
            nbit=nbit+1;
            NbIt=NbIt+1;

            %如果本次循环迭代次数超过最大迭代次数MAXITERATIONS，则显示出错信息：强制终止时迭代次数为k次，终止时的第三个控制
            %条件为：s
            if(nbit==(MAXITERATIONS-1))
                warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
            end
        end
        
        %记录当前内蕴模式函数
        imf(k,:) = m;
        %如果tst为真，则显示提示信息：记录第k个内蕴模式函数
        if tst
            disp(['mode ',int2str(k),' enregistre'])
        end
        %记录当次循环所用的迭代次数
        nbits(k) = nbit;
        %迭代次数计数
        k = k+1;
        %保存当次循环后的残余分量
        r = r - m;
        %调用extr函数，记录极小值点序号，和极大值点序号，过0点序号到向量中
        [indmin,indmax,indzer] = extr(r);
        %极值点的数目
        ner = length(indmin) + length(indmax);
        %过0点的数目
        nzr = length(indzer);
        %迭代次数初始化为1
        nbit=1;

        %如果残余分量中的最大值减去最小值小于原始信号的最大值减去最小值的10的10次方分之一，则显示出错信息：
        %1：如果极值点数大于2，则显示警告信息：强制终止，太小的频率
        %2：如果极值点数小于3，则显示出错信息：强制终止，太小的频率
        if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
            if ner > 2
                warning('forced stop of EMD : too small amplitude')
            else
                disp('forced stop of EMD : too small amplitude')
            end
            break
        end  
    end

    %记录当前分解的残余分量
    imf(k,:) = r;
    %调用函数io(),计算当前分解是否正交
    ort = io(x,imf);

    %如果tst为非0值，则调用close 来Delete specified figure
    if tst
        close
    end