0

我正在使用 python 中的 emd 包将时间序列分解为 IMF 的。然后我想估计边际希尔伯特谱。

在此处输入图像描述

根据 Huang 的原始论文

https://arxiv.org/pdf/1401.4211.pdf

，傅里叶谱的幂律指数应该与边缘希尔伯特谱的一个匹配。边际 Hibert 谱由下式给出：

在此处输入图像描述,

在哪里，在此处输入图像描述

其中 Ai 是模式 i 的幅度和相位时间序列。

我正在使用的代码如下：

import emd


B = Btotal[2]

## Get IMF's
imf = emd.sift.sift(B)


T = len(B)
    
phase, inst_freq, inst_amp = emd.spectra.frequency_transform(imf, sample_rate, 'hilbert')#for i in range(len(IA.T)):    
h = np.zeros(len(inst_freq.T))
for i in range(len(inst_freq.T)):
    ### Estimate marginal Hilbert Spectrum
    h[i] =(1/T)*np.sum(inst_amp.T[i])

freq = []
for i in range(len(inst_freq.T)):
    freq.append(np.nanmean(inst_freq.T[i]))
### Estimate Fourier Power Spectrum
a11 =psd(B, timeu='s')
        
plt.loglog(a11[0], a11[1],label=r'$Fourier \ Transform$')
plt.loglog(freq, h,label=r'$Marginal\ Hilbert \ Transform$')

fit_fourier =fun.curve_fit_log(a11[0][5::], a11[1][5::])
plt.loglog(a11[0][5::], fit_fourier[2],'k-.',lw=0.5,label=r'$k \ = $'+str(round(fit_fourier[0][1],2)))


fit_H =fun.curve_fit_log(freq[2:], h[2:])
plt.loglog(freq[2:], fit_H[2],'k-.',lw=0.5,label=r'$k \ = $'+str(round(fit_H[0][1],2)))
plt.legend(frameon=False)

关于我可能做错了什么的任何想法？

有第二个方程可以用来估计 h 为：

在此处输入图像描述 其中，频率 [ωi] 和幅度 [Ai] 的联合概率密度函数 (pdf) P(ω, A)，它们是从所有模式 i = 1 · · · N 一起提取的。

在这种情况下，有人可以指出我如何在给定矩阵 IA、IF 的情况下执行积分？

4

0 回答 0