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.

  • 2
    $\begingroup$ It's not clear what your actual question is -- in the title you say you want a variable IIR filter, but in the text you don't mention IIR filtering except a vague reference to Chebychev and/or Butterworth. You never say how you want the filter to vary. Could you edit your question, ending it with a clear request for information (I suggest a one or two sentence paragraph, followed by a question mark). If the main thrust of the first 13 paragraphs of your question are "I tried this with FIR filters", then perhaps condense that part -- to that statement even. $\endgroup$
    – TimWescott
    Jul 24 at 16:55
  • $\begingroup$ Fair enough. I don't know how to express what "what's the mathematics?" is in the context of finding but not understanding the mathematics (this was the same problem I had with FIR filtering and phase modulation synthesis; once I got it after a few months of reading stuff, it was like oh, this is just a few multiplications and simple algebra). Probably "this is how I'm implementing some other filter" is too much context about where I am and how I'm thinking at the moment, and so unclear? $\endgroup$
    – John Moser
    Jul 24 at 17:28
  • $\begingroup$ I've rewritten it some and moved the context for how I approached the problem with FIR into a separate question which I self-answered. That should be more newbie-friendly anyway (assuming my math is right) if anyone else is trying to solve the same problem in FIR that I was there. $\endgroup$
    – John Moser
    Jul 24 at 18:17

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.

  • $\begingroup$ I was trying to do it in fixed-point and concluded that this was terrible, so decided to use 24-bit floating point (18-bit significand, bias to 57, so I can go from like 64 down to $7\times 10^{-18}$ I think?). Sinc was easy to implement and I was using it to interpolate (I have wavetable synthesis, too, plus internal lower sample rate than I want to use when dithering output). If you have a better filter that's easier to implement I'm good with that; Chebyshev was just what I saw in the DSP book and elsewhere with other people implementing filters in software. $\endgroup$
    – John Moser
    Jul 24 at 16:06
  • $\begingroup$ I edited the question to reflect some of this and refine it into something more sensible based on this feedback. $\endgroup$
    – John Moser
    Jul 24 at 16:11
  • 1
    $\begingroup$ OK, John, there is the Audio EQ Cookbook that can be used to define each Second-Order Section (SOS) that will go in series. And for Butterworth, if the order is greater than 2, all SOS have the same resonant frequency but different values for Q that are well defined. I can look up the formula if you want. Additionally, if you do this in fixed point, you might need to have a word width of at least 24 bits for data and coefficients. For the coefficients, you will need 2 or maybe 3 bits to the left of the binary point. $\endgroup$ Jul 24 at 17:42
  • 1
    $\begingroup$ also, if this is for time-variant filters for musical use, such as a sweeping resonant frequency, or some other sweeping parameter, you should use the Lattice Filter topology (also this) (whether fixed or floating point) and for fixed point, maybe the normalized ladder topology. Lastly, you might want to learn something about noise shaping in signal quatizers. $\endgroup$ Jul 24 at 17:47
  • 1
    $\begingroup$ Lastly, if FPGA architecting and programming is your forte and DSP math and audio effects is less so, you might want to consider contracting with an audio DSP jock to help you. I can nominate @Hilmar. There are others I can nominate. $\endgroup$ Jul 24 at 17:54

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.

  • $\begingroup$ Interesting. My issue was largely that when I started adding up numbers, many of my coefficients had one figure if that above $10^{-6}$. I can count that up in floating point, because I can accumulate numbers that are a million times smaller, but things that add up to the LSB of an integer would add up to 0 instead. Basically I thought I was supposed to feed back quantization error, but it looked like I'd have what amounts to .4 + .4 + .4 + .4 +.4 +.4 + .4 = 0.0 going on. I considered using TWO fixed point formats, one just starting several digits below zero…for like 30 seconds. $\endgroup$
    – John Moser
    Jul 24 at 19:57
  • $\begingroup$ Also yes, I have one block that does a thing a thousand times. My approach with FIR is essentially stateless (the cutoff and sample to process are given, and it spits out the results for that one sample). IIR gonna need to hold some state for feedback. I'll check for that book; lots of good recommendations today. $\endgroup$
    – John Moser
    Jul 24 at 20:05
  • $\begingroup$ As to coefficient sizing for fixed-point FIR: you normalize everything so that you get the precision you need, do the multiply-accumulate, then adjust the magnitude after the fact. And you generally used fractional arithmetic, or at least fixed-point -- i.e., signed 0x1ffff may represent $2 - 2^{16}$, not $2^{17} - 1$. $\endgroup$
    – TimWescott
    Jul 24 at 20:15
  • $\begingroup$ True. I could just use large accumulators but small inputs. Also my options for multipliers are probably more like $18\times 18$ or $34\times 34$; a $36\times 36$ requires 4 $18\times 18$ instead of the 3 needed for $34\times 34$. $\endgroup$
    – John Moser
    Jul 24 at 20:41
  • $\begingroup$ You may not need more than 18-bit coefficients -- it depends on what frequencies you go down to, and how concerned you are that the ripple of the filter isn't exactly right. For 16-bit audio in, you need more than 18, but 32 would be more than plenty -- the old Analog Devices 56000 DSP chips used 24-bit data paths and they were just fine. $\endgroup$
    – TimWescott
    Jul 24 at 21:03

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?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.