希尔伯特频谱算法Hilbert-Huang spectral analysis（matlab代码）

た入场券 2023-07-02 13:25 14阅读 0赞

前段时间磕盐接触到了希尔伯特频谱，它是一种信号分解方法，1998年提出来的，主旨是把复杂信号分解为简单信号的加权和，就像傅里叶变换小波变换一样，但是他和傅里叶变换等方法的区别是他是纯粹时间域的分解，但是每个子信号却可以表示不同的频率成分，于是可以得到像小波变换那样的时频平面，但是这个方法明显比小波分解冷门的多，而且在我的实验结果里确实远远远远弱于小波分解，不过也算是自己辛苦几天看论文和写代码的成果，特此记录，如果后面有人用到这个，可以快速入门

[文献链接：The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis][The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis]

伪代码，我读论文总结出来的，主要就是经验模式分解EMD和希尔伯特变换两步  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0_size_16_color_FFFFFF_t_70]  
伪代码的latex 代码

\begin{algorithm}
    	\caption{Hilbert-Huang spectral analysis \cite{Branch1998}}
    	\label{hhsa}
    	\begin{algorithmic}[1] 
    		\Require The original signal vector $\boldsymbol x$.
    		\Ensure  The Hilbert-Huang Spectrum, i.e. an energy-time-frequency distribution of $\boldsymbol x$. 
    		\Function{EMD}{$\boldsymbol x, SegLen, ResidueThreshold, SD_T$}		
    		\State $\boldsymbol{IMF} \gets \boldsymbol 0$
    		\State $i \gets 0$
    		\State $ N \gets length(\boldsymbol{x}) / SegLen$
    		\State $residue \gets \infty$
    		\While {$residue > ResidueThreshold$}
    		\State $i \gets i+1$	
    		\State $\boldsymbol{x_i} \gets \boldsymbol{x}-\sum_i \boldsymbol{IMF}$
    		\State $SD \gets \infty$
    		\While{$SD > SD_T$}	
    		\For {$j=1 \to j=N$}	
    		\State $\boldsymbol{seg_{ij}} \gets \boldsymbol{x_i}[(1+(j-1)*SegLen) : (j*SegLen)]$
    		\State $[LocalMax_{ij}, IndMax_{ij}]\gets max(\boldsymbol{seg_{ij}})$
    		\State $IndMax_{ij}\gets IndMax_{ij}+(j-1)*SegLen$
    		\State $[LocalMin_{ij}, IndMin_{ij}]\gets min(\boldsymbol{seg_{ij}})$
    		\State $IndMin_{ij}\gets IndMin_{ij}+(j-1)*SegLen$
    		\EndFor
    		\State $\boldsymbol{UpperEnv_i} \gets spline(\boldsymbol{IndMax_i},\boldsymbol{LocalMax_i},1:length(\boldsymbol{x_i}))$
    		\State $\boldsymbol{LowerEnv_i} \gets spline(\boldsymbol{IndMin_i},\boldsymbol{LocalMin_i},1:length(\boldsymbol{x_i}))$
    		\State $\boldsymbol{LocalMeanApprox_i} \gets (\boldsymbol{UpperEnv_i} + \boldsymbol{LowerEnv_i})/2$
    		\State $\boldsymbol {x_i} \gets \boldsymbol {x_i} - \boldsymbol{LocalMeanApprox_i} $
    		
    		\State $SD\gets \sum[\frac{(\boldsymbol{xTemp}-\boldsymbol{x_i})^2}{\boldsymbol{xTemp}^2}]$
    		\EndWhile
    		\State $\boldsymbol{IMF_i} \gets \boldsymbol x_i$
    		\State $residue \gets mean(\boldsymbol x-\sum_i \boldsymbol{IMF})$
    		\EndWhile
    		\State \Return{$\boldsymbol{IMF}$}
    		\EndFunction
    		\State
    		\Function{HT}{$\boldsymbol{IMF},fs$}
    		\State $\boldsymbol{zt} \gets hilbert(\boldsymbol{IMF})$
    		\State $\boldsymbol{Energy} \gets (imag(\boldsymbol{zt}))^2+(real(\boldsymbol{zt}))^2$
    		\State $\boldsymbol{Phase} \gets arctan\frac{imag(\boldsymbol{zt})}{real(\boldsymbol{zt})+eps}$
    		\State $\boldsymbol{frequency} \gets \boldsymbol{Phase}/(2*\pi*fs)$
    		\State \Return{$\boldsymbol{Energy,frequency}$} 
    		\EndFunction
    	\end{algorithmic}
    \end{algorithm}

导言区需要包含的包

\usepackage{algorithm}
    \usepackage{algorithmicx}
    \usepackage{algpseudocode}  
    \usepackage{amsmath}  % 用于算法伪代码中的数学公式
    \renewcommand{\algorithmicrequire}{ \textbf{Input:}}     %Use Input in the format of Algorithm
    \renewcommand{\algorithmicensure}{ \textbf{Output:}}    %UseOutput in the format of Algorithm

matlab代码后面加，被新冠病毒阻挡了去学校的脚步，代码在教研室电脑上，后面补上

结果图  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0_size_16_color_FFFFFF_t_70 1]

代码：  
这个代码计算了EEG信号的HHS，不能直接用，但是大致实现可以参考，写的很烂，，，大神随意看看就好

% EMD.m
    % HHS hilbert huang spectrum  emd
    function IMF = EMD(eeg)
    IMF = zeros(size(eeg,1), 9);
    for j=1:size(eeg,1)
        ch = eeg(j,:);
        x = ch;
        
        IMF = zeros(100, length(x)); % 每条IMF都和x等长
        seg_len = 4;
        iter_num = length(x)/seg_len;
        local_max = zeros(iter_num,1);local_min = zeros(iter_num,1);
        indmax = zeros(iter_num,1);indmin = zeros(iter_num,1);
        for s=1:100
            res = x - sum(IMF); % IMF第一行是IMF0,初始化为0,res残差
            mean(res)
            if(abs(mean(res))<.017) break; end
            xs = res;
            %                  figure,plot(1:200,res(1:200)),grid on;title('res')
            
            for q=1:100
                for i=1:iter_num
                    seg = xs((1+(i-1)*seg_len):(i*seg_len)); % 对x分段
                    [local_max(i),indma] = max(seg);[local_min(i),indmi] = min(seg); % 局部极值点
                    indmax(i) = indma + (i-1)*seg_len; indmin(i) = indmi + (i-1)*seg_len;% 索引
                end
                % 插值计算上下包络
                up = spline(indmax,local_max,1:length(x));
                lo = spline(indmin,local_min,1:length(x));
                up1 = interp1(indmax,local_max,1:length(x),'pchip');
                lo1 = interp1(indmin,local_min,1:length(x),'pchip');
                %                     up2 = interp1(indmax,local_max,1:length(x),'nearest');
                %                     lo2 = interp1(indmin,local_min,1:length(x),'nearest');
                %                     up3 = interp1(indmax,local_max,1:length(x),'linear');
                %                     lo3 = interp1(indmin,local_min,1:length(x),'linear');
                %                     up4 = interp1(indmax,local_max,1:length(x),'nearest');
                %                     lo4 = interp1(indmin,local_min,1:length(x),'nearest');
                %                     for r=1:length(x)
                %                         if(up(r)<lo(r)) up(r)=lo(r); end
                %                     end
                %                     figure,
                %                     plot(1:200,lo(1:200),'b',1:200,xs(1:200),'k',1:200,up(1:200),'r')
                %                     set(gca,'XLim',[1 200]);
                
                %                     figure,
                %                     plot(1:200,lo1(1:200),'b',1:200,xs(1:200),'k',1:200,up1(1:200),'r')
                %                     figure,
                %                     plot(1:200,lo2(1:200),'b',1:200,xs(1:200),'k',1:200,up2(1:200),'r')
                %                     figure,
                %                     plot(1:200,lo3(1:200),'b',1:200,xs(1:200),'k',1:200,up3(1:200),'r')
                %                     figure,
                %                     plot(1:200,lo4(1:200),'b',1:200,xs(1:200),'k',1:200,up4(1:200),'r')
                
                %                     set(gca,'XLim',[1 200]);
                xtemp = xs;
                xs = xs - (up1+lo1)/2;
                SD = mean(((xtemp-xs).^2)./((xtemp).^2));
                if(SD<.3) break; end % 不满足则继续sifting，最终的x就是IMF % 满足则不再sifting
                
            end
            %
            IMF(s,:) = xs;
            figure,plot(1:200,xs(1:200),'k'),grid on;title('IMF')
            %                 res = x - sum(IMF);
            
            %                 set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',.5);%虚线网格
            %                 if(j<=3)
            %                 subplot(3,1,j)
            %                 plot(1:256,xs(1:256),'k',1:256,up(1:256),'r',1:256,lo(1:256),'b')
            %                 axis([1 256 -2 2])
        end
        % 去除多余的行
        for i=1:100
            if(sum(IMF(i,:))==0)
                IMF(i:end,:) = [];
                break;
            end
        end
    end

[The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis]: https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1998.0193
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/2020013014061087.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20200130142602725.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM2NjA3ODk0,size_16,color_FFFFFF,t_70