%%**************************************************************************************************************************************************************
%   本程序為函數(shù)功能，目的是應(yīng)用Hilbert-Huang變換尋找信號的各階本征模態(tài)和殘余分量
%      
%   中國北方發(fā)動機研究所試驗測試中心  蘇志霄   2006年10月
%**************************************************************************************************************************************************************
%
%function [h,r,t]=hilbertmotai(t0,data,n,err);   
%
%t0為輸入信號的時間序列，，data為輸入信號序列，n為所要求的經(jīng)驗?zāi)B(tài)數(shù)量，err為迭代時的前后兩次均值的絕對最大值(迭代誤差)，
%返回值h為經(jīng)驗?zāi)B(tài)矩陣(第i行為第i個模態(tài))，r為最終殘余分量，t為返回的時間序列(由于Hilbert-huang變換受到邊界效應(yīng)的影響，
%雖然采取了一定措施，實際計算發(fā)現(xiàn)其仍然難以克服，因此該函數(shù)中采取了一個方法，即僅提取中間2/3的時間序列和模態(tài)序列).
%
%**************************************************************************************************************************************************************
%jjj=1000;                                                                         %賦最終殘余分量極值點數(shù)目初始值
j=0;                                                                              %提取模態(tài)階數(shù)的初始值
%panju=2*err;
rrtemp=acc(1,:);                                           
hhtemp=rrtemp;                                                                    %本征模態(tài)的初始賦值
len=max(size(t));
sd=1000;
while (j<5);
    rr=rrtemp;
    hhtemp=rrtemp;
    while (sd>0.3);
        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)=t(k);
                maxpoint(jj)=hh(k);
            end;
            if (hh(k-1)>=hh(k))&(hh(k+1)>=hh(k));                                  %計算極小值點序列
                kk=kk+1;
                tmin(kk)=t(k);
                minpoint(kk)=hh(k);
            end;
        end;
        lenmax=max(size(maxpoint));                                                %賦極大極小端點極值，以降低邊界效應(yīng)影響
        lenmin=max(size(minpoint));
        tmax(lenmax+1)=max(t);
        tmin(lenmin+1)=max(t);
        maxpoint(1)=maxpoint(2);
        maxpoint(lenmax+1)=maxpoint(lenmax);
        minpoint(1)=minpoint(2);
        minpoint(lenmin+1)=minpoint(lenmin);
        intermax=spline(tmax,maxpoint,t);
        intermin=spline(tmin,minpoint,t);
 %       pause;
 for lll=1:len;
     if intermax(lll)<hh(lll);
         intermax(lll)=hh(lll);
     end;
     if intermin(lll)>hh(lll);
         intermin(lll)=hh(lll);
     end;
 end;
 pause;
        average=(intermax+intermin)/2;                                             %計算上下包絡(luò)線的平均值
        hhtemp=hh-average';                                                        %計算本征模態(tài)前后兩次迭代的差值
        sd=sum((average').^2./(hhtemp+eps).^2),                                       %計算迭代中的標準偏差
%        panju=max(abs(average)),
    end;
    j=j+1;
    h(j,:)=hhtemp;                                                                 %迭代成功找到了某階本征模態(tài)
    rrtemp=rr-h(j,:);                                                              %計算原始信號去掉已經(jīng)找到的殘余分量
end;
r=rrtemp; 
%len=round(len/6);                                                                  %提取中間2/3的時間序列和模態(tài)序列
%t=t(1:4*len+1);
%h=h(:,len:5*len);