Arbitrarily high-order IIR filters

My little application allows one to change coefficients at runtime Using shared pointers and atomic operations, including changing the order of the filter. It’s slightly different than a sc-plugin, with 64-bit signal processing and single-sample operation, but that’s easy to change. I’m not sure there is interest in turning it into a plugin. Maybe not.

Well, as I read the documentation for LTI, I see that it takes one buffer for feed-forward coefficients, and a second buffer for feedback coefficients.

So LTI is an IIR.

At this point, I don’t quite understand what is the disconnect here. The LTI documentation is very clear as to the formula being implemented. “Different type of filter” just isn’t true.

This isn’t the case, if buffers are used to supply the coefficients – as LTI already does.

It’s really quite strange to me – “here’s a thing that exists” and then the conversation goes on as if it doesn’t exist. Anyway, I won’t belabor the point; just read the help and the source code.

hjh

James, I think the documentation uses feedback terms involving previous output values. This feedback means the filter’s impulse response is theoretically infinite. At the same time, it’s a specific type of LTI system because it satisfies both linearity and time invariance. However, the documentation is not precise enough to explain not all LTI systems are IIR filters.

my implementation is different; coeffs can change, including their number, and processing is 64bit, single-sample processing. Also, it uses TDF-II for memory efficiency. For example, maybe a different precision might make a difference for that kind of processing.

and it’s ok it something similar already exists, right? it’s not a problem, even if I fail

I’m pretty sure James is right. The transfer equation in the documentation is the same the transfer equation of an IIR, unless I’m missing something. What is the difference?

Also - and I’m happy to be proven wrong - I doubt that whatever minuscule differences can be found in memory efficiency or whatever will actually make a difference for the sound or usability of a multiband compressor, right?

My bad - I skimmed over your post and didn’t actually check out the LTI documentation before commenting, taking Bernardo’s word for it that the filters are different.

Yes, the transfer (the equation precisely, not in the text) is recursive (The SC documentation equation is an IIR filter). Still, as I said before, the documentation could be more complete.

Again, the IIR filter is a specific type of LTI.

TDF-II is significantly more memory efficient, requiring half the delay elements than TDF-I.

My implementation can change coeffs and order at runtime without glitches or data races. Since it’s recursive, I thought 64-bit could be something to TRY.

I almost feel I need to apologize for coding this. It’s a simple filter, friends, why the troubles?

Things can be different and not be a dichotomy. A sc help file is not the canonical reference.

(I’m silencing this thread)

@smoge I apologise, I did not mean to disparage you in any way; my comment was purely a clarification.

I certainly don’t have a problem with you coding filters. I do think that, where there is potential misinformation on the forum, it’s important to clear up misunderstandings to help people learn.

1 Like

Well, so let’s read an actual book. In theory, those words are not the same.

However, if you see the code, you can do similar things with them (it seems that is the case).

However, the implementation is not quite the same regarding real-time safety, seamless coefficients without interrupting the audio stream, and TDF-II, which reduces memory usage and improves numerical stability. One implementation is more optimized for real-time audio processing.

Where is that SC community always at the edge? No more?

Have you considered unrolling the loop to test performance? I think biquad is a good candidate, I would try myself.

// IIR.sc
IIR : UGen {
    *ar { |input, a1, a2, b0, b1, b2|
        ^this.multiNew('audio', input, a1, a2, b0, b1, b2);
    }

    checkInputs {
        ^this.checkValidInputs;
    }
}





// PluginIIR.hpp
// Jordan White (jordanwhitede@gmail.com)

#pragma once

#include "SC_PlugIn.hpp"

namespace IIR {

class IIR : public SCUnit {
public:
    IIR();

private:
    // Calc function
    void next(int nSamples);

    // Member variables 
    float prevSample1;
    float prevSample2;
};

} // namespace IIR




// PluginIIR.cpp
// Jordan White (jordanwhitede@gmail.com)

#include "SC_PlugIn.hpp"
#include "IIR.hpp"

static InterfaceTable* ft;

namespace IIR {

IIR::IIR() {
    mCalcFunc = make_calc_function<IIR, &IIR::next>();
    prevSample1 = 0.0;
    prevSample2 = 0.0;
    next(1);
}

void IIR::next(int nSamples) {

    // args
    const float* input = in(0); 
    const float a1 = in0(1);
    const float a2 = in0(2);
    const float b0 = in0(3);
    const float b1 = in0(4);
    const float b2 = in0(5);

    float* output = out(0);

    int i = 0;


    for (; i < nSamples - 4; i += 4) {
        float filtered0 = ((input[i] - a1 * prevSample1 - a2 * prevSample2) * b0) + (b1 * prevSample1) + (b2 * prevSample2);
        output[i] = zapgremlins(filtered0);
        prevSample2 = prevSample1;
        prevSample1 = output[i];

        float filtered1 = ((input[i + 1] - a1 * prevSample1 - a2 * prevSample2) * b0) + (b1 * prevSample1) + (b2 * prevSample2);
        output[i + 1] = zapgremlins(filtered1);
        prevSample2 = prevSample1;
        prevSample1 = output[i + 1];

        float filtered2 = ((input[i + 2] - a1 * prevSample1 - a2 * prevSample2) * b0) + (b1 * prevSample1) + (b2 * prevSample2);
        output[i + 2] = zapgremlins(filtered2);
        prevSample2 = prevSample1;
        prevSample1 = output[i + 2];

        float filtered3 = ((input[i + 3] - a1 * prevSample1 - a2 * prevSample2) * b0) + (b1 * prevSample1) + (b2 * prevSample2);
        output[i + 3] = zapgremlins(filtered3);
        prevSample2 = prevSample1;
        prevSample1 = output[i + 3];
    }

    // Process remaining samples
    for (; i < nSamples; ++i) {
        float filtered = ((input[i] - a1 * prevSample1 - a2 * prevSample2) * b0) + (b1 * prevSample1) + (b2 * prevSample2);
        output[i] = zapgremlins(filtered);
        prevSample2 = prevSample1;
        prevSample1 = output[i];
    }
}

} // namespace IIR

PluginLoad(IIRUGens) {
    // Plugin magic
    ft = inTable;
    registerUnit<IIR::IIR>(ft, "IIR", false);
}

1 Like

A new implementation, especially if more efficient or precise (since higher-order filters depend on precision), would be very welcome!

It sounded to me like you were saying that LTI is a “different type of filter” from an IIR, hence not an IIR. Probably some confusion over terminology,

In any case, you’re now saying there would be benefit from a different (TDF-II) implementation of higher-order filtering. No dispute there! I was never saying “don’t write a new one” – I was only saying that the LTI UGen as it now exists can do IIR filtering, as it’s a general-case recursive filter. That’s all. I don’t see any room for disagreement about that. When you were saying that LTI as it is today is not an optimal implementation, that is a/ true and b/ a rather different statement from what it seemed like you were saying.

hjh

2 Likes

No Worries. I got confused because LTI is a broader term (often emphasizing one aspect over another, although they are all the same kind), but after seeing the equation and the code, I realized one could do the same things with either of them. I should have checked it out first.

I’m not quite sure what this means sorry. I’m pretty new to c++ and ugen programming…

10 posts were split to a new topic: Implementing various filters, using scipy

Forgive me if I’m missing the point or if this comment is silly or unproductive…
I know there was some mention at some point about Biquad filters, but I’m not finding that in the linked thread, so I don’t really know what else was said or by who about it, but

Whyquad? (sorry, not sorry :smirk:)

  • Stability/Quantization Error Accumulation: Higher-order filters are generally more complex, involving more coefficients in a more complicated structure. This complexity can lead to potential stability issues, especially in higher frequencies or when the filter specifications are stringent. Biquads are simpler and inherently more stable when designed correctly, plus the poles and zeros of each section are decoupled. Higher order filters continually accumulate quantization error in coefficient calculation over a long feedback path.
  • Computational Efficiency: Implementing a single higher-order filter I believe is more computationally intensive due to the larger number of heavy operations required per sample, whereas cascading biquad filters divides the overall task into smaller, more manageable chunks. Again, I’m not certain about this, but I would be inclined to think that this could be handled rather elegantly with SIMD.
  • Frequency Response: My experience has been that it’s more difficult to get the intended frequency response from a higher order filter without oscillation/overshoot/instability/noise.
  • Phase Response: Higher-order filters can introduce significant phase distortion. Biquad cascades also introducing phase shifts, but can be configured (e.g., using Linkwitz-Riley) to maintain a more coherent phase response across the cascade.
  • Modularity/Flexibility: You can add/remove stages in a biquad cascade, but you have to redesign a higher order filter to achieve the same result.

I also can’t help thinking that this is exactly what FOS and SOS are for…
Why not make a pseudo ugen wrapper for these? Or better still, why not instantiate them as individual synths? That’s what the biquad filter suite does (the B prefix ugens). I think you could actually still just make this an SC wrapper that manages an arbitrary number of FOS/SOS Synths.

1 Like

Of course, transfer functions can be rewritten as a chain of simpler filter segments, and the biquad filter, in particular, has many practical advantages. Theoretically and mathematically, any of these filter segments can be implemented in any order without altering the overall filter behavior.

I haven’t encountered any differing opinions on this. :blush:

Some methods, such as the TDF-II (Transposed Direct Form II), can reduce memory consumption by half. My implementation of high-order filters is exceptionally efficient, using minimal memory and CPU resources (0.3% or less).

There are also some clever techniques that can be employed. For instance, while it may not be practical to change the coefficients of more complex filters in real-time due to computational costs just to calculate the new coeffs (!!!), combining them with an all-pass filter allows for very straightforward and cost-effective adjustments, with very drammatical effects. It sounds great.

They are mathematically equivalent, but the calculation is more complicated.

I want to make clear that my intention here is experimentation and exercise in technique and optimization.

Example: my original idea was to create a C++ class template that could be compiled to any order. It is just a question of simplifying the code that generates filters. That’s why I wrote the code in the first place…

BiQuad filters also can benefit from SIMD operation and loop unrolling. I think I posted an example here.

Yes, that’s the theory when we consider just the equations. I believe we’ve mentioned it once or twice here.

You are incorrect regarding performance. A higher-order filter can potentially have fewer coefficients than multiple cascaded biquads, which can reduce computational complexity. Additionally, there are techniques that can optimize memory usage.

Higher-order filters can achieve better numerical stability when implemented correctly with the necessary attention. Also, it’s a very different kind of flexibility compared to complicated cascading Biquad filters that you would need to change at once. It all depends on what you’re doing. For regular commercial use, biquad is a rule.

But let’s put them to the test. Why not? It’s the same c++ template code that generates all of them)))) that was the idea (in the beginning…)

1 Like

Nice resource: Filter Playground

Sorry if my response was redundant, I’m guessing there was a lot more discussion around this in a separate post that I missed.
I’m very interested in your design though.

Thanks for pointing this out!

Filter design is something I’m continually trying to wrap my tiny brain around.

1 Like

Yes, me too. There are so many things when those things are put into practical code… and the hardware at hand…

I learned that software is not “the” platform; there are many details to pay attention to… and none of the things in the software were designed with digital filters in mind, so in the end, you have to figure out everything by yourself…

I posted. It’s the gist link above. There is also a Haskell version, with around 50 lines of code, including imports. Just like C++, one needs to know what works better for real-time dsp. In the case of Haskell, strictness is better than laziness. Stream fusion to eliminate intermediate data structures during list and vector processing. Unlike C++, all those optimizations require less coding. SIMD in C++, you have to rewrite the loops by hand, and the resulting code is not beautiful.

And using a vector module is a kind of rule (Data.Vector.Storable) for audio. This last detail is also helpful for interoperating with C.

C++ has an experiment simd library if you’re looking for one, one day it might make it into the standard.

1 Like

That’s good to know.

Although, as that talk by SPJ I sent you, taking care of all cases is not simple, and it should be compiler work. Not just simple vectors, but more cases including nested/tree-like types.