⬡ Hub
Skip to content

SciPy: Signal Processing

The scipy.signal module provides a comprehensive set of tools for processing analog and digital signals. This includes filtering, convolution, spectral analysis, windowing, and various transformations. It's an essential module for anyone working with time-series data, audio, images, or other forms of signals.

1. Filtering

Filtering is a fundamental operation to remove unwanted components (like noise) or extract desired features from a signal.

a. Basic Filtering (e.g., Savitzky-Golay Filter)

A Savitzky-Golay filter smooths a signal by fitting a local polynomial regression to each data point and replacing the data point with the value of the fitting polynomial.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

# Generate noisy signal
np.random.seed(42)
x = np.linspace(0, 2 * np.pi, 100)
y_true = np.sin(x) + np.cos(2*x)
y_noisy = y_true + np.random.normal(0, 0.5, len(x))

# Apply Savitzky-Golay filter
# window_length: must be odd
# polyorder: order of the polynomial to fit
y_filtered = savgol_filter(y_noisy, window_length=15, polyorder=3)

plt.figure(figsize=(10, 6))
plt.plot(x, y_noisy, label='Noisy Signal', alpha=0.6)
plt.plot(x, y_true, label='True Signal', linestyle='--', color='green')
plt.plot(x, y_filtered, label='Filtered Signal (Savitzky-Golay)', color='red', linewidth=2)
plt.title('Signal Filtering with Savitzky-Golay')
plt.xlabel('Time/Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()

b. Design and Apply Digital Filters (e.g., Butterworth Filter)

scipy.signal allows you to design various types of digital filters (e.g., Butterworth, Chebyshev) and apply them to your data.

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# Generate a signal with two frequencies and noise
fs = 1000 # Sampling frequency
t = np.linspace(0, 1, fs, endpoint=False) # 1 second duration
f1 = 5 # Hz
f2 = 50 # Hz
y = 0.6 * np.sin(2 * np.pi * f1 * t) + 1.0 * np.sin(2 * np.pi * f2 * t)
y_noisy = y + 0.5 * np.random.randn(len(t)) # Add some random noise

# Design a Butterworth low-pass filter
# N: Order of the filter (higher order -> sharper cutoff)
# Wn: Critical frequency (normalized by Nyquist frequency, fs/2)
# b, a are the numerator (b) and denominator (a) polynomials of the IIR filter.
b, a = signal.butter(N=5, Wn=30 / (fs / 2), btype='low', analog=False)

# Apply the filter
y_filtered_butter = signal.lfilter(b, a, y_noisy)

plt.figure(figsize=(12, 6))
plt.plot(t, y_noisy, label='Noisy Signal', alpha=0.7)
plt.plot(t, y, label='Original Signal (without noise)', linestyle='--', color='green')
plt.plot(t, y_filtered_butter, label='Filtered Signal (Butterworth Low-pass)', color='red', linewidth=2)
plt.title('Butterworth Low-pass Filter Application')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.xlim(0, 0.5) # Focus on a part of the signal
plt.show()

2. Convolution

Convolution is a mathematical operation that is fundamental to many signal processing tasks, including filtering, smoothing, and feature extraction.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve

# Generate two 1D signals
signal1 = np.array([1, 2, 3, 4, 5])
kernel = np.array([0.5, 1, 0.5]) # A simple smoothing kernel

# Perform convolution
convolved_signal = convolve(signal1, kernel, mode='valid') # 'valid' returns only parts where kernel fully overlaps
convolved_signal_full = convolve(signal1, kernel, mode='full') # 'full' returns convolution at each point of overlap
convolved_signal_same = convolve(signal1, kernel, mode='same') # 'same' returns output with length equal to signal1

print("Signal 1:", signal1)
print("Kernel:", kernel)
print("Convolved Signal (mode='valid'):", convolved_signal)
print("Convolved Signal (mode='full'):", convolved_signal_full)
print("Convolved Signal (mode='same'):", convolved_signal_same)

plt.figure(figsize=(10, 4))
plt.plot(signal1, 'o-', label='Original Signal')
plt.plot(convolved_signal_same, 'x-', label='Convolved Signal (same mode)')
plt.title('Convolution of 1D Signals')
plt.xlabel('Index')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.show()

3. Spectral Analysis (Fourier Transforms)

Understanding the frequency components of a signal is crucial. scipy.fft (or scipy.fftpack for older versions) provides functions for Fast Fourier Transforms (FFT).

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

# Generate a signal with two distinct frequencies
fs = 100 # Sampling frequency
t = np.linspace(0, 1, fs, endpoint=False) # 1 second duration
f1 = 10 # Hz
f2 = 25 # Hz
amplitude1 = 0.7
amplitude2 = 1.0
signal = amplitude1 * np.sin(2 * np.pi * f1 * t) + amplitude2 * np.sin(2 * np.pi * f2 * t)

# Perform FFT
yf = fft(signal)
xf = fftfreq(fs, 1 / fs) # Frequencies corresponding to FFT output

# Plot the signal in time domain
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t, signal)
plt.title('Time Domain Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

# Plot the magnitude of the FFT (only positive frequencies)
plt.subplot(1, 2, 2)
# Take absolute value and normalize by N/2 (for single-sided spectrum)
plt.plot(xf[:fs//2], 2.0/fs * np.abs(yf[0:fs//2]))
plt.title('Frequency Domain Signal (FFT Magnitude)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid(True)
plt.show()

4. Windowing Functions

Windowing functions are used in spectral analysis to reduce spectral leakage when performing Fourier transforms on finite-length signals.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import windows

# Generate a rectangular window and a Hamming window
window_len = 51
rect_window = windows.boxcar(window_len) # Rectangular window
ham_window = windows.hamming(window_len) # Hamming window

plt.figure(figsize=(10, 5))
plt.plot(rect_window, label='Rectangular Window')
plt.plot(ham_window, label='Hamming Window')
plt.title('Different Window Functions')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()

# Applying a window to a signal before FFT
# signal_windowed = signal * windows.hann(len(signal))

Further Topics:

  • Correlation (correlate)
  • Discrete Wavelet Transform (dwt, idwt)
  • System identification (lsim, impulse)
  • Resampling (resample)
  • Peak finding (find_peaks)

SciPy's signal module is an indispensable toolkit for anyone working with signals and systems, offering highly optimized routines for analysis, filtering, and transformation.