How to Use Fast Fourier Transform (FFT) in Julia 📂Julia

How to Use Fast Fourier Transform (FFT) in Julia

Overview ¹ ²

The Fastest Fourier Transform in the West (FFTW) is a software library developed by Matteo Frigo and Steven G. Johnson at the Massachusetts Institute of Technology (MIT) for computing the Discrete Fourier Transform. While there exists a Julia package named AbstractFFTs.jl for FFT implementation, it is not intended to be used on its own but rather to aid in the implementation of fast Fourier transforms, such as with FFTW.jl.

This package is mainly not intended to be used directly. Instead, developers of packages that implement FFTs (such as FFTW.jl or FastTransforms.jl) extend the types/functions defined in AbstractFFTs. This allows multiple FFT packages to co-exist with the same underlying fft(x) and plan_fft(x) interface.³

Overview

Fourier Transform: fft()
- Transform by column in a $2$-dimensional array: fft( ,[1])
- Transform by row in a $2$-dimensional array: fft( ,[2])
- Transform specific dimensions of an array: fft( ,[n₁, n₂, ...])
Inverse Fourier Transform: ifft()
Frequency centering $0$: fftshift()
- Inverse: ifftshift()
Frequency Sampling: fftfreq(n, fs=1)

Code

Fourier Transform

In Julia, the notations for the Fourier Transform such as $\mathcal{F}[f]$, $\hat{f}$ are used directly in code. Let’s sample sine waves of frequencies $100$, $200$, $350$ at intervals of $1/1000$ and add them together.

using FFTW
using Plots
using LaTeXStrings

Fs = 1000 #진동수
T = 1/1000 #샘플링 간격
L = 1000 #신호의 길이
x = [i for i in 0:L-1].*T #신호의 도메인

f₁ = sin.(2π*100*x) #진동수가 100인 사인파
f₂ = 0.5sin.(2π*200*x) #진동수가 100인 사인파
f₃ = 2sin.(2π*350*x) #진동수가 100인 사인파
f = f₁ + f₂ + f₃

Fourier Transform: fft()
Inverse Fourier Transform: ifft()

According to the definition, the Fourier Transform $f$ of $\mathcal{F}f$ is only nonzero at $50$, $100$, $200$. By the definition of the Discrete Fourier Transform, we obtain symmetric values around the $y$ axis, but the frequency at $0$ is the first value by default. Therefore, if the aim is to check the frequency and amplitude of the signal, plotting only the first half is sufficient.

Fs = 1000 # 샘플링 주파수
ℱf = fft(f) # 푸리에 변환
ξ = Fs*[i for i in 0:L/2-1]/L #주파수 도메인(절반)

plot(ξ, abs.(ℱf[1:Int(L/2)])*2/L, title=L"Fourier transform of ▷eq10◁", label="") 
xlabel!("frequency")
ylabel!("amplitude")
savefig("fft.png")

Frequency Centering $0$

The output of a Fourier Transform inherently places the value at frequency $0$ first. To center the value at frequency $0$, use fftshift(). To reverse this, use ifftshift(), which is not ifft + shift, but rather the inversion + fftshift, meaning it is the inverse operation of fftshift().

p1 = plot(ξ, abs.(ℱf), title=L"▷eq11◁", label="", yticks=[], xticks=[0, 100, 200, 350, 500, 1000]) 
p2 = plot(ξ.-500, abs.(fftshift(ℱf)), title=L"▷eq22◁", label="", yticks=[], xticks=[-500,-350,-200,-100,0,100,200,350,500]) 
plot(p1, p2, size=(800,400))
savefig("fftshift.png")

Multi-dimensional Fourier Transform

To compare with the values of a 2-dimensional Fourier Transform, let’s first calculate the Fourier Transform values of $x = [1\ 2\ 3\ 4]^{T}$.

julia> x = [1.0; 2; 3; 4]
4-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0

julia> fft(x)
4-element Vector{ComplexF64}:
 10.0 + 0.0im
 -2.0 + 2.0im
 -2.0 + 0.0im
 -2.0 - 2.0im

fft() automatically returns a 2-dimensional Fourier Transform when given a 2-dimensional array as input. Alternatively, fft(, [1,2]) means to compute the transform along the first and second dimensions, and it returns the same result.

julia> y = [x x x x]
4×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0
 3.0  3.0  3.0  3.0
 4.0  4.0  4.0  4.0


julia> fft(y)
4×4 Matrix{ComplexF64}:
 40.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0+8.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0-8.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im


julia> fft(y, [1,2])
4×4 Matrix{ComplexF64}:
 40.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0+8.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 -8.0-8.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

Therefore, to transform each column, use fft(, [1]), and to transform each row, use fft(, [2]).

julia> fft(y, [1])
4×4 Matrix{ComplexF64}:
 10.0+0.0im  10.0+0.0im  10.0+0.0im  10.0+0.0im
 -2.0+2.0im  -2.0+2.0im  -2.0+2.0im  -2.0+2.0im
 -2.0+0.0im  -2.0+0.0im  -2.0+0.0im  -2.0+0.0im
 -2.0-2.0im  -2.0-2.0im  -2.0-2.0im  -2.0-2.0im

julia> fft(y, [2])
4×4 Matrix{ComplexF64}:
  4.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  8.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 12.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 16.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

Frequency Sampling

fftfreq(n, fs=1)

Returns the frequency domain of length $n$ with interval $fs/n$. As mentioned earlier, since the Fourier Transform positions the frequency at $0$ at the forefront, the first value sampled by fftfreq() is $0$. The first half is positive frequencies, and the latter half is negative. Hence, using fftshift() arranges them in ascending order according to the index.

julia> fftfreq(4, 1)
4-element Frequencies{Float64}:
  0.0
  0.25
 -0.5
 -0.25

julia> fftfreq(5, 1)
5-element Frequencies{Float64}:
  0.0
  0.2
  0.4
 -0.4
 -0.2

julia> fftshift(fftfreq(4, 1))
-0.5:0.25:0.25

Environment

OS: Windows11
Version: Julia 1.8.2, FFTW 1.5.0