%%**************************************************************************************************************************************************************
%   本程序為函數功能，目的是應用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);