How to Perform Spectral Analysis and Filtering with NumPy

Introduction to Spectral Analysis

Spectral analysis is a technique used in signal processing where signals are analyzed in terms of their frequency content. It involves the decomposition of a time-series signal into its constituent frequencies, which can be visualized using a spectrum. One common way to perform spectral analysis is by using the Fast Fourier Transform (FFT), which efficiently computes the Discrete Fourier Transform (DFT) of a sequence.

In this tutorial, we’ll explore the basics of spectral analysis and filtering using Python’s NumPy library, a powerful package for numerical computing.

Getting Started with NumPy

Before diving into spectral analysis, ensure that you have NumPy installed in your environment. You can install NumPy using pip:

pip install numpy

Next, we import NumPy and other necessary libraries:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter

Generating a Sample Signal

To demonstrate spectral analysis, let’s first generate a synthetic signal composed of multiple sine waves:

fs = 500  # Sampling frequency
T = 1/fs  # Sampling interval
t = np.arange(0, 1, T)  # Time vector

# Signal components
f1, f2 = 5, 45  # Frequencies of the sine waves
A1, A2 = 1, 0.5  # Amplitudes of the sine waves

signal = A1 * np.sin(2 * np.pi * f1 * t) + A2 * np.sin(2 * np.pi * f2 * t)

plt.plot(t, signal)
plt.title('Generated Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.show()

Performing the FFT

The next step is to perform the FFT on our signal:

n = len(signal)  # Length of the signal
F = np.fft.fft(signal)  # FFT of the signal

# Frequency vector
freqs = np.fft.fftfreq(n, T)

magnitude = np.abs(F)  # Magnitude of the FFT

plt.plot(freqs, magnitude)
plt.title('Magnitude Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim([0, 250])  # Only show first half of frequencies
plt.show()

Applying a Low-Pass Filter

Sometimes we are only interested in low frequencies and want to remove high-frequency noise. Here’s how you can apply a low-pass filter using the Butterworth filter design from the SciPy library:

def low_pass_filter(signal, cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    filtered_signal = lfilter(b, a, signal)
    return filtered_signal

# Apply the filter
cutoff = 10  # Cutoff frequency
filtered_signal = low_pass_filter(signal, cutoff, fs)

plt.plot(t, filtered_signal)
plt.title('Filtered Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.show()

Inspecting the filtered signal spectrum

Now let’s see the effect of filtering on the signal’s spectrum:

F_filtered = np.fft.fft(filtered_signal)
magnitude_filtered = np.abs(F_filtered)

plt.plot(freqs, magnitude_filtered)
plt.title('Magnitude Spectrum of Filtered Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim([0, 250])
plt.show()

Conclusion

The spectral analysis and filtering of signals is crucial in various fields of engineering and science. Using NumPy and related Python libraries provides a powerful and accessible means to perform these tasks. From generating and analyzing signals to applying sophisticated filters, the flexibility of these tools enables researchers and engineers to effectively manipulate and interpret signal data.