Fast way to compute integral of type $\int dx f(x) \cos(n \pi x)$ in SciPy

I have an integral of the form $$I(n) = \int_0^1 dx f(x) \cos(n \pi x) ,$$ where $n$ is an integer. In other words, I calculate the cosine Fourier coefficients of function $f$, which is real and continuous on the interval. I need to calculate this for a large number of large values of $n$. Currently I'm just doing this:

iint = array([ integrate.quad(lambda x: f(x)*cos(n*pi*x), 0, 1, limit=1000) for n in range(0,nlim) ])

where $nlim$ is typically of the order of $\sim 10^3$. I have the feeling that this is not the most optimal way of doing this. How do I make my calculation faster?

• what is f(x) ....? Apr 11 '14 at 11:31
• Is there a reason you don't want to use one of the Fourier Transform libraries? Most of them can compute Discrete Cosine Transforms, too. Apr 11 '14 at 11:53
• $f(x)$ can be one of several well-behaved (smooth, no fast oscillations, real valued, positive) functions. The FFT library in scipy requires me to input the function at discrete space points. If I want to get the output for large values of $n$, I'd need to first make a very dense grid, calculate this function on this grid and then give it to the function that calculates FFT. I doubt this can be much faster than what I'm doing now, but correct me if I'm wrong. Apr 11 '14 at 13:02
• It depends, do you need this for only certain values of $n$ or for all $n$ in a range $[1,N]$ where $N$ is big? Apr 11 '14 at 14:10
• i'd build my own quadrature Apr 11 '14 at 19:02

This is part of the (complex-valued) Fourier transform for which there is (provably) no more efficient way than the Fast Fourier Transform (FFT) if you want to compute the integrals for at least a significant fraction of all frequencies $n$ between zero and the largest frequency you care about.

In other words, while your intuition may tell you that evaluating the function at a large number of points is probably in efficient, the intuition is wrong in this case.

• To amplify Wolfgang: At the very least, you will need a large number of function evaluations to capture the oscillatory nature of $cos(n\pi x)$ if $N$ is large ($O(N)$ of them) and an equivalent number of arithmetic operations on the results. The only thing that using quadrature buys you is the ability not to store the $f(x)$ values at the integration points. By switching to FFTs, you get the advantage of getting all the cosine transform values between 1 and $N$. Apr 12 '14 at 14:23
Isn't this the same as the real part of the Fourier transform of f(x) \begin{align} I(n) &= \int^1_0dxf(x)\cos(n\pi x) \\ &=\Re{{\int^1_0dx f(x) e^{-j n \pi x}}} \\ &= \Re\int^1_0dx f(x) (cos(n\pi x) - j\sin(n \pi x)) \\ &= \int^1_0dxf(x)\cos(n\pi x) \end{align}
If $f$ is very nice, then you can approximate $f$ with a sequence of piecewise polynomials and integrate over the resulting intervals exactly. This may be much, much cheaper than using an FFT. This is also true for approximations of $f$ by any functions where you can easily know or determine the exact antiderivative in each interval. If $f$ is a black box, then you'll have more trouble with this method.