How to compute an FFT with SciPy

Sampled sensor, audio, and simulation signals often hide periodic components in time-domain values. SciPy provides scipy.fft routines that convert those samples into frequency bins, making a dominant tone or vibration easier to identify numerically.

The fft() result is a complex spectrum, and fftfreq() builds the matching frequency axis from the sample count and sample spacing. Pairing those arrays keeps each magnitude tied to the frequency bin that produced it instead of treating spectrum indexes as hertz values.

A one-second signal sampled at 200 Hz gives a 1 Hz bin spacing, so 30 Hz and 75 Hz tones fall directly on FFT bins. Real measurement data may need windowing, detrending, or a longer sample window when a component falls between bins, because energy can spread into nearby bins.

1. Create a Python script named fft_compute.py with a signal that contains known frequency components.

   fft_compute.py

   import numpy as np
   from scipy.fft import fft, fftfreq
    
    
   sample_rate = 200.0
   duration = 1.0
   sample_count = int(sample_rate * duration)
   time = np.arange(sample_count) / sample_rate
    
   signal = (
       1.2 * np.sin(2.0 * np.pi * 30.0 * time)
       + 0.4 * np.sin(2.0 * np.pi * 75.0 * time)
   )
    
   spectrum = fft(signal)
   frequencies = fftfreq(sample_count, d=1.0 / sample_rate)
    
   positive = frequencies >= 0
   positive_frequencies = frequencies[positive]
   amplitudes = (2.0 / sample_count) * np.abs(spectrum[positive])
    
   peak_indexes = np.flatnonzero(amplitudes > 0.1)
   detected_peaks = positive_frequencies[peak_indexes]
   expected_peaks = np.array([30.0, 75.0])
    
   print(f"sample_count: {sample_count}")
   print(f"bin_spacing_hz: {sample_rate / sample_count:.1f}")
   print("detected peaks:")
   for index in peak_indexes:
       print(f"  {positive_frequencies[index]:5.1f} Hz  amplitude {amplitudes[index]:.3f}")
   print(f"dominant_frequency_hz: {positive_frequencies[np.argmax(amplitudes)]:.1f}")
   print(f"matches_expected_peaks: {np.allclose(detected_peaks, expected_peaks)}")

   The two sine waves complete an integer number of cycles inside the one-second sample, so their energy appears in the 30 Hz and 75 Hz bins.

2. Run the script to compute the spectrum and print the detected peaks.

   $ python fft_compute.py
   sample_count: 200
   bin_spacing_hz: 1.0
   detected peaks:
      30.0 Hz  amplitude 1.200
      75.0 Hz  amplitude 0.400
   dominant_frequency_hz: 30.0
   matches_expected_peaks: True

3. Pair the FFT output with frequency bins from the same sample spacing.

   spectrum = fft(signal)
   frequencies = fftfreq(sample_count, d=1.0 / sample_rate)

   The d value is the spacing between samples. For a 200 Hz sample rate, each sample is 1 / 200 second apart.

4. Keep the nonnegative frequency bins for a one-sided view of the real-valued signal.

   positive = frequencies >= 0
   positive_frequencies = frequencies[positive]

   A full fft() spectrum includes positive and negative frequency bins. For real-valued input, the negative half mirrors the positive half.

5. Scale the one-sided magnitudes before comparing peak amplitudes.

   amplitudes = (2.0 / sample_count) * np.abs(spectrum[positive])

   This scaling reports the sine-wave amplitudes for this zero-mean sample. Treat DC and Nyquist-bin amplitudes separately when those bins matter.

6. Confirm that the detected peaks match the known signal frequencies.

   peak_indexes = np.flatnonzero(amplitudes > 0.1)
   detected_peaks = positive_frequencies[peak_indexes]
   expected_peaks = np.array([30.0, 75.0])
   print(f"matches_expected_peaks: {np.allclose(detected_peaks, expected_peaks)}")

7. Remove the demo script after adapting the calculation.

   $ rm fft_compute.py