Skip to content

Instantly share code, notes, and snippets.

@jpcima

jpcima/bassFilter.cpp

Last active July 9, 2021 00:24
Show Gist options
  • Save jpcima/9067cae5864de3490860fcc7d73be54d to your computer and use it in GitHub Desktop.
Save jpcima/9067cae5864de3490860fcc7d73be54d to your computer and use it in GitHub Desktop.
Download ZIP
Raw
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));
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment