Faster Fourier - Bessel Coefficient Calculations on Python

Question

So I'm trying to get the fourier-bessel coefficients of a very large array of numbers that are around 1 million points in size, but I'm coming across some speed issues with calculating $J_{o}(x)$ for a large array.

I have my put code below (Python)

import numpy as np
import scipy.integrate as integrate
from scipy.special import j0, j1, jn_zeros

def fourier_bessel_coeffs(x, num_coeffs=None):
   if num_coeffs is None:
       num_coeffs = len(x)

   # Roots of the Zeroth-Order Bessel Function
   zeroth_roots = jn_zeros(0, num_coeffs)

   # Define domain
   t = np.arange(0, len(x))

   # Define end point of domain, which is required as paramter
   a = len(x)

   # Calculate Coefficients
   coeffs = np.array([2*integrate.trapz(t*x*j0(i*t/a)) / (a**2 * j1(i)**2) 
   for i in zeroth_roots])

   return coeffs

The code is based on the following: Any signal $y(t)$ can be represented as a fourier-bessel expansion: \begin{equation} y(t) = \sum_{n=1}^{N} D_{n}J_{0} \left( \frac{\lambda_{n}}{a}t \right) \end{equation} Where $J_{0}(\cdot)$ are zero-order Bessel functions and the $\lambda_{n}; n = 1,2,..., N$ is the ascending order positive roots of $J_{0}(\lambda) = 0$

Using the orthogonality of the set $ \left\{ J_{0} \left( \frac{\lambda_{n}}{a}t \right) \right\}$ , the FB coefficients $D_{n}$ are computed by the following: \begin{equation} D_{n} = \frac{2 \int_{0}^{a} ty(t) J_{0} \left( \frac{\lambda_{n}}{a}t \right) dt}{a^2 |J_{1}(\lambda_{n})|^2} \end{equation} with $1 \leq n \leq N$ and $J_{1}(\lambda_{n})$ are the first order bessel functions.

From what I've been trying out, is that it gets hung up when trying to compute $J_{0}(\frac{\lambda_{n}}{a} t)$

Edit: Note x = $y(t)$ in the code