代码收藏家技术教程 2024-11-01

希尔伯特黄变换（Python-emd-fft-hilbert）

一：希尔伯特黄变换概述

希尔伯特-黄变换（Hilbert-Huang Transform，HHT）是一种处理非线性非平稳信号的时频分析方法，原始信号首先进行EMD经验模态分解得到若干个IMF本征模态函数后，利用希尔伯特Hilbert变换对信号进行进一步分析处理。在上一篇文章中已经对EMD经验模态分解进行了非常详细的概述，这里仅对希尔伯特Hilbert变换进行进一步解析。

针对一个实数信号x(t)，对其进行希尔伯特变换（Hilbert Transform），得到的变换信号记为H(x(t))，由此我们可以得到该原始信号的解析信号z(t)：

$z(t) = x(t) + jH\left(x(t) \right)$

通过将原始信号转变为解析信号，将为我们后续的信号处理带来巨大的便利，这里我们借助如下图像来进行分析：

通过hilbert变换，我们将仅包含实部的实数信号变换为解析信号（复数信号），相当于将一个一维信号扩展到二维复平面。二维的解析信号振幅与原信号相同，相位相差90度。通过计算解析信号的模和幅角，可以求解原信号的瞬时振幅（包络函数）a(t)、瞬时相位φ(t)、瞬时频率f(t)等等。

$a(t) = \sqrt{x^{2}(t) + H^{2}\left(x(t) \right)}$

$\varphi (t) =\arctan \frac{H\left(x(t) \right)}{x(t)}$

$f(t)=\frac{1}{2\pi }\frac{\mathrm{d[\varphi (t)]} }{\mathrm{d} x}$

二：Python-Hilbert

这里我们同样利用从kaggle上下载的一段音频数据集，首先对其进行emd经验模态分解：

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD
import pandas as pd
from scipy.signal import hilbert

file_path = "D:/kaggle/data.txt"
df = pd.DataFrame()

with open(file_path, "r", encoding="ISO-8859-1") as file:
    for line_number, line in enumerate(file, 0):  # 调用enumerate函数返回索引和值的元组
        line = line.strip()  # str.strip()去除字符串两端空白字符(空格、制表符、换行符等)
        if len(line) != 0:
            if line_number < 1000:
                print(f"Line{line_number}：{line}")
                df.loc[line_number, "time"] = line.split(";")[1]  # 时间
                df.loc[line_number, "sample"] = float(line.split(";")[7].split("=")[1])

# 原始信号数据
time = df["time"]
original_signal = df["sample"].values  # <class 'numpy.ndarray'>多维数组

# 进行EMD分解
emd = EMD()
imfs = emd(original_signal)
num_imfs = len(imfs)  # 计算imf数量

# 绘制原始信号和imfs
plt.figure(figsize=(17, 7))

# 原始信号
plt.subplot(num_imfs + 1, 1, 1)
plt.plot(time, original_signal, color="orange")
plt.xlabel('Time', labelpad=-13)  # 手动调整标签位置
plt.ylabel('Amplitude')
plt.xticks([time[0], time.iloc[-1]])  # 首尾时间
plt.title("original signal")

# imfs
for i in range(num_imfs):
    plt.subplot(num_imfs + 1, 1, i + 2)
    plt.plot(time, imfs[i], label=f"imf_{i + 1}")
    plt.xticks([time[0], time.iloc[-1]])  # 首尾时间
    plt.legend()
plt.tight_layout()  # 自动调整布局

emd分解后，继续对分解得到的本征模态函数imf做hilbert希尔伯特变换，得到解析信号，可以计算信号的瞬时振幅（包络）：

# imfs---hilbert---瞬时振幅
plt.figure(figsize=(17, 10))
for imf_num, value in enumerate(imfs):
    # imfs
    plt.subplot(np.shape(imfs)[0], 1, imf_num + 1)
    plt.plot(time, value)
    plt.xticks([time.iloc[0], time.iloc[-1]])  # 首尾时间
    plt.title(f"imf_{imf_num + 1}")
    # imf_hilbert
    # plt.subplot(np.shape(imfs)[0], 1, imf_num + 1)
    value_hilbert = hilbert(value)  # hilbert
    instantaneous_amplitude = abs(value_hilbert)  # 瞬时振幅(振幅包络)
    plt.plot(time, instantaneous_amplitude, linestyle="--")
plt.tight_layout()  # 自动调整布局

同时，通过hilbert变换，我们也可以通过计算瞬时相位进一步得到瞬时频率绘制hilbert谱（时间-频率），与fft傅里叶变换得到的频谱图（频率-幅度）进行对比分析：

# 瞬时频率
plt.figure(figsize=(17, 10))
for imf_num, value in enumerate(imfs):
    # imfs_fft
    plt.subplot(np.shape(imfs)[0], 2, 2 * imf_num + 1)
    N = len(value)  # 采样个数
    fs = 10  # 采样频率
    freq = np.fft.fftfreq(N, d=1 / fs)  # 频率轴
    value_fft = np.fft.fft(value)  # fft
    value_fft = 2 / N * abs(value_fft)  # 归一化
    value_fft[0] = value_fft[0] / 2  # 正负频直流分量在同一轴显示
    plt.plot(freq[:N // 2], value_fft[:N // 2], color='orange')  # 取一半
    plt.title(f"imf_{imf_num + 1}_fft")
    # imf_hilbert
    plt.subplot(np.shape(imfs)[0], 2, 2 * imf_num + 2)
    value_hilbert = hilbert(value)  # hilbert
    instantaneous_phase = np.angle(value_hilbert)  # 瞬时振幅(振幅包络)
    instantaneous_phase = np.unwrap(instantaneous_phase)  # 消除相位跳跃
    instantaneous_fre = np.diff(instantaneous_phase) / (2.0 * np.pi)
    plt.plot(time[:-1], abs(instantaneous_fre), color='#6AA84F')  # 截取掉最后一个时间元素
    plt.title(f"imf_{imf_num + 1}_hilbert")
    plt.xticks([time.iloc[0], time.iloc[-1]])  # 首尾时间
    # 图像紧凑不对过多图像信息进行标注
plt.tight_layout()  # 自动调整布局

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD
import pandas as pd
from scipy.signal import hilbert

file_path = "D:/kaggle/data.txt"
df = pd.DataFrame()

with open(file_path, "r", encoding="ISO-8859-1") as file:
    for line_number, line in enumerate(file, 0):  # 调用enumerate函数返回索引和值的元组
        line = line.strip()  # str.strip()去除字符串两端空白字符(空格、制表符、换行符等)
        if len(line) != 0:
            if line_number < 1000:
                print(f"Line{line_number}：{line}")
                df.loc[line_number, "time"] = line.split(";")[1]  # 时间
                df.loc[line_number, "sample"] = float(line.split(";")[7].split("=")[1])

# 原始信号数据
time = df["time"]
original_signal = df["sample"].values  # <class 'numpy.ndarray'>多维数组

# 进行EMD分解
emd = EMD()
imfs = emd(original_signal)
num_imfs = len(imfs)  # 计算imf数量

# 绘制原始信号和imfs
plt.figure(figsize=(17, 7))

# 原始信号
plt.subplot(num_imfs + 1, 1, 1)
plt.plot(time, original_signal, color="orange")
plt.xlabel('Time', labelpad=-13)  # 手动调整标签位置
plt.ylabel('Amplitude')
plt.xticks([time[0], time.iloc[-1]])  # 首尾时间
plt.title("original signal")

# imfs
for i in range(num_imfs):
    plt.subplot(num_imfs + 1, 1, i + 2)
    plt.plot(time, imfs[i], label=f"imf_{i + 1}")
    plt.xticks([time[0], time.iloc[-1]])  # 首尾时间
    plt.legend()
plt.tight_layout()  # 自动调整布局

# imfs---hilbert---瞬时振幅
plt.figure(figsize=(17, 10))
for imf_num, value in enumerate(imfs):
    # imfs
    plt.subplot(np.shape(imfs)[0], 1, imf_num + 1)
    plt.plot(time, value)
    plt.xticks([time.iloc[0], time.iloc[-1]])  # 首尾时间
    plt.title(f"imf_{imf_num + 1}")
    # imf_hilbert
    # plt.subplot(np.shape(imfs)[0], 1, imf_num + 1)
    value_hilbert = hilbert(value)  # hilbert
    instantaneous_amplitude = abs(value_hilbert)  # 瞬时振幅(振幅包络)
    plt.plot(time, instantaneous_amplitude, linestyle="--")
plt.tight_layout()  # 自动调整布局

# 瞬时频率
plt.figure(figsize=(17, 10))
for imf_num, value in enumerate(imfs):
    # imfs_fft
    plt.subplot(np.shape(imfs)[0], 2, 2 * imf_num + 1)
    N = len(value)  # 采样个数
    fs = 10  # 采样频率
    freq = np.fft.fftfreq(N, d=1 / fs)  # 频率轴
    value_fft = np.fft.fft(value)  # fft
    value_fft = 2 / N * abs(value_fft)  # 归一化
    value_fft[0] = value_fft[0] / 2  # 正负频直流分量在同一轴显示
    plt.plot(freq[:N // 2], value_fft[:N // 2], color='orange')  # 取一半
    plt.title(f"imf_{imf_num + 1}_fft")
    # imf_hilbert
    plt.subplot(np.shape(imfs)[0], 2, 2 * imf_num + 2)
    value_hilbert = hilbert(value)  # hilbert
    instantaneous_phase = np.angle(value_hilbert)  # 瞬时振幅(振幅包络)
    instantaneous_phase = np.unwrap(instantaneous_phase)  # 消除相位跳跃
    instantaneous_fre = np.diff(instantaneous_phase) / (2.0 * np.pi)
    plt.plot(time[:-1], abs(instantaneous_fre), color='#6AA84F')  # 截取掉最后一个时间元素
    plt.title(f"imf_{imf_num + 1}_hilbert")
    plt.xticks([time.iloc[0], time.iloc[-1]])  # 首尾时间
    # 图像紧凑不对过多图像信息进行标注
plt.tight_layout()  # 自动调整布局

plt.show()

作者：默寫無言