Fourier expansion package

[baizhan]

I have a periodic function, say f(x+T)=f(x), then it can be decomposed into a sum of sin and cos functions. For example,

f(x)=\sum_{s=-\infty}^{\infty}A_{s}e^{isx},

where

A_{s}=\frac{1}{T}\int_{-T/2}^{T/2}f(x)e^{-isx}dx.

Given the f(x), I want to calculate the coefficients A_s. Of course this can be done by numerical integration, but it would be slow, especially for large s when the integration kernel oscillates fast.

I want to know if there is any package which can do this calculation both fast and accurate.
I believe there should be since this is a quite mature field.
However, I can’t find one. The closest I know is the FFTW, but that is only for discrete Fourier transform which is not my case.

Can anybody give me some suggestions?

[Bruno_Amorim]

Why can’t you approximate the Fourier Transform by a Discrete Fourier Transform and just use FFTW?

[stevengj]

@baizhan

The closest I know is the FFTW, but that is only for discrete Fourier transform which is not my case.

The discrete Fourier transform (DFT) is equivalent to a trapezoidal-rule approximation for the Fourier-series integrals. For a smooth periodic function, this converges exponentially fast to the exact integral as you increase the number of points.

So, fast Fourier transforms are, in fact, a fast and accurate way to compute Fourier-series coefficients, especially if you want all the Fourier-series coefficients up to some cutoff (i.e., until the series converges to some tolerance).

If you don’t want all the coefficients, but only want a few high-frequency terms, then there are specialized quadrature methods for highly oscillatory integrals that you might want to look into, e.g. the OscillatoryIntegralsODE.jl package. If your function isn’t smooth and periodic, then you may also need to be more careful.

[mike.ingold]

@baizhan

Given the f(x), I want to calculate the coefficients A_s.

How is your f(x) defined? Do you have a symbolic definition of the function? Is it coming from some source that you need to sample/measure?

[baizhan]

Thank you!

I think what I need is an integration over oscillatory function.
The package you mentioned, however, cannot run on my computer.
Luckily it is not a long program so I think I can read the source to see what I can do.

If I fail to install the package, maybe I can try FFTW as you say.

[baizhan]

I have a formula for my function

f(x)=e^{i\alpha\sin(x)+i\beta\sin(2x)}

I think this function has a quite good behavior.

[stevengj]

@baizhan

I think this function has a quite good behavior.

It has excellent behavior — it’s an entire function, so both its Fourier series and the trapezoidal-rule approximation thereof (the DFT) should converge faster than exponentially.

For example, with \alpha = \beta = 1 the Fourier series converges to machine precision after only about 30 terms (at both positive and negative frequencies), so a DFT with 100 terms should give a Fourier series that is accurate to nearly machine precision everywhere in [0,2\pi].

using FFTW
f(x,α,β) = cis(α*sin(x) + β*sin(2x))
n = 100
x = 2π/n * (0:n-1)
F = fft(f.(x, 1,1) .* (1/n)) # ≈ Fourier series coefficients

using Plots
scatter(abs.(F), yscale=:log10, xlabel="index", ylabel="|Fourier coefficient|", legend=false)

(The cis function is convenient here.) Note that the frequency ordering of the DFT takes some getting used to: F[1] corresponds to zero frequency, and F[2:(n+2)÷2] correspond to “positive” frequencies while F[n:-1:(n+3)÷2] correspond to “negative” frequencies (depending on your sign convention). (The middle n÷2+1 “Nyquist” element for even n is ambiguous, but for sufficiently large n this term is negligible anyway.) This is one context where zero-based indexing makes a lot of sense, which you could do (if you want) by wrapping F with the OffsetArrays.jl package.

If we compare to adaptive numerical integration via QuadGK.jl, we see that they agree to machine precision as expected:

using QuadGK

fourier(f, n) = quadgk(x -> f(x) * cis(-n*x), 0, 2π, rtol=1e-14)[1] / 2π

fourier.(x -> f(x, 1,1), 0:3) ./ F[1:4] # compare the first 4 coefficients

gives

4-element Vector{ComplexF64}:
 0.9999999999999999 + 2.4523045117069327e-17im
 0.9999999999999993 - 2.91266466383833e-16im
 0.9999999999999999 - 2.969943128507037e-16im
                1.0 - 1.6216552943397132e-16im

(i.e. about 1.0 ± ε where ε is close to machine precision).

From the plot, you will also see that the accuracy hits a “noise floor” once the coefficients get so small that they reach machine precision. In this regime, any direct summation approach is vulnerable to cancellation errors that limit accuracy. If you actually want to calculate such tiny coefficients accurately (e.g. you are interested in high-frequency asymptotics), then you need to turn to more specialized methods (either analytical or computational) as described here: Package to evaluate oscillating integrals with increasing frequency - #10 by stevengj

PS. Note also that your function f(x) = \overline{f(-x)} is conjugate-symmetric (assuming \alpha and \beta are real), so the Fourier-series coefficients F are purely real. You could therefore get another factor of ≈2 in performance by using the brfft function in FFTW. Of course if you care that much about performance, you should also precompute your FFT plans, maybe work with a pre-allocated output array using mul!, etcetera etcetera. I didn’t want to get into the complexities of performance optimization here.

[stevengj]

@stevengj

So, fast Fourier transforms are, in fact, a fast and accurate way to compute Fourier-series coefficients, especially if you want all the Fourier-series coefficients up to some cutoff (i.e., until the series converges to some tolerance).

Note also that this process is automated by the ApproxFun.jl package, which implements Fourier-series representation of functions on periodic domains using FFTs internally, as well as lots of other amazing features.

[baizhan]

Thank you so much!

It is really helpful!