Integration packages to calculate the first N Fourier coefficients

[ccmejia]

I am trying to calculate the first N Fourier coefficients for sampled data (harmonic + gaussian like). The Riemann sum was my first approach, however I think the results are not enough accurate. Therefore, I am trying to choose a more suitable method. The QuadGK package seems to be extremely accurate, but it only works (I think…) for analytical functions. Among the others integration packages, which would be your recommendation?

[stevengj]

@ccmejia

I am trying to calculate the first N Fourier coefficients for sampled data (harmonic + gaussian like).

More infomation is needed.

Sampled data = experimental data? From a simulation? Noisy? Smooth function? If it is from a computer program, then can you evaluate it on demand for any desired x?

When you say “N Fourier coefficients” do you mean that the data is periodic and you want a Fourier series? If it’s not periodic, then what do you mean?

@ccmejia

The Riemann sum was my first approach, however I think the results are not enough accurate.

If you have equally-spaced samples of a periodic function, then a Riemann sum actually converges exponentially fast (with the number of samples) for smooth functions.

(However, adaptive methods like QuadGK can still be advantageous for functions with localized features like sharp peaks.)

And, of course, a Riemann sum of equally-spaced samples is exactly equivalent (with suitable scaling) to a discrete Fourier transform (DFT), so you can compute it very efficiently with FFT algorithms.

@ccmejia

The QuadGK package seems to be extremely accurate, but it only works (I think…) for analytical functions.

QuadGK and similar methods assume that you can evaluate f(x) at any requested x in the integration domain. It doesn’t require you to have an analytical formula for f(x), but you need a way to compute it on demand (e.g. via a computer program).

[IlianPihlajamaa]

In addition to the above answer, if your integrant is not periodic and you have a discrete set of fixed (possibly irregular) sampling points, you can try Filon’s method or related methods. I implemented it once in Julia, but I’m not sure it is available in a package.

Edit: a quick search gives

GitHub - dkm2/FCCQuad.jl: Filon-Clenshaw-Curtis Quadrature

And

[RFC/ANN] OscillatoryIntegralsODE.jl: Levin method + OrdinaryDiffEq

I’m requesting feedback on OscillatoryIntegralsODE.jl, which numerically evaluates integrals I = \int_a^b f(x) S(rx) \, dx where f(x) is smooth and not oscillatory, but S(x) is highly oscillatory for large frequencies r. Examples of oscillatory functions include Bessel functions J_{\nu}(rx), spherical Bessel j_{\nu}(rx), \cos(\omega x), and the harmonic transform e^{i \omega x}. Gaussian quadrature is inefficient for these integrals, since the number of nodes required for accurate integration …

[ccmejia]

Sampled data comes from a function H(V, u) that depends (sum) of transversal harmonic complex potential V(x) and nonlinear terms that depends of a localized complex function u(x), where u(x) is a mode for a nonlinear wave equation in periodical media.

Here is the Riemann approach:

N  # Number of modes
M # Grid size
L # Transversal length
k0 = 2*pi/L
dx = L/M 
x = -L/2:dx=L/2-dx
c = zeros(N)
for n in 1:N
    c[n]=dx*sum(H.*exp.(-im*k0*n*x))/L
end