希尔伯特变换获取频率数组时，如何确定各IMF对应的正确频率

边际希尔伯特谱的正确计算与绘图方法 #

1. 现有代码的问题 #

你当前的实现是对每个IMF的全时间域振幅平方求和，本质是计算每个IMF的总平均能量，不符合边际希尔伯特谱的定义。边际谱的核心是对希尔伯特时频谱在时间维度积分，输出是频率的函数，而非IMF序号的函数。

注意你提供的示例代码中emd.spectra.frequency_transform的sample_rate参数需要传入1，对应dt=1s的采样设置。

2. 正确计算步骤 #

步骤1：定义频率分箱

由于采样间隔dt=1s，采样率为1Hz，奈奎斯特频率为0.5Hz，首先将需要分析的频率范围划分为若干等宽分箱，分箱数量可根据所需分辨率调整：

import numpy as np
# 示例：将0~0.5Hz划分为100个等宽分箱
freq_bins = np.linspace(0, 0.5, 100)
h = np.zeros_like(freq_bins)
T = len(B)

步骤2：按频率维度累加能量

遍历所有IMF的所有时刻点，将每个时刻的能量累加到其瞬时频率对应的分箱中，最后除以总时长得到边际谱：

for imf_idx in range(IF.shape[1]):
    # 取当前IMF的瞬时频率、振幅序列
    cur_if = IF[:, imf_idx]
    cur_a2 = A[:, imf_idx] ** 2
    # 过滤掉超出频率范围的异常值（希尔伯特变换偶尔会产生负频率/超奈奎斯特频率的噪点）
    valid_mask = (cur_if >= freq_bins.min()) & (cur_if <= freq_bins.max())
    valid_if = cur_if[valid_mask]
    valid_a2 = cur_a2[valid_mask]
    # 映射频率到对应分箱索引
    bin_idxs = np.digitize(valid_if, freq_bins) - 1
    # 累加能量
    for idx, a2 in zip(bin_idxs, valid_a2):
        if 0 <= idx < len(h):
            h[idx] += a2
# 归一化到总时长
h = h / T

步骤3：绘图

直接以分箱的频率值为x轴，边际谱能量为y轴绘图即可：

import matplotlib.pyplot as plt
plt.plot(freq_bins, h)
plt.xlabel('频率 (Hz)')
plt.ylabel('边际谱能量')
plt.show()

常见疑问说明 #

不需要为每个IMF指定单一频率值：IMF的瞬时频率是随时间变化的，每个IMF覆盖的是一段频率范围，若强行用平均频率/中心频率对应单个能量值，会丢失频率分辨率，不符合边际谱的物理意义。

内容的提问来源于stack exchange，提问作者jokerp