Python Lab Reference

Every library, method and pattern you need for the Signals and Systems labs.

88 entries across 8 libraries

🔢NumPy

A sampled signal is a NumPy array: a time axis, the values on it, and complex numbers whenever frequency is involved. Angles are in radians.

import numpy as np

Signals on a grid

np.arange / np.linspace

Time axes: n = 0 … N−1 for discrete time, t = n / fs for samples of a continuous signal.

fs = 1000
n = np.arange(1000)
t = n / fs                      # 1 s at 1 kHz
tau = np.linspace(-1, 1, 2001)  # includes both ends

arange excludes its end point; linspace includes it. Build time axes from integers (n / fs), not from a float step.

np.cos / np.sin / np.exp

Sinusoids and exponentials: x(t) = A cos(2π f t + φ), the decaying e^(−at), the complex exponential e^(jωn).

x = 2 * np.cos(2 * np.pi * 50 * t + np.pi / 4)
y = np.exp(-3 * t) * (t >= 0)
z = np.exp(1j * 0.2 * np.pi * n)

np.where / boolean masks

Piecewise signals: steps, rectangles, one-sided signals.

u = np.where(t >= 0, 1.0, 0.0)          # unit step
rect = (np.abs(t) <= 0.5).astype(float)

np.sinc

np.sinc(x)

The normalised sinc, sin(πx)/(πx): the transform of a rectangle and the ideal interpolator.

h = 0.25 * np.sinc(0.25 * (n - 50))    # ideal low-pass, cutoff π/4

np.sinc already contains the π: sinc(1) = 0. For sin(x)/x write np.sinc(x / np.pi).

np.zeros / np.ones / np.r_ / np.concatenate

Impulses, padding and signals built from pieces.

delta = np.zeros(32); delta[0] = 1.0
h = np.r_[np.ones(4), np.zeros(4)]
xp = np.concatenate([x, np.zeros(len(x))])   # zero padding

np.roll / x[::-1]

Circular shifts and time reversal.

xs = np.roll(x, 3)     # x[(n − 3) mod N]
xr = x[::-1]           # x[−n] on a symmetric grid

np.cumsum

A running sum: the accumulator y[n] = Σ x[k], a discrete integral.

step = np.cumsum(delta)

Complex numbers

1j and complex arrays

Python writes j as 1j; every arithmetic operation works on complex arrays.

z = 3 + 4j
H = 1 / (1 + 1j * w * RC)

np.abs / np.angle

Magnitude and phase of a phasor or a frequency response.

mag = np.abs(H)
phase = np.angle(H)              # radians, in (−π, π]
phase_deg = np.degrees(np.angle(H))

np.angle(z, deg=True) gives degrees directly. Use np.unwrap(np.angle(H)) to remove the 2π jumps.

np.real / np.imag / np.conj

Parts of a complex signal; a real signal’s spectrum is conjugate-symmetric.

x = np.real(np.fft.ifft(X))   # drop the round-off imaginary part

decibels

20 log10 of a magnitude, 10 log10 of a power.

mag_db = 20 * np.log10(np.abs(H))
psd_db = 10 * np.log10(Pxx)

Sums, peaks and checks

np.sum / np.mean

Energy, average power and inner products of sampled signals.

E = np.sum(np.abs(x) ** 2)            # discrete-time energy
P = np.mean(np.abs(x) ** 2)           # average power
Ec = np.sum(np.abs(x) ** 2) / fs      # approximates ∫|x(t)|² dt

np.argmax / np.argmin / np.max

Where a spectrum peaks, where a response is smallest.

k = np.argmax(np.abs(X))
print(f[k], np.abs(X[k]))

np.allclose

np.allclose(a, b, rtol=1e-5, atol=1e-8)

Check two results agree to rounding — the labs’ favourite test.

assert np.allclose(np.convolve(x, h), lfilter(h, 1, np.r_[x, np.zeros(len(h) - 1)]))

np.round / np.sort

Tidy printed results; order poles or frequencies.

print(np.round(np.sort(np.abs(p)), 4))

Convolution and polynomials

np.convolve

np.convolve(x, h, mode='full')

The convolution sum y[n] = Σ x[k] h[n − k].

y = np.convolve(x, h)       # length len(x) + len(h) − 1

mode='same' keeps the length of x (centred); 'valid' keeps only the full overlaps.

np.polymul

Multiply polynomials given as coefficient arrays, highest power first: cascaded systems.

den = np.polymul([1, 2], [1, 2, 5])   # (s + 2)(s² + 2s + 5)

np.roots

Roots of a polynomial: poles from a denominator, zeros from a numerator.

poles = np.roots([1, 3, 2])     # s² + 3s + 2 → −2, −1
stable = np.all(poles.real < 0)

The FFT (numpy.fft)

np.fft.fft / np.fft.ifft

The DFT of N samples and its inverse, computed by the FFT.

X = np.fft.fft(x)          # X[k], k = 0 … N−1
x2 = np.fft.ifft(X)

Bin k is the frequency k·fs/N; bins above N/2 are the negative frequencies.

np.fft.rfft / np.fft.rfftfreq

One-sided spectrum of a real signal, with its frequency axis in Hz.

X = np.fft.rfft(x)
f = np.fft.rfftfreq(len(x), d=1 / fs)
amp = 2 * np.abs(X) / len(x)       # amplitude of each sinusoid

The DC and Nyquist bins are not doubled in the amplitude scaling.

np.fft.fftshift / np.fft.fftfreq

Put zero frequency in the middle for a two-sided plot or table.

Xs = np.fft.fftshift(np.fft.fft(x))
fs_axis = np.fft.fftshift(np.fft.fftfreq(len(x), d=1 / fs))

zero padding

Pass n > len(x) to sample the DTFT more finely — it adds no resolution.

X = np.fft.rfft(x, n=8 * len(x))

Random signals

np.random.default_rng

rng = np.random.default_rng(seed)

Reproducible noise, random bits and codes.

rng = np.random.default_rng(7)
noise = rng.normal(0, 0.5, size=len(t))
bits = rng.integers(0, 2, size=100)

Seed it so the lab’s printed numbers are the same every run.

📡scipy.signal

The working library of the course: run a system on a signal, compute its frequency response, design filters, estimate spectra, change sample rates.

from scipy import signal   # or: from scipy.signal import lfilter, freqz, …

Running systems

lfilter

y = lfilter(b, a, x)

Run a difference equation a[0]y[n] + a[1]y[n−1] + … = b[0]x[n] + b[1]x[n−1] + … on x.

b, a = [0.1], [1, -0.9]       # y[n] = 0.9 y[n−1] + 0.1 x[n]
y = lfilter(b, a, x)

b is the numerator (feed-forward) and a the denominator (feedback); for an FIR filter a = 1.

filtfilt

y = filtfilt(b, a, x)

Filter forward and backward: zero phase, squared magnitude. Offline only.

y = filtfilt(b, a, noisy)

sosfilt

y = sosfilt(sos, x)

Run a filter as cascaded second-order sections — the safe way for high orders.

sos = butter(8, 100, fs=fs, output='sos')
y = sosfilt(sos, x)

fftconvolve

Convolution by FFT: the same answer as np.convolve, much faster for long signals.

y = fftconvolve(x, h)

lti / lsim / impulse / step

Continuous-time systems from numerator and denominator: responses to any input, an impulse or a step.

sys = lti([1], [1, 2, 5])
t, y = step(sys, T=np.linspace(0, 6, 601))
t, h = impulse(sys, T=t)
tout, yout, _ = lsim(sys, U=np.sin(3 * t), T=t)

Frequency response

freqz

w, H = freqz(b, a=1, worN=512, fs=2*np.pi)

H(e^(jω)) of a discrete-time system.

f, H = freqz(b, a, worN=2048, fs=fs)   # f in Hz because fs is given
mag_db = 20 * np.log10(np.abs(H))

Without fs, w is in radians per sample from 0 to π.

freqs

w, H = freqs(b, a, worN)

H(jω) of a continuous-time system (analog filters, RC circuits).

w = np.logspace(0, 4, 400)
w, H = freqs([1], [1e-3, 1], worN=w)

group_delay

w, gd = group_delay((b, a), w=512, fs=…)

Delay of each frequency in samples — constant for a linear-phase FIR filter.

f, gd = group_delay((h, 1), w=1024, fs=fs)

bode

Magnitude (dB) and phase (degrees) of a continuous-time system in one call.

w, mag, phase = bode(lti([100], [1, 10]), w=np.logspace(0, 3, 300))

Filter design

butter / cheby1 / ellip

butter(N, Wn, btype='low', analog=False, output='ba', fs=None)

Classical IIR designs: maximally flat, equiripple pass band, equiripple both bands.

b, a = butter(4, 200, btype='low', fs=fs)
b, a = cheby1(4, 1, 200, fs=fs)          # 1 dB ripple
b, a = ellip(4, 1, 60, 200, fs=fs)       # 1 dB ripple, 60 dB stop band

Without fs, Wn is a fraction of the Nyquist frequency (0 to 1). With analog=True, Wn is in rad/s.

buttord / cheb1ord

N, Wn = buttord(wp, ws, gpass, gstop, fs=…)

The lowest order that meets a pass-band and stop-band specification.

N, Wn = buttord(100, 150, 1, 40, fs=fs)
b, a = butter(N, Wn, fs=fs)

iirnotch

b, a = iirnotch(w0, Q, fs=…)

A second-order notch: remove one frequency (mains hum) and keep the rest.

b, a = iirnotch(50, 30, fs=fs)

firwin / firwin2

firwin(numtaps, cutoff, window='hamming', pass_zero=True, fs=…)

FIR design by the window method: from a cutoff, or from any gain curve.

h = firwin(101, 200, fs=fs)                         # low-pass
h = firwin(101, [100, 300], pass_zero=False, fs=fs) # band-pass
h2 = firwin2(101, [0, 200, 250, fs / 2], [1, 1, 0, 0], fs=fs)

kaiserord

N, beta = kaiserord(ripple_db, width)

Kaiser’s formula: taps and window shape for a given attenuation and transition width.

N, beta = kaiserord(60, (fst - fp) / (fs / 2))
h = firwin(N | 1, fp, window=('kaiser', beta), fs=fs)

width is a fraction of the Nyquist frequency, not of fs. N | 1 makes the length odd (a type I filter).

remez

remez(numtaps, bands, desired, weight=None, fs=…)

Parks–McClellan equiripple FIR design.

h = remez(61, [0, fp, fst, fs / 2], [1, 0], weight=[1, dp / ds], fs=fs)

tf2zpk

Zeros, poles and gain from numerator and denominator.

z, p, k = tf2zpk(b, a)
print(np.abs(p).max() < 1)     # stable if every pole is inside the unit circle

Partial fractions

residue

r, p, k = residue(b, a)

Partial fractions of H(s) = Σ r_i/(s − p_i) + k(s): inverse Laplace transforms by table.

r, p, k = residue([1, 3], [1, 3, 2])   # (s + 3)/((s + 1)(s + 2))

residuez

r, p, k = residuez(b, a)

The same for H(z) in powers of z^(−1): each term r/(1 − p z^(−1)) gives r·pⁿ.

r, p, k = residuez([1], [1, -1.5, 0.56])
h = sum(rk * pk ** n for rk, pk in zip(r, p)).real

Spectra

get_window

get_window('hann', N)

Window samples by name: hann, hamming, blackman, (kaiser, β).

w = get_window('hann', len(x))
X = np.fft.rfft(x * w)

get_window returns the periodic (DFT-even) window by default; pass fftbins=False for the symmetric one used in filter design.

periodogram / welch

Power spectral density: one FFT, or the average of many windowed segments (less variance).

f, P = periodogram(x, fs=fs)
f, Pw = welch(x, fs=fs, nperseg=256)

spectrogram

f, t, Sxx = spectrogram(x, fs, nperseg=…, noverlap=…)

The short-time Fourier transform as a power grid: frequency against time.

f, tt, S = spectrogram(x, fs, nperseg=256, noverlap=192)
ridge = f[np.argmax(S, axis=0)]     # the strongest frequency at each time

find_peaks

peaks, props = find_peaks(x, height=…, distance=…)

Indices of local maxima above a height: tones in a spectrum, pulses in a signal.

pk, _ = find_peaks(np.abs(X), height=0.1 * np.abs(X).max())
print(f[pk])

Rates, analytic signals and test signals

resample_poly / decimate

Change the sample rate by a ratio with built-in anti-aliasing filters.

y = resample_poly(x, up=3, down=2)    # 48 kHz → 72 kHz
d = decimate(x, 4)                    # low-pass, then keep every 4th sample

hilbert

The analytic signal x + j·H{x}: its magnitude is the envelope, its phase the instantaneous phase.

z = hilbert(x)
env = np.abs(z)
f_inst = np.diff(np.unwrap(np.angle(z))) * fs / (2 * np.pi)

chirp / square

Test signals: a swept sine, a square wave.

c = chirp(t, f0=10, t1=1, f1=200, method='linear')
sq = square(2 * np.pi * 5 * t)        # ±1, 50 % duty

max_len_seq

seq, state = max_len_seq(nbits)

A maximal-length pseudo-random sequence of 2^nbits − 1 chips (0 and 1).

seq, _ = max_len_seq(10)
code = 2.0 * seq - 1     # ±1 chips, length 1023

correlate

correlate(a, b, mode='full')

Cross-correlation: find a known code or pulse inside a noisy signal.

r = correlate(rx, code, mode='full')
lag = np.argmax(np.abs(r)) - (len(code) - 1)

∫SciPy: integrals, ODEs, special functions

Numerical checks on the formulas the lessons derive: an integral done by machine, a differential equation stepped in time, a Bessel function for FM.

from scipy.integrate import quad, solve_ivp
from scipy.special import jv, sici

Integration

quad

value, err = quad(f, a, b)

A definite integral to near machine precision: Fourier coefficients, energies, convolution integrals.

T = 2.0
ak = quad(lambda t: np.cos(2 * np.pi * 3 * t / T) * x(t), 0, T)[0] * 2 / T

quad returns a pair; take [0] for the value. Pass points=[…] at known jumps.

solve_ivp

solve_ivp(f, (t0, t1), y0, t_eval=…, max_step=…)

Step a differential equation in time — check a transform answer against direct simulation.

sol = solve_ivp(lambda t, y: [y[1], -4 * y[0] - 0.4 * y[1] + u(t)], (0, 10), [0, 0],
                t_eval=np.linspace(0, 10, 1001), max_step=0.01)
y = sol.y[0]

Special functions

jv

jv(n, beta)

Bessel functions of the first kind: the amplitudes of FM sidebands.

beta = 2.4
side = [jv(n, beta) for n in range(6)]   # carrier vanishes near β = 2.405

sici

Si, Ci = sici(x)

The sine integral: the Gibbs overshoot is Si(π)/π − ½ ≈ 9 % of the jump.

Si, _ = sici(np.pi)
print(Si / np.pi - 0.5)

∑SymPy

Exact transforms, integrals and partial fractions: derive by machine what the lesson derives by hand, then compare.

import sympy as sp

Symbols and exact numbers

sp.symbols

Declare symbols, with assumptions that help simplification.

t, f = sp.symbols('t f', real=True)
s = sp.symbols('s')
a, T = sp.symbols('a T', positive=True)

sp.Rational / sp.nsimplify / sp.pi / sp.I / sp.oo

Exact fractions, constants and infinity.

half = sp.Rational(1, 2)
print(sp.nsimplify(0.375))      # 3/8

sp.simplify / sp.expand / sp.factor

Tidy a result before comparing it with the lesson.

print(sp.factor(s**2 + 3*s + 2))     # (s + 1)*(s + 2)

sp.lambdify

Turn an expression into a NumPy function to evaluate on arrays.

Xf = sp.lambdify(f, X, 'numpy')
vals = Xf(np.linspace(-5, 5, 101))

Calculus

sp.integrate / sp.limit / sp.diff

Exact integrals (energies, coefficients), limits (initial and final values) and derivatives.

E = sp.integrate(sp.exp(-2*a*t), (t, 0, sp.oo))     # 1/(2a)
print(sp.limit(s * F, s, 0))                         # final value

sp.summation

Closed forms of series: geometric sums and z-transforms.

n = sp.symbols('n', integer=True, nonnegative=True)
z = sp.symbols('z')
X = sp.summation((sp.Rational(1, 2) / z)**n, (n, 0, sp.oo))

sp.Heaviside / sp.DiracDelta / sp.Piecewise

Steps, impulses and piecewise signals.

x = sp.exp(-a*t) * sp.Heaviside(t)
rect = sp.Piecewise((1, sp.Abs(t) <= sp.Rational(1, 2)), (0, True))

sp.dsolve

Solve a linear differential equation with initial conditions.

y = sp.Function('y')
sol = sp.dsolve(sp.Eq(y(t).diff(t, 2) + 3*y(t).diff(t) + 2*y(t), 0), y(t),
                ics={y(0): 1, y(t).diff(t).subs(t, 0): 0})

Transforms and partial fractions

sp.laplace_transform / sp.inverse_laplace_transform

Laplace pairs, both ways.

F = sp.laplace_transform(sp.exp(-2*t) * sp.cos(3*t), t, s, noconds=True)
f = sp.inverse_laplace_transform(1 / (s**2 + 3*s + 2), s, t)

The inverse comes back multiplied by Heaviside(t): the one-sided signal.

sp.apart

Partial fractions of a rational function, before reading terms off a table.

print(sp.apart((s + 3) / ((s + 1) * (s + 2)), s))

sp.fourier_transform

sp.fourier_transform(x, t, f)

The Fourier transform with frequency in hertz: X(f) = ∫ x(t) e^(−j2πft) dt.

X = sp.simplify(sp.fourier_transform(sp.exp(-a * sp.Abs(t)), t, f))

SymPy’s convention uses f in Hz (e^(−j2πft)), not ω: substitute f = ω/(2π) to compare with a table in ω.

sp.fourier_series

The Fourier series of a function on an interval, truncated to a few terms.

saw = sp.fourier_series(t, (t, -sp.pi, sp.pi))
print(saw.truncate(4))

sp.Poly

Work with a polynomial’s coefficients and degree.

den = sp.Poly(s**3 + 2*s**2 + 3*s + 1, s)
print(den.all_coeffs())

🎛️python-control

Transfer functions as objects: poles, gains, bandwidths, step specifications, feedback loops and stability margins — from Unit 4 on, and the core of the Laplace and feedback lessons.

import control as ct

Transfer functions

ct.tf

ct.tf(num, den)

A transfer function from coefficient lists, or algebra on the symbol s.

G = ct.tf([10], [1, 2, 10])
s = ct.tf('s')
H = 100 / (s * (s + 10))

ct.poles / ct.zeros

Poles and zeros: stable if every pole has a negative real part.

print(ct.poles(G))
print(all(p.real < 0 for p in ct.poles(G)))

ct.dcgain / ct.bandwidth

The gain at ω = 0 and the frequency where the gain falls 3 dB below it.

print(ct.dcgain(G), ct.bandwidth(G))     # bandwidth in rad/s

Responses and feedback

ct.step_info

ct.step_info(sys, SettlingTimeThreshold=0.02)

Rise time, settling time, overshoot (%) and peak of the step response.

info = ct.step_info(G)
print(info['Overshoot'], info['SettlingTime'])

ct.feedback

ct.feedback(G, H=1, sign=-1)

The closed loop G/(1 + GH), negative feedback by default.

T = ct.feedback(K * G)
print(ct.poles(T))

ct.margin

gm, pm, wcg, wcp = ct.margin(L)

Gain margin (as a ratio), phase margin (degrees) and the frequencies where they are read.

gm, pm, wg, wp = ct.margin(L)
print(20 * np.log10(gm), pm)

The gain margin comes back as a ratio: 20·log10(gm) gives decibels. An infinite margin means the phase never reaches −180°.

〰️PyWavelets

The continuous wavelet transform for the time–frequency lesson: short wavelets for high frequencies, long ones for low.

import pywt

Scalograms

pywt.frequency2scale

Choose scales from the frequencies you want to see.

wavelet = 'cmor1.5-1.0'
scales = pywt.frequency2scale(wavelet, np.geomspace(10, 400, 80) / fs)

pywt.cwt

C, freqs = pywt.cwt(x, scales, wavelet, sampling_period=1 / fs)

The continuous wavelet transform: a complex coefficient per scale and time, with each scale’s frequency in Hz.

C, freqs = pywt.cwt(x, scales, wavelet, sampling_period=1 / fs)
ridge = freqs[np.argmax(np.abs(C), axis=0)]

🐼pandas

Measured signals with real timestamps: daily readings, averages and the moving average as a filter.

import pandas as pd

Series with a time index

pd.Series / pd.date_range

Values on a calendar: one reading per day or per hour.

days = pd.date_range('2024-01-01', periods=365, freq='D')
temp = pd.Series(values, index=days)

.rolling(N).mean()

The N-point moving average — an FIR low-pass filter on a table.

smooth = temp.rolling(7, center=True).mean()

.mean() / .std() / pd.DataFrame

Summary statistics, and several signals side by side.

df = pd.DataFrame({'raw': temp, 'smooth': smooth})
print(df.std())

📈Signals cookbook

The small pieces of code the labs build again and again, and the traps that catch most first attempts.

Recipes

frequency of the strongest tone

One-sided spectrum, then its largest bin.

X = np.fft.rfft(x * np.hanning(len(x)))
f = np.fft.rfftfreq(len(x), d=1 / fs)
print(f[np.argmax(np.abs(X))])

the aliased frequency

Where a tone at f lands after sampling at fs: fold into 0 … fs/2.

def alias(f, fs):
    r = f % fs
    return min(r, fs - r)

a frequency response at one frequency

Evaluate H(e^(jΩ)) or H(jω) directly from the coefficients.

z = np.exp(1j * 2 * np.pi * 50 / fs)
H50 = np.polyval(b[::-1], 1 / z) / np.polyval(a[::-1], 1 / z)

b and a are in powers of z^(−1): reverse them for np.polyval, which wants the highest power first.

impulse response of a difference equation

Feed an impulse through lfilter.

d = np.zeros(50); d[0] = 1
h = lfilter(b, a, d)

energy in a band

Parseval: sum |X[k]|² over the bins of the band.

band = (f >= 40) & (f <= 60)
E_band = np.sum(np.abs(X[band]) ** 2)

Common traps

normalised frequency

scipy.signal designs take cutoffs as a fraction of Nyquist (fs/2) unless you pass fs=.

butter(4, 0.2)            # 0.2 × fs/2
butter(4, 100, fs=1000)   # 100 Hz — the same filter

convolution length

np.convolve returns len(x) + len(h) − 1 samples, not len(x).

y = np.convolve(x, h)[:len(x)]     # or mode='same' for a centred filter

two-sided against one-sided

np.fft.fft of a real cosine of amplitude A has two peaks of height A·N/2, at +f and −f; rfft keeps one.

A_est = 2 * np.abs(np.fft.rfft(x)[k]) / len(x)

radians and hertz

ω = 2πf. freqs and python-control work in rad/s; freqz with fs= and the FFT axes in Hz.

f_hz = w / (2 * np.pi)

circular against linear convolution

Multiplying two DFTs gives circular convolution. Pad both to len(x) + len(h) − 1 first.

N = len(x) + len(h) - 1
y = np.real(np.fft.ifft(np.fft.fft(x, N) * np.fft.fft(h, N)))