Computing Power Spectral Density with FFT in Python: From Theory to Practice

Introduction and Problem Context

In signal processing, power spectral density analysis is a crucial technique for revealing the energy distribution of signals in the frequency domain. This article is based on a practical case: a user collected video data with 301 frames at a sampling rate of 30 frames per second, totaling 10 seconds in duration. The user aims to compute the PSD of this signal and display the frequency axis correctly in Hertz (Hz). Initial attempts involved directly applying FFT to obtain magnitude spectra but failed to set the frequency axis properly, making results difficult to interpret.

Theoretical Foundation: FFT and Power Spectral Density

Fast Fourier Transform (FFT) is an efficient algorithm for the Discrete Fourier Transform (DFT), converting time-domain signals to frequency-domain representations. For a discrete signal x[n] of length N, its FFT result X[k] contains complex components, where k=0,1,...,N-1. Power spectral density (PSD) is defined as the energy distribution of a signal in the frequency domain, typically obtained by squaring the magnitude of the FFT result: PSD[k] = |X[k]|². This calculation removes phase information, emphasizing the energy contributions of each frequency component.

In Python, the NumPy library offers robust FFT support. The core function np.fft.fft() computes the FFT, while np.fft.fftfreq() generates the corresponding frequency array. The latter takes two parameters: the number of data points N and the sampling time interval time_step. The frequency array is calculated as: f[k] = k / (N * time_step), where k ranges from 0 to N-1, adjusted to include positive and negative frequencies.

Practical Implementation: Code Analysis and Optimization

Referring to the best answer, we first import necessary libraries and generate sample data. Assume data is an array of 301 values representing the signal extracted from video. The sampling frequency is 30 Hz, so the sampling time interval time_step = 1 / 30 seconds.

import numpy as np
import matplotlib.pyplot as plt

data = np.random.rand(301) - 0.5  # Sample data; replace with real signal in practice
time_step = 1 / 30  # Sampling time interval in seconds

Next, compute the FFT and PSD. Use np.fft.fft(data) for transformation, then calculate the squared magnitude:

fft_result = np.fft.fft(data)
psd = np.abs(fft_result)**2  # Power spectral density

To correctly set the frequency axis, call np.fft.fftfreq(data.size, time_step) to generate the frequency array. This function returns an array of the same length as fft_result, containing corresponding frequency values in Hz. Since the FFT result includes positive and negative frequencies, the frequency array is ordered accordingly, with negative frequencies in the latter half.

freqs = np.fft.fftfreq(data.size, time_step)

For visualization, sort the frequencies and PSD. Use np.argsort(freqs) to get sorting indices, then plot:

idx = np.argsort(freqs)
plt.plot(freqs[idx], psd[idx])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density')
plt.title('Power Spectrum of Signal')
plt.grid(True)
plt.show()

This code produces a power spectrum plot with Hz units. Note that the maximum visible frequency is not the sampling frequency of 30 Hz, but the Nyquist frequency. For 301 data points (odd length), the maximum frequency is calculated as: f_max = (N-1)/(2 * N * time_step) ≈ 14.95 Hz. This aligns with the sampling theorem, where the highest resolvable frequency is half the sampling frequency.

In-Depth Analysis: Nyquist Frequency and Spectral Symmetry

In FFT analysis, the Nyquist frequency (f_Nyquist = sampling frequency / 2) is a key concept. It represents the highest frequency that can be resolved without aliasing. For real-valued signals, the FFT result exhibits conjugate symmetry, i.e., X[k] = X*[N-k], where * denotes complex conjugate. This implies that the PSD is symmetric between positive and negative frequencies, so typically only the positive frequency portion (0 to f_Nyquist) is plotted.

In the user's case, with a sampling frequency of 30 Hz, the Nyquist frequency is 15 Hz. Since the number of data points is 301 (odd), the maximum positive frequency in the FFT frequency array is slightly below 15 Hz, as shown in the code output: max(freqs) ≈ 14.95 Hz. If the data point count were even, the maximum frequency would exactly equal the Nyquist frequency (e.g., for 300 points, f_max = 15 Hz).

Supplementary Methods: Using rfft and Welch Method

Beyond standard FFT, NumPy provides the np.fft.rfft() function, specifically for real-valued input signals. It returns only the positive frequency portion (including zero frequency), reducing computational and storage requirements. Correspondingly, np.fft.rfftfreq() generates the positive frequency array. Example code:

rfft_result = np.fft.rfft(data)
rpsd = np.abs(rfft_result)**2
rfreqs = np.fft.rfftfreq(data.size, time_step)
plt.plot(rfreqs, rpsd)
plt.show()

Another advanced approach is using SciPy's signal.welch function, which estimates PSD via Welch's method. This technique reduces variance by segmenting data and averaging, making it suitable for noisy signals. Example:

from scipy import signal
f, Pxx = signal.welch(data, fs=1/time_step, nperseg=256)
plt.semilogy(f, Pxx)
plt.show()

Welch's method offers smoother spectral estimates but requires tuning parameters like segment length and overlap to balance resolution and variance.

Common Issues and Debugging Tips

In practice, users may encounter issues such as incorrect frequency axis display, abnormal power values, or unclear plots. Solutions include: ensuring time_step is set correctly (e.g., for 30 fps, time_step = 1/30); checking for DC offset in data (preprocess by subtracting the mean); using logarithmic scales (e.g., plt.semilogy()) to better display dynamic range.

For the user's specific case, it is advisable to verify input data quality. If the signal originates from video frames, consider effects like motion blur or compression artifacts. Additionally, 301 data points provide a frequency resolution of approximately 0.033 Hz (Δf = 1/(N * time_step)), which may suffice for analyzing low-frequency components but could be limited for high-frequency details.

Conclusion

This article details methods for computing power spectral density using FFT in Python, focusing on correctly setting frequency axes. By combining NumPy's FFT functions with matplotlib visualization, users can effectively analyze signal frequency-domain characteristics. Key steps include: computing FFT, generating frequency arrays, sorting, and plotting. Understanding the Nyquist frequency and spectral symmetry is essential for interpreting results. For advanced applications, consider using rfft for optimization or Welch's method for improved estimation stability. These techniques provide powerful tools for fields such as signal processing, audio analysis, and video processing.