Scipy-Numpy计算联合概率密度积分求解边际希尔伯特谱方法问询

我使用Python的emd包将时间序列分解为IMF，接下来需要估算边际希尔伯特谱（Marginal-Hilbert Spectrum）。
根据Huang发表的原始论文，边际希尔伯特谱的计算公式如下：

式中Ai为第i个模态的振幅时间序列，P(ω, A)是从所有i = 1 · · · N个模态中提取的频率[ωi]和振幅[Ai]的联合概率密度函数（pdf）。

问题1：现有瞬时振幅IA、瞬时频率IF两个一维数组，如何正确实现上述积分计算？ #

最小可复现代码 #

import numpy as np
import emd
from matplotlib import pyplot as plt
from scipy.signal import hilbert

### 定义函数
def hilb(signal,fs, unwrap=False):
    analytic_signal = hilbert(signal)
    A = np.abs(analytic_signal)
    i_p = np.unwrap(np.angle(analytic_signal))
    i_f = np.diff(i_p) /(2.0*np.pi) * fs
    return  A,i_f

### 定义参数
dt = 1
fs = 1/dt
B = np.sin(np.linspace(0,2*np.pi,1000))+0.05*np.random.rand(1000)
T = dt*len(B)     # 信号长度
nbins = 200

## 分解得到IMF
imf = emd.sift.sift(B)

inst_freq = np.zeros((len(imf)-1,len(imf.T)))
inst_amp = np.zeros((len(imf)-1,len(imf.T)))
for i in range(len(imf.T)):
    # 计算瞬时频率(IF)和瞬时振幅(IA)
    inst_amp1, inst_f1 = hilb(imf.T[i],fs, unwrap=True)
    inst_amp[:,i] = inst_amp1[1:]
    inst_freq[:,i] = inst_f1

### 展开为1维数组方便后续处理 
IA_flat = np.ravel(inst_amp)
IF_flat1 = np.ravel(inst_freq)
# 过滤掉负频率
IF_flat = IF_flat1[IF_flat1>0]
IA_flat = IA_flat[IF_flat1>0]

问题2：当前尝试的积分估算方法是否存在错误？ #

现有尝试代码 #

# 第一步：计算联合概率密度函数

# 用对数坐标分箱大概率是错的？因为梯形积分假设dA是常数？
#IA_bins = np.logspace(np.log10(min(IA_flat)),np.log10(max(IA_flat)),nbins)
#IF_bins = np.logspace(np.log10(min(IF_flat[IF_flat>0])),np.log10(max(IF_flat[IF_flat>0])),nbins)
IA_b = np.linspace((min(IA_flat[IF_flat>0])),(max(IA_flat[IF_flat>0])),nbins)
IF_b = np.linspace((min(IF_flat[IF_flat>0])),(max(IF_flat[IF_flat>0])),nbins)
p_out = plt.hist2d(IF_flat[IF_flat>0],IA_flat[IF_flat>0],bins=[IA_bins,IF_bins],density=True)
plt.clf()

### 提取结果：1）联合概率密度 2）x轴边界 3）y轴边界
Pjoint,f_edges, A_edges =p_out[0],p_out[1],p_out[2]

### 计算频率中心用于后续绘图，积分不需要
f_centers = (0.5)*np.diff(f_edges) + f_edges[:-1]

h =[]
for i in range(len(f_centers)):
    index = np.where((IF_flat>f_edges[i]) & (IF_flat<f_edges[i+1]))[0].astype(int)
    H = np.nanmean(IA_flat[index]**2)
    P = Pjoint[i]
    P = P[P>0]
    h.append(np.trapz(P*H))

h = np.array(h)

plt.loglog(f_centers[h>0],h[h>0],label=r'$边际希尔伯特谱$')

错误说明与正确实现 #

现有代码的核心错误 #

1. hist2d参数顺序错误：plt.hist2d第一个入参对应x轴（频率），第二个对应y轴（振幅），你传入的bins顺序是[振幅分箱, 频率分箱]，和数据顺序完全相反，得到的联合密度矩阵维度错位，后续索引全部错误。
2. 积分逻辑不符合公式要求：原公式是对每个固定频率ω，计算∫ A² * P(ω,A) dA，你提前对每个频率区间的A平方取均值，再和密度相乘积分，完全偏离了公式定义。
3. 冗余计算：已经通过hist2d得到了离散化的联合密度，不需要再循环从原始数组里捞对应频率的样本，重复计算且容易出错。

正确实现代码 #

import numpy as np
from matplotlib import pyplot as plt

nbins = 200
# 分别对频率、振幅生成线性分箱
IF_bins = np.linspace(np.min(IF_flat), np.max(IF_flat), nbins)
IA_bins = np.linspace(np.min(IA_flat), np.max(IA_flat), nbins)

# 计算联合概率密度，注意bins顺序和数据顺序对应
p_out = plt.hist2d(IF_flat, IA_flat, bins=[IF_bins, IA_bins], density=True)
plt.clf()
P_joint, f_edges, a_edges = p_out[0], p_out[1], p_out[2]

# 计算分箱中心
f_centers = (f_edges[:-1] + f_edges[1:]) / 2
a_centers = (a_edges[:-1] + a_edges[1:]) / 2

# 对每个频率维度，计算A²*P(ω,A)对A的积分
marginal_hs = np.trapz((a_centers**2)[np.newaxis, :] * P_joint, x=a_centers, axis=1)

# 绘图
plt.loglog(f_centers[marginal_hs>0], marginal_hs[marginal_hs>0], label='边际希尔伯特谱')
plt.xlabel('频率')
plt.ylabel('能量')
plt.legend()
plt.show()

补充说明 #

• 如果你的频率/振幅动态范围很大，必须用对数分箱的话，积分的时候不需要额外调整，np.trapz会根据传入的a_centers自动适配不等步长的积分。
• 也可以直接用emd包自带的希尔伯特谱计算接口，不需要手动实现。

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