Fourier Analysis with Python

Published

February 2026

This page is a reference for Fourier analysis in PHYS 4430: what the Fourier Transform computes, how to use it in Python, and how your acquisition parameters affect the spectrum you see. All examples use synthetic data and require only NumPy and Matplotlib—no hardware needed.

Use this page as a reference for topics like windowing and aliasing as they come up during your analysis.

For a quick-reference function and parameter table, see the FFT section of Data Analysis.

1 What Is the Discrete Fourier Transform?

The discrete Fourier Transform (DFT) of a set of data \(\{y_0, y_1, \ldots, y_{N-1}\}\) is given by

\[Y_m = \sum_{n=0}^{N-1} y_n \cdot e^{-2\pi i \frac{m}{N} n}\]

The basic idea: the DFT decomposes your data into a set of different frequency components. The amplitude of \(Y_m\) tells you how much of your signal was formed by an oscillation at the \(m\)-th frequency. If your signal is a pure 100 Hz sine wave, the DFT will have a large value at the bin corresponding to 100 Hz and small values everywhere else.

Key facts:

  • The DFT output \(Y_m\) is complex-valued (it encodes both amplitude and phase). To see how much power is at each frequency, compute \(|Y_m|^2\).
  • The DFT of \(N\) data points produces \(N\) complex values.
  • For real-valued input (which is always the case for experimental data), the DFT is symmetric: \(|Y_m| = |Y_{N-m}|\). This means the second half of the output is redundant, so we typically only plot the first half up to the Nyquist frequency.

1.1 Conceptual Questions

Before computing any FFTs, think through these questions. They will help you interpret your results correctly.

  1. How do the units of the Fourier Transform array \(Y_m\) relate to the units of the data \(y_n\)?
  2. Does the data \(y_n\) have to be taken at equally spaced intervals for the DFT formula above to apply?
  3. Is it possible for two different sets of data to have the same Fourier Transform?
  4. If a data set has \(N\) elements, how many elements does its DFT have?

2 Computing the Power Spectrum in Python

NumPy provides efficient FFT (Fast Fourier Transform) functions for spectral analysis. For real-valued signals—which is what you always have in the lab—use np.fft.rfft() and np.fft.rfftfreq(). These automatically return only the positive-frequency half of the spectrum.

import numpy as np

def compute_power_spectrum(signal, sample_rate):
    """
    Compute the one-sided power spectrum of a real-valued signal.

    Parameters:
        signal: 1D array of signal values
        sample_rate: Sampling rate in Hz

    Returns:
        frequencies: Array of positive frequency values (Hz)
        power: Power spectrum (magnitude squared, normalized)
    """
    n = len(signal)
    frequencies = np.fft.rfftfreq(n, d=1 / sample_rate)
    fft_result = np.fft.rfft(signal)

    # Power spectrum: magnitude squared, normalized by number of points
    power = (np.abs(fft_result) / n) ** 2

    # Double the power for frequencies that have both +/- contributions
    # (all except DC at index 0 and Nyquist at index -1)
    power[1:-1] *= 2

    return frequencies, power

2.1 Basic example

import numpy as np
import matplotlib.pyplot as plt

# Generate a test signal: 50 Hz + 120 Hz + noise
sample_rate = 1000  # Hz
duration = 1.0      # seconds
t = np.arange(0, duration, 1 / sample_rate)

signal = (np.sin(2 * np.pi * 50 * t)
          + 0.5 * np.sin(2 * np.pi * 120 * t)
          + 0.2 * np.random.randn(len(t)))

# Compute power spectrum
frequencies, power = compute_power_spectrum(signal, sample_rate)

# Plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))

ax1.plot(t * 1000, signal, 'b-', linewidth=0.5)
ax1.set_xlabel('Time (ms)')
ax1.set_ylabel('Amplitude')
ax1.set_title('Time Domain')
ax1.set_xlim(0, 100)
ax1.grid(True, alpha=0.3)

ax2.plot(frequencies, power, 'r-')
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Power')
ax2.set_title('Power Spectrum')
ax2.set_xlim(0, 200)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

You should see two clear peaks: a large one at 50 Hz and a smaller one at 120 Hz, sitting on a low noise floor.

2.2 rfft vs. fft

The functions above use np.fft.rfft and np.fft.rfftfreq, which are designed for real-valued input. NumPy also provides np.fft.fft and np.fft.fftfreq for the general case (complex-valued signals). For real signals, rfft is simpler—it returns only the positive-frequency half directly, so you don’t have to slice the output manually. Both give the same result for the positive frequencies.

3 Key Relationships: Resolution and Nyquist

Two parameters control what you can see in the spectrum:

Parameter Formula Meaning
Frequency resolution \(\Delta f = f_s / N = 1/T\) Smallest frequency difference you can distinguish
Maximum frequency (Nyquist) \(f_N = f_s / 2\) Highest frequency you can measure
Total duration \(T = N / f_s\) Longer duration = finer resolution

where \(f_s\) is the sample rate, \(N\) is the number of samples, and \(T\) is the total acquisition duration.

# Example: compute frequency resolution and Nyquist frequency
sample_rate = 1000     # Hz
num_samples = 2000

freq_resolution = sample_rate / num_samples  # = 1 / duration
max_frequency = sample_rate / 2              # Nyquist frequency

print(f"Duration: {num_samples / sample_rate} s")
print(f"Frequency resolution: {freq_resolution} Hz")
print(f"Maximum frequency: {max_frequency} Hz")
# Duration: 2.0 s
# Frequency resolution: 0.5 Hz
# Maximum frequency: 500.0 Hz

A common misconception: Increasing the sample rate does not improve frequency resolution. What actually controls resolution is the total acquisition duration \(T = N/f_s\). Increasing the sample rate raises the Nyquist frequency (letting you see higher frequencies), but to get finer frequency resolution you need to record for a longer time.

4 Windowing and Spectral Leakage

The DFT assumes your signal is periodic—that the first and last samples connect smoothly. In practice, your signal has finite duration and the edges usually don’t match up. This mismatch creates artifacts called spectral leakage: the FFT spreads power from a narrow peak into broad “skirts” that extend across many frequency bins.

4.1 Why it matters

Spectral leakage makes peaks appear broader than they really are and raises the noise floor around strong peaks. In Fourier Transform spectroscopy, leakage directly degrades your spectral resolution.

4.2 The fix: window functions

Multiplying your signal by a window function that tapers smoothly to zero at both ends eliminates the edge discontinuity. The most common choice is the Hanning window (np.hanning). The following example demonstrates the problem and the fix. The signal contains a strong 47 Hz component and a weak 58 Hz component (only 8% of the strong signal’s amplitude). Without windowing, the leakage from the strong peak buries the weak one; with windowing, both peaks are clearly visible.

import numpy as np
import matplotlib.pyplot as plt

# Two-frequency signal: strong at 47 Hz, weak at 58 Hz
# The weak signal is only 8% of the strong one — can you find it in the spectrum?
sample_rate = 1000
t = np.arange(0, 0.5, 1 / sample_rate)  # 0.5 seconds (resolution = 2 Hz)
signal = np.sin(2 * np.pi * 47 * t) + 0.08 * np.sin(2 * np.pi * 58 * t)

# Apply Hanning window
window = np.hanning(len(signal))
signal_windowed = signal * window

# Compute power spectra
freqs = np.fft.rfftfreq(len(signal), d=1 / sample_rate)
power_raw = np.abs(np.fft.rfft(signal)) ** 2
power_windowed = np.abs(np.fft.rfft(signal_windowed)) ** 2

# 2x2 comparison: time domain (top) and frequency domain (bottom)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# Top left: raw signal (note abrupt start and end)
axes[0, 0].plot(t * 1000, signal, 'b-', linewidth=0.5)
axes[0, 0].set_xlabel('Time (ms)')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Signal (no window)')
axes[0, 0].grid(True, alpha=0.3)

# Top right: windowed signal (tapers to zero at edges)
axes[0, 1].plot(t * 1000, window, 'r-', alpha=0.4, linewidth=2,
                label='Hanning window')
axes[0, 1].plot(t * 1000, signal_windowed, 'b-', linewidth=0.5,
                label='Windowed signal')
axes[0, 1].set_xlabel('Time (ms)')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_title('Signal × Hanning window')
axes[0, 1].legend(loc='upper right')
axes[0, 1].grid(True, alpha=0.3)

# Bottom row: zoom into the region around both peaks
axes[1, 0].plot(freqs, power_raw)
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Power')
axes[1, 0].set_title('Spectrum (no window)')
axes[1, 0].set_xlim(30, 75)
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(freqs, power_windowed)
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('Power')
axes[1, 1].set_title('Spectrum (Hanning window)')
axes[1, 1].set_xlim(30, 75)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Top row (inputs): The raw signal (left) starts and ends at arbitrary values — when the DFT wraps this signal to repeat it, those edges create discontinuities. The windowed signal (right) tapers smoothly to zero at both ends, so the periodicity assumption is satisfied. The red curve shows the shape of the Hanning window.

Bottom row (spectra): Without windowing (left), the leakage skirts from the strong 47 Hz peak spread across the spectrum and bury the weak 58 Hz signal. With windowing (right), the leakage is suppressed and the weak peak emerges clearly. This is the practical payoff: windowing lets you detect weak signals near strong ones. The window does reduce overall peak height (it attenuates the signal near the edges), but the improvement in peak isolation is well worth the tradeoff.

Tip

When to window: Apply a window function whenever your signal doesn’t naturally start and end at zero. This is almost always the case for experimental data.

5 Aliasing

If your signal contains frequencies above the Nyquist frequency (\(f_s/2\)), those frequencies “fold back” and appear at incorrect positions in the spectrum. This is called aliasing, and it cannot be undone after the fact.

Example: A 400 Hz signal sampled at 500 Hz appears at 100 Hz in the spectrum (it folds around the 250 Hz Nyquist frequency).

Rule of thumb: Sample at least 2x your highest frequency of interest (preferably 5-10x for clean spectra).

For an interactive aliasing demonstration with code, see the Data Analysis page.

Note

Debugging unexpected peaks: If your FFT shows a peak you didn’t expect (e.g., at 120 Hz), consider physical sources: power line interference (60 Hz and harmonics at 120, 180 Hz), fluorescent light flicker (100-120 Hz), or mechanical vibrations. To identify the source, try shielding the detector, turning off lights, or isolating the optics from the table.

These exercises let you verify your understanding using only Python. No hardware is needed.

5.1 Exercise 1: Single frequency

Generate a 100 Hz sine wave sampled at 1000 Hz for 1 second. Compute and plot its power spectrum.

sample_rate = 1000
duration = 1.0
t = np.arange(0, duration, 1 / sample_rate)
signal = np.sin(2 * np.pi * 100 * t)

freqs, power = compute_power_spectrum(signal, sample_rate)

plt.plot(freqs, power)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.xlim(0, 200)
plt.grid(True, alpha=0.3)
plt.show()

Verify:

  • The peak appears at exactly 100 Hz.
  • The frequency resolution is \(\Delta f = 1000 / 1000 = 1\) Hz.
  • The Nyquist frequency is \(f_N = 1000 / 2 = 500\) Hz.

5.2 Exercise 2: Two frequencies and resolution

Generate a signal containing 95 Hz and 105 Hz components (10 Hz apart). Try two different durations:

a) Duration = 1.0 s (frequency resolution = 1 Hz). Can you resolve the two peaks?

b) Duration = 0.05 s (frequency resolution = 20 Hz). Can you still resolve them? Why or why not?

sample_rate = 1000

for duration in [1.0, 0.05]:
    t = np.arange(0, duration, 1 / sample_rate)
    signal = np.sin(2 * np.pi * 95 * t) + np.sin(2 * np.pi * 105 * t)

    freqs, power = compute_power_spectrum(signal, sample_rate)

    plt.figure(figsize=(8, 3))
    plt.plot(freqs, power)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power')
    plt.title(f'Duration = {duration} s (resolution = {sample_rate / len(t):.0f} Hz)')
    plt.xlim(50, 150)
    plt.grid(True, alpha=0.3)
    plt.show()

This exercise demonstrates concretely why longer recordings give better frequency resolution.

5.3 Exercise 3: Aliasing

Generate a 400 Hz sine wave sampled at only 500 Hz. Before running the code, predict where the peak will appear. Then verify.

# Undersampled: 400 Hz signal at 500 Hz sample rate
sample_rate = 500
t = np.arange(0, 1.0, 1 / sample_rate)
signal = np.sin(2 * np.pi * 400 * t)

freqs, power = compute_power_spectrum(signal, sample_rate)

plt.plot(freqs, power)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.title(f'400 Hz signal sampled at {sample_rate} Hz')
plt.grid(True, alpha=0.3)
plt.show()

Now increase the sample rate to 1000 Hz and confirm the peak appears at the correct 400 Hz.

Back to Python Resources