# Fourier expansion package

**URL:** <https://discourse.julialang.org/t/fourier-expansion-package/108621>\
**Category:** New to Julia\
**Tags:** fft\
**Created:** [January 10, 2024, 1:04pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621 "2024-01-10T13:04:25Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![baizhan](https://avatars.discourse-cdn.com/v4/letter/b/ecae2f/32.png) [@baizhan](https://discourse.julialang.org/u/baizhan)\
**Post date:** [January 10, 2024, 1:04pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/1 "2024-01-10T13:04:25Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Bruno\_Amorim](https://avatars.discourse-cdn.com/v4/letter/b/f6c823/32.png) [@Bruno\_Amorim](https://discourse.julialang.org/u/Bruno_Amorim)\
**Post date:** [January 10, 2024, 1:14pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/2 "2024-01-10T13:14:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2024, 1:21pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/3 "2024-01-10T13:21:14Z")

</div>

> [@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](https://github.com/xzackli/OscillatoryIntegralsODE.jl) package. If your function isn’t smooth and periodic, then you may also need to be more careful.

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [January 10, 2024, 1:55pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/4 "2024-01-10T13:55:04Z")

</div>

> [@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?

---

<div class="post-metadata">

**Author:** ![baizhan](https://avatars.discourse-cdn.com/v4/letter/b/ecae2f/32.png) [@baizhan](https://discourse.julialang.org/u/baizhan)\
**Post date:** [January 10, 2024, 2:19pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/5 "2024-01-10T14:19:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baizhan](https://avatars.discourse-cdn.com/v4/letter/b/ecae2f/32.png) [@baizhan](https://discourse.julialang.org/u/baizhan)\
**Post date:** [January 10, 2024, 2:25pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/6 "2024-01-10T14:25:31Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2024, 3:11pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/7 "2024-01-10T15:11:12Z")

</div>

> [@baizhan](#):
>
> I think this function has a quite good behavior.

It has excellent behavior — it’s an [entire function](https://en.wikipedia.org/wiki/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].

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/4/0495d934577b2f27210552e18234ae0c81b57a36.png)  
(The [`cis` function](https://docs.julialang.org/en/v1/base/math/#Base.cis) 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:

```julia
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

```julia
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](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/10)

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2024, 3:29pm UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/8 "2024-01-10T15:29:03Z")

</div>

> [@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](https://juliaapproximation.github.io/ApproxFun.jl/stable/generated/Periodic/#Periodic-Domains) using FFTs internally, as well as lots of other amazing features.

---

<div class="post-metadata">

**Author:** ![baizhan](https://avatars.discourse-cdn.com/v4/letter/b/ecae2f/32.png) [@baizhan](https://discourse.julialang.org/u/baizhan)\
**Post date:** [January 11, 2024, 1:33am UTC](https://discourse.julialang.org/t/fourier-expansion-package/108621/9 "2024-01-11T01:33:42Z")

</div>

Thank you so much!

It is really helpful!