Python Tutorial for DSP

1. Getting Started with DSP in Python

Installing Required Libraries

# Install via pip
pip install numpy scipy matplotlib

# Or using conda
conda install numpy scipy matplotlib

Essential Imports

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft, ifft, fftshift, fftfreq
from numpy.polynomial import polynomial as P

# Set plot style
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = [10, 6]
plt.rcParams['font.size'] = 12

Basic DSP Setup

import numpy as np
import matplotlib.pyplot as plt

# Sampling parameters
Fs = 1000           # Sampling frequency (Hz)
T = 1/Fs            # Sampling period
N = 1024            # Number of samples
n = np.arange(N)    # Sample indices
t = n * T           # Time vector

# Create a test signal: two sinusoids
f1, f2 = 50, 120    # Frequencies (Hz)
x = np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t)

# Plot the signal
plt.figure(figsize=(10, 4))
plt.plot(t, x, 'b-', linewidth=1.5)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Test Signal: 50 Hz + 120 Hz Sinusoids')
plt.grid(True)
plt.tight_layout()
plt.show()

2. Discrete-Time Signals

Basic Signal Generation

import numpy as np
import matplotlib.pyplot as plt

# Sample indices
n = np.arange(-10, 11)

# Unit Impulse (Delta Function)
delta = (n == 0).astype(float)

# Unit Step Function
u = (n >= 0).astype(float)

# Exponential Sequence
a = 0.9
x_exp = (a ** n) * (n >= 0)

# Create subplots
fig, axes = plt.subplots(3, 1, figsize=(10, 8))

axes[0].stem(n, delta, basefmt=' ')
axes[0].set_xlabel('n')
axes[0].set_ylabel(r'$\delta[n]$')
axes[0].set_title('Unit Impulse')
axes[0].grid(True)

axes[1].stem(n, u, basefmt=' ', linefmt='r-', markerfmt='ro')
axes[1].set_xlabel('n')
axes[1].set_ylabel('u[n]')
axes[1].set_title('Unit Step')
axes[1].grid(True)

axes[2].stem(n, x_exp, basefmt=' ', linefmt='g-', markerfmt='go')
axes[2].set_xlabel('n')
axes[2].set_ylabel('x[n]')
axes[2].set_title(f'Exponential Sequence: x[n] = ({a})^n u[n]')
axes[2].grid(True)

plt.tight_layout()
plt.show()

Sinusoidal Sequences

import numpy as np
import matplotlib.pyplot as plt

n = np.arange(51)
omega0 = 0.2 * np.pi    # Digital frequency (rad/sample)
phi = 0                  # Phase

x = np.cos(omega0 * n + phi)

fig, axes = plt.subplots(2, 1, figsize=(10, 6))

# Stem plot
axes[0].stem(n, x, basefmt=' ')
axes[0].set_xlabel('n')
axes[0].set_ylabel('x[n]')
axes[0].set_title(f'Discrete Sinusoid: ω₀ = {omega0/np.pi:.2f}π rad/sample')
axes[0].grid(True)

# With continuous envelope
n_fine = np.linspace(0, 50, 500)
x_envelope = np.cos(omega0 * n_fine + phi)

axes[1].plot(n_fine, x_envelope, 'r--', linewidth=1, label='Continuous')
axes[1].stem(n, x, basefmt=' ', label='Discrete')
axes[1].set_xlabel('n')
axes[1].set_ylabel('x[n]')
axes[1].set_title('Discrete Samples of Continuous Sinusoid')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

3. Convolution & Correlation

Linear Convolution

import numpy as np
import matplotlib.pyplot as plt

# Define two sequences
x = np.array([1, 2, 3, 4])      # Input signal
h = np.array([1, 1, 1])          # Impulse response

# Compute convolution
y = np.convolve(x, h)

# Print lengths
print(f"Length of x: {len(x)}")
print(f"Length of h: {len(h)}")
print(f"Length of y: {len(y)}")  # len(x) + len(h) - 1
print(f"y = {y}")

# Plot results
fig, axes = plt.subplots(3, 1, figsize=(10, 8))

axes[0].stem(np.arange(len(x)), x, basefmt=' ')
axes[0].set_xlabel('n')
axes[0].set_ylabel('x[n]')
axes[0].set_title('Input Signal x[n]')
axes[0].grid(True)

axes[1].stem(np.arange(len(h)), h, basefmt=' ', linefmt='r-', markerfmt='ro')
axes[1].set_xlabel('n')
axes[1].set_ylabel('h[n]')
axes[1].set_title('Impulse Response h[n]')
axes[1].grid(True)

axes[2].stem(np.arange(len(y)), y, basefmt=' ', linefmt='g-', markerfmt='go')
axes[2].set_xlabel('n')
axes[2].set_ylabel('y[n]')
axes[2].set_title('Output y[n] = x[n] * h[n] (Convolution)')
axes[2].grid(True)

plt.tight_layout()
plt.show()

Correlation

import numpy as np
import matplotlib.pyplot as plt

x = np.array([1, 2, 3, 4, 3, 2, 1])
y = np.array([0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0])  # Shifted version

# Cross-correlation using numpy
r_xy = np.correlate(x, y, mode='full')
lags = np.arange(-(len(y)-1), len(x))

# Find lag of maximum correlation
max_idx = np.argmax(r_xy)
lag_at_max = lags[max_idx]
print(f"Maximum correlation at lag = {lag_at_max}")

# Autocorrelation
r_xx = np.correlate(x, x, mode='full')
lags_auto = np.arange(-(len(x)-1), len(x))

fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].plot(lags, r_xy, 'b-', linewidth=2)
axes[0].axvline(lag_at_max, color='r', linestyle='--', label=f'Max at lag={lag_at_max}')
axes[0].set_xlabel('Lag')
axes[0].set_ylabel(r'$R_{xy}[k]$')
axes[0].set_title('Cross-Correlation')
axes[0].legend()
axes[0].grid(True)

axes[1].stem(lags_auto, r_xx, basefmt=' ', linefmt='g-', markerfmt='go')
axes[1].set_xlabel('Lag')
axes[1].set_ylabel(r'$R_{xx}[k]$')
axes[1].set_title('Autocorrelation of x[n]')
axes[1].grid(True)

plt.tight_layout()
plt.show()

4. Discrete Fourier Transform (DFT) & FFT

Basic FFT Usage

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

# Sampling parameters
Fs = 1000               # Sampling frequency
N = 1024                # Number of samples
t = np.arange(N) / Fs   # Time vector

# Create signal with two frequencies + noise
f1, f2 = 50, 120
x = np.sin(2*np.pi*f1*t) + 0.6*np.sin(2*np.pi*f2*t)
x_noisy = x + 0.5*np.random.randn(N)

# Compute FFT
X = fft(x_noisy)

# Frequency axis
f = fftfreq(N, 1/Fs)
f_shifted = fftshift(f)

# Magnitude spectrum
X_mag = np.abs(X) / N

# One-sided spectrum (for real signals)
X_onesided = X_mag[:N//2+1].copy()
X_onesided[1:-1] *= 2
f_onesided = f[:N//2+1]

fig, axes = plt.subplots(3, 1, figsize=(10, 9))

axes[0].plot(t, x_noisy, 'b-', linewidth=0.5)
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Noisy Signal in Time Domain')
axes[0].grid(True)

axes[1].plot(f_onesided, X_onesided, 'r-', linewidth=1.5)
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('|X(f)|')
axes[1].set_title('One-Sided Amplitude Spectrum')
axes[1].set_xlim([0, 200])
axes[1].grid(True)

axes[2].plot(f_shifted, fftshift(X_mag), 'g-', linewidth=1.5)
axes[2].set_xlabel('Frequency (Hz)')
axes[2].set_ylabel('|X(f)|')
axes[2].set_title('Two-Sided Amplitude Spectrum (Centered)')
axes[2].set_xlim([-200, 200])
axes[2].grid(True)

plt.tight_layout()
plt.show()

Zero-Padding for Better Resolution

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

# Short sequence
x = np.array([1, 1, 1, 1, 1])
N_original = len(x)
N_padded = 64

# FFT without padding
X1 = fft(x)
f1 = np.arange(N_original) / N_original

# FFT with zero-padding
X2 = fft(x, N_padded)
f2 = np.arange(N_padded) / N_padded

fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].stem(f1, np.abs(X1), basefmt=' ')
axes[0].set_xlabel('ω/2π')
axes[0].set_ylabel('|X[k]|')
axes[0].set_title(f'DFT without Zero-Padding (N = {N_original})')
axes[0].grid(True)

axes[1].stem(f2, np.abs(X2), basefmt=' ', linefmt='r-', markerfmt='ro')
axes[1].set_xlabel('ω/2π')
axes[1].set_ylabel('|X[k]|')
axes[1].set_title(f'DFT with Zero-Padding (N = {N_padded})')
axes[1].grid(True)

plt.tight_layout()
plt.show()

5. z-Transform Analysis

Pole-Zero Analysis

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

# Define system: H(z) = (z^2 - 0.5z) / (z^2 - 1.2z + 0.72)
# Coefficients in descending powers of z
b = [1, -0.5, 0]            # Numerator
a = [1, -1.2, 0.72]         # Denominator

# Find zeros and poles
zeros = np.roots(b)
poles = np.roots(a)

print("Zeros:", zeros)
print("Poles:", poles)
print("Pole magnitudes:", np.abs(poles))

# Check stability
if np.all(np.abs(poles) < 1):
    print("System is STABLE (all poles inside unit circle)")
else:
    print("System is UNSTABLE")

# Pole-Zero Plot
fig, ax = plt.subplots(figsize=(8, 8))

# Draw unit circle
theta = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(theta), np.sin(theta), 'k--', linewidth=1, alpha=0.5)

# Plot zeros (o) and poles (x)
ax.scatter(np.real(zeros), np.imag(zeros), s=100, marker='o', 
           facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(np.real(poles), np.imag(poles), s=100, marker='x', 
           c='red', linewidths=2, label='Poles')

ax.axhline(0, color='k', linewidth=0.5)
ax.axvline(0, color='k', linewidth=0.5)
ax.set_xlabel('Real Part')
ax.set_ylabel('Imaginary Part')
ax.set_title('Pole-Zero Plot')
ax.legend()
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.set_xlim([-2, 2])
ax.set_ylim([-2, 2])

plt.tight_layout()
plt.show()

Frequency Response

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

# Define system
b = [1, 0.5]
a = [1, -0.8, 0.15]

# Compute frequency response
w, H = signal.freqz(b, a, worN=1024)

# Convert to different representations
mag_db = 20 * np.log10(np.abs(H))
phase_deg = np.unwrap(np.angle(H)) * 180 / np.pi

fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].plot(w/np.pi, mag_db, 'b-', linewidth=2)
axes[0].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Frequency Response')
axes[0].grid(True)

axes[1].plot(w/np.pi, phase_deg, 'r-', linewidth=2)
axes[1].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[1].set_ylabel('Phase (degrees)')
axes[1].set_title('Phase Response')
axes[1].grid(True)

plt.tight_layout()
plt.show()

Impulse and Step Response

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

# Define system
b = [1, 0.5]
a = [1, -0.8, 0.15]

# Create discrete-time system
system = signal.dlti(b, a)

# Impulse response
n_imp, h = signal.dimpulse(system, n=50)

# Step response  
n_step, s = signal.dstep(system, n=50)

fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].stem(n_imp, h[0].flatten(), basefmt=' ')
axes[0].set_xlabel('n')
axes[0].set_ylabel('h[n]')
axes[0].set_title('Impulse Response')
axes[0].grid(True)

axes[1].stem(n_step, s[0].flatten(), basefmt=' ', linefmt='r-', markerfmt='ro')
axes[1].set_xlabel('n')
axes[1].set_ylabel('s[n]')
axes[1].set_title('Step Response')
axes[1].grid(True)

plt.tight_layout()
plt.show()

6. Digital Filter Design

FIR Filter Design

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

# Filter specifications
Fs = 1000               # Sampling frequency
Fc = 100                # Cutoff frequency
N = 50                  # Filter order
Wn = Fc / (Fs/2)        # Normalized frequency (0 to 1)

# Design FIR filters with different windows
b_rect = signal.firwin(N+1, Wn, window='boxcar')
b_hamm = signal.firwin(N+1, Wn, window='hamming')
b_black = signal.firwin(N+1, Wn, window='blackman')
b_kaiser = signal.firwin(N+1, Wn, window=('kaiser', 5))

# Compute frequency responses
w, H_rect = signal.freqz(b_rect, 1, worN=1024, fs=Fs)
_, H_hamm = signal.freqz(b_hamm, 1, worN=1024, fs=Fs)
_, H_black = signal.freqz(b_black, 1, worN=1024, fs=Fs)
_, H_kaiser = signal.freqz(b_kaiser, 1, worN=1024, fs=Fs)

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# Magnitude response comparison
axes[0].plot(w, 20*np.log10(np.abs(H_rect)), label='Rectangular')
axes[0].plot(w, 20*np.log10(np.abs(H_hamm)), label='Hamming')
axes[0].plot(w, 20*np.log10(np.abs(H_black)), label='Blackman')
axes[0].plot(w, 20*np.log10(np.abs(H_kaiser)), label='Kaiser')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('FIR Lowpass Filter - Window Comparison')
axes[0].legend()
axes[0].set_xlim([0, Fs/2])
axes[0].set_ylim([-100, 5])
axes[0].grid(True)

# Filter coefficients
axes[1].stem(np.arange(N+1), b_hamm, basefmt=' ')
axes[1].set_xlabel('n')
axes[1].set_ylabel('h[n]')
axes[1].set_title('Filter Coefficients (Hamming Window)')
axes[1].grid(True)

plt.tight_layout()
plt.show()

IIR Filter Design

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

# Filter specifications
Fs = 1000
Fc = 100
Wn = Fc / (Fs/2)
order = 4

# Design IIR filters
b_butter, a_butter = signal.butter(order, Wn)
b_cheby1, a_cheby1 = signal.cheby1(order, 0.5, Wn)  # 0.5 dB ripple
b_cheby2, a_cheby2 = signal.cheby2(order, 40, Wn)   # 40 dB stopband
b_ellip, a_ellip = signal.ellip(order, 1, 40, Wn)   # 1 dB pass, 40 dB stop

# Compute frequency responses
w, H_butter = signal.freqz(b_butter, a_butter, worN=1024, fs=Fs)
_, H_cheby1 = signal.freqz(b_cheby1, a_cheby1, worN=1024, fs=Fs)
_, H_cheby2 = signal.freqz(b_cheby2, a_cheby2, worN=1024, fs=Fs)
_, H_ellip = signal.freqz(b_ellip, a_ellip, worN=1024, fs=Fs)

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(w, 20*np.log10(np.abs(H_butter)), linewidth=2, label='Butterworth')
axes[0].plot(w, 20*np.log10(np.abs(H_cheby1)), linewidth=2, label='Chebyshev I')
axes[0].plot(w, 20*np.log10(np.abs(H_cheby2)), linewidth=2, label='Chebyshev II')
axes[0].plot(w, 20*np.log10(np.abs(H_ellip)), linewidth=2, label='Elliptic')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('IIR Lowpass Filters Comparison (Order 4)')
axes[0].legend()
axes[0].set_xlim([0, Fs/2])
axes[0].set_ylim([-80, 5])
axes[0].grid(True)

# Pole-Zero plot for Butterworth
zeros_b = np.roots(b_butter)
poles_b = np.roots(a_butter)

theta = np.linspace(0, 2*np.pi, 100)
axes[1].plot(np.cos(theta), np.sin(theta), 'k--', linewidth=1, alpha=0.5)
axes[1].scatter(np.real(zeros_b), np.imag(zeros_b), s=100, marker='o', 
               facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
axes[1].scatter(np.real(poles_b), np.imag(poles_b), s=100, marker='x', 
               c='red', linewidths=2, label='Poles')
axes[1].set_xlabel('Real Part')
axes[1].set_ylabel('Imaginary Part')
axes[1].set_title('Butterworth Filter Pole-Zero Plot')
axes[1].legend()
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Applying Filters

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

# Create noisy signal
Fs = 1000
t = np.arange(0, 1, 1/Fs)
signal_clean = np.sin(2*np.pi*50*t)
noise = 0.5*np.sin(2*np.pi*300*t) + 0.3*np.random.randn(len(t))
x = signal_clean + noise

# Design lowpass filter
b, a = signal.butter(6, 100/(Fs/2))

# Apply filter
y_filter = signal.lfilter(b, a, x)           # Causal filtering
y_filtfilt = signal.filtfilt(b, a, x)        # Zero-phase filtering

fig, axes = plt.subplots(4, 1, figsize=(10, 10))

axes[0].plot(t, signal_clean, 'g-', linewidth=1.5)
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Original Signal (50 Hz)')
axes[0].grid(True)

axes[1].plot(t, x, 'b-', linewidth=0.5)
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Noisy Signal')
axes[1].grid(True)

axes[2].plot(t, y_filter, 'r-', linewidth=1.5)
axes[2].set_ylabel('Amplitude')
axes[2].set_title('After lfilter() - Note Phase Delay')
axes[2].grid(True)

axes[3].plot(t, y_filtfilt, 'm-', linewidth=1.5)
axes[3].set_xlabel('Time (s)')
axes[3].set_ylabel('Amplitude')
axes[3].set_title('After filtfilt() - Zero-Phase Filtering')
axes[3].grid(True)

plt.tight_layout()
plt.show()

7. MATLAB to Python Equivalents

Common DSP Functions

Array Indexing Differences

% MATLAB                    # Python
x(1)                        x[0]
x(end)                      x[-1]
x(2:5)                      x[1:5]  # Note: Python excludes end
x(1:2:end)                  x[::2]

8. Tips & Best Practices

Useful Libraries Quick Reference

Need Help?

Contact Dr. Al Bataineh during office hours (Tuesday/Thursday 11:00 AM - 12:00 PM)

Email: mffbataineh@uaeu.ac.ae | Office: F1-1175

← Back to ECOM 451 Course Page