What is a simple method, if any, to calculate a variable IIR filter for audio effects in real time?

Question

Question edited based on first answer which raised a lot of good points about what problem I'm trying to solve. Also reordered some things to stick context of how I did this with FIR so far at the end, to make it more coherent.

Sorry if I'm not entirely sure what, directly, is the right question to ask. The cone is a bit wide.

I've been trying to implement variable audio filters on FPGA for an (unreasonable) synthesizer chip, allowing low-pass, high-pass, and band-pass with movable cutoff frequencies. The cutoff frequency might move around constantly with an LFO, so I recalculate the filter for every sample output; at the clock rate $48000Hz\times 288\times 12$ and the sample rate $40960Hz$, I have roughly 4,000 clock ticks to work with.

I've figured out how to do it with FIR 83-tap sinc filters (detailed here) with the coefficients calculated in real time, able to filter about 500 samples per $40960Hz$ sample-tick via a pipeline that requires 3 $18\times 18$ multipliers, running 6 of these pipelines in parallel. That was fun.

Now I want to figure out how to do it with an IIR filter like a Type-1 Chebyshev (based on some material I'd read) or a Butterworth. I've not yet been able to penetrate the math.

Internally I'm using 24-bit floating point with an 18-bit significand (17 bits plus 1 hidden bit), because I have 18×18 multipliers. I'm looking at a bias of (if this is correct—still working out how IEEE754 works) 57 to give a range of about 64 to about $7\times 10^{-18}$ positive and negative. I looked at doing this in fixed point, and started having trouble with not having enough precision to avoid piling up massive amounts of error over many additions and multiplications.

From what I've gathered, 6-pole seems (more than?) fine, or 3-pole and then run it back through itself, so this shouldn't be very intensive to use. I'm not sure about calculating it because I haven't been able to grasp precisely how to do that. Yet. Also not sure if I can combine IIR filters the same way as FIR sinc filters, which is just done by adding the kernels; inverting a low-pass to a high-pass or a low-cut to a band-cut is done by flipping the signs and adding 1 to the midpoint. All nice and easy stuff for the FIR approach, but I don't know if these principles translate to recursive filters.

Besides not understanding the math (I didn't understand the math behind a sinc filter until I read some code that does it and some other basic material), I'm not entirely sure how to do imaginary calculations on a real FPGA (or on paper). Based on some responses here so far I've found more stuff related to Butterworth filters and have started trying to unpack what that means.

Hilmar mentioned doing all that ahead of time and then moving the cutoff frequency by frequency warping; no idea how to do that, still asking Google. It sounds like a lower-cost approach (doesn't need 8,000 multipliers, doesn't take 5 years to finish calculating…) that might fit my needs…if I can find resources explaining that, or figure out how to actually apply the principle.

I guess what I need is:

1. An approach to generating a filter appropriate for audio (synthesis) whose cutoff I can modify in a few calculations, or at least in a direct enough process that I can pipeline the calculations.
2. The math behind modifying it so I can work from there to design said pipeline.
3. A math degree I don't have!

Did it one way, now want to do it another way.

Answer 1

If I understand correctly you want to implement a low pass filter with a programmable cutoff frequency, but I'm not entirely sure.

A few things to consider:

1. FIR filters are not a great fit for audio. With a window length of 82 you will not be able do anything meaningful at low frequencies.
2. IIR is indeed the better choice for typical audio frequencies.
3. Why Chebyshev ? Butterworth is mathematically far easier and has nicer phase properties. If you really need Chebyshev: which one ? There are Type I and Type II and they are rather different animals.
4. Are you implementing this in fixed point or in floating point? IIR in fixed point are VERY complicated: you need to figure out section scaling and gains, section topology, section order, clipping prevention, transfer functions to all state variables, rounding behavior, limit cycles, etc. If you are lucky you can just throw extra bits at the problem, but in general it requires heavy math and analysis especially for a variable filter.

My recommendation would be to a Butterworth filter of the right order (which is probably fairly low) in floating point math. The poles are equally distributed on half circle on the left side of the s-plane and all the zeros are at $s=\infty$ which maps to $z=-1$ in the z-plane. You can simply design a prototype filter for your desired order and than warp it to the desired cutoff frequency.

Come to think of it, that approach actually works for any IIR low-pass filter, you can store a prototype filter and than simply apply frequency warping to adjust the cut off frequency.

Answer 2

Internally I'm using 24-bit floating point ... I looked at doing this in fixed point, and started having trouble with not having enough precision ...

If you don't have enough precision with 18-bit fixed point, floating point with an 18-bit mantissa won't help. Floating point would help if you had some small number -- call it $\epsilon$ -- that you needed to know with sufficient precision by itsel. But for filtering, you essentially need to distinguish $1$ from $1 \pm \epsilon$. For that, floating point with an $n$-bit mantissa gains you nothing over fixed-point math with careful scaling.

I suggest you do a compare/contrast between the amount of logic you need to implement your floating point vs. the amount of logic you'd need to implement 36-bit wide data paths inside the filters (note that you don't need 36 bits -- 24 would probably do -- but you have 18 or 36, so...).

For generating the filter, I'd suggest that you implement a microprocessor in your FPGA (or move the coefficient calculations to an external micro, if you already have one). You probably don't need to be changing the filter coefficients each step, and if you do you can do some linear interpolation between points if the cutoff frequencies aren't changing much.

You don't mention it, but you want to use filters separated into blocks of no more than 2nd-order. This cuts down on the size of the data path you need.

Note that IIR filters suffer a speed hit in an FPGA, because they are little feedback loops. To get the maximal speed out of the thing you need to pipeline, but that messes up the filter behavior. For a beginner, you probably just want to hold your clock rate down -- but if you really want to spend time dotting 'i's and crossing 't's, you can probably make one IIR filter engine that loops through all your channels and all your filter stages, to essentially take in a vector of inputs and spit out a vector of outputs at each sample time of your audio system.

For modifying the filter, as long as you're not changing the filter parameters too fast, or in too large of jumps, you can just store them in registers and write them as necessary. You probably want to change all of the parameters in a filter in one go.

You may want to look for a book on DSP. My recommendation is "Understanding Digital Signal Processing" by Richard Lyons. It may still be over your head, but it's the most accessible DSP book I know of.

Answer 3

Think I have an answer, from this answer on StackOverflow found after digesting some of the answers and comments here and getting a better idea of what I'm looking for. Also helpful was PySDR and stuff on DSP Related.

As with the FIR approach, use $c=\dfrac{\text{Cutoff Frequency}}{\text{Sample Rate}}$. Hard-code $\sqrt{2}$ and a $tan()$ look-up table.

$\begin{align} ita &= \dfrac{1}{\tan\left(\pi c\right)} \\ b_0 &= \dfrac{1}{1+ita\sqrt{2}+ita^2} \\ b_1 &= 2b_0 \\ b_2 &= b_0 \\ a_1 &= 2b_0\times\left(ita^2-1\right) \\ a_2 &= -b_0\times \left(1-ita\sqrt{2}+ita^2\right) \end{align} $

To make this a bit more amenable to a hardware implementation, ensure the $\tan()$ table is indexed in terms of multiples of $\pi$ (skipping the $\pi\times c$ calculation) and do something like below:

$\begin{align} ita &= \dfrac{1}{\tan\left(\pi c\right)} \\ itsq &= ita^2 \\ itsr &= ita\sqrt{2} \\ b_0 &= \dfrac{1}{1+itsr+itsq} \\ b_1 &= b_0 \ll 1\\ b_2 &= b_1 \\ a_1 &= \left(b_0\times itsq - b_0\right) \ll 1\\ a_2 &= b_0\times itsr - b_0\times itsq - b_0 \end{align} $

The partial calculation $\left(b_0\times itsq\right)$ for $a_1$ can be stored, shifted right one place, and subtracted when producing $a_2$. DOI 10.1109/SECON.1993.465659 also gives a $O\left(log N\right)$ algorithm to calculate the $N$ significant bits of $\dfrac{1}{x}$ as a special case, much faster than a divider. The bit shifts can be carried out by assigning to registers simultaneously, just offset; and there's no reason to actually store $b_1$ or $b_2$…except to support high-pass.

To produce a high-pass, as per this answer and stuff in the audio EQ cookbook, further perform:

$ \begin{align} b_0 &= b_0\times ita^2 \\ b_1 &= -b_1\times ita^2 \\ b_2 &= b_2\times ita^2 \end{align} $

Now $b_1$ must actually be computed and stored separately.

If I've got this right, then all I need to do is implement a biquad direct 1 filter topology or other appropriate topology—the $z^{-1}$ blocks are just registers to delay until the next sample, and biquad direct 1 can't overflow if using fixed-point integer instead of floating point—and we're done here. For high order, send the output to the input of another identical filter.

With the single-cycle multipliers in FPGA DSPs, $ita^2$ and $ita\sqrt{2}$ can be calculated easily, parallel or serial. The reciprocals take a little more time (for 32 bit integers, apparently 5 cycles). $b_1$ and $b_2$ take literally no calculations. $a_1$ and $a_2$ share a multiply, saving time. Can be done in 15 clock cycles, which is just under what I have available to perform 256 of these per sample output. It's also possible to only keep $ita$ and $b_0$ and recompute those every 4 samples (staggered) (0.1ms) regardless of constant changes to the filter cutoff (e.g. from an LFO driving vibrato and other effects).

I believe this is faster to calculate than a FIR windowed sinc filter with a narrow pass-band, cheaper to calculate (I need a lot of multipliers to get the coefficients for FIR in parallel or I can't keep up), and cheaper to apply (same). It's more flexible. I can use bigger multipliers and intermediate bit depth. It requires a little more state, but not much. The delay shrinks…a lot, actually…wow.

Should be right. Maybe?