bassFilter.cpp

/**
   Calculate the analog filter and the bilinear transform
   of the Bass and EQ section.

   SPDX-License-Identifier: MIT

   Compile
       c++ -O2 -g -o bassFilter bassFilter.cpp

   Run
       ./bassFilter [bass] [treble] > bassFilter.dat
           with bass and treble values in [0;1]

    Plot
       gnuplot
       plot "bassFilter.dat" u 1:(20*log10($2)) w lines t "Analog", "" u 1:(20*log10($4)) w lines t "Digital"
*/

#include <vector>
#include <complex>
#include <cmath>
#include <cstdio>
#include <cstdlib>

using R = double;
using C = std::complex<R>;
static constexpr C j{0, 1};

struct TF {
    std::vector<R> num;
    std::vector<R> den;
};

TF bassEQfilterAnalog(R bass, R treble)
{
    R B0 =
        0.0;
    R B1 =
        R(10000.0L)*bass - R(11000.0L);
    R B2 =
        (R(4722366482869644988655.0L/147573952589676412928.0L)*(bass*bass))
        - (R(3149228148263694434375.0L/147573952589676412928.0L)*bass)
        - R(598633738680022361533287.0L/36893488147419103232000.0L);
    R B3 =
        - (R(1.0L/100.0L)*(bass*bass)*treble)
        + (R(95513923305516444051222397770391.0L/2028240960365167042394725128601600.0L)*(bass*bass))
        - (R(1.0L/100.0L)*bass*(treble*treble))
        + (R(1.0L/50.0L)*bass*treble)
        - (R(11855575473574492365904017309744059.0L/253530120045645880299340641075200000.0L)*bass)
        + (R(11.0L/1000.0L)*(treble*treble))
        - (R(87.0L/10000.0L)*treble)
        - R(14370594264427299590777634582379451.0L/2535301200456458802993406410752000000.0L);
    R B4 =
        - (R(7486212072260646522045135498046875.0L/340282366920938463463374607431768211456.0L)*(bass*bass)*(treble*treble))
        + (R(4491727243356387913227081298828125.0L/340282366920938463463374607431768211456.0L)*(bass*bass)*treble)
        + (R(43930290233957032614896735927734375.0L/2722258935367507707706996859454145691648.0L)*(bass*bass))
        + (R(7486212072260646522045135498046875.0L/340282366920938463463374607431768211456.0L)*bass*(treble*treble))
        - (R(4491727243356387913227081298828125.0L/340282366920938463463374607431768211456.0L)*bass*treble)
        - (R(43930290233957032614896735927734375.0L/2722258935367507707706996859454145691648.0L)*bass)
        + (R(1497242414452129304409027099609375.0L/680564733841876926926749214863536422912.0L)*(treble*treble))
        - (R(2096139380232981026172637939453125.0L/1361129467683753853853498429727072845824.0L)*treble)
        - R(3988653792100472466945648193359375.0L/5444517870735015415413993718908291383296.0L);
    R B5 =
        (R(55263030661574758589267730712890625.0L/11417981541647679048466287755595961091061972992.0L)*bass*(bass - R(1.0L))*(- R(1000.0L)*(treble*treble) + R(700.0L)*treble + R(333.0L)));
    R B6 =
        0.0;

    R A0 =
        R(1120000.0L) - R(100000.0L)*bass;
    R A1 =
        (R(827159382962765777100875.0L/73786976294838206464.0L)*bass)
        - (R(97398808709186431407275.0L/73786976294838206464.0L)*(bass*bass))
        + R(22486138303994174387839511.0L/3689348814741910323200.0L);
    R A2 =
        (R(1.0L/10.0L)*(bass*bass)*treble)
        - (R(554005154673064069875538171729.0L/12379400392853802748991242240.0L)*(bass*bass))
        + (R(1.0L/10.0L)*bass*(treble*treble))
        - (R(11.0L/5.0L)*bass*treble)
        + (R(11172233007643856170213158767360603.0L/198070406285660843983859875840000.0L)*bass)
        - (R(28.0L/25.0L)*(treble*treble))
        + (R(2097.0L/1000.0L)*treble)
        + R(4803318945294176732167021742106863.0L/495176015714152109959649689600000.0L);
    R A3 =
        (R(395912635463280607188582045966881.0L/324518553658426726783156020576256000.0L)*(bass*bass)*(treble*treble))
        - (R(924536224793224579754221815612658771.0L/83076749736557242056487941267521536000.0L)*(bass*bass)*treble)
        - (R(507040844370222621749497226472328135711.0L/10384593717069655257060992658440192000000.0L)*(bass*bass))
        - (R(931846978770041617107737440612658771.0L/83076749736557242056487941267521536000.0L)*bass*(treble*treble))
        + (R(877514784442333180710841126228898003.0L/41538374868278621028243970633760768000.0L)*bass*treble)
        + (R(4083733035316798552790018918194514614443.0L/83076749736557242056487941267521536000000.0L)*bass)
        - (R(3611720156421957836512892203189315491.0L/1038459371706965525706099265844019200000.0L)*(treble*treble))
        + (R(47940277510652695890303360485452680333.0L/8307674973655724205648794126752153600000.0L)*treble)
        + R(1936978606034868499295094264265095385289.0L/519229685853482762853049632922009600000000.0L);
    R A4 =
        (R(38827074815011218430459499359130859375.0L/1393796574908163946345982392040522594123776.0L)*(bass*bass)*(treble*treble))
        - (R(51068199162681861845266819000244140625.0L/1393796574908163946345982392040522594123776.0L)*(bass*bass)*treble)
        - (R(99107653864029446716585253626748046875.0L/11150372599265311570767859136324180752990208.0L)*(bass*bass))
        - (R(38827074815011218430459499359130859375.0L/1393796574908163946345982392040522594123776.0L)*bass*(treble*treble))
        + (R(51068199162681861845266819000244140625.0L/1393796574908163946345982392040522594123776.0L)*bass*treble)
        + (R(99107653864029446716585253626748046875.0L/11150372599265311570767859136324180752990208.0L)*bass)
        - (R(1585472873686109630191326141357421875.0L/696898287454081973172991196020261297061888.0L)*(treble*treble))
        + (R(16362191671670370027673244476318359375.0L/5575186299632655785383929568162090376495104.0L)*treble)
        + R(54095420370500304200094670353834375.0L/696898287454081973172991196020261297061888.0L);
    R A5 =
        (R(1833247147700098857916891574859619140625.0L/365375409332725729550921208179070754913983135744.0L)*(bass*bass)*(treble*treble))
        - (R(9454388675080322685949504375457763671875.0L/1461501637330902918203684832716283019655932542976.0L)*(bass*bass)*treble)
        - (R(500354302003399815320782081071474609375.0L/2923003274661805836407369665432566039311865085952.0L)*(bass*bass))
        - (R(1833247147700098857916891574859619140625.0L/365375409332725729550921208179070754913983135744.0L)*bass*(treble*treble))
        + (R(9454388675080322685949504375457763671875.0L/1461501637330902918203684832716283019655932542976.0L)*bass*treble)
        + (R(500354302003399815320782081071474609375.0L/2923003274661805836407369665432566039311865085952.0L)*bass)
        - (R(3183150566106706671416759490966796875.0L/1461501637330902918203684832716283019655932542976.0L)*(treble*treble))
        + (R(3183150566106706671416759490966796875.0L/1461501637330902918203684832716283019655932542976.0L)*treble)
        + R(840351749452170561254024505615234375.0L/11692013098647223345629478661730264157247460343808.0L);
    R A6 =
        -(R(117489690137807889841496944427490234375.0L/24519928653854221733733552434404946937899825954937634816.0L)*bass*(bass - R(1.0L))*(- R(1000.0L)*(treble*treble) + R(1000.0L)*treble + R(33.0L)));

    TF tf;
    tf.num = {B0, B1, B2, B3, B4, B5, B6};
    tf.den = {A0, A1, A2, A3, A4, A5, A6};
    return tf;
}

C calcTF(const TF &tf, C s)
{
    size_t nn = tf.num.size();
    size_t nd = tf.den.size();

    C hn = 0;
    C hd = 0;

    C spowi = 1;
    for (size_t i = 0; i < nn; ++i) {
        hn += tf.num[i] * spowi;
        spowi *= s;
    }
    spowi = 1;
    for (size_t i = 0; i < nd; ++i) {
        hd += tf.den[i] * spowi;
        spowi *= s;
    }

    return hn / hd;
}

TF normalizeTF(TF tf)
{
    const size_t nn = tf.num.size();
    const size_t nd = tf.den.size();

    if (nd > 0) {
        const R a0 = tf.den[0];
        for (size_t i = 0; i < nn; ++i)
            tf.num[i] /= a0;
        for (size_t i = 0; i < nd; ++i)
            tf.den[i] /= a0;
    }

    return tf;
}

void printTF(const TF &tf, FILE *out)
{
    fprintf(out, "B = [");
    for (size_t i = 0, n = tf.num.size(); i < n; ++i)
        fprintf(out, "%s%e", i ? " " : "", (double)tf.num[i]);
    fprintf(out, "];\n");

    fprintf(out, "A = [");
    for (size_t i = 0, n = tf.den.size(); i < n; ++i)
        fprintf(out, "%s%e", i ? " " : "", (double)tf.den[i]);
    fprintf(out, "];\n");
}

TF bilinearTransformOrd6(const TF &inTF, R fs)
{
    const R pi = M_PI;
    //const R k = 1/std::tan(pi*fc/fs);
    const R k = 2*fs;

    R
        b0 = inTF.num[0], b1 = inTF.num[1], b2 = inTF.num[2], b3 = inTF.num[3],
        b4 = inTF.num[4], b5 = inTF.num[5], b6 = inTF.num[6];
    R
        a0 = inTF.den[0], a1 = inTF.den[1], a2 = inTF.den[2], a3 = inTF.den[3],
        a4 = inTF.den[4], a5 = inTF.den[5], a6 = inTF.den[6];

    TF outTF;
    outTF.num = {
        b6*(k*k*k*k*k*k) + b5*(k*k*k*k*k) + b4*(k*k*k*k) + b3*(k*k*k) + b2*(k*k) + b1*k + b0,
        R(-6.0)*b6*(k*k*k*k*k*k) - R(4.0)*b5*(k*k*k*k*k) - R(2.0)*b4*(k*k*k*k) + R(2.0)*b2*(k*k) + R(4.0)*b1*k + R(6.0)*b0,
        R(15.0)*b6*(k*k*k*k*k*k) + R(5.0)*b5*(k*k*k*k*k) - b4*(k*k*k*k) - R(3.0)*b3*(k*k*k) - b2*(k*k) + R(5.0)*b1*k + R(15.0)*b0,
        R(-20.0)*b6*(k*k*k*k*k*k) + R(4.0)*b4*(k*k*k*k) - R(4.0)*b2*(k*k) + R(20.0)*b0,
        R(15.0)*b6*(k*k*k*k*k*k) - R(5.0)*b5*(k*k*k*k*k) - b4*(k*k*k*k) + R(3.0)*b3*(k*k*k) - b2*(k*k) - R(5.0)*b1*k + R(15.0)*b0,
        R(-6.0)*b6*(k*k*k*k*k*k) + R(4.0)*b5*(k*k*k*k*k) - R(2.0)*b4*(k*k*k*k) + R(2.0)*b2*(k*k) - R(4.0)*b1*k + R(6.0)*b0,
        b6*(k*k*k*k*k*k) - b5*(k*k*k*k*k) + b4*(k*k*k*k) - b3*(k*k*k) + b2*(k*k) - b1*k + b0,
    };
    outTF.den = {
        a6*(k*k*k*k*k*k) + a5*(k*k*k*k*k) + a4*(k*k*k*k) + a3*(k*k*k) + a2*(k*k) + a1*k + a0,
        R(-6.0)*a6*(k*k*k*k*k*k) - R(4.0)*a5*(k*k*k*k*k) - R(2.0)*a4*(k*k*k*k) + R(2.0)*a2*(k*k) + R(4.0)*a1*k + R(6.0)*a0,
        R(15.0)*a6*(k*k*k*k*k*k) + R(5.0)*a5*(k*k*k*k*k) - a4*(k*k*k*k) - R(3.0)*a3*(k*k*k) - a2*(k*k) + R(5.0)*a1*k + R(15.0)*a0,
        R(-20.0)*a6*(k*k*k*k*k*k) + R(4.0)*a4*(k*k*k*k) - R(4.0)*a2*(k*k) + R(20.0)*a0,
        R(15.0)*a6*(k*k*k*k*k*k) - R(5.0)*a5*(k*k*k*k*k) - a4*(k*k*k*k) + R(3.0)*a3*(k*k*k) - a2*(k*k) - R(5.0)*a1*k + R(15.0)*a0,
        R(-6.0)*a6*(k*k*k*k*k*k) + R(4.0)*a5*(k*k*k*k*k) - R(2.0)*a4*(k*k*k*k) + R(2.0)*a2*(k*k) - R(4.0)*a1*k + R(6.0)*a0,
        a6*(k*k*k*k*k*k) - a5*(k*k*k*k*k) + a4*(k*k*k*k) - a3*(k*k*k) + a2*(k*k) - a1*k + a0,
    };

    return outTF;
}

int main(int argc, char *argv[])
{
    R param1 = 0;
    R param2 = 0;
    if (argc > 1)
        param1 = atof(argv[1]);
    if (argc > 2)
        param2 = atof(argv[2]);

    const R Fs = 44100;

    TF tfA = bassEQfilterAnalog(param1, param2);
    TF tfD = bilinearTransformOrd6(tfA, Fs);

    fprintf(stderr, "---- Analog\n");
    printTF(tfA, stderr);
    fprintf(stderr, "---- Normalized Analog\n");
    printTF(normalizeTF(tfA), stderr);
    fprintf(stderr, "---- Digital\n");
    printTF(tfD, stderr);
    fprintf(stderr, "---- Normalized Digital\n");
    printTF(normalizeTF(tfD), stderr);
    fprintf(stderr, "----\n");

    for (int key = 0; key < 140; ++key) {
        R f = R(440) * std::pow(R(2), (key - R(69)) / R(12));
        R w = R(2 * M_PI) * f;
        C s = j * w;
        C z = std::polar(R(1), w / Fs);
        C hA = calcTF(tfA, s);
        C hD = calcTF(tfD, R(1)/z);
        std::printf("%e %e %e %e %e\n", f, std::abs(hA), std::arg(hA), std::abs(hD), std::arg(hD));
    }
}