hilbertmotai.m

自己编写的比较全面的故障诊断matlab函数程序

%%**************************************************************************************************************************************************************
%   本程序为函数功能，目的是应用Hilbert-Huang变换寻找信号的各阶本征模态和残余分量
%      
%   中国北方发动机研究所试验测试中心  苏志霄   2006年10月
%**************************************************************************************************************************************************************
%
function [h,r,t]=hilbertmotai(t0,data,n,err);   
%
%t0为输入信号的时间序列，，data为输入信号序列，n为所要求的经验模态数量，err为迭代时的前后两次均值的绝对最大值(迭代误差)，
%返回值h为经验模态矩阵(第i行为第i个模态)，r为最终残余分量，t为返回的时间序列(由于Hilbert-huang变换受到边界效应的影响，
%虽然采取了一定措施，实际计算发现其仍然难以克服，因此该函数中采取了一个方法，即仅提取中间2/3的时间序列和模态序列).
%
%**************************************************************************************************************************************************************
%jjj=1000;                                                                         %赋最终残余分量极值点数目初始值
j=0;                                                                              %提取模态阶数的初始值
panju=2*err;
rrtemp=data;                                           
hhtemp=rrtemp;                                                                    %本征模态的初始赋值
len=max(size(t0));
%sd=1000;
while (j<n);
    rr=rrtemp;
    hhtemp=rrtemp;
    while (panju>err);
        hh=hhtemp;
        jj=1;
        kk=1;
        tmax(1)=0;
        maxpoint(1)=0;
        tmin(1)=0;
        minpoint(1)=0;
        for k=2:len-1;
            if (hh(k-1)<=hh(k))&(hh(k+1)<=hh(k));                                  %计算极大值点序列
                jj=jj+1;
                tmax(jj)=t0(k);
                maxpoint(jj)=hh(k);
            end;
            if (hh(k-1)>=hh(k))&(hh(k+1)>=hh(k));                                  %计算极小值点序列
                kk=kk+1;
                tmin(kk)=t0(k);
                minpoint(kk)=hh(k);
            end;
        end;
        lenmax=max(size(maxpoint));                                                %赋极大极小端点极值，以降低边界效应影响
        lenmin=max(size(minpoint));
        tmax(lenmax+1)=max(t0);
        tmin(lenmin+1)=max(t0);
        maxpoint(1)=maxpoint(2);
        maxpoint(lenmax+1)=maxpoint(lenmax);
        minpoint(1)=minpoint(2);
        minpoint(lenmin+1)=minpoint(lenmin);
        intermax=spline(tmax,maxpoint,t0);
        intermin=spline(tmin,minpoint,t0);
        average=(intermax+intermin)/2;                                             %计算上下包络线的平均值
        hhtemp=hh-average';                                                        %计算本征模态前后两次迭代的差值
%        sd=sum((average'./(hhtemp+eps)).^2),                                       %计算迭代中的标准偏差
        panju=max(abs(average));
    end;
    j=j+1;
    h(j,:)=hhtemp;                                                                 %迭代成功找到了某阶本征模态
    rrtemp=rr-h(j,:);                                                              %计算原始信号去掉已经找到的残余分量
end;
r=rrtemp; 
len=round(len/6);                                                                  %提取中间2/3的时间序列和模态序列
t=t0(1:4*len+1);
h=h(:,len:5*len);