首页 > 解决方案 > 如何使用 Scipy 在 Python 中执行离散数据的数值积分(有限制)?

问题描述

我有一个实验时间信号,我需要计算一些积分。具体来说,我需要计算 PSD,然后计算某些频段的功率。因此,这似乎scipy是计算积分的最佳方法。但是整个阵列的算法simpsontrapezoid计算。没有集成限制

我可以编写一个函数来执行数组的搜索并获取积分限制的索引,然后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')

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

标签: pythonscipysimpsons-rule

解决方案


推荐阅读