我有一个实验时间信号,我需要计算一些积分。具体来说,我需要计算 PSD,然后计算某些频段的功率。因此,这似乎scipy
是计算积分的最佳方法。但是整个阵列的算法simpson
和trapezoid
计算。没有集成限制
我可以编写一个函数来执行数组的搜索并获取积分限制的索引,然后simpson
将其应用于原始数组的一个切片。但我想知道是否还有其他方法。
谢谢
MWE
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from scipy.integrate import simpson
t_0 = 0
t_N = 5
N = 10000
freq = N/(t_N-t_0)
w_1 = 1
w_2 = 5
x = np.linspace(t_0,t_N,N)
y = np.cos(2*np.pi*w_1*x) + np.sin(2*np.pi*w_2*x) + \
0.05*np.random.randn(x.shape[0])
这是时间信号
plt.plot(x,y)
plt.xlabel('t')
plt.ylabel('signal')
PSD
f,p = welch(y,fs=freq,nperseg=2**13)
plt.plot(f,p, '-*')
plt.xlim([0,10])
plt.xlabel('f [Hz]')
plt.ylabel('PSD')
我的整合功能
def psd_integrate(f, p, f0, fN):
k0, = np.where(f >= f0) # tuple unpack
k0 = k0[0] # get 1st index
kN, = np.where(f <= fN) # tuple unpack
kN = kN[-1] # get last index
return simpson(p[k0:kN], f[k0:kN])
编辑:忘记添加。我正在使用seaborn
,这就是为什么绘图看起来与默认值不同的原因matplotlib