How to perform a numerical integration (with limits) of discrete data in Python using Scipy?

Question

I have an experimental time signal, and I need to compute some integral out of it. In detail, I need to compute the PSD, and then compute the power in some bands of frequencies. So, it seems that scipy is the best way to compute integrals. But the algorithms simpson and trapezoid compute for the whole array. There are no integration limits

I can write a function to perform the search of the arrays and get the index of the integration limits, applying then simpson to an slice of the original array. But I was wondering if there is any other way.

Thanks

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])

This is the time signal

plt.plot(x,y)
plt.xlabel('t')
plt.ylabel('signal')

The 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')

My integration function

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])

EDIT: Forgot to add. I'm using seaborn, that's why the plots looks different from default matplotlib