HHT代码

HHT（Hilbert-Huang 变换）是一种用于非线性和非平稳信号分析的方法。

HHT 主要由两部分组成：经验模态分解（Empirical Mode Decomposition，EMD）和 Hilbert 变换。EMD 可以将非线性和非平稳信号分解为一系列本征模态函数（Intrinsic Mode Functions，IMFs），每个 IMF 代表信号中的一个局部特征。然后，对每个 IMF 应用 Hilbert 变换，得到该模态函数的即时频率。最后，将所有 IMFs 的即时频率叠加，得到信号的时频分布图。

下面是 HHT 的 Python 实现代码：

import numpy as np
import scipy.signal as signal

def emd(s):
    """Empirical Mode Decomposition"""
    # Calculate the mean value of the envelope
    mean_env = np.mean(signal.hilbert(s).real)
    # Initialize the residual as the original signal
    r = s.copy()
    imfs = []
    while True:
        # Calculate the upper and lower envelopes
        u = r.copy()
        d = r.copy()
        for i in range(10):
            # Calculate the mean value of the envelope
            mean_env = (signal.hilbert(u).real + signal.hilbert(d).real) / 2
            # Extract the local maxima and minima
            u = u - mean_env
            d = d - mean_env
        # Calculate the IMF
        imf = r - (u + d) / 2
        # Check if the IMF is valid
        if np.sum(np.abs(imf)) < np.sum(np.abs(r)) * 0.2:
            imfs.append(imf)
            # Update the residual
            r = r - imf
            if np.sum(np.abs(r)) < 1e-10:
                break
    return imfs

def hht(s, fs):
    """Hilbert-Huang Transform"""
    # EMD
    imfs = emd(s)
    # Hilbert Transform
    freqs = []
    for imf in imfs:
        analytic_signal = signal.hilbert(imf)
        instantaneous_phase = np.unwrap(np.angle(analytic_signal))
        instantaneous_frequency = np.diff(instantaneous_phase) / (2 * np.pi) * fs
        freqs.append(instantaneous_frequency)
    return freqs

其中，`emd(s)` 函数实现 EMD，输入信号 `s`，输出 IMFs。`hht(s, fs)` 函数实现 HHT，输入信号 `s` 和采样率 `fs`，输出 IMFs 的即时频率。