希尔伯特黄变换(hht)-编程知识

提示：主要对希尔伯特黄变换进行简略的介绍

一、希尔伯特黄变换是什么？

1.1 定义

为了纪念故事中两位老先生（Hilbert和Huang）的突出贡献，人们决定把“经过EMD分解出的IMF分量再经过Hilbert变换，最终得到信号瞬时频率和瞬时幅值”的方法叫做希尔伯特黄变换
来源链接

1.2 公式

a(t):瞬时幅值
在这里插入图片描述
图片链接

二、代码

2.1.代码

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from pyhht import EMD
from scipy.signal import hilbert
import tftb.processingmatplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
matplotlib.rcParams['axes.unicode_minus'] = False  # 显示负号def picture(x, y, N):'''画信号的时域图和频谱输入：x: 0-1时间序列y: 信号N: 1s内采样点数输出：信号的时域图和频谱'''ax1 = plt.subplot(2, 1, 1)plt.plot(x, y)plt.xlabel('时间/s')plt.ylabel('幅值')plt.title('合成信号时域曲线')yf = np.fft.fft(y)xf = np.linspace(0.0, N / 2, N // 2)ax2 = plt.subplot(2, 1, 2)plt.plot(xf, 2.0 / N * np.abs(yf[:N // 2]))  # 频谱幅值归一化，需要*2/Nplt.xlabel('频率/Hz')plt.ylabel('幅值')plt.title('合成信号频谱')plt.show()def HHTAnalysis(t, signal, N):'''进行HHT分析，画出每一个IMF的时域图和频谱输入：t: 0-1时间序列signal: 信号N: 1s内采样点数输出：信号每一个IMF的时域图和频谱且返回IMFs，维度为(IMFs个数，信号点数N)'''# 进行EMD分解decomposer = EMD(signal)# 获取EMD分解后的IMF成分imfs = decomposer.decompose()# 分解后的IMF个数n_components = imfs.shape[0]# 画出每一个IMF的时域图和频谱fig1, axes = plt.subplots(n_components, 2, figsize=(10, 7), sharex='col', sharey=False)for i in range(n_components):# 画时域图axes[i][0].plot(t, imfs[i])axes[i][0].set_title('IMF{}'.format(i + 1))# 画fft图yf = np.fft.fft(imfs[i])xf = np.linspace(0.0, N / 2, N // 2)axes[i][1].plot(xf, 2.0 / N * np.abs(yf[:N // 2]))  # 频谱幅值归一化，需要*2/Naxes[i][1].set_title('IMF{}'.format(i + 1))plt.show()return imfsdef HHTPicture(t, imfs, N, n):'''画出指定个数的IMFs的时域图和时频图输入：t: 0-1时间序列imfs: IMFs成分N: 1s内采样点数n: 指定画前几个IMFs成分输出：前n个IMFs的时域图和时频图'''fig2, axes = plt.subplots(n, 2, figsize=(10, 7), sharex='col', sharey=False)# 计算并绘制各个组分for iter in range(n):# 绘制分解后的IMF时域图axes[iter][0].plot(t, imfs[iter])axes[iter][0].set_xlabel('时间/s')axes[iter][0].set_ylabel('幅值')# 计算各组分的Hilbert变换imfsHT = hilbert(imfs[iter])# 计算各组分Hilbert变换后的瞬时频率# 使用梯形积分法计算解析信号在特定时间瞬时的瞬时频率。instf, timestamps = tftb.processing.inst_freq(imfsHT)# 绘制瞬时频率，这里乘以fs是正则化频率到真实频率的转换# plt.figure(figsize=(10, 7))# plt.plot(timestamps, instf)# # 绘制标题# plt.title('hilbert')fs = Naxes[iter][1].plot(timestamps / fs, instf * fs)axes[iter][1].set_xlabel('时间/s')axes[iter][1].set_ylabel('频率/Hz')# 计算瞬时频率的均值和中位数axes[iter][1].set_title('Freq_Mean{:.2f}----Freq_Median{:.2f}'.format(np.mean(instf * fs), np.median(instf * fs)))def HHTFilter(signal, componentsRetain):'''定义HHT的滤波函数，提取部分EMD组分输入：signol: 信号componentsRetain: IMF的列表 []输出：一幅图，同时包含原始信号和合成信号'''# 进行EMD分解decomposer = EMD(signal)# 获取EMD分解后的IMF成分imfs = decomposer.decompose()# 选取需要保留的EMD组分，并且将其合成信号signalRetain = np.sum(imfs[componentsRetain], axis=0)# 绘图plt.figure(figsize=(10, 7))# 绘制原始数据plt.plot(signal, label='RawData')# 绘制保留组分合成的数据plt.plot(signalRetain, label='HHTData')# 绘制标题plt.title('RawData-----HHTData')# 绘制图例plt.legend()plt.show()return signalRetain# 生成0-1时间序列，共2048个点
N = 2000
t = np.linspace(0, 1, N)
# 生成信号
# signal = (2 + np.cos(8 * np.pi * t)) * np.cos(40 * np.pi * (t + 1) ** 2) + np.cos(
#     20 * np.pi * t + 5 * np.sin(200 * np.pi * t))
signal = 0.8 * np.cos(20 * np.pi * t*2) + 0.5 * np.cos(40 * np.pi * (t + 1) * 2) + 0.2*np.cos(200 * np.pi * t*2) + 0.1 * np.cos(400 * np.pi * t*2)
# 画出原始信号的时域图和频谱
picture(t, signal, N)
# 进行HHT分析，画出所有IMFs的时域图和频谱
imfs = HHTAnalysis(t, signal, N)
# 画出前2个IMFs的时域图和时频图
HHTPicture(t, imfs, N, 4)
# 为了纪念故事中两位老先生（Hilbert和Huang）的突出贡献，人们决定把“经过EMD分解出的IMF分量再经过Hilbert变换，最终得到信号瞬时频率和瞬时幅值”的方法叫做希尔伯特黄变换
# 进行验证，判断与原始信号的差异
signalRetain = HHTFilter(signal, [0, 1])