1

一点背景

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

根据 Huang 的原始论文

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

边际 Hibert 谱由下式给出：

在此处输入图像描述

，其中 Ai 是模式 i 的幅度和相位时间序列。频率 [ωi] 和幅度 [Ai] 的联合概率密度函数 (pdf) P(ω, A)，它们是从所有模式 i = 1 · · · N 一起提取的。

在这种情况下，有人可以指出我如何在给定两个一维数组 IA、IF 的情况下执行积分？

### For minimum reproducible example  
import emd
### Define functions
def hilb(signal,fs, unwrap=False):
    from scipy.signal import hilbert
    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

### Define parameters
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)     ### Length of signal
nbins = 200

## Get IMF's
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)):
    ### Calculate instantaneous frequencies (IF) 
    ### and instantaneous amplitudes (IA)
    
    inst_amp1, inst_f1 = hilb(imf.T[i],fs, unwrap=True)
    inst_amp[:,i]    = inst_amp1[1:]
    inst_freq[:,i]   = inst_f1

### End of code required for minimum reproducible example####   

### Collapse the arrays into 1d, make easier to use     
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]

我试图估计积分的方式。大概是错的？

# 1st step create a Joint probability density function 

# Doing it with log space is probably wrong? Becaus trapezoid integration assumes dA = constant?

#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()

### Get 1) Probability density 2) Edges in x axis 3) Edges in y axis
Pjoint,f_edges, A_edges =p_out[0],p_out[1],p_out[2]


### So you can plot later, not needed for integration
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'$Marginal\ Hilbert \ Transform$')

关于如何估计开头显示的积分的任何想法？

4

0 回答 0