Implementing various filters, using scipy

I edited your files as biquad and did just that. Look above.

That’s the part where a loop interacts less and performs four operations on each one. It may not work; it depends on a series of things. But I think a biquad filter could work in this case.

The “advanced” version of this technique would do precisely that, but it uses SIMD operations provided by the hardware and has been known for a long time. If you have a recent processor, you could eventually try that.

It works best with the SAME OPERATION on different data in parallel, always vectors, never something nested like trees, and always the same function for all.

It doesn’t need to be different lists; often, the data is sliced in chinks and processed like that.

C++ is a mess; it doesn’t do it for you (unless you find a library for your use case). Other languages, especially purely functional languages, have that as part of the compiler.

BTW, maybe this is more correct, i did a mistake I think. use this:

// 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) {

    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 tempPrevSample1 = prevSample1;
        float tempPrevSample2 = prevSample2;

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

        tempPrevSample1 = prevSample1;
        tempPrevSample2 = prevSample2;

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

        tempPrevSample1 = prevSample1;
        tempPrevSample2 = prevSample2;

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

        tempPrevSample1 = prevSample1;
        tempPrevSample2 = prevSample2;

        float filtered3 = ((input[i + 3] - a1 * tempPrevSample1 - a2 * tempPrevSample2) * b0) + (b1 * tempPrevSample1) + (b2 * tempPrevSample2);
        output[i + 3] = zapgremlins(filtered3);
        prevSample2 = tempPrevSample1;
        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) {
    ft = inTable;
    registerUnit<IIR::IIR>(ft, "IIR", false);
}

C++ is a thing of the devil

1 Like

Or maybe try a library? Maybe save this stuff for later. Just try the first one I sent.

Cool, thanks for the suggestions! I think they are currently not necessarily relevant to my own work, and I don’t think I need to go that deep down the c++ rabbit hole right now. But certainly super interesting.

To conclude, an uncomplicated utility to use Scipy. It has an outstanding module for digital signals, many things not easy to find


import sys
import json
import numpy as np
from scipy import signal

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def calculate_iir_coeffs(freqData, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order)
    return b, a

if __name__ == "__main__":
    input_json = sys.argv[1]
    params = json.loads(input_json)
    freqData = params["freqData"]
    cutoff = params["cutoff"]
    fs = params["fs"]
    order = params["order"]

    b, a = calculate_iir_coeffs(freqData, cutoff, fs, order)
    result = {
        "b": b.tolist(),
        "a": a.tolist()
    }
    print(json.dumps(result))

(
var output;
var result;
var jsonString;

var freqData = [0.1, 0.2, 0.3, 0.4, 0.5]; // example frequency data
var cutoff = 0.3; // normalized cutoff frequency (0.0 to 1.0, 1.0 being Nyquist frequency)
var fs = 50; // sampling frequency
var order = 5; // filter order
var pythonScriptPath = "~/scwork/python/iir_filter_coeffs.py".standardizePath;

var params = (
    freqData: freqData,
    cutoff: cutoff,
    fs: fs,
    order: order
);

jsonString = params.asJSON;

"JSON string: ".postln;
jsonString.postln;

p = Pipe.new("python " ++ pythonScriptPath ++ " '" ++ jsonString ++ "'", "r");

if (p.isOpen) {
    output = p.getLine;
    while ({output.notNil}) {
        output.postln;
        output = p.getLine;
    };
    p.close;
} {
    "Error: Could not open pipe".postln;
}
)
{ "cutoff": 0.3, "freqData": [ 0.1, 0.2, 0.3, 0.4, 0.5 ], "fs": 50, "order": 5 }

{"b": [2.240111257258062e-09, 1.120055628629031e-08, 2.240111257258062e-08, 2.240111257258062e-08, 1.120055628629031e-08, 2.240111257258062e-09], "a": [1.0, -4.878005218496593, 9.519429131062548, -9.289979563992423, 4.533698940535282, -0.8851432174252531]}

@smoge, fearing I may be heading off topic here…

atk-sc3’s Near-field Effect pseudo-UGens are built upon Ambisonic degree utilities, which are implemented as cascades of SOS and FOS UGens.

I have entertained the thought of putting together a simple separate quark that assembles the SOS (+ FOS) cascades from dictionaries of coefficients.

Other(?) additional quarks could provide the various recipes for higher order Butterworth (and other types of) filters, using this mechanism.

2 Likes

DSP Macros, I love it already!

EDIT: “requires suitably pre-conditioned input to avoid overflow”, that’s the part I forgot to do. There is an algorithm with the “forbidden” coefficients, and also “denormalization” to avoid tiny floating-point numbers and turn them to zero.

Create a new thread if you want to talk more about it, it’s quite interesting (or DSP chat room???)

EDIT2: That’s the sketch of the “checks” to avoid " Bad Coefficients" (this one I’m adapting from scipy, if you prefer python look for “bad coefficients”) . I think it’s kind of nice not to explode the sound))). Take a look, if you get those maybe you can use this too ( I also need to ‘denormalize’ the tiny floats too)

Correct, actually it’s been off the topic of multiband compressors for 20+ posts already.

High-order filters are a valid topic in themselves – good to continue the topic in its own thread.

hjh

1 Like

Thank you for organizing the topics, James. good call

I’ve tried out the “move posts” functionality. Feel free to edit!

2 Likes

Friendly reminder: if you realize that the thread has diverted significantly from its original topic, consider “Reply as linked topic” (the arrow icon on the top left in the reply text box), so that the moderators don’t have to do the work and split up the threads after the fact.

1 Like