Vtyjkyuk

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:
beginequation
y(t) = sum_n=1^N D_nJ_0 left( fraclambda_nat right)
endequation
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( fraclambda_nat right) right$ , the FB coefficients $D_n$ are computed by the following:
beginequation
D_n = frac2 int_0^a ty(t) J_0 left( fraclambda_nat right) dta^2
endequation
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(fraclambda_na t)$

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