Level 03 · Signals, Systems & Control

Frequency

Signals in frequency with the DFT, and P and PI control.

Level 3 · Signals, Systems & Control: Frequency

Level 1 measured a frequency by counting zero crossings; Level 2 described a system by its answer to one short tap. This level sees both in frequency. Any recording can be taken apart into waves, and the discrete Fourier transform (DFT) says how much of each frequency it holds. An LTI system treats each wave separately, and its frequency response says how much of each gets through. Behind both is the complex number, used as a turning arrow.

The second half builds the first real controllers. Level 1's proportional controller stopped short; integral action gives it a memory that closes the gap, and PI control holds most of the world's ovens, motors and tanks at their setpoints. The memory brings its own problems, windup and bumps, and their standard cures.

The other paths meet these ideas. The DFT is two nested loops, O(N²) in the Algorithms path; the fast Fourier transform gets the same numbers in O(N log N) by divide and conquer. A bin's amplitude is the least-squares fit of one cosine, the Statistics path's line-fitting with a wave for a line. And in the Machine Learning path a spectrum is a favourite input for a model that recognises words or a failing bearing.

Complex numbers as turning arrows

A complex number is a pair, best seen as an arrow: the real part reaches right, the imaginary part up. Python writes the imaginary unit 1j (a bare j is a variable name). The arrow of length 1 at angle \theta (radians, anticlockwise from the right) is cmath.exp(1j * angle) (math.exp refuses complex numbers), with parts \cos\theta and \sin\theta:

e^{j\theta} = \cos\theta + j\sin\theta,

Euler's formula, with j the imaginary unit. abs(z) is the length, cmath.phase(z) the angle. Multiplying arrows multiplies lengths and adds angles, so a unit arrow turns without stretching. A pod on an observation wheel:

import cmath, math
z = 60 * cmath.exp(1j * math.radians(30))   # 60 m from the hub, 30 degrees above the right
print(round(z.real, 2), round(z.imag, 2), "  length", abs(z), "  angle", round(math.degrees(cmath.phase(z)), 1))
turn = cmath.exp(1j * math.radians(12))     # the wheel turns 12 degrees a minute
later = z * turn ** 5                       # five minutes later
print("after 5 minutes:", round(math.degrees(cmath.phase(later)), 1), "degrees,", round(abs(later), 2), "m out")

51.96 30.0   length 60.0   angle 30.0
after 5 minutes: 90.0 degrees, 60.0 m out

Five turns of 12 degrees reach the top; 1j * 1j is -1, two quarter turns. An arrow e^{2\pi j f t} goes round f times a second: a wave seen end on, whose real part is Level 1's cosine.

Practice: Filming a Reflector on a Bicycle Wheel

The DFT from its formula

For N samples x_0, \dots, x_{N-1} and each whole number k, the DFT makes one complex number, bin k:

X_k = \sum_{n=0}^{N-1} x_n \, e^{-2\pi j k n / N}

The arrow turns backwards k times over the recording, and each sample is multiplied by it. In Python the sum over n is a generator and the bins a list comprehension over k: [sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N)) for k in range(N)]. A wave making k cycles undoes the turning, so its products point one way and pile up; other whole numbers of cycles, or a constant, spread round the circle and cancel. A wave a\cos(2\pi k n / N + p) gives X_k = \frac{aN}{2} e^{jp}: the amplitude is a = 2|X_k| / N and the angle is the phase p. Bin 0's arrow never turns, so X_0 / N is the mean, the offset. An outdoor thermometer, read hourly for two days:

import cmath, math
N = 48                                       # hourly readings from midnight
temps = [12 + 5 * math.cos(2 * math.pi * 2 * n / N + math.radians(135))
         + 1.5 * math.cos(2 * math.pi * 4 * n / N) for n in range(N)]
X = [sum(temps[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N)) for k in range(N)]
print("offset:", round(X[0].real / N, 2), "  amplitudes, bins 1-5:", [round(2 * abs(X[k]) / N, 2) for k in range(1, 6)])
print("phase of bin 2:", round(math.degrees(cmath.phase(X[2]))), "degrees")
plot([2 * abs(X[k]) / N for k in range(N // 2 + 1)], kind="stem", xlabel="bin k", ylabel="amplitude (°C)")

offset: 12.0   amplitudes, bins 1-5: [0.0, 5.0, 0.0, 1.5, 0.0]
phase of bin 2: 135 degrees

Press Run: two stems rise from a flat floor. Bin 2, two cycles in 48 hours, is the daily rhythm, 5 °C either side of 12 °C; bin 4 a smaller twice-daily wave. The daily wave is 135 degrees past its peak at midnight, so it peaks 225 degrees later, 225 / 360 * 24 = 15 hours after midnight. Each bin is one pass over the data, so the whole DFT costs N * N products: two billion for one second of audio.

Common mistakes: +2j in the forward transform. Doubling bin 0. Dividing by N inside the sum and again outside.

Practice: How Much Hum Is in the Recording?

Bins and resolution

The list knows nothing of seconds; fs does. Bin k is k cycles in N / f_s seconds, so it stands for

f_k = \frac{k f_s}{N} \text{ hertz.}

The spacing f_s / N is the frequency resolution, one over the recording time: only a longer recording sharpens it. For a real signal bins 0 to N // 2 (the Nyquist frequency) carry the information; the rest mirror them. A piano tuner hears two strings, at 220 and 222 Hz:

import cmath, math
fs = 1000
for N in (250, 500):
    x = [math.cos(2 * math.pi * 220 * n / fs) + math.cos(2 * math.pi * 222 * n / fs) for n in range(N)]
    ks = range(round(216 * N / fs), round(226 * N / fs) + 1)       # only the bins near 220 Hz
    amps = [2 * abs(sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N))) / N for k in ks]
    print(f"{fs / N} Hz apart:", "  ".join(f"{k * fs / N:g}:{a:.2f}" for k, a in zip(ks, amps)))
    plot([k * fs / N for k in ks], amps, kind="stem", label=f"{N / fs} s", xlabel="Hz")

4.0 Hz apart: 216:0.21  220:1.19  224:0.64
2.0 Hz apart: 216:0.00  218:0.00  220:1.00  222:1.00  224:0.00  226:0.00

A quarter of a second shows one lump; half a second shows both strings. (Only the bins of interest were computed, one pass each.) A faster sample rate would not have helped.

Practice: The Frequency Ruler of a Spectrum App, The Spectrum of a Motor Current

The dominant frequency

"What is the main rhythm?" is the strongest bin, leaving out bin 0, which holds the offset and is often the biggest, turned into hertz. A chest strap stretches as its wearer breathes:

import cmath, math, random
random.seed(7)
fs, N = 4, 120                               # 4 readings a second for 30 s
strap = [50 + 3 * math.sin(2 * math.pi * 0.27 * n / fs) + random.gauss(0, 1) for n in range(N)]
amps = [2 * abs(sum(strap[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N))) / N for k in range(N // 2 + 1)]
best = max(range(1, N // 2 + 1), key=lambda k: amps[k])
print("largest bin:", amps.index(max(amps)), "  largest wave: bin", best, "=", round(60 * best * fs / N, 1), "breaths a minute")
plot([k * fs / N for k in range(1, N // 2 + 1)], amps[1:], kind="stem", xlabel="Hz")

largest bin: 0   largest wave: bin 8 = 16.0 breaths a minute

Noise spreads thinly over all bins while the breathing piles up in one. The true rate was 16.2; bins are 2 breaths a minute apart, so the answer is good to about one.

Practice: Why Does the Footbridge Sway?

Aliasing and the sampling theorem

Level 1 showed a 9 Hz wave sampled 10 times a second posing as 1 Hz. In a spectrum, a tone above fs / 2 turns up at a lower apparent frequency, for two reasons: adding a whole multiple of fs changes no sample, since each moves on by whole cycles; and a tone at fs - f gives the samples of one at f (a cosine exactly, a sine with its sign flipped). So frequencies fold into 0 to fs / 2, like paper folded at every multiple of fs / 2:

import cmath, math
fs, N = 400, 40                              # bins 10 Hz apart
def peak_hz(f):
    x = [math.cos(2 * math.pi * f * n / fs) for n in range(N)]
    size = [abs(sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N))) for k in range(N // 2 + 1)]
    return size.index(max(size)) * fs / N
for f in (130, 270, 330, 730):
    print(f, "Hz shows up at", peak_hz(f), "Hz")

130 Hz shows up at 130.0 Hz
270 Hz shows up at 130.0 Hz
330 Hz shows up at 70.0 Hz
730 Hz shows up at 70.0 Hz

270 Hz is 70 above the 200 Hz fold and lands 70 below it; 730 Hz loses a whole 400, then folds like 330. The sampling theorem says when no harm is done: if every frequency is strictly below fs / 2, the samples hold everything, so fs must exceed twice the highest frequency. Exactly fs / 2 fails: sampled at its zero crossings, a wave reads as silence. Converters put an analogue anti-aliasing filter in front and leave room to spare.

Common mistakes: <= where the theorem needs <. Folding only once. Expecting the spectrum to warn you: an aliased peak looks like any other.

Practice: The Propeller Hum the Logger Gets Wrong, Is the Gearbox Logger Fast Enough?

Leakage and the Hann window

The DFT treats a recording as one period of an endless repeat. A tone with a whole number of cycles sits in one bin; otherwise the end does not meet the start, and the jump smears the tone over every bin: leakage, a skirt that can bury a weak tone. The cure is a window that fades the recording in and out. The Hann window is one cycle of a raised cosine, 0 at the ends and 1 in the middle:

w_n = 0.5 - 0.5\cos\left(\frac{2\pi n}{N}\right), \qquad n = 0, \dots, N - 1

A flute note of 117 Hz makes 7.49 cycles in 64 samples:

import cmath, math
fs, N = 1000, 64
note = [math.cos(2 * math.pi * 117 * n / fs) for n in range(N)]
hann = [0.5 - 0.5 * math.cos(2 * math.pi * n / N) for n in range(N)]
for name, x in (("plain", note), ("Hann", [a * w for a, w in zip(note, hann)])):
    size = [abs(sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N))) for k in range(N // 2 + 1)]
    db = [20 * math.log10(s / max(size)) for s in size]          # dB below the strongest bin
    print(f"{name:5}  bin 15: {db[15]:6.1f} dB   bin 30: {db[30]:6.1f} dB")
    plot(db, kind="stem", label=name, xlabel="bin", ylabel="dB below the peak")

plain  bin 15:  -21.0 dB   bin 30:  -26.3 dB
Hann   bin 15:  -60.7 dB   bin 30:  -96.6 dB

Press Run. Without the window, far bins are only 21 to 26 dB down, enough to hide a rattle 30 dB quieter; with it, the skirt falls below -60 dB. The price: the peak widens to about three bins, and amplitudes read half size (the window's mean is 0.5). Multiply the samples, not the spectrum.

Practice: Leakage in an Engine Recording

The frequency response of an FIR filter

Feed a sine into an LTI system and the same frequency comes out, scaled and shifted. For an FIR filter with weights b (Level 2), both are one complex number, the frequency response:

\omega = \frac{2\pi f}{f_s} \ \text{(radians per sample)}, \qquad H(\omega) = \sum_{k} b_k \, e^{-j\omega k}

Weight k meets the wave k samples late, an arrow turned back by \omega k; where the arrows cancel, the wave is removed. The gain |H| (abs(H)) multiplies the amplitude, 20\log_{10}|H| in decibels, and the angle of H is the phase shift. A 1-2-1 smoother for speech at 8 kHz:

import cmath, math
fs, b = 8000, [0.25, 0.5, 0.25]
def response(f):
    w = 2 * math.pi * f / fs
    return sum(b[k] * cmath.exp(-1j * w * k) for k in range(len(b)))
for f in (1000, 2000, 3000):
    H = response(f)
    print(f"{f:4} Hz: gain {abs(H):.4f} = {20 * math.log10(abs(H)):6.2f} dB, phase {math.degrees(cmath.phase(H)):6.1f} degrees")
x = [math.sin(2 * math.pi * 2000 * n / fs) for n in range(10)]                 # 0, 1, 0, -1, ...
print("2000 Hz out:", [round(sum(b[k] * x[n - k] for k in range(3)), 2) + 0 for n in range(2, 10)])
freqs = [20 * 1.05 ** i for i in range(106)]                                    # 20 Hz to 3.4 kHz
plot(freqs, [20 * math.log10(abs(response(f))) for f in freqs], logx=True, xlabel="Hz", ylabel="gain (dB)")

1000 Hz: gain 0.8536 =  -1.38 dB, phase  -45.0 degrees
2000 Hz: gain 0.5000 =  -6.02 dB, phase  -90.0 degrees
3000 Hz: gain 0.1464 = -16.69 dB, phase -135.0 degrees
2000 Hz out: [0.5, 0.0, -0.5, 0.0, 0.5, 0.0, -0.5, 0.0]

The 2000 Hz tone, 0, 1, 0, -1 repeating, comes out at half size and one sample late, a quarter cycle: -90 degrees, as H said. This low-pass filter's gain at fs / 2 is exactly 0, minus infinity decibels, so code printing decibels needs a floor for a true null. On the logarithmic axis every octave has the same width; the gain crosses -3 dB (half the power) near 1460 Hz, the cut-off.

Common mistakes: 10 * log10 for a gain. Forgetting / fs. Taking log10 of 0.

Practice: What the Drone's Sensor Filter Lets Through, Choosing the Filter for a Weighing Scale

The inverse DFT and negative frequencies

The DFT loses nothing. The inverse DFT adds every bin's arrow, now turning forwards, and divides by N:

x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k \, e^{2\pi j k n / N}

Why do bins above N / 2 exist? A cosine is two arrows, \cos\theta = \frac{e^{j\theta} + e^{-j\theta}}{2}, turning opposite ways so their imaginary parts cancel. At the sampling moments, turning N - k times forwards is the same as k times backwards, so bin N - k holds the negative frequency -k * fs / N. A lighthouse lamp, read 8 times per turn:

import cmath, math
N = 8
lamp = [3 + 2 * math.cos(2 * math.pi * n / N) for n in range(N)]
X = [sum(lamp[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N)) for k in range(N)]
idft = lambda X: [sum(X[k] * cmath.exp(2j * math.pi * k * n / N) for k in range(N)) / N for n in range(N)]
print("bins:", [round(abs(v), 6) + 0 for v in X])
print("round trip within 1e-12:", max(abs(a - b) for a, b in zip(idft(X), lamp)) < 1e-12)
half = idft([X[0], X[1], 0, 0, 0, 0, 0, 0])              # keep bin 1, drop its mirror, bin 7
print("without bin 7:", [complex(round(v.real, 2), round(v.imag, 2) + 0) for v in half[:3]])

bins: [24.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0]
round trip within 1e-12: True
without bin 7: [(4+0j), (3.71+0.71j), (3+1j)]

The wave is shared between bins 1 and 7, which is why the amplitude rule doubles one side. Keep only one and the result is not real: half a wheel turning one way. So removing a frequency from a spectrum means removing both arrows of the pair. Bin k above N / 2 has frequency (k - N) * fs / N.

Practice: Cleaning a Heart-Rate Signal in the Frequency Domain

P control and its steady-state error

A 3D printer's nozzle must melt plastic at 210 °C. Its heater adds 0.1 °C a second per watt, and the nozzle loses heat to the 25 °C room with a time constant of 40 s. (For now the heater has no limit.) Setting the rate of change to zero, as in Level 2, with u = K_p e and T = 210 - e, gives P control's error:

e_\infty = \frac{210 - 25}{1 + K_p \cdot 0.1 \cdot 40}
dt, r = 0.1, 210
for Kp in (1, 4, 16):                        # watts per °C of error
    T = 25.0
    for step in range(3000):                 # 5 minutes
        T = T + dt * ((25 - T) / 40 + 0.1 * Kp * (r - T))
    print(f"Kp = {Kp:2}: {r - T:5.2f} °C short; formula {(r - 25) / (1 + Kp * 0.1 * 40):5.2f} °C")

Kp =  1: 37.00 °C short; formula 37.00 °C
Kp =  4: 10.88 °C short; formula 10.88 °C
Kp = 16:  2.85 °C short; formula  2.85 °C

K_p \cdot 0.1 \cdot 40 is the loop gain, the controller's gain times the plant's steady gain (4 °C per watt). The error is the step over 1 + loop gain: smaller with more gain, never zero.

Practice: The Gap a P Controller Leaves

Integral action

The fix is a memory. Integral action keeps a running total of the error, each sample weighted by how long it lasted:

I_k = I_{k-1} + e_k \, \Delta t = \sum_{i=0}^{k} e_i \, \Delta t, \qquad u_k = K_p e_k + K_i I_k

In code the running sum is one line, I = I + e * dt, and the command u = Kp * e + Ki * I. I is the area under the error curve (°C·s); leaving out * dt would make Ki depend on the step size. I stops changing only when e = 0, so a PI loop that settles, settles at the setpoint. The order within a step changes the exact numbers; this loop integrates first:

dt, r, Kp = 0.1, 210, 4
for Ki in (0, 0.1):
    T, I, temps = 25.0, 0.0, [25.0]
    for step in range(1200):                 # 2 minutes
        e = r - T
        I = I + e * dt                       # integrate first,
        u = Kp * e + Ki * I                  # then command
        T = T + dt * ((25 - T) / 40 + 0.1 * u)
        temps.append(T)
    print(f"Ki = {Ki}: {T:.2f} °C after 2 minutes; P term {Kp * e:.1f} W, I term {Ki * I:.1f} W")
    plot([k * dt for k in range(len(temps))], temps, label=f"Ki = {Ki}", xlabel="seconds", ylabel="nozzle (°C)")

Ki = 0: 199.12 °C after 2 minutes; P term 43.5 W, I term 0.0 W
Ki = 0.1: 210.00 °C after 2 minutes; P term 0.0 W, I term 46.2 W

The P controller needs an error to make its 43.5 W. With PI the P term does the early work and fades, and the I term ends supplying all 46.2 W, exactly what holding 210 °C needs, (210 - 25) / 40 / 0.1.

Practice: A Controller With a Memory, Spinning Up a Motor With PI

Choosing Ki: overshoot and settling

Too little integral and the last degrees take minutes; too much and the controller remembers the cold start too long and sails past. The integral time T_i = K_p / K_i is how long a constant error takes to make the I term equal the P term. Measure overshoot and 2 % settling against the setpoint, where a PI loop ends:

dt, r, Kp = 0.1, 210, 4
for Ki in (0.02, 0.1, 0.3, 3):
    T, I, temps = 25.0, 0.0, [25.0]
    for _ in range(3000):
        I = I + (r - T) * dt
        T = T + dt * ((25 - T) / 40 + 0.1 * (Kp * (r - T) + Ki * I))
        temps.append(T)
    over = max(0.0, 100 * (max(temps) - r) / (r - 25))
    outside = [k for k, T in enumerate(temps) if abs(T - r) > 0.02 * (r - 25)]
    print(f"Ki = {Ki:4}: Ti = {Kp / Ki:5.1f} s  overshoot {over:4.1f} %  settled after {(outside[-1] + 1) * dt:5.1f} s")
    plot([k * dt for k in range(600)], temps[:600], label=f"Ki = {Ki}", xlabel="seconds")

Ki = 0.02: Ti = 200.0 s  overshoot  0.0 %  settled after 187.1 s
Ki =  0.1: Ti =  40.0 s  overshoot  0.0 %  settled after   9.6 s
Ki =  0.3: Ti =  13.3 s  overshoot  6.7 %  settled after  28.6 s
Ki =    3: Ti =   1.3 s  overshoot 36.4 %  settled after  18.3 s

The best row has Ti equal to the plant's time constant, 40 s: for a first-order plant that cancels its sluggishness, a good start for tuning. Below it the response creeps; above it the overshoot grows, and the settling time no longer improves. A plant that must never overshoot rules out the lively rows however fast they settle.

Practice: How Much Integral Is Enough?

Windup and clamping

The integral does not know the actuator's limits. A compressor fills a workshop's air tank: each per cent of motor power (0 to 100) adds 0.0012 bar a second, and air leaks at the pressure over 100 s. Asked for 8 bar from empty, it runs flat out while the integral piles up error; when the pressure arrives, the I term alone demands more than full power, and unwinding it takes an overshoot: integrator windup. The usual cure, clamping: keep the new integral only if the command it gives, before limiting, is one the actuator can obey.

dt, r, Kp, Ki = 0.5, 8, 20, 0.4              # s; bar; % per bar; % per bar-second
for clamp in (False, True):
    P, I, reached, pressures = 0.0, 0.0, None, [0.0]
    for step in range(800):
        e = r - P
        I_try = I + e * dt
        if not clamp or 0 <= Kp * e + Ki * I_try <= 100:
            I = I_try                        # clamped: integrate only while the motor can obey
        u = max(0, min(100, Kp * e + Ki * I))
        P = P + dt * (0.0012 * u - P / 100)
        pressures.append(P)
        if reached is None and P >= r:
            reached = (step + 1) * dt, Ki * I
    name = "clamped" if clamp else "plain"
    print(f"{name:7}: 8 bar at {reached[0]:5.1f} s, I term then {reached[1]:5.1f} %, peak {max(pressures):.2f} bar")
    plot([k * dt for k in range(len(pressures))], pressures, label=name, xlabel="seconds", ylabel="bar")

plain  : 8 bar at 110.0 s, I term then 144.7 %, peak 9.73 bar
clamped: 8 bar at 140.0 s, I term then  75.1 %, peak 8.20 bar

The plain controller arrives asking 145 % from its integral alone and overshoots by 1.7 bar; the clamped one arrives with 75 %, near the 67 % that holding 8 bar needs. Back-calculation, another cure, bleeds the integral by the gap between asked and delivered command.

Common mistakes: limiting the integral instead of the command. Testing the limit with the old integral. Forgetting to limit u afterwards.

Practice: The Oven That Would Not Stop Heating

Bumpless transfer and setpoint weighting

A helmsman has steered by hand and hands over to the autopilot. With an empty integral, the first command is just the P term, and the rudder jumps. Bumpless transfer loads the integral beforehand so the first automatic command equals what the actuator is doing: solve K_p e + K_i I = u_{\text{man}} for I, giving I = (u_{\text{man}} - K_p e) / K_i. With Kp = 1.5, Ki = 0.05, the ship 2 degrees off course and the helmsman holding 5 degrees of rudder, an empty integral would command 3 degrees, a jump of -2; loading I = (5 - 1.5 * 2) / 0.05 = 40 makes the first command exactly 5. The loaded integral hands over the helmsman's knowledge: this current needs 5 degrees. A controller that integrates before commanding adds e * dt first, so the loaded value must leave room for that update.

When the setpoint jumps, the P term jumps with the error: a setpoint kick. Setpoint weighting lets the P term see only a share b of the setpoint r, while the integral keeps the full error:

u_k = K_p\,(b\,r - y_k) + K_i I_k, \qquad I_k = I_{k-1} + (r - y_k)\,\Delta t

The integral still decides where the output settles, so it ends at r; b changes only the path. A wind tunnel's fan, settled at 10 m/s, is asked for 15:

dt, Kp, Ki = 0.1, 8, 1.6
for b in (1, 0):
    v, u_prev, kick = 10.0, 40.0, 0.0
    I = (40 - Kp * (b * 10 - v)) / Ki        # the integral that holds 10 m/s on 40 %
    for step in range(400):
        I = I + (15 - v) * dt
        u = Kp * (b * 15 - v) + Ki * I
        kick, u_prev = max(kick, abs(u - u_prev)), u
        v = v + dt * (0.05 * u - v / 5)
        if abs(v - 15) > 0.1:
            late = step + 1
    print(f"b = {b}: largest change of command {kick:4.1f} %, within 0.1 m/s after {late * dt:4.1f} s")

b = 1: largest change of command 40.8 %, within 0.1 m/s after  9.8 s
b = 0: largest change of command  0.8 %, within 0.1 m/s after 23.2 s

With b = 1 the command leaps 41 % in one step; with b = 0 it never moves more than 0.8 %, for a slower climb. (With b below 1 the P term does not vanish at the setpoint, so a settled loop holds a different I.) A disturbance still meets the full P term.

Practice: Handing Over Without a Bump, A Treadmill That Does Not Jerk

Rejecting a constant disturbance

The real test is a disturbance nobody planned: after eight minutes a sandblaster starts drawing 0.01 bar a second from the workshop tank:

dt, r, Kp = 0.5, 8, 20
for Ki in (0, 0.4):
    P, I, pressures = 0.0, 0.0, []
    for step in range(1920):                 # 16 minutes
        e = r - P
        I_try = I + e * dt
        if 0 <= Kp * e + Ki * I_try <= 100:
            I = I_try
        u = max(0, min(100, Kp * e + Ki * I))
        P = P + dt * (0.0012 * u - P / 100 - (0.01 if step >= 960 else 0))
        pressures.append(P)
        if step == 959:
            before = P
    print(f"Ki = {Ki}: {before:.2f} bar before, {P:.2f} bar after; I term {Ki * I:.1f} %")
    plot([k * dt for k in range(len(pressures))], pressures, label=f"Ki = {Ki}", xlabel="seconds", ylabel="bar")

Ki = 0: 5.65 bar before, 5.35 bar after; I term 0.0 %
Ki = 0.4: 8.00 bar before, 8.00 bar after; I term 75.0 %

The P controller was 2.35 bar short, (8 / 100) / (20 * 0.0012 + 1 / 100), and loses 0.3 bar more for good. Press Run: the PI controller dips to 7.8 bar and climbs back, because its integral grows from the 67 % that held 8 bar until it adds the 8.3 % the sandblaster takes, 0.01 / 0.0012. The integral learns the disturbance without anyone measuring it, and a settled PI loop cancels any constant one exactly: a tank's outflow, an open lid, a hill.

Practice: The Hot Tub With the Lid Off

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

Checklist

You can now:

  • Use a complex number as an arrow: build it with cmath.exp(1j * angle), read its length and angle, turn it by multiplying.
  • Compute DFT bins from the formula, read offset, amplitude and phase, and find the dominant frequency.
  • Give a bin's frequency and a recording's resolution.
  • Fold a frequency above Nyquist, and choose a sample rate by the sampling theorem.
  • Explain leakage, apply a Hann window and compare spectra in decibels.
  • Compute an FIR filter's gain, decibels and phase at any frequency.
  • Rebuild a signal with the inverse DFT and explain negative frequencies.
  • Predict P control's steady-state error from the loop gain.
  • Simulate PI control step by step, explain why it settles at the setpoint and rejects a constant disturbance, and choose Ki from overshoot and settling time.
  • Prevent windup by clamping, load the integral for a bumpless transfer, and soften setpoint kicks by weighting.