Signal Processing 101: The Math Behind the Fourier Transform

Connect an oscilloscope to a network cable or a radio antenna and you won’t see clean 1s and 0s. You’ll see a continuous, noisy voltage waveform — multiple frequencies mixed together, corrupted by thermal noise, smeared by the channel. The hardware’s job is to pull structured information out of that mess.

The tool for doing that isn’t time-domain analysis. It’s a change of representation: transforming the signal from the time domain into the frequency domain, where individual frequency components become visible and separable.

Most engineers reach for fftw_execute() or numpy.fft.fft() and treat what comes out as oracle output. That’s fine for a lot of applications. But if you’re building signal processors, writing audio codecs, or debugging a spectrum analyzer, you need to know what’s actually happening inside — and it turns out to be something conceptually clean: a continuous, infinite-dimensional dot product.

The Rotor: Euler’s Complex Exponential

The mathematical engine of the Fourier Transform is Euler’s formula:

$e^{i\theta} = \cos(\theta) + i\sin(\theta)$

This isn’t just a trigonometric identity. It describes rotation in the complex plane — a 2D space where the horizontal axis is the real numbers and the vertical axis is the imaginary numbers ($i\mathbb{R}$). As $\theta$ increases from 0 to $2\pi$, the point $e^{i\theta}$ traces one full counterclockwise circle of radius 1 around the origin.

Adding a frequency parameter gives us control over how fast that circle is traced:

$e^{-i 2\pi f t}$

• $2\pi$ normalizes the rotation so that $t = 1$ corresponds to exactly one full revolution
• $f$ is the frequency in Hz — how many revolutions per second
• The negative sign makes the rotation clockwise (convention in signal processing; some fields use the positive sign)

This term is the probe. For any target frequency $f$ we care about, $e^{-i 2\pi f t}$ is a reference oscillator spinning at exactly that rate. The Fourier Transform uses it to interrogate a signal.

The Continuous Fourier Transform

The Continuous Fourier Transform (CFT) takes a time-domain signal $x(t)$ and produces its frequency-domain representation $X(f)$:

$X(f) = \int_{-\infty}^{\infty} x(t)\, e^{-i 2\pi ft}\, dt$

Breaking this into two operations makes the mechanics concrete.

Wrapping the Signal

The product $x(t) \cdot e^{-i 2\pi ft}$ takes the signal’s amplitude at each moment in time and uses it to modulate the complex rotor. Geometrically, this “wraps” the signal around the complex plane at frequency $f$: think of winding a wire around a cylinder, where the wire’s distance from the center at each point is proportional to the signal’s amplitude at that time.

The Center of Mass

The integral sums this wrapped shape over all time — it computes the center of mass of the wound-up signal.

If $x(t)$ does not contain frequency $f$, the wrapping is incoherent. Peaks fall on different sides of the circle at different times, and they cancel. The center of mass stays at the origin: $X(f) = 0 + 0i$.

If $x(t)$ does contain a sinusoid at frequency $f$, the peaks always land on the same side of the circle — because the signal and the rotor are spinning at the same rate, keeping them phase-aligned. The wound shape is lopsided, and the center of mass shifts away from the origin.

The resulting complex number $X(f)$ carries two pieces of information:

• Magnitude $|X(f)|$: how far the center of mass shifted — the strength of frequency $f$ in the signal
• Phase $\angle X(f)$: the angle of the shift — the time offset of that frequency component relative to the reference

The Linear Algebra View: Orthogonal Projection

There’s a clean algebraic reason why the integral cancels for non-matching frequencies. It comes from orthogonality.

In standard linear algebra, projecting vector $\vec{A}$ onto vector $\vec{B}$ uses the dot product. If $\vec{A}$ and $\vec{B}$ are orthogonal, the projection is zero. The Fourier Transform is the same operation for continuous functions — the dot product becomes an integral:

$\int_{-\infty}^{\infty} e^{i 2\pi f_1 t} \cdot e^{-i 2\pi f_2 t}\, dt = 0 \quad \text{when } f_1 \neq f_2$

The set of all complex exponentials $\{e^{-i 2\pi f t}\}$ forms an orthogonal basis for the space of square-integrable functions — the function-space analogue of the unit vectors $\hat{i}, \hat{j}, \hat{k}$ in 3D space. When you apply the Fourier Transform, you’re projecting your time-domain signal onto each of these basis vectors and reading off the coefficient. That coefficient is $X(f)$.

The change of basis framing also makes the inverse Fourier Transform intuitive: you’re just reconstructing the original vector from its basis coefficients.

From Calculus to Silicon: The DFT

A processor can’t integrate to infinity over continuous time. Real hardware works with finite buffers of discrete samples. Discretizing the CFT gives the Discrete Fourier Transform (DFT):

• Continuous time $t \to$ sample index $n \in \{0, 1, \ldots, N-1\}$
• Continuous frequency $f \to$ frequency bin $k \in \{0, 1, \ldots, N-1\}$
• Infinite integral $\to$ finite sum over $N$ samples

$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i \frac{2\pi}{N} k n}$

Expanding via Euler’s formula separates the real and imaginary components into plain floating-point arithmetic:

$\text{Re}[X_k] = \sum_{n=0}^{N-1} x_n \cos\!\left(\frac{2\pi k n}{N}\right)$

$\text{Im}[X_k] = -\sum_{n=0}^{N-1} x_n \sin\!\left(\frac{2\pi k n}{N}\right)$

The complexity is $O(N^2)$: for each of the $N$ output bins $k$, you perform $N$ multiply-accumulate operations over all input samples. That’s the nested loop structure directly.

Bin Resolution and the Nyquist Limit

Two practical constraints shape how you interpret DFT output.

Bin resolution is the frequency spacing between adjacent bins: $\Delta f = \frac{f_s}{N}$, where $f_s$ is the sample rate in Hz. With a 1000 Hz sample rate and $N = 256$ samples, each bin spans approximately 3.9 Hz. A signal at exactly 15 Hz will appear in the nearest bin — 15.62 Hz in this case. To improve resolution, increase $N$ (longer capture window) or lower $f_s$.

The Nyquist limit is more fundamental. The DFT produces $N$ output bins, but only the first $N/2$ correspond to physically meaningful frequencies. The upper half is the mirror image of the lower half (for real-valued input signals). The highest frequency the DFT can represent is $f_s / 2$, known as the Nyquist frequency. Any sinusoid in the signal with frequency above $f_s / 2$ cannot be distinguished from a lower-frequency alias — it folds back and appears at $f_s - f_{actual}$. This is aliasing, and it’s the reason audio hardware sampling at 44.1 kHz uses an analog anti-aliasing filter to remove components above 22.05 kHz before the ADC.

Spectral Leakage and Windowing

One more practical issue worth knowing before implementing anything: the DFT implicitly assumes the input buffer is one period of a signal that repeats indefinitely. If the signal frequencies don’t align exactly with bin centers — which they usually don’t — there’s a discontinuity at the buffer boundaries. This smears energy from one frequency bin into its neighbors, a phenomenon called spectral leakage.

The fix is a window function: multiply the signal by a smooth taper that goes to zero at both ends of the buffer before computing the DFT. Common choices are the Hann window ($w_n = 0.5(1 - \cos(2\pi n / N))$) and the Hamming window. The trade-off is frequency resolution — windowing widens the main lobe of each frequency peak. For the demonstration below, we skip windowing to keep the code focused on the DFT itself, but any production spectrum analyzer applies one.

Implementation in C

The following generates a composite signal of 15 Hz and 40 Hz sine waves with added noise, then runs a raw DFT to recover the frequency components. No complex number libraries — the real and imaginary accumulators are explicit.

A fixed random seed makes the output deterministic. The magnitude is normalized by $N/2$ for non-DC bins (DC bin $k=0$ would normalize by $N$, since it has no negative-frequency mirror).

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N           256
#define SAMPLE_RATE 1000
#define PI          3.14159265358979323846

int main(void) {
    double x[N];
    double X_re[N], X_im[N], mag[N];

    /* Reproducible noise across runs */
    srand(42);

    /* 1. Build signal: 15 Hz (amplitude 1.0) + 40 Hz (amplitude 0.5) + noise */
    for (int n = 0; n < N; n++) {
        double t     = (double)n / SAMPLE_RATE;
        double noise = ((double)rand() / RAND_MAX) - 0.5;
        x[n] = 1.0 * sin(2.0 * PI * 15.0 * t)
             + 0.5 * sin(2.0 * PI * 40.0 * t)
             + noise;
    }

    /* 2. DFT — O(N^2): for each bin k, correlate against the reference rotor */
    for (int k = 0; k < N; k++) {
        X_re[k] = 0.0;
        X_im[k] = 0.0;
        for (int n = 0; n < N; n++) {
            double angle = (2.0 * PI * k * n) / N;
            X_re[k] += x[n] *  cos(angle);
            X_im[k] -= x[n] *  sin(angle);
        }
        /* Normalize: divide by N/2 for bins 1..N/2-1 (bin 0 would use N) */
        mag[k] = sqrt(X_re[k]*X_re[k] + X_im[k]*X_im[k]) / (N / 2.0);
    }

    /* 3. Report peaks up to the Nyquist limit */
    printf("Frequency bins above threshold (magnitude > 0.3):\n");
    printf("--------------------------------------------------\n");
    for (int k = 1; k <= N / 2; k++) {
        double freq_hz = (double)k * SAMPLE_RATE / N;
        if (mag[k] > 0.3) {
            printf("  %6.2f Hz  |  magnitude: %.3f\n", freq_hz, mag[k]);
        }
    }

    return 0;
}

Frequency bins above threshold (magnitude > 0.3):
--------------------------------------------------
   15.62 Hz  |  magnitude: 0.981
   39.06 Hz  |  magnitude: 0.497

The detected frequencies are 15.62 Hz and 39.06 Hz rather than exactly 15 and 40 — this is the bin resolution at work: $1000 / 256 \approx 3.9$ Hz per bin, so neither target frequency falls exactly on a bin center. The magnitudes come out close to 1.0 and 0.5 (the original amplitudes), with slight deviation from noise and the bin-center mismatch. “Close” is the honest word here, not “perfect” — the noise does affect the output, and the frequency displacement is a real measurement artifact, not an implementation bug.

The FFT: Exploiting Symmetry

This $O(N^2)$ implementation is correct but slow. For $N = 1{,}048{,}576$ samples (one megasample buffer), the naive DFT requires about $10^{12}$ operations. At 1 GHz that’s 1000 seconds per buffer — useless for real-time processing.

The Fast Fourier Transform (FFT), formalized by Cooley and Tukey in 1965, reduces this to $O(N \log N)$ by exploiting a structural property of the DFT matrix: its twiddle factors ($e^{-i 2\pi kn / N}$) are periodic and symmetric. The most common variant, decimation-in-time (DIT), splits the $N$-point DFT recursively into two $N/2$-point DFTs — one over the even-indexed samples and one over the odd-indexed samples. Each level of recursion halves the problem size and contributes one $\log_2 N$ factor to the total work:

$X_k = \underbrace{\sum_{n \text{ even}} x_n \, e^{-i2\pi kn/N}}_{\text{N/2-point DFT}} + e^{-i2\pi k/N} \underbrace{\sum_{n \text{ odd}} x_n \, e^{-i2\pi kn/N}}_{\text{N/2-point DFT}}$

The $e^{-i2\pi k/N}$ term in front of the odd sum is the twiddle factor. The “butterfly” operation at each stage combines two values using one complex multiply and two additions, and the entire computation reorganizes into $\log_2 N$ stages of $N/2$ butterflies each.

For the same megasample buffer: $N \log_2 N \approx 2 \times 10^7$ operations. At 1 GHz that’s 20 milliseconds — real-time capable. That gap between $O(N^2)$ and $O(N \log N)$ is why FFT hardware exists in every modern ADC, SDR, and audio DSP chip.

The mathematics is identical to the DFT above. The FFT is purely an algorithmic optimization — a smarter order of operations that avoids redundant computation by reusing intermediate results. Understanding the DFT first means the FFT is just an implementation detail, not a separate concept.