HHT matlab实现

HHT（Hilbert-Huang Transform）希尔伯特黄变换，HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。

function hht=HHT（x）
    clc;
    N=length（x）;
    t=0:0.01:8;
    deta=t（2）-t（1）;  %采样间隔
    fs=1/deta;       %采样频率
%fft默认计算的信号是从0开始的
    z=x;
    c=emd（z）;
%计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性
    [mn]=size（c）;
    for i=1:m
        a=corrcoef（c（i:）z）;
        xg（i）=a（12）;
    end
        xg;
    for i=1:m-1
    %--------------------------------------------------------------------
    %计算各IMF的方差贡献率
    %定义：方差为平方的均值减去均值的平方
    %均值的平方
    %imfp2=mean（c（i:）2）.^2
    %平方的均值
    %imf2p=mean（c（i:）.^22）
    %各个IMF的方差
        mse（i）=mean（c（i:）.^22）-mean（c（i:）2）.^2;
    end;
      mmse=sum（mse）;
    for i=1:m-1
        mse（i）=mean（c（i:）.^22）-mean（c（i:）2）.^2; 
    %方差百分比，也就是方差贡献率
        mseb（i）=mse（i）/mmse*100;
    %显示各个IMF的方差和贡献率
    end;
%画出每个IMF分量及最后一个剩余分量residual的图形
        figure（1）
    for i=1:m-1
        disp（[‘imf‘int2str（i）]） ;disp（[mse（i） mseb（i）]）;   %disp函数直接将内容输出在Matlab命令窗口中
    end;
         subplot（m+111）
         plot（tz）
         set（gca‘fontname‘‘times New Roman‘）
         set（gca‘fontsize‘14.0）
         ylabel（[‘signal‘‘Amplitude‘]）

    for i=1:m-1
        subplot（m+11i+1）;
        set（gcf‘color‘‘w‘）
        plot（tc（i:）‘k‘）
        set（gca‘fontname‘‘times New Roman‘）
        set（gca‘fontsize‘14.0）
        ylabel（[‘imf‘int2str（i）]）
    end
        subplot（m+11m+1）;
        set（gcf‘color‘‘w‘）
        plot（tc（m:）‘k‘）
        set（gca‘fontname‘‘times New Roman‘）
        set（gca‘fontsize‘14.0）
        ylabel（[‘r‘int2str（m-1）]）

%画出每个IMF分量及剩余分量residual的幅频曲线
%     figure（2）
%     subplot（m+111）
%     set（gcf‘color‘‘w‘）
%     [fz]=fftfenxi（tz）;
%     plot（fz‘k‘）
%     set（gca‘fontname‘‘times New Roman‘）
%     set（gca‘fontsize‘14.0）
%     ylabel（[‘initial signal‘int2str（m-1）‘Amplitude‘]）

%     for i=1:m-1
%         subplot（m+11i+1）;
%         set（gcf‘color‘‘w‘）
%         [fz]=fftfenxi（tc（i:））;
%         plot（fz‘k‘）
%         set（gca‘fontname‘‘times New Roman‘）
%         set（gca‘fontsize‘14.0）
%         ylabel（[‘imf‘int2str（i）‘Amplitude‘]）
%     end
%         subplot（m+11m+1）;
%         set（gcf‘color‘‘w‘）
%         [fz]=fftfenxi（tc（m:））;
%         plot（fz‘k‘）
%         set（gca‘fontname‘‘times New Roman‘）
%         set（gca‘fontsize‘14.0）
%         ylabel（[‘r‘int2str（m-1）‘Amplitude‘]）

        hx=hilbert（z）;
        xr=real（hx）;xi=imag（hx）;
    %计算瞬时振幅
        sz=sqrt（xr.^2+xi.^2）;
    %计算瞬时相位
        sx=angle（hx）;
    %计算瞬时频率
        dt=diff（t）;
        dx=diff（sx）;
        sp=dx./dt;
        figure（6）
        plot（t（1:N-1）sp）
        title（‘瞬时频率‘）

    %计算HHT时频谱和边际谱
        [Afatt]=hhspectrum（c）;