Understanding the Fourier Transform with Julia

Summary
A comprehensive guide to the Fourier Transform, covering its mathematical foundations, the relationship between time and frequency domains, and practical implementations in Julia. Includes examples of spectral analysis, Gaussian pulses, cosine waves, and the Gibbs phenomenon when building square waves from harmonics.

Quick Overview

TransformDomainUse Case
CFTContinuousTheoretical analysis
DFTDiscreteDigital signal processing
FFTDiscreteFast DFT computation

What is Fourier Transform (FT)?

The Fourier Transform acts as a mathematical “prism” for data. Just as a glass prism splits pure white light into a rainbow of colors (frequencies), the Fourier Transform decomposes a complex signal into the individual sine waves that comprise it.

From a mathematical perspective, the Fourier Transform is essentially a change of basis. It uses integration to compute the inner product (equivalent to a “dot product” for infinite-dimensional vectors), projecting a complex time-domain function onto an orthogonal basis consisting of sine and cosine functions (or complex exponentials).

The core value of this operation lies in demonstrating that “arbitrary” functions (meeting certain conditions) can be disassembled into a linear combination of simple, well-studied sine functions. Through this decomposition, difficult or irregular signals are transformed into clear coefficients of frequency, amplitude, and phase, allowing us to analyze them within a much simpler mathematical framework: the frequency domain.

Why sine functions?

Sine and cosine functions (or the complex exponential $e^{j\omega t}$) are used because they are the eigenfunctions of Linear Time-Invariant (LTI) systems. Sine waves are the only signals that maintain their fundamental sinusoidal shape—changing only in amplitude and phase—when passing through an LTI system. Since the purpose of the Fourier Transform is often to analyze how a signal behaves inside these systems, using sine waves as basis functions is the most mathematically efficient and stable approach. It allows the system’s response to be analyzed one frequency at a time.

Time domain vs. Frequency domain

Time and frequency domain analyses represent two distinct facets of signal observation. The time domain plots dynamic signals against a time axis, offering an intuitive and visual perspective. In contrast, the frequency domain maps signals onto a frequency axis, providing a more concise framework for deep analysis. While modern signal processing increasingly favors the frequency domain, the two approaches are complementary and indispensable, bridged directly by the Fourier Transform.

Continuous Fourier Transform (CFT)

The Continuous Fourier Transform (CFT) utilizes integrals to bridge the time and frequency domains, serving as a fundamental tool for the theoretical analysis of continuous analog signals. By decomposing a complex signal—such as a musical track mixing instruments and vocals—into its constituent frequencies, the CFT transforms difficult-to-analyze time-domain data into a manageable frequency spectrum. This enables precise processing, such as isolating vocals or filtering noise, before the modified signal is converted back to its original form via the Inverse Fourier Transform.

Forward Transform (Time $\to$ Frequency):

$$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-j\omega t} , dt$$

Inverse Transform (Frequency $\to$ Time):

$$f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{j\omega t} , d\omega$$

Discrete Fourier Transform (DFT)

Since computers are fundamentally limited to processing finite, discrete data (binary 0s and 1s) rather than continuous or infinite signals, they cannot utilize the integrals found in continuous mathematics. Consequently, the Discrete Fourier Transform (DFT) replaces the continuous integral with a finite summation ($\sum$). This adaptation allows real-world analog signals (like audio) to be discretized and processed digitally, making the DFT the essential mathematical engine for computer-based signal processing.

Forward Transform:

$$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}nk}$$

Inverse Transform:

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j\frac{2\pi}{N}nk}$$

Key Variables

  • $t$: Continuous time
  • $\omega$: Angular frequency ($2\pi f$)
  • $n$: Discrete time sample index
  • $k$: Discrete frequency bin index
  • $N$: Total number of samples
  • $j$: Imaginary unit ($\sqrt{-1}$)

Fast Fourier Transform (FFT)

The Fast Fourier Transform (FFT) is not a distinct transformation but rather a highly efficient algorithm used to compute the Discrete Fourier Transform (DFT). By drastically reducing computational complexity from $O(N^2)$ to $O(N \log N)$, the FFT—combined with advancements in digital hardware—has made the DFT a practical and indispensable tool in modern signal processing.

Implementation in Julia

FFTW.jl is the Julia interface to the “Fastest Fourier Transform in the West” (FFTW) library, a widely used and highly optimized C library for computing Discrete Fourier Transforms (DFTs). It serves as the standard, high-performance FFT engine within the Julia ecosystem.

Signal with noise

Here is a specific example of using the Fourier Transform to extract frequency components from a signal hidden in noise (Spectral Analysis).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
using FFTW
using Plots
using Random

# Sampling parameters
Fs = 1000.0
T = 1/Fs
L = 1500
t = (0:L-1) * T

# Signal: 0.7*sin(100Hz) + sin(300Hz)
S = 0.7 * sin.(2π * 100 * t) .+ sin.(2π * 300 * t)

# Add Gaussian noise (std = 2)
X = S .+ 2 * randn(length(t))

# Plot noisy signal (first 50 samples)
plot(1000 * t[1:50], X[1:50],
    title = "Signal with Noise",
    xlabel = "t (ms)",
    ylabel = "X(t)",
    legend = false)

By observing the time-domain signal $X(t)$, it is difficult to determine the frequency components. Therefore, we use the Fourier Transform to perform spectral analysis. Here, we use the Fast Fourier Transform: Y = fft(X).

First, we compute the normalized two-sided spectrum ($P_2$). Exploiting the symmetry of real-valued signals, we then derive the single-sided spectrum ($P_1$) by discarding the redundant negative frequencies. The abs function is used here to calculate the modulus (magnitude) of the complex FFT coefficients, which corresponds to the signal amplitude. Note that while complex numbers also possess a phase attribute (representing direction), this analysis focuses strictly on magnitude. Finally, we multiply the positive frequency amplitudes by 2 to conserve the total energy from the discarded negative side.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# Compute FFT
Y = fft(X)

# Two-sided spectrum -> Single-sided spectrum
P2 = abs.(Y / L)
P1 = P2[1:L÷2+1]
P1[2:end-1] .= 2 * P1[2:end-1]

# Frequency axis
f = Fs * (0:(L÷2)) / L

# Plot single-sided amplitude spectrum
p2 = plot(f, P1,
    title = "Single-Sided Amplitude Spectrum",
    xlabel = "f (Hz)",
    ylabel = "|P1(f)|",
    legend = false)

# Find two largest peaks
sorted_idx = sortperm(P1, rev=true)
peak1_idx, peak2_idx = sorted_idx[1], sorted_idx[2]
peak1_f, peak1_amp = f[peak1_idx], P1[peak1_idx]
peak2_f, peak2_amp = f[peak2_idx], P1[peak2_idx]

# Annotate peaks
scatter!([peak1_f, peak2_f], [peak1_amp, peak2_amp], 
    markersize = 6, color = :red)
annotate!(peak1_f, peak1_amp + 0.05, 
    text("$(round(peak1_f, digits=1)) Hz", 8, :center))
annotate!(peak2_f, peak2_amp + 0.05, 
    text("$(round(peak2_f, digits=1)) Hz", 8, :center))

Analysis: As seen in the plot (referencing the original image), the amplitudes at 100Hz and 300Hz are close to 0.7 and 1, respectively. They are not exact because noise was added to the signal; without noise, they would correspond exactly to 0.7 and 1. You can verify this by removing the noise.

This demonstrates that while it is difficult to view frequency components in the time domain, analyzing the signal becomes much simpler in the frequency domain after applying the Fourier Transform.

Gaussian pulses

A “pulse” essentially refers to a signal transmission formed by momentarily switching a connection on and off. This mechanism is standard in modern stepper and servo motor drivers, which are controlled via pulses.

A Gaussian pulse is a signal defined by a time-intensity curve that follows a Gaussian shape, given by the function $f(x)=ae^{-\frac{(x - b)^2}{2c^2}}$. Here, $a$ determines the peak height, $b$ sets the center position (horizontal offset), and $c$ controls the width (standard deviation). Because the Gaussian pulse exhibits the characteristics of a normal distribution, it offers excellent performance in both time and frequency domains, leading to its widespread application in many fields.

Let’s examine the transformation of a Gaussian pulse from the time domain to the frequency domain:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
using FFTW, Plots

# Sampling parameters
Fs = 100.0
t = -0.5:1/Fs:0.5
L = length(t)

# Gaussian pulse: X(t) = A * exp(-t² / 2σ²)
σ² = 0.01
A = 1 / (4 * sqrt(2π * σ²))
X = A * exp.(-t.^2 / (2σ²))

# Plot in time domain
plot(t, X,
    title = "Gaussian Pulse in Time Domain",
    xlabel = "Time (s)",
    ylabel = "X(t)",
    linewidth = 2,
    legend = false)

As shown in the graph, the signal forms a Gaussian distribution centered at $0$ with a peak of approximately $0.9974$. The formula indicates a width parameter of $0.01$. This function is a classic eigenfunction in Fourier analysis—meaning a Gaussian in the time domain transforms into a Gaussian in the frequency domain.

Before converting to the frequency domain, we determine a new input length equal to the next power of 2 relative to the original signal length. Using a power-of-2 length improves FFT performance (nextpow2) by padding the signal with trailing zeros.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
# Zero-padding to next power of 2 for efficient FFT
n = nextpow(2, L)
X_padded = [X; zeros(n - L)]

# Compute FFT
Y = fft(X_padded)

# Frequency axis and normalized magnitude
f = Fs * (0:n÷2) / n
P = abs.(Y[1:n÷2+1]) / n

# Plot frequency domain
p2 = plot(f, P,
    title = "Gaussian Pulse - Frequency Domain",
    xlabel = "Frequency (Hz)",
    ylabel = "|P(f)|",
    linewidth = 2,
    legend = false)

Cosine wave

Let’s examine one of the most common waveforms and compare how it behaves at different frequencies:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
using FFTW, Plots

# Sampling parameters
Fs = 1000.0
T = 1 / Fs
L = 1000
t = (0:L-1) * T

# Cosine waves at different frequencies
frequencies = [50, 150, 300]
X = [cos.(2π * f * t) for f in frequencies]

# Plot first 100 samples of each cosine wave
plots = [plot(t[1:100], x[1:100],
    title = "Cosine Wave $(f) Hz",
    xlabel = "Time (s)",
    ylabel = "Amplitude",
    linewidth = 2,
    legend = false) for (x, f) in zip(X, frequencies)]

plot(plots..., layout = (3, 1), size = (800, 600))

Here, we observe three cosine waves. While we can say that the lower waves have higher frequencies than the ones above, it is still difficult to clearly discern their essential characteristics. Let’s convert them into the frequency domain to see the effect.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Zero-padding to next power of 2 for efficient FFT
n = nextpow(2, L)

# Compute FFT for each signal
Y = [fft([x; zeros(n - L)]) for x in X]

# Two-sided spectrum -> Single-sided spectrum
P2 = [abs.(y) / L for y in Y]
P1 = [p[1:n÷2+1] for p in P2]
for p in P1
    p[2:end-1] .= 2 * p[2:end-1]
end

# Frequency axis
f = 0:(Fs/n):(Fs/2)

# Plot single-sided amplitude spectrum (Frequency Domain)
p_freq = [plot(f[1:n÷2], p[1:n÷2],
    title = "Cosine Wave $(freq) Hz - Frequency Domain",
    xlabel = "Frequency (Hz)",
    ylabel = "|P(f)|",
    linewidth = 2,
    legend = false) for (p, freq) in zip(P1, frequencies)]

plot(p_freq..., layout = (3, 1), size = (800, 600))

As we can see, once transformed from the time domain to the frequency domain, the problem becomes simple and clear. Of course, for the final three cosine waves, we can use the method from the earlier noise example: add the three together, transform from time to frequency, and compare. This confirms that while the signal is relatively difficult to process in the time domain, it becomes much simpler to analyze once transformed into the frequency domain.

Odd harmonics

Harmonics are obtained by decomposing a periodic AC quantity using a Fourier Series. The component with the same frequency as the power supply (e.g., 50Hz) is called the fundamental frequency. Components with frequencies that are integer multiples of the fundamental (greater than 1) are called harmonics.

  • A wave with a frequency 3 times the fundamental is called the 3rd harmonic; 5 times is the 5th harmonic, and so on. These odd-numbered multiples are known as odd harmonics.
  • Harmonics are typically caused by applying a sinusoidal voltage to a non-linear load, causing the fundamental current to distort.

Common non-linear loads include UPS systems, Switching Mode Power Supplies (SMPS), rectifiers, Variable Frequency Drives (VFDs), and inverters.

Below, we will see how to construct a square wave by superimposing sine waves.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
using Plots

# Time vector (0 to 2π)
t = 0:0.1:2π

# Fundamental frequency: sin(t)
y1 = sin.(t)

# Add 3rd harmonic: sin(t) + sin(3t)/3
y2 = sin.(t) .+ sin.(3t) / 3

# Add odd harmonics: sin(t) + sin(3t)/3 + sin(5t)/5 + sin(7t)/7 + sin(9t)/9
y3 = sum(sin.(n * t) / n for n in 1:2:9)

# Plot fundamental wave
p1 = plot(t, y1,
    title = "Fundamental Frequency",
    xlabel = "t (rad)",
    ylabel = "Amplitude",
    linewidth = 2,
    legend = false)

# Plot with 3rd harmonic
p2 = plot(t, y2,
    title = "Fundamental + 3rd Harmonic",
    xlabel = "t (rad)",
    ylabel = "Amplitude",
    linewidth = 2,
    legend = false)

# Plot with odd harmonics (1, 3, 5, 7, 9)
p3 = plot(t, y3,
    title = "Odd Harmonics (1st-9th)",
    xlabel = "t (rad)",
    ylabel = "Amplitude",
    linewidth = 2,
    legend = false)

plot(p1, p2, p3, layout = (3, 1), size = (800, 700))

We can see the sine wave beginning to morph into a square wave shape. Now, let’s add even more harmonics:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# Build square wave using Fourier series (Gibbs phenomenon)
harmonics = 1:2:29
y = zeros(length(harmonics), length(t))
x = zeros(length(t))

for (i, k) in enumerate(harmonics)
    x .+= sin.(k * t) / k
    y[i, :] = x
end

# Plot progressive harmonic addition
plot(t, y',
    title = "Building Square Wave: Gibbs Phenomenon",
    xlabel = "t (rad)",
    ylabel = "Amplitude",
    linewidth = 1.5,
    legend = false)

Here, we can clearly see that despite superimposing many sine waves, there is essentially permanent, violent oscillation/jitter at the edges of the waveform. This is the Gibbs phenomenon. In other words, no matter how many sine waves you add, the peak amplitude of this overshoot remains constant.

The Gibbs phenomenon demonstrates that a Fourier Series cannot perfectly represent discontinuous signals (signals with jumps). Therefore, the “arbitrary” functions Fourier referred to in his transform theory must technically be continuous functions. Looking back, when Fourier submitted his original paper to Lagrange for review, Lagrange’s critique—that “sine curves cannot be combined to form a signal with sharp corners”—was mathematically correct. It can only result in an approximation.

Let’s look at a 3D view for a more intuitive perspective:

1
2
3
4
5
6
7
8
# 3D surface plot: sine wave to square wave transition
surface(1:length(t), 1:length(harmonics), y,
    title = "Sine to Square Wave Transition",
    xlabel = "Time",
    ylabel = "Harmonics",
    zlabel = "Amplitude",
    colorbar = false,
    yflip = true)