A build log on a 5-band parametric equalizer written entirely in MATLAB: the filter design in detail, the pole-zero reading that tells you what a filter does before you plot it, and an honest audit of where my implementation diverges from its own UI. Five cascaded biquads, RBJ cookbook coefficients, a live response curve, a rolling spectrogram and a pole-zero view. One classdef, no App Designer, no toolboxes beyond base MATLAB. The design: twenty lines that contain a bilinear transform Here is the entire filter design, verbatim: function [b, a] = biquadCoeffs(app, i) wc = 2 * pi * app.Freqs(i) / app.Fs; wc = max(0.005, min(pi - 0.01, wc)); alpha = sin(wc) / (2 * app.Qs(i)); cosw = cos(wc); switch app.Tipos{i} case 'Peaking' A = 10^(app.Gains(i) / 40); b = [1+alpha*A, -2*cosw, 1-alpha*A]; a = [1+alpha/A, -2*cosw, 1-alpha/A]; case 'High-Pass' b = [(1+cosw)/2, -(1+cosw), (1+cosw)/2]; a = [1+alpha, -2*cosw, 1-alpha ]; case 'Low-Pass' b = [(1-cosw)/2, (1-cosw), (1-cosw)/2]; a = [1+alpha, -2*cosw, 1-alpha ]; end end Almost everything interesting about this project is compressed into those coefficients, so it is worth unpacking properly. Where the bilinear transform went These are digital filters, but nobody designs a biquad directly in the z-domain. You start from an analog prototype, in the peaking case s² + (A·ω0/Q)·s + ω0² H(s) = ------------------------ s² + (ω0/(A·Q))·s + ω0² and map the s-plane onto the z-plane with the bilinear transform, s = (2/T)·(1−z⁻¹)/(1+z⁻¹). That map is not linear in frequency. It squeezes the entire infinite analog axis into the finite interval [0, π], according to ω_digital = 2·arctan(ω_analog·T / 2) which is why you prewarp: you deliberately design the analog prototype at the distorted frequency that maps back onto the digital frequency you actually wanted. Prewarping at ω0 is what makes the center frequency come out exact. And that is the reason cos(wc) and sin(wc) appear in the code at all. If you carry the algebra through with prewarping at ω0, every tan(ω0/2) collapses into trigonometric functions of ω0 itself, and the coefficients come out in the closed form above. The trig is not a shortcut or an approximation. It is the bilinear transform, already solved. Which also tells you what the design does not promise. Prewarping is exact at one frequency and only one. Everything else is warped, and the warping compresses hardest near Nyquist. That is measurable, and I measured it further down. Alpha, and what Q means here α = sin(ω0) / (2Q) α is the only place the bandwidth enters. For the low-pass and high-pass this is the classical resonance Q: Q = 1/√2 ≈ 0.707 gives the maximally flat Butterworth response, higher Q gives a resonant peak at the corner. One caveat worth stating precisely, because it is a common misreading: for the peaking filter, RBJ defines the bandwidth between the half-gain points (half the boost in dB), not the conventional −3 dB points. A 12 dB bell at Q = 2 has its stated bandwidth measured at 6 dB. So Q numbers are not directly comparable between a bell and a shelf, and they are not comparable between EQs unless both tell you which convention they use. Why 10^(G/40) and not 10^(G/20) The obvious way to turn decibels into a linear gain is 10^(G/20). The code uses forty, and the reason is a small piece of algebra that is worth doing once. Evaluate the peaking filter on the unit circle at the center frequency, z = e^{jω0}, and factor e^{−jω0} out of the numerator: N(ω0) = e^{−jω0} · [ (1+αA)·e^{jω0} − 2cos ω0 + (1−αA)·e^{−jω0} ] = e^{−jω0} · [ 2cos ω0 − 2cos ω0 + j·2αA·sin ω0 ] = e^{−jω0} · j·2αA·sin ω0 The cosine terms annihilate and only the imaginary part survives. The denominator is the same computation with α/A in place of αA, so almost everything cancels: H(ω0) = (2αA·sin ω0) / (2(α/A)·sin ω0) = A² The filter delivers A², not A, because A multiplies the numerator and divides the denominator. So to get G decibels of boost you need A² = 10^(G/20), hence A = 10^(G/40). Verified numerically on the shipping coefficients, which is the only reason I trust the derivation rather than my memory of it: | Slider | 20·log10|H(ω0)| | |---|---| | +6 dB | +6.0000 | | −9 dB | −9.0000 | | +12 dB | +12.0000 | Reading the poles and zeros This is the part I would want a student to take away, because it turns filter design from coefficient soup into something you can predict on paper. The zeros are the filter's identity Look at the numerators, not the denominators. High-pass: b = (1+cos ω0)/2 · [1, −2, 1]. The polynomial 1 − 2z⁻¹ + z⁻² factors as (1 − z⁻¹)², a double zero at z = +1. On the unit circle z = 1 is DC, so the response is forced to exactly zero at 0 Hz. Not attenuated. Zero. Low-pass: b = (1−cos ω0)/2 · [1, 2, 1], which is (1 + z⁻¹)², a double zero at z = −1, and z = −1 is Nyquist. Both are second-order nulls, which is where the 12 dB/octave rolloff comes from: each zero contributes 6 dB/octave. The denominators of both are identical. A low-pass and a high-pass at the same frequency and Q have exactly the same poles and differ only in where they put their zeros. That is the whole distinction. The pole radius has a closed form For the LP/HP denominator a = [1+α, −2cos ω0, 1−α], normalise by a(1) and use the fact that for a monic quadratic the product of the roots is the constant term. The poles are a complex conjugate pair, so their product is |p|²: |p|² = (1 − α)/(1 + α) ⟹ r = sqrt( (1−α)/(1+α) ) Two consequences fall straight out. The filter is unconditionally stable. α > 0 for any positive Q, so (1−α)/(1+α) < 1, so r < 1 always. There is no Q, no frequency and no sample rate at which this design puts a pole outside the unit circle. That is a property of the form, not of your parameter clamping. Q is pole radius. As Q → ∞, α → 0 and r → 1: the pole walks toward the unit circle and the response sharpens into a resonance. Checked against roots() on the actual coefficients at 1 kHz: Q α sqrt((1−α)/(1+α)) abs(roots(a)) 0.5 0.14199 0.866788 0.866788 0.7071 0.10041 0.904163 0.904163 2.0 0.03550 0.965110 0.965110 10.0 0.00710 0.992925 0.992925 Six decimal places. The formula is not an approximation. Boost sharpens the poles, cut sharpens the zeros The peaking filter is where this gets elegant. Its numerator carries αA and its denominator α/A, so by the same argument: r_pole = sqrt( (1 − α/A)/(1 + α/A) ) r_zero = sqrt( (1 − αA)/(1 + αA) ) At 1 kHz, Q = 2: Gain pole radius zero radius +12 dB 0.982364 0.931511 −12 dB 0.931511 0.982364 They swap. Exactly. A boost pulls the poles toward the unit circle to build the peak; a cut pulls the zeros toward the circle to dig the notch, and leaves the poles further in. And that symmetry is not a coincidence of these numbers. Substituting A → 1/A maps the numerator coefficients onto the denominator coefficients and vice versa, so H_{+G}(z) · H_{−G}(z) = 1 exactly, at every frequency The RBJ peaking filter is its own inverse. I checked it on 512 log-spaced points: the largest deviation of the product from unity was 4.6 × 10⁻¹⁶, which is floating point noise and nothing else. Boost 9 dB at 1 kHz with Q = 2, then cut 9 dB at 1 kHz with Q = 2, and you have mathematically reconstructed the input signal. Not approximately undone it. The cascade Five biquads in series. Two consequences, one for the audio and one for the picture. For the audio, applyEQ runs them sequentially: for i = 1:app.NumBands [b, a] = app.biquadCoeffs(i); b = b / a(1); a = a / a(1); out = filter(b, a, out); end (The explicit a(1) normalisation is redundant, since filter normalises internally, but it makes the difference equation the code is implementing unambiguous. And it is worth knowing that MATLAB's filter is a transposed direct form II, which has better floating-point behaviour than the plain DF-II you would write by hand.) For the picture, cascading LTI systems multiplies transfer functions, so updateEQ evaluates the product directly on the unit circle: f = logspace(log10(20), log10(min(app.Fs/2 - 1, 20000)), 512); z = exp(1j * 2*pi*f / app.Fs); H = ones(size(f)); for i = 1:app.NumBands [b, a] = app.biquadCoeffs(i); H = H .* (b(1)+b(2).*z.^-1+b(3).*z.^-2) ./ ... (a(1)+a(2).*z.^-1+a(3).*z.^-2); end magdB = 20*log10(abs(H) + eps); z = exp(1j·2πf/Fs) is sampling the DTFT: walking the unit circle and reading H(z) at each stop. Points are log-spaced, because the display is logarithmic in frequency and linear spacing would put hundreds of points in the top octave where nothing changes and a handful below 200 Hz where everything does. The important structural property is that both loops call the same biquadCoeffs. The curve cannot drift away from the audio, because there is exactly one place the coefficients exist. The tempting alternative, drawing an idealised bell shape and separately writing the DSP, is how EQ displays end up lying to their users. The POLES button does the third version of the same idea. Cascading in the coefficient domain is polynomial multiplication: b_tot = conv(b_tot, b); a_tot = conv(a_tot, a); Five biquads convolve into a single 10th-order transfer function with 11 coefficients each side, and zplane(b_tot, a_tot) draws the composite. So the plot above is not five overlaid biquads, it is the one filter the audio actually passes through. Three places the implementation does not match its own UI Deriving the maths carefully is also how you find out where your code stops obeying it. All three of these came out of writing this post. The frequency slider lies below 35 Hz. That clamp on line 3: wc = max(0.005, min(pi - 0.01, wc)); is there to keep cos(ω0) away from the degenerate values at DC and Nyquist. Reasonable. But solve it for frequency at 44.1 kHz: f_min = 0.005 · 44100 / 2π = 35.09 Hz The slider goes down to 20 Hz, and 20 Hz maps to ω = 0.00285, well under the clamp. So every band set between 20 and 35 Hz produces an identical filter at 35.09 Hz, while the label happily reads 20 Hz. The upper clamp is harmless by comparison: π − 0.01 corresponds to 21,980 Hz, above the 20 kHz the UI allows. The warping is visible at the top of the band, as predicted. A 12 dB bell at Q = 2, measured between its half-gain points: f0 half-gain band below f0 above f0 1 kHz 781 Hz to 1279 Hz 219 Hz 279 Hz 16 kHz 14579 Hz to 17212 Hz 1421 Hz 1212 Hz At 1 kHz the bell is almost perfectly log-symmetric: f0²/f_lo = 1280 Hz against a measured 1279 Hz. At 16 kHz log symmetry would put the upper edge at 17,561 Hz and it actually lands at 17,212 Hz. The upper skirt has been pulled in by 350 Hz, squashed against Nyquist. This is the bilinear transform doing exactly what the theory says it does, and it is not a bug; it is the cost of the design, and now it has a number. The spectrum behind the curve shares an axis but not a scale. The overlay is computed like this: dB = dB - max(dB); % normalise to own peak dB = max(dB * 0.35, -25); % compress and floor That 0.35 is a display factor with no acoustic meaning. So the FFT fill sitting behind the response curve is on the same dB axis as the curve but is neither calibrated nor absolute; you cannot read a level off it. It is there to show spectral shape. You can watch the floor doing its work in the first screenshot: the blue fill flattens along the bottom of the plot at exactly −25 dB, and that flat line is the clamp, not the signal. The spectrogram gets the honest version, normalised to peak with a −60 dB floor and no compression. There is a fourth issue in the same function that is subtler. The FFT is 8192 points, giving a bin spacing of Fs/N = 5.38 Hz, and it is resampled onto a 300-point log-spaced display grid with interp1(..., 'linear') in linear magnitude. That resampling is doing two opposite and equally wrong things at once: Between 20 and 100 Hz there are about 69 display points drawn from only 14 FFT bins, so the bottom of the display is interpolation, not measurement. Between 10 and 20 kHz there are about 30 display points point-sampling 1857 bins, so 98% of the bins are simply never read. A narrow peak that lands between two sampled bins does not appear at all. The fix is standard and I have not done it: bin the FFT into the display bands and take the max or the RMS per band, instead of point-sampling. Point-sampling a spectrum is how analysers under-report high-frequency content. The rendering strategy, which came out of a performance bug The first version redrew the response the obvious way: cla, then re-plot. Every parameter change destroyed every graphics object and built new ones, and dragging a band was visibly laggy. Every graphics object is now created once in initPlot, and afterwards only ever updated: set(app.HEQCurve, 'YData', magdB) % not plot(...) again set(app.HSpecFill, 'YData', specdB) set(app.HBandDot{i}, 'XData', f, 'YData', g) Creating a MATLAB graphics object is expensive; updating one's YData is cheap. Same lesson as the DOM, and as ctx.fillRect versus rebuilding a canvas: the drawing was never the problem, the allocation was. The spectrogram uses the buffer form of the same trick, a 300 × 100 circular buffer with a head index, refreshed into a single imagesc handle. New columns overwrite the oldest, nothing is reallocated and nothing scrolls. There is also an offline path: computeFullSpectrogram() runs a proper STFT over the whole file and opens it as a 3-D surface. That one is allowed to be slow, because you asked for it explicitly and it renders once. Three timers, and the constraint they work around 80 ms, spectrum. An 8192-point Hann-windowed FFT of the current playback window, feeding both the overlay and one new spectrogram column. One transform, two consumers. 50 ms, debounce poll. Fires restartPlayback() 150 ms after the last parameter change. 10-second chunks, playback. startChunk() processes and plays ten seconds, chaining the next from the audioplayer object's StopFcn. The last two are one mechanism, and they exist because of a hard limit: MATLAB's audioplayer plays a buffer you hand it, and gives you no way to swap coefficients mid-buffer the way a real-time audio callback would. Changing the EQ therefore means re-filtering from the current position and restarting the player. Chunking makes that affordable. Re-filtering a five-minute file on every knob movement is unusable; re-filtering ten seconds is instant. The 150 ms debounce stops a knob drag, which emits dozens of callbacks a second, from queuing dozens of restarts. Drag freely, and the audio catches up 150 ms after you let go. It is scaffolding around a missing feature, but naming the constraint honestly is what produced the right structure: the buffer you can afford to recompute sets the latency you can offer. The mouse maths Grab a band directly on the response plot and move it. onAxesButtonDown picks the nearest band, onMouseMove maps cursor position to frequency and gain, onMouseUp clears it. Nearest is measured in normalised log-frequency, not in pixels or hertz: logDist = abs(log10(app.Freqs) - log10(clickF)) / (log10(20000) - log10(20)); [minD, i] = min(logDist); if minD > 0.2, return; end The denominator is 3 decades, so the 0.2 threshold is 0.6 decades, a factor of 3.98, almost exactly two octaves. You can grab a band from up to two octaves away and no further. The drag mapping is logarithmic in X and linear in dB in Y, so a centimetre of mouse travel near 100 Hz changes the frequency far less than the same centimetre near 10 kHz, which is what makes it feel right. The scroll wheel adjusts Q multiplicatively, ×1.2 or ÷1.2 per tick, clamped to [0.1, 10]. Additive steps would crawl at high Q and leap at low Q; a constant ratio feels identical everywhere in the range. Every gesture ends by calling updateEQ(), scheduleRestart() and markCustomPreset(), the last of which flips the preset dropdown to "Custom" the moment you touch anything, so the label never lies about what you are hearing. Presets Six built-ins plus user presets saved as .mat files into a presets/ folder created on first launch. Each preset is a struct of Freqs, Gains, Qs and Tipos. savePreset() validates the name against the built-ins, sanitises it for the filesystem, warns on overwrite and rebuilds the dropdown. Unremarkable, and it is what makes the thing usable rather than a demo. What MATLAB was good at, and what it was not I write most of my audio code in C++ with JUCE, so the comparison is worth making. The maths is the code. filter(b,a,x), conv, roots, zplane: the gap between the equation on the page and the line in the editor is nearly zero. Verifying the pole-radius formula against roots() took one line. Inspection is free. Being able to stop, print the actual coefficient vectors, convolve the cascade by hand and check H(+G)·H(−G) = 1 in the command window is a different debugging experience from attaching a debugger to a plugin host. Every verified number in this post came out of that. And what it was worse at is the thing the timers work around: there is no real-time audio callback. The chunking, the debouncing, the restarting are all scaffolding around that one absence. In JUCE the same feature is a coefficient update inside processBlock and nothing else. Which is the honest summary. MATLAB is where I would design and verify a filter. It is not where I would ship one. Takeaways The trig in the RBJ coefficients is a solved bilinear transform. cos(ω0) and sin(ω0) are what prewarping collapses into. Knowing that tells you the center frequency is exact and everything else is warped. A = 10^(G/40) because the peaking filter delivers A². The gain multiplies the numerator and divides the denominator. r = sqrt((1−α)/(1+α)) is the pole radius, exactly. It proves the form is unconditionally stable and it makes Q concrete: Q is how close the pole gets to the unit circle. Read the zeros. z = +1 is DC, z = −1 is Nyquist, and a double zero is 12 dB/octave. You can predict a filter's shape before plotting anything. Derive your own constraints and then check them numerically. The 35 Hz clamp had been in the code for months and I found it by solving the clamp for frequency, not by listening. Create graphics objects once, set() them forever. cla in a live update path is the performance bug. Source is on GitHub. I'm an audio DSP student at UPC. Most of my work is C++/JUCE (six VST3 plugins here), but this one started in MATLAB and stayed there.