From ad8c3fbe741b0250bf8cfc8049459afa2365ad9a Mon Sep 17 00:00:00 2001 From: Christopher Herb Date: Fri, 7 Jul 2023 10:07:53 +0200 Subject: [PATCH] initial commit --- clip/aw_cliponly2.h | 100 +++++++ clip/aw_clipsoftly.h | 97 +++++++ clip/aw_tube2.h | 193 ++++++++++++ companding/mulaw.h | 47 +++ companding/ulaw.h | 93 ++++++ dynamics/aw_pop2.h | 312 ++++++++++++++++++++ filter/aw_eq.h | 677 +++++++++++++++++++++++++++++++++++++++++++ filter/chebyshev.h | 80 +++++ filter/ybandpass.h | 313 ++++++++++++++++++++ filter/yhighpass.h | 313 ++++++++++++++++++++ filter/ylowpass.h | 313 ++++++++++++++++++++ filter/ynotch.h | 313 ++++++++++++++++++++ filter/ysvf.h | 112 +++++++ synth/tx_envelope.h | 280 ++++++++++++++++++ synth/tx_operator.h | 39 +++ synth/tx_sineosc.h | 95 ++++++ synth/tx_voice.h | 166 +++++++++++ util/audio_math.h | 11 + 18 files changed, 3554 insertions(+) create mode 100644 clip/aw_cliponly2.h create mode 100644 clip/aw_clipsoftly.h create mode 100644 clip/aw_tube2.h create mode 100644 companding/mulaw.h create mode 100644 companding/ulaw.h create mode 100644 dynamics/aw_pop2.h create mode 100644 filter/aw_eq.h create mode 100644 filter/chebyshev.h create mode 100644 filter/ybandpass.h create mode 100644 filter/yhighpass.h create mode 100644 filter/ylowpass.h create mode 100644 filter/ynotch.h create mode 100644 filter/ysvf.h create mode 100644 synth/tx_envelope.h create mode 100644 synth/tx_operator.h create mode 100644 synth/tx_sineosc.h create mode 100644 synth/tx_voice.h create mode 100644 util/audio_math.h diff --git a/clip/aw_cliponly2.h b/clip/aw_cliponly2.h new file mode 100644 index 0000000..f7901b3 --- /dev/null +++ b/clip/aw_cliponly2.h @@ -0,0 +1,100 @@ +#pragma once +#include + +namespace trnr::lib::clip { +// Clipper based on ClipOnly2 by Chris Johnson +class aw_cliponly2 { +public: + aw_cliponly2() { + samplerate = 44100; + + lastSampleL = 0.0; + wasPosClipL = false; + wasNegClipL = false; + lastSampleR = 0.0; + wasPosClipR = false; + wasNegClipR = false; + for (int x = 0; x < 16; x++) {intermediateL[x] = 0.0; intermediateR[x] = 0.0;} + //this is reset: values being initialized only once. Startup values, whatever they are. + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_block(double** inputs, double** outputs, long sample_frames) { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + int spacing = floor(overallscale); //should give us working basic scaling, usually 2 or 4 + if (spacing < 1) spacing = 1; if (spacing > 16) spacing = 16; + + while (--sample_frames >= 0) + { + double inputSampleL = *in1; + double inputSampleR = *in2; + + //begin ClipOnly2 stereo as a little, compressed chunk that can be dropped into code + if (inputSampleL > 4.0) inputSampleL = 4.0; if (inputSampleL < -4.0) inputSampleL = -4.0; + if (wasPosClipL == true) { //current will be over + if (inputSampleL0.9549925859) {wasPosClipL=true;inputSampleL=0.7058208+(lastSampleL*0.2609148);} + if (wasNegClipL == true) { //current will be -over + if (inputSampleL > lastSampleL) lastSampleL=-0.7058208+(inputSampleL*0.2609148); + else lastSampleL=-0.2491717+(lastSampleL*0.7390851); + } wasNegClipL = false; + if (inputSampleL<-0.9549925859) {wasNegClipL=true;inputSampleL=-0.7058208+(lastSampleL*0.2609148);} + intermediateL[spacing] = inputSampleL; + inputSampleL = lastSampleL; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateL[x-1] = intermediateL[x]; + lastSampleL = intermediateL[0]; //run a little buffer to handle this + + if (inputSampleR > 4.0) inputSampleR = 4.0; if (inputSampleR < -4.0) inputSampleR = -4.0; + if (wasPosClipR == true) { //current will be over + if (inputSampleR0.9549925859) {wasPosClipR=true;inputSampleR=0.7058208+(lastSampleR*0.2609148);} + if (wasNegClipR == true) { //current will be -over + if (inputSampleR > lastSampleR) lastSampleR=-0.7058208+(inputSampleR*0.2609148); + else lastSampleR=-0.2491717+(lastSampleR*0.7390851); + } wasNegClipR = false; + if (inputSampleR<-0.9549925859) {wasNegClipR=true;inputSampleR=-0.7058208+(lastSampleR*0.2609148);} + intermediateR[spacing] = inputSampleR; + inputSampleR = lastSampleR; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateR[x-1] = intermediateR[x]; + lastSampleR = intermediateR[0]; //run a little buffer to handle this + //end ClipOnly2 stereo as a little, compressed chunk that can be dropped into code + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + + double lastSampleL; + double intermediateL[16]; + bool wasPosClipL; + bool wasNegClipL; + double lastSampleR; + double intermediateR[16]; + bool wasPosClipR; + bool wasNegClipR; //Stereo ClipOnly2 + //default stuff +}; +} \ No newline at end of file diff --git a/clip/aw_clipsoftly.h b/clip/aw_clipsoftly.h new file mode 100644 index 0000000..c80e7c2 --- /dev/null +++ b/clip/aw_clipsoftly.h @@ -0,0 +1,97 @@ +#pragma once +#include +#include + +namespace trnr::lib::clip { +// soft clipper based on ClipSoftly by Chris Johnson +class aw_clipsoftly { +public: + aw_clipsoftly() { + samplerate = 44100; + + lastSampleL = 0.0; + lastSampleR = 0.0; + for (int x = 0; x < 16; x++) {intermediateL[x] = 0.0; intermediateR[x] = 0.0;} + fpdL = 1.0; while (fpdL < 16386) fpdL = rand()*UINT32_MAX; + fpdR = 1.0; while (fpdR < 16386) fpdR = rand()*UINT32_MAX; + //this is reset: values being initialized only once. Startup values, whatever they are. + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_block(double** inputs, double** outputs, long sample_frames) { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + int spacing = floor(overallscale); //should give us working basic scaling, usually 2 or 4 + if (spacing < 1) spacing = 1; if (spacing > 16) spacing = 16; + + while (--sample_frames >= 0) + { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL)<1.18e-23) inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR)<1.18e-23) inputSampleR = fpdR * 1.18e-17; + + double softSpeed = fabs(inputSampleL); + if (softSpeed < 1.0) softSpeed = 1.0; else softSpeed = 1.0/softSpeed; + if (inputSampleL > 1.57079633) inputSampleL = 1.57079633; + if (inputSampleL < -1.57079633) inputSampleL = -1.57079633; + inputSampleL = sin(inputSampleL)*0.9549925859; //scale to what cliponly uses + inputSampleL = (inputSampleL*softSpeed)+(lastSampleL*(1.0-softSpeed)); + + softSpeed = fabs(inputSampleR); + if (softSpeed < 1.0) softSpeed = 1.0; else softSpeed = 1.0/softSpeed; + if (inputSampleR > 1.57079633) inputSampleR = 1.57079633; + if (inputSampleR < -1.57079633) inputSampleR = -1.57079633; + inputSampleR = sin(inputSampleR)*0.9549925859; //scale to what cliponly uses + inputSampleR = (inputSampleR*softSpeed)+(lastSampleR*(1.0-softSpeed)); + + intermediateL[spacing] = inputSampleL; + inputSampleL = lastSampleL; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateL[x-1] = intermediateL[x]; + lastSampleL = intermediateL[0]; //run a little buffer to handle this + + intermediateR[spacing] = inputSampleR; + inputSampleR = lastSampleR; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateR[x-1] = intermediateR[x]; + lastSampleR = intermediateR[0]; //run a little buffer to handle this + + //begin 64 bit stereo floating point dither + //int expon; frexp((double)inputSampleL, &expon); + fpdL ^= fpdL << 13; fpdL ^= fpdL >> 17; fpdL ^= fpdL << 5; + //inputSampleL += ((double(fpdL)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //frexp((double)inputSampleR, &expon); + fpdR ^= fpdR << 13; fpdR ^= fpdR >> 17; fpdR ^= fpdR << 5; + //inputSampleR += ((double(fpdR)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //end 64 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + + double lastSampleL; + double intermediateL[16]; + double lastSampleR; + double intermediateR[16]; + uint32_t fpdL; + uint32_t fpdR; + //default stuff +}; +} \ No newline at end of file diff --git a/clip/aw_tube2.h b/clip/aw_tube2.h new file mode 100644 index 0000000..51f4212 --- /dev/null +++ b/clip/aw_tube2.h @@ -0,0 +1,193 @@ +#pragma once +#include +#include + +namespace trnr::lib::clip { +// modeled tube preamp based on tube2 by Chris Johnson +class aw_tube2 { +public: + aw_tube2() { + samplerate = 44100; + + A = 0.5; + B = 0.5; + previousSampleA = 0.0; + previousSampleB = 0.0; + previousSampleC = 0.0; + previousSampleD = 0.0; + previousSampleE = 0.0; + previousSampleF = 0.0; + fpdL = 1.0; while (fpdL < 16386) fpdL = rand()*UINT32_MAX; + fpdR = 1.0; while (fpdR < 16386) fpdR = rand()*UINT32_MAX; + //this is reset: values being initialized only once. Startup values, whatever they are. + } + + void set_input(double value) { + A = clamp(value); + } + + void set_tube(double value) { + B = clamp(value); + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_block(double **inputs, double **outputs, long sampleframes) { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + double inputPad = A; + double iterations = 1.0-B; + int powerfactor = (9.0*iterations)+1; + double asymPad = (double)powerfactor; + double gainscaling = 1.0/(double)(powerfactor+1); + double outputscaling = 1.0 + (1.0/(double)(powerfactor)); + + while (--sampleframes >= 0) + { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL)<1.18e-23) inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR)<1.18e-23) inputSampleR = fpdR * 1.18e-17; + + if (inputPad < 1.0) { + inputSampleL *= inputPad; + inputSampleR *= inputPad; + } + + if (overallscale > 1.9) { + double stored = inputSampleL; + inputSampleL += previousSampleA; previousSampleA = stored; inputSampleL *= 0.5; + stored = inputSampleR; + inputSampleR += previousSampleB; previousSampleB = stored; inputSampleR *= 0.5; + } //for high sample rates on this plugin we are going to do a simple average + + if (inputSampleL > 1.0) inputSampleL = 1.0; + if (inputSampleL < -1.0) inputSampleL = -1.0; + if (inputSampleR > 1.0) inputSampleR = 1.0; + if (inputSampleR < -1.0) inputSampleR = -1.0; + + //flatten bottom, point top of sine waveshaper L + inputSampleL /= asymPad; + double sharpen = -inputSampleL; + if (sharpen > 0.0) sharpen = 1.0+sqrt(sharpen); + else sharpen = 1.0-sqrt(-sharpen); + inputSampleL -= inputSampleL*fabs(inputSampleL)*sharpen*0.25; + //this will take input from exactly -1.0 to 1.0 max + inputSampleL *= asymPad; + //flatten bottom, point top of sine waveshaper R + inputSampleR /= asymPad; + sharpen = -inputSampleR; + if (sharpen > 0.0) sharpen = 1.0+sqrt(sharpen); + else sharpen = 1.0-sqrt(-sharpen); + inputSampleR -= inputSampleR*fabs(inputSampleR)*sharpen*0.25; + //this will take input from exactly -1.0 to 1.0 max + inputSampleR *= asymPad; + //end first asym section: later boosting can mitigate the extreme + //softclipping of one side of the wave + //and we are asym clipping more when Tube is cranked, to compensate + + //original Tube algorithm: powerfactor widens the more linear region of the wave + double factor = inputSampleL; //Left channel + for (int x = 0; x < powerfactor; x++) factor *= inputSampleL; + if ((powerfactor % 2 == 1) && (inputSampleL != 0.0)) factor = (factor/inputSampleL)*fabs(inputSampleL); + factor *= gainscaling; + inputSampleL -= factor; + inputSampleL *= outputscaling; + factor = inputSampleR; //Right channel + for (int x = 0; x < powerfactor; x++) factor *= inputSampleR; + if ((powerfactor % 2 == 1) && (inputSampleR != 0.0)) factor = (factor/inputSampleR)*fabs(inputSampleR); + factor *= gainscaling; + inputSampleR -= factor; + inputSampleR *= outputscaling; + + if (overallscale > 1.9) { + double stored = inputSampleL; + inputSampleL += previousSampleC; previousSampleC = stored; inputSampleL *= 0.5; + stored = inputSampleR; + inputSampleR += previousSampleD; previousSampleD = stored; inputSampleR *= 0.5; + } //for high sample rates on this plugin we are going to do a simple average + //end original Tube. Now we have a boosted fat sound peaking at 0dB exactly + + //hysteresis and spiky fuzz L + double slew = previousSampleE - inputSampleL; + if (overallscale > 1.9) { + double stored = inputSampleL; + inputSampleL += previousSampleE; previousSampleE = stored; inputSampleL *= 0.5; + } else previousSampleE = inputSampleL; //for this, need previousSampleC always + if (slew > 0.0) slew = 1.0+(sqrt(slew)*0.5); + else slew = 1.0-(sqrt(-slew)*0.5); + inputSampleL -= inputSampleL*fabs(inputSampleL)*slew*gainscaling; + //reusing gainscaling that's part of another algorithm + if (inputSampleL > 0.52) inputSampleL = 0.52; + if (inputSampleL < -0.52) inputSampleL = -0.52; + inputSampleL *= 1.923076923076923; + //hysteresis and spiky fuzz R + slew = previousSampleF - inputSampleR; + if (overallscale > 1.9) { + double stored = inputSampleR; + inputSampleR += previousSampleF; previousSampleF = stored; inputSampleR *= 0.5; + } else previousSampleF = inputSampleR; //for this, need previousSampleC always + if (slew > 0.0) slew = 1.0+(sqrt(slew)*0.5); + else slew = 1.0-(sqrt(-slew)*0.5); + inputSampleR -= inputSampleR*fabs(inputSampleR)*slew*gainscaling; + //reusing gainscaling that's part of another algorithm + if (inputSampleR > 0.52) inputSampleR = 0.52; + if (inputSampleR < -0.52) inputSampleR = -0.52; + inputSampleR *= 1.923076923076923; + //end hysteresis and spiky fuzz section + + //begin 64 bit stereo floating point dither + //int expon; frexp((double)inputSampleL, &expon); + fpdL ^= fpdL << 13; fpdL ^= fpdL >> 17; fpdL ^= fpdL << 5; + //inputSampleL += ((double(fpdL)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //frexp((double)inputSampleR, &expon); + fpdR ^= fpdR << 13; fpdR ^= fpdR >> 17; fpdR ^= fpdR << 5; + //inputSampleR += ((double(fpdR)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //end 64 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + + double previousSampleA; + double previousSampleB; + double previousSampleC; + double previousSampleD; + double previousSampleE; + double previousSampleF; + + uint32_t fpdL; + uint32_t fpdR; + //default stuff + + float A; + float B; + + double clamp(double& value) { + if (value > 1) { + value = 1; + } else if (value < 0) { + value = 0; + } + return value; + } +}; +} \ No newline at end of file diff --git a/companding/mulaw.h b/companding/mulaw.h new file mode 100644 index 0000000..95306b1 --- /dev/null +++ b/companding/mulaw.h @@ -0,0 +1,47 @@ +#pragma once +#include + +namespace trnr::lib::companding { +// mulaw companding based on code by Emilie Gillet / Mutable Instruments +class mulaw { +public: + int8_t encode_samples(int16_t pcm_val) { + int16_t mask; + int16_t seg; + uint8_t uval; + pcm_val = pcm_val >> 2; + if (pcm_val < 0) { + pcm_val = -pcm_val; + mask = 0x7f; + } else { + mask = 0xff; + } + if (pcm_val > 8159) pcm_val = 8159; + pcm_val += (0x84 >> 2); + + if (pcm_val <= 0x3f) seg = 0; + else if (pcm_val <= 0x7f) seg = 1; + else if (pcm_val <= 0xff) seg = 2; + else if (pcm_val <= 0x1ff) seg = 3; + else if (pcm_val <= 0x3ff) seg = 4; + else if (pcm_val <= 0x7ff) seg = 5; + else if (pcm_val <= 0xfff) seg = 6; + else if (pcm_val <= 0x1fff) seg = 7; + else seg = 8; + if (seg >= 8) + return static_cast(0x7f ^ mask); + else { + uval = static_cast((seg << 4) | ((pcm_val >> (seg + 1)) & 0x0f)); + return (uval ^ mask); + } + } + + int16_t decode_samples(uint8_t u_val) { + int16_t t; + u_val = ~u_val; + t = ((u_val & 0xf) << 3) + 0x84; + t <<= ((unsigned)u_val & 0x70) >> 4; + return ((u_val & 0x80) ? (0x84 - t) : (t - 0x84)); + } +}; +} \ No newline at end of file diff --git a/companding/ulaw.h b/companding/ulaw.h new file mode 100644 index 0000000..83ff6c3 --- /dev/null +++ b/companding/ulaw.h @@ -0,0 +1,93 @@ +#pragma once +#include +#include +#include + +namespace trnr::lib::companding { +// ulaw compansion based on code by Chris Johnson +class ulaw { +public: + ulaw() { + fpd_l = 1.0; while (fpd_l < 16386) fpd_l = rand()*UINT32_MAX; + fpd_r = 1.0; while (fpd_r < 16386) fpd_r = rand()*UINT32_MAX; + } + + void encode_samples(double& input_sample_l, double& input_sample_r) { + + // ulaw encoding + static int noisesource_l = 0; + static int noisesource_r = 850010; + int residue; + double applyresidue; + + noisesource_l = noisesource_l % 1700021; noisesource_l++; + residue = noisesource_l * noisesource_l; + residue = residue % 170003; residue *= residue; + residue = residue % 17011; residue *= residue; + residue = residue % 1709; residue *= residue; + residue = residue % 173; residue *= residue; + residue = residue % 17; + applyresidue = residue; + applyresidue *= 0.00000001; + applyresidue *= 0.00000001; + input_sample_l += applyresidue; + if (input_sample_l<1.2e-38 && -input_sample_l<1.2e-38) { + input_sample_l -= applyresidue; + } + + noisesource_r = noisesource_r % 1700021; noisesource_r++; + residue = noisesource_r * noisesource_r; + residue = residue % 170003; residue *= residue; + residue = residue % 17011; residue *= residue; + residue = residue % 1709; residue *= residue; + residue = residue % 173; residue *= residue; + residue = residue % 17; + applyresidue = residue; + applyresidue *= 0.00000001; + applyresidue *= 0.00000001; + input_sample_r += applyresidue; + if (input_sample_r<1.2e-38 && -input_sample_r<1.2e-38) { + input_sample_r -= applyresidue; + } + + if (input_sample_l > 1.0) input_sample_l = 1.0; + if (input_sample_l < -1.0) input_sample_l = -1.0; + + if (input_sample_r > 1.0) input_sample_r = 1.0; + if (input_sample_r < -1.0) input_sample_r = -1.0; + + if (input_sample_l > 0) input_sample_l = log(1.0+(255*fabs(input_sample_l))) / log(256); + if (input_sample_l < 0) input_sample_l = -log(1.0+(255*fabs(input_sample_l))) / log(256); + + if (input_sample_r > 0) input_sample_r = log(1.0+(255*fabs(input_sample_r))) / log(256); + if (input_sample_r < 0) input_sample_r = -log(1.0+(255*fabs(input_sample_r))) / log(256); + } + + void decode_samples(double& input_sample_l, double& input_sample_r) { + + // ulaw decoding + if (fabs(input_sample_l)<1.18e-23) input_sample_l = fpd_l * 1.18e-17; + if (fabs(input_sample_r)<1.18e-23) input_sample_r = fpd_r * 1.18e-17; + + if (input_sample_l > 1.0) input_sample_l = 1.0; + if (input_sample_l < -1.0) input_sample_l = -1.0; + + if (input_sample_r > 1.0) input_sample_r = 1.0; + if (input_sample_r < -1.0) input_sample_r = -1.0; + + if (input_sample_l > 0) input_sample_l = (pow(256,fabs(input_sample_l))-1.0) / 255; + if (input_sample_l < 0) input_sample_l = -(pow(256,fabs(input_sample_l))-1.0) / 255; + + if (input_sample_r > 0) input_sample_r = (pow(256,fabs(input_sample_r))-1.0) / 255; + if (input_sample_r < 0) input_sample_r = -(pow(256,fabs(input_sample_r))-1.0) / 255; + + // 64 bit stereo floating point dither + fpd_l ^= fpd_l << 13; fpd_l ^= fpd_l >> 17; fpd_l ^= fpd_l << 5; + fpd_r ^= fpd_r << 13; fpd_r ^= fpd_r >> 17; fpd_r ^= fpd_r << 5; + } + +private: + uint32_t fpd_l; + uint32_t fpd_r; +}; +} \ No newline at end of file diff --git a/dynamics/aw_pop2.h b/dynamics/aw_pop2.h new file mode 100644 index 0000000..8c05f60 --- /dev/null +++ b/dynamics/aw_pop2.h @@ -0,0 +1,312 @@ +#pragma once +#include +#include + +namespace trnr::lib::dynamics { +// compressor based on pop2 by Chris Johnson +class aw_pop2 { +public: + aw_pop2() { + samplerate = 44100; + + A = 0.5; + B = 0.5; + C = 0.5; + D = 0.5; + E = 1.0; + fpdL = 1.0; while (fpdL < 16386) fpdL = rand()*UINT32_MAX; + fpdR = 1.0; while (fpdR < 16386) fpdR = rand()*UINT32_MAX; + + lastSampleL = 0.0; + wasPosClipL = false; + wasNegClipL = false; + lastSampleR = 0.0; + wasPosClipR = false; + wasNegClipR = false; + for (int x = 0; x < 16; x++) {intermediateL[x] = 0.0; intermediateR[x] = 0.0;} + + muVaryL = 0.0; + muAttackL = 0.0; + muNewSpeedL = 1000.0; + muSpeedAL = 1000.0; + muSpeedBL = 1000.0; + muCoefficientAL = 1.0; + muCoefficientBL = 1.0; + + muVaryR = 0.0; + muAttackR = 0.0; + muNewSpeedR = 1000.0; + muSpeedAR = 1000.0; + muSpeedBR = 1000.0; + muCoefficientAR = 1.0; + muCoefficientBR = 1.0; + + flip = false; + //this is reset: values being initialized only once. Startup values, whatever they are. + } + + void set_compression(double value) { + A = clamp(value); + } + + void set_attack(double value) { + B = clamp(value); + } + + void set_release(double value) { + C = clamp(value); + } + + void set_drive(double value) { + D = clamp(value); + } + + void set_drywet(double value) { + E = clamp(value); + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_block(double **inputs, double **outputs, long sampleframes) { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + int spacing = floor(overallscale); //should give us working basic scaling, usually 2 or 4 + if (spacing < 1) spacing = 1; if (spacing > 16) spacing = 16; + + double threshold = 1.0 - ((1.0-pow(1.0-A,2))*0.9); + double attack = ((pow(B,4)*100000.0)+10.0)*overallscale; + double release = ((pow(C,5)*2000000.0)+20.0)*overallscale; + double maxRelease = release * 4.0; + double muPreGain = 1.0/threshold; + double muMakeupGain = sqrt(1.0 / threshold)*D; + double wet = E; + //compressor section + + while (--sampleframes >= 0) + { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL)<1.18e-23) inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR)<1.18e-23) inputSampleR = fpdR * 1.18e-17; + double drySampleL = inputSampleL; + double drySampleR = inputSampleR; + + //begin compressor section + inputSampleL *= muPreGain; + inputSampleR *= muPreGain; + //adjust coefficients for L + if (flip) { + if (fabs(inputSampleL) > threshold) { + muVaryL = threshold / fabs(inputSampleL); + muAttackL = sqrt(fabs(muSpeedAL)); + muCoefficientAL = muCoefficientAL * (muAttackL-1.0); + if (muVaryL < threshold) muCoefficientAL = muCoefficientAL + threshold; + else muCoefficientAL = muCoefficientAL + muVaryL; + muCoefficientAL = muCoefficientAL / muAttackL; + muNewSpeedL = muSpeedAL * (muSpeedAL-1.0); + muNewSpeedL = muNewSpeedL + release; + muSpeedAL = muNewSpeedL / muSpeedAL; + if (muSpeedAL > maxRelease) muSpeedAL = maxRelease; + } else { + muCoefficientAL = muCoefficientAL * ((muSpeedAL * muSpeedAL)-1.0); + muCoefficientAL = muCoefficientAL + 1.0; + muCoefficientAL = muCoefficientAL / (muSpeedAL * muSpeedAL); + muNewSpeedL = muSpeedAL * (muSpeedAL-1.0); + muNewSpeedL = muNewSpeedL + attack; + muSpeedAL = muNewSpeedL / muSpeedAL;} + } else { + if (fabs(inputSampleL) > threshold) { + muVaryL = threshold / fabs(inputSampleL); + muAttackL = sqrt(fabs(muSpeedBL)); + muCoefficientBL = muCoefficientBL * (muAttackL-1); + if (muVaryL < threshold) muCoefficientBL = muCoefficientBL + threshold; + else muCoefficientBL = muCoefficientBL + muVaryL; + muCoefficientBL = muCoefficientBL / muAttackL; + muNewSpeedL = muSpeedBL * (muSpeedBL-1.0); + muNewSpeedL = muNewSpeedL + release; + muSpeedBL = muNewSpeedL / muSpeedBL; + if (muSpeedBL > maxRelease) muSpeedBL = maxRelease; + } else { + muCoefficientBL = muCoefficientBL * ((muSpeedBL * muSpeedBL)-1.0); + muCoefficientBL = muCoefficientBL + 1.0; + muCoefficientBL = muCoefficientBL / (muSpeedBL * muSpeedBL); + muNewSpeedL = muSpeedBL * (muSpeedBL-1.0); + muNewSpeedL = muNewSpeedL + attack; + muSpeedBL = muNewSpeedL / muSpeedBL; + } + } + //got coefficients, adjusted speeds for L + + //adjust coefficients for R + if (flip) { + if (fabs(inputSampleR) > threshold) { + muVaryR = threshold / fabs(inputSampleR); + muAttackR = sqrt(fabs(muSpeedAR)); + muCoefficientAR = muCoefficientAR * (muAttackR-1.0); + if (muVaryR < threshold) muCoefficientAR = muCoefficientAR + threshold; + else muCoefficientAR = muCoefficientAR + muVaryR; + muCoefficientAR = muCoefficientAR / muAttackR; + muNewSpeedR = muSpeedAR * (muSpeedAR-1.0); + muNewSpeedR = muNewSpeedR + release; + muSpeedAR = muNewSpeedR / muSpeedAR; + if (muSpeedAR > maxRelease) muSpeedAR = maxRelease; + } else { + muCoefficientAR = muCoefficientAR * ((muSpeedAR * muSpeedAR)-1.0); + muCoefficientAR = muCoefficientAR + 1.0; + muCoefficientAR = muCoefficientAR / (muSpeedAR * muSpeedAR); + muNewSpeedR = muSpeedAR * (muSpeedAR-1.0); + muNewSpeedR = muNewSpeedR + attack; + muSpeedAR = muNewSpeedR / muSpeedAR; + } + } else { + if (fabs(inputSampleR) > threshold) { + muVaryR = threshold / fabs(inputSampleR); + muAttackR = sqrt(fabs(muSpeedBR)); + muCoefficientBR = muCoefficientBR * (muAttackR-1); + if (muVaryR < threshold) muCoefficientBR = muCoefficientBR + threshold; + else muCoefficientBR = muCoefficientBR + muVaryR; + muCoefficientBR = muCoefficientBR / muAttackR; + muNewSpeedR = muSpeedBR * (muSpeedBR-1.0); + muNewSpeedR = muNewSpeedR + release; + muSpeedBR = muNewSpeedR / muSpeedBR; + if (muSpeedBR > maxRelease) muSpeedBR = maxRelease; + } else { + muCoefficientBR = muCoefficientBR * ((muSpeedBR * muSpeedBR)-1.0); + muCoefficientBR = muCoefficientBR + 1.0; + muCoefficientBR = muCoefficientBR / (muSpeedBR * muSpeedBR); + muNewSpeedR = muSpeedBR * (muSpeedBR-1.0); + muNewSpeedR = muNewSpeedR + attack; + muSpeedBR = muNewSpeedR / muSpeedBR; + } + } + //got coefficients, adjusted speeds for R + + if (flip) { + inputSampleL *= pow(muCoefficientAL,2); + inputSampleR *= pow(muCoefficientAR,2); + } else { + inputSampleL *= pow(muCoefficientBL,2); + inputSampleR *= pow(muCoefficientBR,2); + } + inputSampleL *= muMakeupGain; + inputSampleR *= muMakeupGain; + flip = !flip; + //end compressor section + + //begin ClipOnly2 stereo as a little, compressed chunk that can be dropped into code + if (inputSampleL > 4.0) inputSampleL = 4.0; if (inputSampleL < -4.0) inputSampleL = -4.0; + if (wasPosClipL == true) { //current will be over + if (inputSampleL0.9549925859) {wasPosClipL=true;inputSampleL=0.7058208+(lastSampleL*0.2609148);} + if (wasNegClipL == true) { //current will be -over + if (inputSampleL > lastSampleL) lastSampleL=-0.7058208+(inputSampleL*0.2609148); + else lastSampleL=-0.2491717+(lastSampleL*0.7390851); + } wasNegClipL = false; + if (inputSampleL<-0.9549925859) {wasNegClipL=true;inputSampleL=-0.7058208+(lastSampleL*0.2609148);} + intermediateL[spacing] = inputSampleL; + inputSampleL = lastSampleL; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateL[x-1] = intermediateL[x]; + lastSampleL = intermediateL[0]; //run a little buffer to handle this + + if (inputSampleR > 4.0) inputSampleR = 4.0; if (inputSampleR < -4.0) inputSampleR = -4.0; + if (wasPosClipR == true) { //current will be over + if (inputSampleR0.9549925859) {wasPosClipR=true;inputSampleR=0.7058208+(lastSampleR*0.2609148);} + if (wasNegClipR == true) { //current will be -over + if (inputSampleR > lastSampleR) lastSampleR=-0.7058208+(inputSampleR*0.2609148); + else lastSampleR=-0.2491717+(lastSampleR*0.7390851); + } wasNegClipR = false; + if (inputSampleR<-0.9549925859) {wasNegClipR=true;inputSampleR=-0.7058208+(lastSampleR*0.2609148);} + intermediateR[spacing] = inputSampleR; + inputSampleR = lastSampleR; //Latency is however many samples equals one 44.1k sample + for (int x = spacing; x > 0; x--) intermediateR[x-1] = intermediateR[x]; + lastSampleR = intermediateR[0]; //run a little buffer to handle this + //end ClipOnly2 stereo as a little, compressed chunk that can be dropped into code + + if (wet<1.0) { + inputSampleL = (drySampleL*(1.0-wet))+(inputSampleL*wet); + inputSampleR = (drySampleR*(1.0-wet))+(inputSampleR*wet); + } + + //begin 64 bit stereo floating point dither + //int expon; frexp((double)inputSampleL, &expon); + fpdL ^= fpdL << 13; fpdL ^= fpdL >> 17; fpdL ^= fpdL << 5; + //inputSampleL += ((double(fpdL)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //frexp((double)inputSampleR, &expon); + fpdR ^= fpdR << 13; fpdR ^= fpdR >> 17; fpdR ^= fpdR << 5; + //inputSampleR += ((double(fpdR)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //end 64 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + + uint32_t fpdL; + uint32_t fpdR; + //default stuff + + double muVaryL; + double muAttackL; + double muNewSpeedL; + double muSpeedAL; + double muSpeedBL; + double muCoefficientAL; + double muCoefficientBL; + + double muVaryR; + double muAttackR; + double muNewSpeedR; + double muSpeedAR; + double muSpeedBR; + double muCoefficientAR; + double muCoefficientBR; + + bool flip; + + double lastSampleL; + double intermediateL[16]; + bool wasPosClipL; + bool wasNegClipL; + double lastSampleR; + double intermediateR[16]; + bool wasPosClipR; + bool wasNegClipR; //Stereo ClipOnly2 + + float A; + float B; + float C; + float D; + float E; //parameters. Always 0-1, and we scale/alter them elsewhere. + + double clamp(double& value) { + if (value > 1) { + value = 1; + } else if (value < 0) { + value = 0; + } + return value; + } +}; +} \ No newline at end of file diff --git a/filter/aw_eq.h b/filter/aw_eq.h new file mode 100644 index 0000000..3212ac2 --- /dev/null +++ b/filter/aw_eq.h @@ -0,0 +1,677 @@ +#pragma once +#include +#include + +namespace trnr::lib::filter { +// 3 band equalizer with high/lowpass filters based on EQ by Chris Johnson. +class aw_eq { +public: + aw_eq() { + samplerate = 44100; + + A = 0.5; //Treble -12 to 12 + B = 0.5; //Mid -12 to 12 + C = 0.5; //Bass -12 to 12 + D = 1.0; //Lowpass 16.0K log 1 to 16 defaulting to 16K + E = 0.4; //TrebFrq 6.0 log 1 to 16 defaulting to 6K + F = 0.4; //BassFrq 100.0 log 30 to 1600 defaulting to 100 hz + G = 0.0; //Hipass 30.0 log 30 to 1600 defaulting to 30 + H = 0.5; //OutGain -18 to 18 + + lastSampleL = 0.0; + last2SampleL = 0.0; + lastSampleR = 0.0; + last2SampleR = 0.0; + + iirHighSampleLA = 0.0; + iirHighSampleLB = 0.0; + iirHighSampleLC = 0.0; + iirHighSampleLD = 0.0; + iirHighSampleLE = 0.0; + iirLowSampleLA = 0.0; + iirLowSampleLB = 0.0; + iirLowSampleLC = 0.0; + iirLowSampleLD = 0.0; + iirLowSampleLE = 0.0; + iirHighSampleL = 0.0; + iirLowSampleL = 0.0; + + iirHighSampleRA = 0.0; + iirHighSampleRB = 0.0; + iirHighSampleRC = 0.0; + iirHighSampleRD = 0.0; + iirHighSampleRE = 0.0; + iirLowSampleRA = 0.0; + iirLowSampleRB = 0.0; + iirLowSampleRC = 0.0; + iirLowSampleRD = 0.0; + iirLowSampleRE = 0.0; + iirHighSampleR = 0.0; + iirLowSampleR = 0.0; + + tripletLA = 0.0; + tripletLB = 0.0; + tripletLC = 0.0; + tripletFactorL = 0.0; + + tripletRA = 0.0; + tripletRB = 0.0; + tripletRC = 0.0; + tripletFactorR = 0.0; + + lowpassSampleLAA = 0.0; + lowpassSampleLAB = 0.0; + lowpassSampleLBA = 0.0; + lowpassSampleLBB = 0.0; + lowpassSampleLCA = 0.0; + lowpassSampleLCB = 0.0; + lowpassSampleLDA = 0.0; + lowpassSampleLDB = 0.0; + lowpassSampleLE = 0.0; + lowpassSampleLF = 0.0; + lowpassSampleLG = 0.0; + + lowpassSampleRAA = 0.0; + lowpassSampleRAB = 0.0; + lowpassSampleRBA = 0.0; + lowpassSampleRBB = 0.0; + lowpassSampleRCA = 0.0; + lowpassSampleRCB = 0.0; + lowpassSampleRDA = 0.0; + lowpassSampleRDB = 0.0; + lowpassSampleRE = 0.0; + lowpassSampleRF = 0.0; + lowpassSampleRG = 0.0; + + highpassSampleLAA = 0.0; + highpassSampleLAB = 0.0; + highpassSampleLBA = 0.0; + highpassSampleLBB = 0.0; + highpassSampleLCA = 0.0; + highpassSampleLCB = 0.0; + highpassSampleLDA = 0.0; + highpassSampleLDB = 0.0; + highpassSampleLE = 0.0; + highpassSampleLF = 0.0; + + highpassSampleRAA = 0.0; + highpassSampleRAB = 0.0; + highpassSampleRBA = 0.0; + highpassSampleRBB = 0.0; + highpassSampleRCA = 0.0; + highpassSampleRCB = 0.0; + highpassSampleRDA = 0.0; + highpassSampleRDB = 0.0; + highpassSampleRE = 0.0; + highpassSampleRF = 0.0; + + flip = false; + flipthree = 0; + + fpdL = 1.0; while (fpdL < 16386) fpdL = rand()*UINT32_MAX; + fpdR = 1.0; while (fpdR < 16386) fpdR = rand()*UINT32_MAX; + //this is reset: values being initialized only once. Startup values, whatever they are. + } + + void set_treble(double value) { + A = clamp(value); + } + + void set_mid(double value) { + B = clamp(value); + } + + void set_bass(double value) { + C = clamp(value); + } + + void set_lowpass(double value) { + D = clamp(value); + } + + void set_treble_frq(double value) { + E = clamp(value); + } + + void set_bass_frq(double value) { + F = clamp(value); + } + + void set_hipass(double value) { + G = clamp(value); + } + + void set_out_gain(double value) { + H = clamp(value); + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_block(double **inputs, double **outputs, long sampleframes) { + + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + double overallscale = 1.0; + overallscale /= 44100.0; + double compscale = overallscale; + overallscale = samplerate; + compscale = compscale * overallscale; + //compscale is the one that's 1 or something like 2.2 for 96K rates + + double inputSampleL; + double inputSampleR; + + double highSampleL = 0.0; + double midSampleL = 0.0; + double bassSampleL = 0.0; + + double highSampleR = 0.0; + double midSampleR = 0.0; + double bassSampleR = 0.0; + + double densityA = (A*12.0)-6.0; + double densityB = (B*12.0)-6.0; + double densityC = (C*12.0)-6.0; + bool engageEQ = true; + if ( (0.0 == densityA) && (0.0 == densityB) && (0.0 == densityC) ) engageEQ = false; + + densityA = pow(10.0,densityA/20.0)-1.0; + densityB = pow(10.0,densityB/20.0)-1.0; + densityC = pow(10.0,densityC/20.0)-1.0; + //convert to 0 to X multiplier with 1.0 being O db + //minus one gives nearly -1 to ? (should top out at 1) + //calibrate so that X db roughly equals X db with maximum topping out at 1 internally + + double tripletIntensity = -densityA; + + double iirAmountC = (((D*D*15.0)+1.0)*0.0188) + 0.7; + if (iirAmountC > 1.0) iirAmountC = 1.0; + bool engageLowpass = false; + if (((D*D*15.0)+1.0) < 15.99) engageLowpass = true; + + double iirAmountA = (((E*E*15.0)+1.0)*1000)/overallscale; + double iirAmountB = (((F*F*1570.0)+30.0)*10)/overallscale; + double iirAmountD = (((G*G*1570.0)+30.0)*1.0)/overallscale; + bool engageHighpass = false; + if (((G*G*1570.0)+30.0) > 30.01) engageHighpass = true; + //bypass the highpass and lowpass if set to extremes + double bridgerectifier; + double outA = fabs(densityA); + double outB = fabs(densityB); + double outC = fabs(densityC); + //end EQ + double outputgain = pow(10.0,((H*36.0)-18.0)/20.0); + + while (--sampleframes >= 0) + { + inputSampleL = *in1; + inputSampleR = *in2; + if (fabs(inputSampleL)<1.18e-23) inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR)<1.18e-23) inputSampleR = fpdR * 1.18e-17; + + last2SampleL = lastSampleL; + lastSampleL = inputSampleL; + + last2SampleR = lastSampleR; + lastSampleR = inputSampleR; + + flip = !flip; + flipthree++; + if (flipthree < 1 || flipthree > 3) flipthree = 1; + //counters + + //begin highpass + if (engageHighpass) + { + if (flip) + { + highpassSampleLAA = (highpassSampleLAA * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLAA; + highpassSampleLBA = (highpassSampleLBA * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLBA; + highpassSampleLCA = (highpassSampleLCA * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLCA; + highpassSampleLDA = (highpassSampleLDA * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLDA; + } + else + { + highpassSampleLAB = (highpassSampleLAB * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLAB; + highpassSampleLBB = (highpassSampleLBB * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLBB; + highpassSampleLCB = (highpassSampleLCB * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLCB; + highpassSampleLDB = (highpassSampleLDB * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLDB; + } + highpassSampleLE = (highpassSampleLE * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLE; + highpassSampleLF = (highpassSampleLF * (1.0 - iirAmountD)) + (inputSampleL * iirAmountD); + inputSampleL -= highpassSampleLF; + + if (flip) + { + highpassSampleRAA = (highpassSampleRAA * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRAA; + highpassSampleRBA = (highpassSampleRBA * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRBA; + highpassSampleRCA = (highpassSampleRCA * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRCA; + highpassSampleRDA = (highpassSampleRDA * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRDA; + } + else + { + highpassSampleRAB = (highpassSampleRAB * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRAB; + highpassSampleRBB = (highpassSampleRBB * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRBB; + highpassSampleRCB = (highpassSampleRCB * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRCB; + highpassSampleRDB = (highpassSampleRDB * (1.0 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRDB; + } + highpassSampleRE = (highpassSampleRE * (1 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRE; + highpassSampleRF = (highpassSampleRF * (1 - iirAmountD)) + (inputSampleR * iirAmountD); + inputSampleR -= highpassSampleRF; + + } + //end highpass + + //begin EQ + if (engageEQ) + { + switch (flipthree) + { + case 1: + tripletFactorL = last2SampleL - inputSampleL; + tripletLA += tripletFactorL; + tripletLC -= tripletFactorL; + tripletFactorL = tripletLA * tripletIntensity; + iirHighSampleLC = (iirHighSampleLC * (1.0 - iirAmountA)) + (inputSampleL * iirAmountA); + highSampleL = inputSampleL - iirHighSampleLC; + iirLowSampleLC = (iirLowSampleLC * (1.0 - iirAmountB)) + (inputSampleL * iirAmountB); + bassSampleL = iirLowSampleLC; + + tripletFactorR = last2SampleR - inputSampleR; + tripletRA += tripletFactorR; + tripletRC -= tripletFactorR; + tripletFactorR = tripletRA * tripletIntensity; + iirHighSampleRC = (iirHighSampleRC * (1.0 - iirAmountA)) + (inputSampleR * iirAmountA); + highSampleR = inputSampleR - iirHighSampleRC; + iirLowSampleRC = (iirLowSampleRC * (1.0 - iirAmountB)) + (inputSampleR * iirAmountB); + bassSampleR = iirLowSampleRC; + break; + case 2: + tripletFactorL = last2SampleL - inputSampleL; + tripletLB += tripletFactorL; + tripletLA -= tripletFactorL; + tripletFactorL = tripletLB * tripletIntensity; + iirHighSampleLD = (iirHighSampleLD * (1.0 - iirAmountA)) + (inputSampleL * iirAmountA); + highSampleL = inputSampleL - iirHighSampleLD; + iirLowSampleLD = (iirLowSampleLD * (1.0 - iirAmountB)) + (inputSampleL * iirAmountB); + bassSampleL = iirLowSampleLD; + + tripletFactorR = last2SampleR - inputSampleR; + tripletRB += tripletFactorR; + tripletRA -= tripletFactorR; + tripletFactorR = tripletRB * tripletIntensity; + iirHighSampleRD = (iirHighSampleRD * (1.0 - iirAmountA)) + (inputSampleR * iirAmountA); + highSampleR = inputSampleR - iirHighSampleRD; + iirLowSampleRD = (iirLowSampleRD * (1.0 - iirAmountB)) + (inputSampleR * iirAmountB); + bassSampleR = iirLowSampleRD; + break; + case 3: + tripletFactorL = last2SampleL - inputSampleL; + tripletLC += tripletFactorL; + tripletLB -= tripletFactorL; + tripletFactorL = tripletLC * tripletIntensity; + iirHighSampleLE = (iirHighSampleLE * (1.0 - iirAmountA)) + (inputSampleL * iirAmountA); + highSampleL = inputSampleL - iirHighSampleLE; + iirLowSampleLE = (iirLowSampleLE * (1.0 - iirAmountB)) + (inputSampleL * iirAmountB); + bassSampleL = iirLowSampleLE; + + tripletFactorR = last2SampleR - inputSampleR; + tripletRC += tripletFactorR; + tripletRB -= tripletFactorR; + tripletFactorR = tripletRC * tripletIntensity; + iirHighSampleRE = (iirHighSampleRE * (1.0 - iirAmountA)) + (inputSampleR * iirAmountA); + highSampleR = inputSampleR - iirHighSampleRE; + iirLowSampleRE = (iirLowSampleRE * (1.0 - iirAmountB)) + (inputSampleR * iirAmountB); + bassSampleR = iirLowSampleRE; + break; + } + tripletLA /= 2.0; + tripletLB /= 2.0; + tripletLC /= 2.0; + highSampleL = highSampleL + tripletFactorL; + + tripletRA /= 2.0; + tripletRB /= 2.0; + tripletRC /= 2.0; + highSampleR = highSampleR + tripletFactorR; + + if (flip) + { + iirHighSampleLA = (iirHighSampleLA * (1.0 - iirAmountA)) + (highSampleL * iirAmountA); + highSampleL -= iirHighSampleLA; + iirLowSampleLA = (iirLowSampleLA * (1.0 - iirAmountB)) + (bassSampleL * iirAmountB); + bassSampleL = iirLowSampleLA; + + iirHighSampleRA = (iirHighSampleRA * (1.0 - iirAmountA)) + (highSampleR * iirAmountA); + highSampleR -= iirHighSampleRA; + iirLowSampleRA = (iirLowSampleRA * (1.0 - iirAmountB)) + (bassSampleR * iirAmountB); + bassSampleR = iirLowSampleRA; + } + else + { + iirHighSampleLB = (iirHighSampleLB * (1.0 - iirAmountA)) + (highSampleL * iirAmountA); + highSampleL -= iirHighSampleLB; + iirLowSampleLB = (iirLowSampleLB * (1.0 - iirAmountB)) + (bassSampleL * iirAmountB); + bassSampleL = iirLowSampleLB; + + iirHighSampleRB = (iirHighSampleRB * (1.0 - iirAmountA)) + (highSampleR * iirAmountA); + highSampleR -= iirHighSampleRB; + iirLowSampleRB = (iirLowSampleRB * (1.0 - iirAmountB)) + (bassSampleR * iirAmountB); + bassSampleR = iirLowSampleRB; + } + + iirHighSampleL = (iirHighSampleL * (1.0 - iirAmountA)) + (highSampleL * iirAmountA); + highSampleL -= iirHighSampleL; + iirLowSampleL = (iirLowSampleL * (1.0 - iirAmountB)) + (bassSampleL * iirAmountB); + bassSampleL = iirLowSampleL; + + iirHighSampleR = (iirHighSampleR * (1.0 - iirAmountA)) + (highSampleR * iirAmountA); + highSampleR -= iirHighSampleR; + iirLowSampleR = (iirLowSampleR * (1.0 - iirAmountB)) + (bassSampleR * iirAmountB); + bassSampleR = iirLowSampleR; + + midSampleL = (inputSampleL-bassSampleL)-highSampleL; + midSampleR = (inputSampleR-bassSampleR)-highSampleR; + + //drive section + highSampleL *= (densityA+1.0); + bridgerectifier = fabs(highSampleL)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityA > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (highSampleL > 0) highSampleL = (highSampleL*(1-outA))+(bridgerectifier*outA); + else highSampleL = (highSampleL*(1-outA))-(bridgerectifier*outA); + //blend according to densityA control + + highSampleR *= (densityA+1.0); + bridgerectifier = fabs(highSampleR)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityA > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (highSampleR > 0) highSampleR = (highSampleR*(1-outA))+(bridgerectifier*outA); + else highSampleR = (highSampleR*(1-outA))-(bridgerectifier*outA); + //blend according to densityA control + + midSampleL *= (densityB+1.0); + bridgerectifier = fabs(midSampleL)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityB > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (midSampleL > 0) midSampleL = (midSampleL*(1-outB))+(bridgerectifier*outB); + else midSampleL = (midSampleL*(1-outB))-(bridgerectifier*outB); + //blend according to densityB control + + midSampleR *= (densityB+1.0); + bridgerectifier = fabs(midSampleR)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityB > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (midSampleR > 0) midSampleR = (midSampleR*(1-outB))+(bridgerectifier*outB); + else midSampleR = (midSampleR*(1-outB))-(bridgerectifier*outB); + //blend according to densityB control + + bassSampleL *= (densityC+1.0); + bridgerectifier = fabs(bassSampleL)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityC > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (bassSampleL > 0) bassSampleL = (bassSampleL*(1-outC))+(bridgerectifier*outC); + else bassSampleL = (bassSampleL*(1-outC))-(bridgerectifier*outC); + //blend according to densityC control + + bassSampleR *= (densityC+1.0); + bridgerectifier = fabs(bassSampleR)*1.57079633; + if (bridgerectifier > 1.57079633) bridgerectifier = 1.57079633; + //max value for sine function + if (densityC > 0) bridgerectifier = sin(bridgerectifier); + else bridgerectifier = 1-cos(bridgerectifier); + //produce either boosted or starved version + if (bassSampleR > 0) bassSampleR = (bassSampleR*(1-outC))+(bridgerectifier*outC); + else bassSampleR = (bassSampleR*(1-outC))-(bridgerectifier*outC); + //blend according to densityC control + + inputSampleL = midSampleL; + inputSampleL += highSampleL; + inputSampleL += bassSampleL; + + inputSampleR = midSampleR; + inputSampleR += highSampleR; + inputSampleR += bassSampleR; + } + //end EQ + + //EQ lowpass is after all processing like the compressor that might produce hash + if (engageLowpass) + { + if (flip) + { + lowpassSampleLAA = (lowpassSampleLAA * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLAA; + lowpassSampleLBA = (lowpassSampleLBA * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLBA; + lowpassSampleLCA = (lowpassSampleLCA * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLCA; + lowpassSampleLDA = (lowpassSampleLDA * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLDA; + lowpassSampleLE = (lowpassSampleLE * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLE; + + lowpassSampleRAA = (lowpassSampleRAA * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRAA; + lowpassSampleRBA = (lowpassSampleRBA * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRBA; + lowpassSampleRCA = (lowpassSampleRCA * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRCA; + lowpassSampleRDA = (lowpassSampleRDA * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRDA; + lowpassSampleRE = (lowpassSampleRE * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRE; + } + else + { + lowpassSampleLAB = (lowpassSampleLAB * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLAB; + lowpassSampleLBB = (lowpassSampleLBB * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLBB; + lowpassSampleLCB = (lowpassSampleLCB * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLCB; + lowpassSampleLDB = (lowpassSampleLDB * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLDB; + lowpassSampleLF = (lowpassSampleLF * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleL = lowpassSampleLF; + + lowpassSampleRAB = (lowpassSampleRAB * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRAB; + lowpassSampleRBB = (lowpassSampleRBB * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRBB; + lowpassSampleRCB = (lowpassSampleRCB * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRCB; + lowpassSampleRDB = (lowpassSampleRDB * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRDB; + lowpassSampleRF = (lowpassSampleRF * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + inputSampleR = lowpassSampleRF; + } + lowpassSampleLG = (lowpassSampleLG * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + lowpassSampleRG = (lowpassSampleRG * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + + inputSampleL = (lowpassSampleLG * (1.0 - iirAmountC)) + (inputSampleL * iirAmountC); + inputSampleR = (lowpassSampleRG * (1.0 - iirAmountC)) + (inputSampleR * iirAmountC); + } + + //built in output trim and dry/wet if desired + if (outputgain != 1.0) { + inputSampleL *= outputgain; + inputSampleR *= outputgain; + } + + //begin 64 bit stereo floating point dither + //int expon; frexp((double)inputSampleL, &expon); + fpdL ^= fpdL << 13; fpdL ^= fpdL >> 17; fpdL ^= fpdL << 5; + //inputSampleL += ((double(fpdL)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //frexp((double)inputSampleR, &expon); + fpdR ^= fpdR << 13; fpdR ^= fpdR >> 17; fpdR ^= fpdR << 5; + //inputSampleR += ((double(fpdR)-uint32_t(0x7fffffff)) * 1.1e-44l * pow(2,expon+62)); + //end 64 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + *in1++; + *in2++; + *out1++; + *out2++; + } + } + +private: + double samplerate; + + uint32_t fpdL; + uint32_t fpdR; + //default stuff + + double lastSampleL; + double last2SampleL; + double lastSampleR; + double last2SampleR; + + //begin EQ + double iirHighSampleLA; + double iirHighSampleLB; + double iirHighSampleLC; + double iirHighSampleLD; + double iirHighSampleLE; + double iirLowSampleLA; + double iirLowSampleLB; + double iirLowSampleLC; + double iirLowSampleLD; + double iirLowSampleLE; + double iirHighSampleL; + double iirLowSampleL; + + double iirHighSampleRA; + double iirHighSampleRB; + double iirHighSampleRC; + double iirHighSampleRD; + double iirHighSampleRE; + double iirLowSampleRA; + double iirLowSampleRB; + double iirLowSampleRC; + double iirLowSampleRD; + double iirLowSampleRE; + double iirHighSampleR; + double iirLowSampleR; + + double tripletLA; + double tripletLB; + double tripletLC; + double tripletFactorL; + + double tripletRA; + double tripletRB; + double tripletRC; + double tripletFactorR; + + double lowpassSampleLAA; + double lowpassSampleLAB; + double lowpassSampleLBA; + double lowpassSampleLBB; + double lowpassSampleLCA; + double lowpassSampleLCB; + double lowpassSampleLDA; + double lowpassSampleLDB; + double lowpassSampleLE; + double lowpassSampleLF; + double lowpassSampleLG; + + double lowpassSampleRAA; + double lowpassSampleRAB; + double lowpassSampleRBA; + double lowpassSampleRBB; + double lowpassSampleRCA; + double lowpassSampleRCB; + double lowpassSampleRDA; + double lowpassSampleRDB; + double lowpassSampleRE; + double lowpassSampleRF; + double lowpassSampleRG; + + double highpassSampleLAA; + double highpassSampleLAB; + double highpassSampleLBA; + double highpassSampleLBB; + double highpassSampleLCA; + double highpassSampleLCB; + double highpassSampleLDA; + double highpassSampleLDB; + double highpassSampleLE; + double highpassSampleLF; + + double highpassSampleRAA; + double highpassSampleRAB; + double highpassSampleRBA; + double highpassSampleRBB; + double highpassSampleRCA; + double highpassSampleRCB; + double highpassSampleRDA; + double highpassSampleRDB; + double highpassSampleRE; + double highpassSampleRF; + + bool flip; + int flipthree; + //end EQ + + + float A; + float B; + float C; + float D; + float E; + float F; + float G; + float H; + + double clamp(double& value) { + if (value > 1) { + value = 1; + } else if (value < 0) { + value = 0; + } + return value; + } +}; +} \ No newline at end of file diff --git a/filter/chebyshev.h b/filter/chebyshev.h new file mode 100644 index 0000000..d18ce5a --- /dev/null +++ b/filter/chebyshev.h @@ -0,0 +1,80 @@ +#pragma once +#define _USE_MATH_DEFINES +#include +#include + +namespace trnr::lib::filter { +class chebyshev { +public: + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void process_sample(double& input, double frequency) { + + if (frequency >= 20000.f) { + frequency = 20000.f; + } + + // First calculate the prewarped digital frequency : + auto K = tanf(M_PI * frequency / samplerate); + + // Now we calc some Coefficients : + auto sg = sinh(passband_ripple); + auto cg = cosh(passband_ripple); + cg *= cg; + + std::array coeff; + coeff[0] = 1 / (cg - 0.85355339059327376220042218105097); + coeff[1] = K * coeff[0] * sg * 1.847759065022573512256366378792; + coeff[2] = 1 / (cg - 0.14644660940672623779957781894758); + coeff[3] = K * coeff[2] * sg * 0.76536686473017954345691996806; + + K *= K; // (just to optimize it a little bit) + + // Calculate the first biquad: + a0 = 1 / (coeff[1] + K + coeff[0]); + a1 = 2 * (coeff[0] - K) * a0; + a2 = (coeff[1] - K - coeff[0]) * a0; + b0 = a0 * K; + b1 = 2 * b0; + b2 = b0; + + // Calculate the second biquad: + a3 = 1 / (coeff[3] + K + coeff[2]); + a4 = 2 * (coeff[2] - K) * a3; + a5 = (coeff[3] - K - coeff[2]) * a3; + b3 = a3 * K; + b4 = 2 * b3; + b5 = b3; + + // Then calculate the output as follows: + auto Stage1 = b0 * input + state0; + state0 = b1 * input + a1 * Stage1 + state1; + state1 = b2 * input + a2 * Stage1; + input = b3 * Stage1 + state2; + state2 = b4 * Stage1 + a4 * input + state3; + state3 = b5 * Stage1 + a5 * input; + } + +private: + double samplerate = 0; + double a0 = 0; + double a1 = 0; + double a2 = 0; + double a3 = 0; + double a4 = 0; + double a5 = 0; + double b0 = 0; + double b1 = 0; + double b2 = 0; + double b3 = 0; + double b4 = 0; + double b5 = 0; + double state0 = 0; + double state1 = 0; + double state2 = 0; + double state3 = 0; + double passband_ripple = 1; +}; +} \ No newline at end of file diff --git a/filter/ybandpass.h b/filter/ybandpass.h new file mode 100644 index 0000000..bcb4f51 --- /dev/null +++ b/filter/ybandpass.h @@ -0,0 +1,313 @@ +#pragma once +#define _USE_MATH_DEFINES +#include +#include +#include + +namespace trnr::lib::filter { +// Bandpass filter based on YBandpass by Chris Johnson +class ybandpass { +public: + ybandpass(double _samplerate) + : samplerate { _samplerate } + , A { 0.1f } + , B { 1.0f } + , C { 0.0f } + , D { 0.1f } + , E { 0.9f } + , F { 1.0f } + , fpdL { 0 } + , fpdR { 0 } + , biquad { 0 } + { + for (int x = 0; x < biq_total; x++) { + biquad[x] = 0.0; + } + powFactorA = 1.0; + powFactorB = 1.0; + inTrimA = 0.1; + inTrimB = 0.1; + outTrimA = 1.0; + outTrimB = 1.0; + for (int x = 0; x < fix_total; x++) { + fixA[x] = 0.0; + fixB[x] = 0.0; + } + + fpdL = 1.0; + while (fpdL < 16386) + fpdL = rand() * UINT32_MAX; + fpdR = 1.0; + while (fpdR < 16386) + fpdR = rand() * UINT32_MAX; + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void set_drive(float value) + { + A = value * 0.9 + 0.1; + } + void set_frequency(float value) + { + B = value; + } + void set_resonance(float value) + { + C = value; + } + void set_edge(float value) + { + D = value; + } + void set_output(float value) + { + E = value; + } + void set_mix(float value) + { + F = value; + } + void processblock(double** inputs, double** outputs, int blockSize) + { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + int inFramesToProcess = blockSize; + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + inTrimA = inTrimB; + inTrimB = A * 10.0; + + biquad[biq_freq] = pow(B, 3) * 20000.0; + if (biquad[biq_freq] < 15.0) + biquad[biq_freq] = 15.0; + biquad[biq_freq] /= samplerate; + biquad[biq_reso] = (pow(C, 2) * 15.0) + 0.5571; + biquad[biq_aA0] = biquad[biq_aB0]; + // biquad[biq_aA1] = biquad[biq_aB1]; + biquad[biq_aA2] = biquad[biq_aB2]; + biquad[biq_bA1] = biquad[biq_bB1]; + biquad[biq_bA2] = biquad[biq_bB2]; + // previous run through the buffer is still in the filter, so we move it + // to the A section and now it's the new starting point. + double K = tan(M_PI * biquad[biq_freq]); + double norm = 1.0 / (1.0 + K / biquad[biq_reso] + K * K); + biquad[biq_aB0] = K / biquad[biq_reso] * norm; + // biquad[biq_aB1] = 0.0; //bandpass can simplify the biquad kernel: leave out this multiply + biquad[biq_aB2] = -biquad[biq_aB0]; + biquad[biq_bB1] = 2.0 * (K * K - 1.0) * norm; + biquad[biq_bB2] = (1.0 - K / biquad[biq_reso] + K * K) * norm; + // for the coefficient-interpolated biquad filter + + powFactorA = powFactorB; + powFactorB = pow(D + 0.9, 4); + + // 1.0 == target neutral + + outTrimA = outTrimB; + outTrimB = E; + + double wet = F; + + fixA[fix_freq] = fixB[fix_freq] = 20000.0 / samplerate; + fixA[fix_reso] = fixB[fix_reso] = 0.7071; // butterworth Q + + K = tan(M_PI * fixA[fix_freq]); + norm = 1.0 / (1.0 + K / fixA[fix_reso] + K * K); + fixA[fix_a0] = fixB[fix_a0] = K * K * norm; + fixA[fix_a1] = fixB[fix_a1] = 2.0 * fixA[fix_a0]; + fixA[fix_a2] = fixB[fix_a2] = fixA[fix_a0]; + fixA[fix_b1] = fixB[fix_b1] = 2.0 * (K * K - 1.0) * norm; + fixA[fix_b2] = fixB[fix_b2] = (1.0 - K / fixA[fix_reso] + K * K) * norm; + // for the fixed-position biquad filter + + for (int s = 0; s < blockSize; s++) { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL) < 1.18e-23) + inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR) < 1.18e-23) + inputSampleR = fpdR * 1.18e-17; + double drySampleL = inputSampleL; + double drySampleR = inputSampleR; + + double temp = (double)s / inFramesToProcess; + biquad[biq_a0] = (biquad[biq_aA0] * temp) + (biquad[biq_aB0] * (1.0 - temp)); + // biquad[biq_a1] = (biquad[biq_aA1]*temp)+(biquad[biq_aB1]*(1.0-temp)); + biquad[biq_a2] = (biquad[biq_aA2] * temp) + (biquad[biq_aB2] * (1.0 - temp)); + biquad[biq_b1] = (biquad[biq_bA1] * temp) + (biquad[biq_bB1] * (1.0 - temp)); + biquad[biq_b2] = (biquad[biq_bA2] * temp) + (biquad[biq_bB2] * (1.0 - temp)); + // this is the interpolation code for the biquad + double powFactor = (powFactorA * temp) + (powFactorB * (1.0 - temp)); + double inTrim = (inTrimA * temp) + (inTrimB * (1.0 - temp)); + double outTrim = (outTrimA * temp) + (outTrimB * (1.0 - temp)); + + inputSampleL *= inTrim; + inputSampleR *= inTrim; + + temp = (inputSampleL * fixA[fix_a0]) + fixA[fix_sL1]; + fixA[fix_sL1] = (inputSampleL * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sL2]; + fixA[fix_sL2] = (inputSampleL * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixA[fix_a0]) + fixA[fix_sR1]; + fixA[fix_sR1] = (inputSampleR * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sR2]; + fixA[fix_sR2] = (inputSampleR * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, powFactor); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, powFactor); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, powFactor); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, powFactor); + + temp = (inputSampleL * biquad[biq_a0]) + biquad[biq_sL1]; + biquad[biq_sL1] = -(temp * biquad[biq_b1]) + biquad[biq_sL2]; + biquad[biq_sL2] = (inputSampleL * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleL = temp; // coefficient interpolating biquad filter + temp = (inputSampleR * biquad[biq_a0]) + biquad[biq_sR1]; + biquad[biq_sR1] = -(temp * biquad[biq_b1]) + biquad[biq_sR2]; + biquad[biq_sR2] = (inputSampleR * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleR = temp; // coefficient interpolating biquad filter + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, (1.0 / powFactor)); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, (1.0 / powFactor)); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, (1.0 / powFactor)); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, (1.0 / powFactor)); + + inputSampleL *= outTrim; + inputSampleR *= outTrim; + + temp = (inputSampleL * fixB[fix_a0]) + fixB[fix_sL1]; + fixB[fix_sL1] = (inputSampleL * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sL2]; + fixB[fix_sL2] = (inputSampleL * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixB[fix_a0]) + fixB[fix_sR1]; + fixB[fix_sR1] = (inputSampleR * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sR2]; + fixB[fix_sR2] = (inputSampleR * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + if (wet < 1.0) { + inputSampleL = (inputSampleL * wet) + (drySampleL * (1.0 - wet)); + inputSampleR = (inputSampleR * wet) + (drySampleR * (1.0 - wet)); + } + + // begin 32 bit stereo floating point dither + int expon; + frexpf((float)inputSampleL, &expon); + fpdL ^= fpdL << 13; + fpdL ^= fpdL >> 17; + fpdL ^= fpdL << 5; + inputSampleL += ((double(fpdL) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + frexpf((float)inputSampleR, &expon); + fpdR ^= fpdR << 13; + fpdR ^= fpdR >> 17; + fpdR ^= fpdR << 5; + inputSampleR += ((double(fpdR) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + // end 32 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + enum { + biq_freq, + biq_reso, + biq_a0, + biq_a1, + biq_a2, + biq_b1, + biq_b2, + biq_aA0, + biq_aA1, + biq_aA2, + biq_bA1, + biq_bA2, + biq_aB0, + biq_aB1, + biq_aB2, + biq_bB1, + biq_bB2, + biq_sL1, + biq_sL2, + biq_sR1, + biq_sR2, + biq_total + }; // coefficient interpolating biquad filter, stereo + std::array biquad; + + double powFactorA; + double powFactorB; + double inTrimA; + double inTrimB; + double outTrimA; + double outTrimB; + + enum { + fix_freq, + fix_reso, + fix_a0, + fix_a1, + fix_a2, + fix_b1, + fix_b2, + fix_sL1, + fix_sL2, + fix_sR1, + fix_sR2, + fix_total + }; // fixed frequency biquad filter for ultrasonics, stereo + std::array fixA; + std::array fixB; + + uint32_t fpdL; + uint32_t fpdR; + // default stuff + + float A; + float B; + float C; + float D; + float E; + float F; // parameters. Always 0-1, and we scale/alter them elsewhere. +}; +} \ No newline at end of file diff --git a/filter/yhighpass.h b/filter/yhighpass.h new file mode 100644 index 0000000..063394e --- /dev/null +++ b/filter/yhighpass.h @@ -0,0 +1,313 @@ +#pragma once +#define _USE_MATH_DEFINES +#include +#include +#include + +namespace trnr::lib::filter { +// Highpass filter based on YHighpass by Chris Johnson +class yhighpass { +public: + yhighpass(double _samplerate) + : samplerate { _samplerate } + , A { 0.1f } + , B { 1.0f } + , C { 0.0f } + , D { 0.1f } + , E { 0.9f } + , F { 1.0f } + , fpdL { 0 } + , fpdR { 0 } + , biquad { 0 } + { + for (int x = 0; x < biq_total; x++) { + biquad[x] = 0.0; + } + powFactorA = 1.0; + powFactorB = 1.0; + inTrimA = 0.1; + inTrimB = 0.1; + outTrimA = 1.0; + outTrimB = 1.0; + for (int x = 0; x < fix_total; x++) { + fixA[x] = 0.0; + fixB[x] = 0.0; + } + + fpdL = 1.0; + while (fpdL < 16386) + fpdL = rand() * UINT32_MAX; + fpdR = 1.0; + while (fpdR < 16386) + fpdR = rand() * UINT32_MAX; + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void set_drive(float value) + { + A = value * 0.9 + 0.1; + } + void set_frequency(float value) + { + B = value; + } + void set_resonance(float value) + { + C = value; + } + void set_edge(float value) + { + D = value; + } + void set_output(float value) + { + E = value; + } + void set_mix(float value) + { + F = value; + } + void processblock(double** inputs, double** outputs, int blockSize) + { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + int inFramesToProcess = blockSize; + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + inTrimA = inTrimB; + inTrimB = A * 10.0; + + biquad[biq_freq] = pow(B, 3) * 20000.0; + if (biquad[biq_freq] < 15.0) + biquad[biq_freq] = 15.0; + biquad[biq_freq] /= samplerate; + biquad[biq_reso] = (pow(C, 2) * 15.0) + 0.5571; + biquad[biq_aA0] = biquad[biq_aB0]; + biquad[biq_aA1] = biquad[biq_aB1]; + biquad[biq_aA2] = biquad[biq_aB2]; + biquad[biq_bA1] = biquad[biq_bB1]; + biquad[biq_bA2] = biquad[biq_bB2]; + // previous run through the buffer is still in the filter, so we move it + // to the A section and now it's the new starting point. + double K = tan(M_PI * biquad[biq_freq]); + double norm = 1.0 / (1.0 + K / biquad[biq_reso] + K * K); + biquad[biq_aB0] = norm; + biquad[biq_aB1] = -2.0 * biquad[biq_aB0]; + biquad[biq_aB2] = biquad[biq_aB0]; + biquad[biq_bB1] = 2.0 * (K * K - 1.0) * norm; + biquad[biq_bB2] = (1.0 - K / biquad[biq_reso] + K * K) * norm; + // for the coefficient-interpolated biquad filter + + powFactorA = powFactorB; + powFactorB = pow(D + 0.9, 4); + + // 1.0 == target neutral + + outTrimA = outTrimB; + outTrimB = E; + + double wet = F; + + fixA[fix_freq] = fixB[fix_freq] = 20000.0 / samplerate; + fixA[fix_reso] = fixB[fix_reso] = 0.7071; // butterworth Q + + K = tan(M_PI * fixA[fix_freq]); + norm = 1.0 / (1.0 + K / fixA[fix_reso] + K * K); + fixA[fix_a0] = fixB[fix_a0] = K * K * norm; + fixA[fix_a1] = fixB[fix_a1] = 2.0 * fixA[fix_a0]; + fixA[fix_a2] = fixB[fix_a2] = fixA[fix_a0]; + fixA[fix_b1] = fixB[fix_b1] = 2.0 * (K * K - 1.0) * norm; + fixA[fix_b2] = fixB[fix_b2] = (1.0 - K / fixA[fix_reso] + K * K) * norm; + // for the fixed-position biquad filter + + for (int s = 0; s < blockSize; s++) { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL) < 1.18e-23) + inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR) < 1.18e-23) + inputSampleR = fpdR * 1.18e-17; + double drySampleL = inputSampleL; + double drySampleR = inputSampleR; + + double temp = (double)s / inFramesToProcess; + biquad[biq_a0] = (biquad[biq_aA0] * temp) + (biquad[biq_aB0] * (1.0 - temp)); + biquad[biq_a1] = (biquad[biq_aA1] * temp) + (biquad[biq_aB1] * (1.0 - temp)); + biquad[biq_a2] = (biquad[biq_aA2] * temp) + (biquad[biq_aB2] * (1.0 - temp)); + biquad[biq_b1] = (biquad[biq_bA1] * temp) + (biquad[biq_bB1] * (1.0 - temp)); + biquad[biq_b2] = (biquad[biq_bA2] * temp) + (biquad[biq_bB2] * (1.0 - temp)); + // this is the interpolation code for the biquad + double powFactor = (powFactorA * temp) + (powFactorB * (1.0 - temp)); + double inTrim = (inTrimA * temp) + (inTrimB * (1.0 - temp)); + double outTrim = (outTrimA * temp) + (outTrimB * (1.0 - temp)); + + inputSampleL *= inTrim; + inputSampleR *= inTrim; + + temp = (inputSampleL * fixA[fix_a0]) + fixA[fix_sL1]; + fixA[fix_sL1] = (inputSampleL * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sL2]; + fixA[fix_sL2] = (inputSampleL * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixA[fix_a0]) + fixA[fix_sR1]; + fixA[fix_sR1] = (inputSampleR * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sR2]; + fixA[fix_sR2] = (inputSampleR * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, powFactor); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, powFactor); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, powFactor); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, powFactor); + + temp = (inputSampleL * biquad[biq_a0]) + biquad[biq_sL1]; + biquad[biq_sL1] = (inputSampleL * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sL2]; + biquad[biq_sL2] = (inputSampleL * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleL = temp; // coefficient interpolating biquad filter + temp = (inputSampleR * biquad[biq_a0]) + biquad[biq_sR1]; + biquad[biq_sR1] = (inputSampleR * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sR2]; + biquad[biq_sR2] = (inputSampleR * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleR = temp; // coefficient interpolating biquad filter + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, (1.0 / powFactor)); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, (1.0 / powFactor)); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, (1.0 / powFactor)); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, (1.0 / powFactor)); + + inputSampleL *= outTrim; + inputSampleR *= outTrim; + + temp = (inputSampleL * fixB[fix_a0]) + fixB[fix_sL1]; + fixB[fix_sL1] = (inputSampleL * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sL2]; + fixB[fix_sL2] = (inputSampleL * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixB[fix_a0]) + fixB[fix_sR1]; + fixB[fix_sR1] = (inputSampleR * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sR2]; + fixB[fix_sR2] = (inputSampleR * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + if (wet < 1.0) { + inputSampleL = (inputSampleL * wet) + (drySampleL * (1.0 - wet)); + inputSampleR = (inputSampleR * wet) + (drySampleR * (1.0 - wet)); + } + + // begin 32 bit stereo floating point dither + int expon; + frexpf((float)inputSampleL, &expon); + fpdL ^= fpdL << 13; + fpdL ^= fpdL >> 17; + fpdL ^= fpdL << 5; + inputSampleL += ((double(fpdL) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + frexpf((float)inputSampleR, &expon); + fpdR ^= fpdR << 13; + fpdR ^= fpdR >> 17; + fpdR ^= fpdR << 5; + inputSampleR += ((double(fpdR) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + // end 32 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + enum { + biq_freq, + biq_reso, + biq_a0, + biq_a1, + biq_a2, + biq_b1, + biq_b2, + biq_aA0, + biq_aA1, + biq_aA2, + biq_bA1, + biq_bA2, + biq_aB0, + biq_aB1, + biq_aB2, + biq_bB1, + biq_bB2, + biq_sL1, + biq_sL2, + biq_sR1, + biq_sR2, + biq_total + }; // coefficient interpolating biquad filter, stereo + std::array biquad; + + double powFactorA; + double powFactorB; + double inTrimA; + double inTrimB; + double outTrimA; + double outTrimB; + + enum { + fix_freq, + fix_reso, + fix_a0, + fix_a1, + fix_a2, + fix_b1, + fix_b2, + fix_sL1, + fix_sL2, + fix_sR1, + fix_sR2, + fix_total + }; // fixed frequency biquad filter for ultrasonics, stereo + std::array fixA; + std::array fixB; + + uint32_t fpdL; + uint32_t fpdR; + // default stuff + + float A; + float B; + float C; + float D; + float E; + float F; // parameters. Always 0-1, and we scale/alter them elsewhere. +}; +} \ No newline at end of file diff --git a/filter/ylowpass.h b/filter/ylowpass.h new file mode 100644 index 0000000..30cc750 --- /dev/null +++ b/filter/ylowpass.h @@ -0,0 +1,313 @@ +#pragma once +#define _USE_MATH_DEFINES +#include +#include +#include + +namespace trnr::lib::filter { +// Lowpass filter based on YLowpass by Chris Johnson +class ylowpass { +public: + ylowpass(double _samplerate) + : samplerate { _samplerate } + , A { 0.1f } + , B { 1.0f } + , C { 0.0f } + , D { 0.1f } + , E { 0.9f } + , F { 1.0f } + , fpdL { 0 } + , fpdR { 0 } + , biquad { 0 } + { + for (int x = 0; x < biq_total; x++) { + biquad[x] = 0.0; + } + powFactorA = 1.0; + powFactorB = 1.0; + inTrimA = 0.1; + inTrimB = 0.1; + outTrimA = 1.0; + outTrimB = 1.0; + for (int x = 0; x < fix_total; x++) { + fixA[x] = 0.0; + fixB[x] = 0.0; + } + + fpdL = 1.0; + while (fpdL < 16386) + fpdL = rand() * UINT32_MAX; + fpdR = 1.0; + while (fpdR < 16386) + fpdR = rand() * UINT32_MAX; + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void set_drive(float value) + { + A = value * 0.9 + 0.1; + } + void set_frequency(float value) + { + B = value; + } + void set_resonance(float value) + { + C = value; + } + void set_edge(float value) + { + D = value; + } + void set_output(float value) + { + E = value; + } + void set_mix(float value) + { + F = value; + } + void processblock(double** inputs, double** outputs, int blockSize) + { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + int inFramesToProcess = blockSize; + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + inTrimA = inTrimB; + inTrimB = A * 10.0; + + biquad[biq_freq] = pow(B, 3) * 20000.0; + if (biquad[biq_freq] < 15.0) + biquad[biq_freq] = 15.0; + biquad[biq_freq] /= samplerate; + biquad[biq_reso] = (pow(C, 2) * 15.0) + 0.5571; + biquad[biq_aA0] = biquad[biq_aB0]; + biquad[biq_aA1] = biquad[biq_aB1]; + biquad[biq_aA2] = biquad[biq_aB2]; + biquad[biq_bA1] = biquad[biq_bB1]; + biquad[biq_bA2] = biquad[biq_bB2]; + // previous run through the buffer is still in the filter, so we move it + // to the A section and now it's the new starting point. + double K = tan(M_PI * biquad[biq_freq]); + double norm = 1.0 / (1.0 + K / biquad[biq_reso] + K * K); + biquad[biq_aB0] = K * K * norm; + biquad[biq_aB1] = 2.0 * biquad[biq_aB0]; + biquad[biq_aB2] = biquad[biq_aB0]; + biquad[biq_bB1] = 2.0 * (K * K - 1.0) * norm; + biquad[biq_bB2] = (1.0 - K / biquad[biq_reso] + K * K) * norm; + // for the coefficient-interpolated biquad filter + + powFactorA = powFactorB; + powFactorB = pow(D + 0.9, 4); + + // 1.0 == target neutral + + outTrimA = outTrimB; + outTrimB = E; + + double wet = F; + + fixA[fix_freq] = fixB[fix_freq] = 20000.0 / samplerate; + fixA[fix_reso] = fixB[fix_reso] = 0.7071; // butterworth Q + + K = tan(M_PI * fixA[fix_freq]); + norm = 1.0 / (1.0 + K / fixA[fix_reso] + K * K); + fixA[fix_a0] = fixB[fix_a0] = K * K * norm; + fixA[fix_a1] = fixB[fix_a1] = 2.0 * fixA[fix_a0]; + fixA[fix_a2] = fixB[fix_a2] = fixA[fix_a0]; + fixA[fix_b1] = fixB[fix_b1] = 2.0 * (K * K - 1.0) * norm; + fixA[fix_b2] = fixB[fix_b2] = (1.0 - K / fixA[fix_reso] + K * K) * norm; + // for the fixed-position biquad filter + + for (int s = 0; s < blockSize; s++) { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL) < 1.18e-23) + inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR) < 1.18e-23) + inputSampleR = fpdR * 1.18e-17; + double drySampleL = inputSampleL; + double drySampleR = inputSampleR; + + double temp = (double)s / inFramesToProcess; + biquad[biq_a0] = (biquad[biq_aA0] * temp) + (biquad[biq_aB0] * (1.0 - temp)); + biquad[biq_a1] = (biquad[biq_aA1] * temp) + (biquad[biq_aB1] * (1.0 - temp)); + biquad[biq_a2] = (biquad[biq_aA2] * temp) + (biquad[biq_aB2] * (1.0 - temp)); + biquad[biq_b1] = (biquad[biq_bA1] * temp) + (biquad[biq_bB1] * (1.0 - temp)); + biquad[biq_b2] = (biquad[biq_bA2] * temp) + (biquad[biq_bB2] * (1.0 - temp)); + // this is the interpolation code for the biquad + double powFactor = (powFactorA * temp) + (powFactorB * (1.0 - temp)); + double inTrim = (inTrimA * temp) + (inTrimB * (1.0 - temp)); + double outTrim = (outTrimA * temp) + (outTrimB * (1.0 - temp)); + + inputSampleL *= inTrim; + inputSampleR *= inTrim; + + temp = (inputSampleL * fixA[fix_a0]) + fixA[fix_sL1]; + fixA[fix_sL1] = (inputSampleL * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sL2]; + fixA[fix_sL2] = (inputSampleL * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixA[fix_a0]) + fixA[fix_sR1]; + fixA[fix_sR1] = (inputSampleR * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sR2]; + fixA[fix_sR2] = (inputSampleR * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, powFactor); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, powFactor); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, powFactor); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, powFactor); + + temp = (inputSampleL * biquad[biq_a0]) + biquad[biq_sL1]; + biquad[biq_sL1] = (inputSampleL * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sL2]; + biquad[biq_sL2] = (inputSampleL * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleL = temp; // coefficient interpolating biquad filter + temp = (inputSampleR * biquad[biq_a0]) + biquad[biq_sR1]; + biquad[biq_sR1] = (inputSampleR * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sR2]; + biquad[biq_sR2] = (inputSampleR * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleR = temp; // coefficient interpolating biquad filter + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, (1.0 / powFactor)); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, (1.0 / powFactor)); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, (1.0 / powFactor)); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, (1.0 / powFactor)); + + inputSampleL *= outTrim; + inputSampleR *= outTrim; + + temp = (inputSampleL * fixB[fix_a0]) + fixB[fix_sL1]; + fixB[fix_sL1] = (inputSampleL * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sL2]; + fixB[fix_sL2] = (inputSampleL * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixB[fix_a0]) + fixB[fix_sR1]; + fixB[fix_sR1] = (inputSampleR * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sR2]; + fixB[fix_sR2] = (inputSampleR * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + if (wet < 1.0) { + inputSampleL = (inputSampleL * wet) + (drySampleL * (1.0 - wet)); + inputSampleR = (inputSampleR * wet) + (drySampleR * (1.0 - wet)); + } + + // begin 32 bit stereo floating point dither + int expon; + frexpf((float)inputSampleL, &expon); + fpdL ^= fpdL << 13; + fpdL ^= fpdL >> 17; + fpdL ^= fpdL << 5; + inputSampleL += ((double(fpdL) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + frexpf((float)inputSampleR, &expon); + fpdR ^= fpdR << 13; + fpdR ^= fpdR >> 17; + fpdR ^= fpdR << 5; + inputSampleR += ((double(fpdR) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + // end 32 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + enum { + biq_freq, + biq_reso, + biq_a0, + biq_a1, + biq_a2, + biq_b1, + biq_b2, + biq_aA0, + biq_aA1, + biq_aA2, + biq_bA1, + biq_bA2, + biq_aB0, + biq_aB1, + biq_aB2, + biq_bB1, + biq_bB2, + biq_sL1, + biq_sL2, + biq_sR1, + biq_sR2, + biq_total + }; // coefficient interpolating biquad filter, stereo + std::array biquad; + + double powFactorA; + double powFactorB; + double inTrimA; + double inTrimB; + double outTrimA; + double outTrimB; + + enum { + fix_freq, + fix_reso, + fix_a0, + fix_a1, + fix_a2, + fix_b1, + fix_b2, + fix_sL1, + fix_sL2, + fix_sR1, + fix_sR2, + fix_total + }; // fixed frequency biquad filter for ultrasonics, stereo + std::array fixA; + std::array fixB; + + uint32_t fpdL; + uint32_t fpdR; + // default stuff + + float A; + float B; + float C; + float D; + float E; + float F; // parameters. Always 0-1, and we scale/alter them elsewhere. +}; +} \ No newline at end of file diff --git a/filter/ynotch.h b/filter/ynotch.h new file mode 100644 index 0000000..64c4a20 --- /dev/null +++ b/filter/ynotch.h @@ -0,0 +1,313 @@ +#pragma once +#define _USE_MATH_DEFINES +#include +#include +#include + +namespace trnr::lib::filter { +// Notch filter based on YNotch by Chris Johnson +class ynotch { +public: + ynotch(double _samplerate) + : samplerate { _samplerate } + , A { 0.1f } + , B { 1.0f } + , C { 0.0f } + , D { 0.1f } + , E { 0.9f } + , F { 1.0f } + , fpdL { 0 } + , fpdR { 0 } + , biquad { 0 } + { + for (int x = 0; x < biq_total; x++) { + biquad[x] = 0.0; + } + powFactorA = 1.0; + powFactorB = 1.0; + inTrimA = 0.1; + inTrimB = 0.1; + outTrimA = 1.0; + outTrimB = 1.0; + for (int x = 0; x < fix_total; x++) { + fixA[x] = 0.0; + fixB[x] = 0.0; + } + + fpdL = 1.0; + while (fpdL < 16386) + fpdL = rand() * UINT32_MAX; + fpdR = 1.0; + while (fpdR < 16386) + fpdR = rand() * UINT32_MAX; + } + + void set_samplerate(double _samplerate) { + samplerate = _samplerate; + } + + void set_drive(float value) + { + A = value * 0.9 + 0.1; + } + void set_frequency(float value) + { + B = value; + } + void set_resonance(float value) + { + C = value; + } + void set_edge(float value) + { + D = value; + } + void set_output(float value) + { + E = value; + } + void set_mix(float value) + { + F = value; + } + void processblock(double** inputs, double** outputs, int blockSize) + { + double* in1 = inputs[0]; + double* in2 = inputs[1]; + double* out1 = outputs[0]; + double* out2 = outputs[1]; + + int inFramesToProcess = blockSize; + double overallscale = 1.0; + overallscale /= 44100.0; + overallscale *= samplerate; + + inTrimA = inTrimB; + inTrimB = A * 10.0; + + biquad[biq_freq] = pow(B, 3) * 20000.0; + if (biquad[biq_freq] < 15.0) + biquad[biq_freq] = 15.0; + biquad[biq_freq] /= samplerate; + biquad[biq_reso] = (pow(C, 2) * 15.0) + 0.0001; + biquad[biq_aA0] = biquad[biq_aB0]; + biquad[biq_aA1] = biquad[biq_aB1]; + biquad[biq_aA2] = biquad[biq_aB2]; + biquad[biq_bA1] = biquad[biq_bB1]; + biquad[biq_bA2] = biquad[biq_bB2]; + // previous run through the buffer is still in the filter, so we move it + // to the A section and now it's the new starting point. + double K = tan(M_PI * biquad[biq_freq]); + double norm = 1.0 / (1.0 + K / biquad[biq_reso] + K * K); + biquad[biq_aB0] = (1.0 + K * K) * norm; + biquad[biq_aB1] = 2.0 * (K * K - 1) * norm; + biquad[biq_aB2] = biquad[biq_aB0]; + biquad[biq_bB1] = biquad[biq_aB1]; + biquad[biq_bB2] = (1.0 - K / biquad[biq_reso] + K * K) * norm; + // for the coefficient-interpolated biquad filter + + powFactorA = powFactorB; + powFactorB = pow(D + 0.9, 4); + + // 1.0 == target neutral + + outTrimA = outTrimB; + outTrimB = E; + + double wet = F; + + fixA[fix_freq] = fixB[fix_freq] = 20000.0 / samplerate; + fixA[fix_reso] = fixB[fix_reso] = 0.7071; // butterworth Q + + K = tan(M_PI * fixA[fix_freq]); + norm = 1.0 / (1.0 + K / fixA[fix_reso] + K * K); + fixA[fix_a0] = fixB[fix_a0] = K * K * norm; + fixA[fix_a1] = fixB[fix_a1] = 2.0 * fixA[fix_a0]; + fixA[fix_a2] = fixB[fix_a2] = fixA[fix_a0]; + fixA[fix_b1] = fixB[fix_b1] = 2.0 * (K * K - 1.0) * norm; + fixA[fix_b2] = fixB[fix_b2] = (1.0 - K / fixA[fix_reso] + K * K) * norm; + // for the fixed-position biquad filter + + for (int s = 0; s < blockSize; s++) { + double inputSampleL = *in1; + double inputSampleR = *in2; + if (fabs(inputSampleL) < 1.18e-23) + inputSampleL = fpdL * 1.18e-17; + if (fabs(inputSampleR) < 1.18e-23) + inputSampleR = fpdR * 1.18e-17; + double drySampleL = inputSampleL; + double drySampleR = inputSampleR; + + double temp = (double)s / inFramesToProcess; + biquad[biq_a0] = (biquad[biq_aA0] * temp) + (biquad[biq_aB0] * (1.0 - temp)); + biquad[biq_a1] = (biquad[biq_aA1] * temp) + (biquad[biq_aB1] * (1.0 - temp)); + biquad[biq_a2] = (biquad[biq_aA2] * temp) + (biquad[biq_aB2] * (1.0 - temp)); + biquad[biq_b1] = (biquad[biq_bA1] * temp) + (biquad[biq_bB1] * (1.0 - temp)); + biquad[biq_b2] = (biquad[biq_bA2] * temp) + (biquad[biq_bB2] * (1.0 - temp)); + // this is the interpolation code for the biquad + double powFactor = (powFactorA * temp) + (powFactorB * (1.0 - temp)); + double inTrim = (inTrimA * temp) + (inTrimB * (1.0 - temp)); + double outTrim = (outTrimA * temp) + (outTrimB * (1.0 - temp)); + + inputSampleL *= inTrim; + inputSampleR *= inTrim; + + temp = (inputSampleL * fixA[fix_a0]) + fixA[fix_sL1]; + fixA[fix_sL1] = (inputSampleL * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sL2]; + fixA[fix_sL2] = (inputSampleL * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixA[fix_a0]) + fixA[fix_sR1]; + fixA[fix_sR1] = (inputSampleR * fixA[fix_a1]) - (temp * fixA[fix_b1]) + fixA[fix_sR2]; + fixA[fix_sR2] = (inputSampleR * fixA[fix_a2]) - (temp * fixA[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, powFactor); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, powFactor); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, powFactor); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, powFactor); + + temp = (inputSampleL * biquad[biq_a0]) + biquad[biq_sL1]; + biquad[biq_sL1] = (inputSampleL * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sL2]; + biquad[biq_sL2] = (inputSampleL * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleL = temp; // coefficient interpolating biquad filter + temp = (inputSampleR * biquad[biq_a0]) + biquad[biq_sR1]; + biquad[biq_sR1] = (inputSampleR * biquad[biq_a1]) - (temp * biquad[biq_b1]) + biquad[biq_sR2]; + biquad[biq_sR2] = (inputSampleR * biquad[biq_a2]) - (temp * biquad[biq_b2]); + inputSampleR = temp; // coefficient interpolating biquad filter + + // encode/decode courtesy of torridgristle under the MIT license + if (inputSampleL > 1.0) + inputSampleL = 1.0; + else if (inputSampleL > 0.0) + inputSampleL = 1.0 - pow(1.0 - inputSampleL, (1.0 / powFactor)); + if (inputSampleL < -1.0) + inputSampleL = -1.0; + else if (inputSampleL < 0.0) + inputSampleL = -1.0 + pow(1.0 + inputSampleL, (1.0 / powFactor)); + if (inputSampleR > 1.0) + inputSampleR = 1.0; + else if (inputSampleR > 0.0) + inputSampleR = 1.0 - pow(1.0 - inputSampleR, (1.0 / powFactor)); + if (inputSampleR < -1.0) + inputSampleR = -1.0; + else if (inputSampleR < 0.0) + inputSampleR = -1.0 + pow(1.0 + inputSampleR, (1.0 / powFactor)); + + inputSampleL *= outTrim; + inputSampleR *= outTrim; + + temp = (inputSampleL * fixB[fix_a0]) + fixB[fix_sL1]; + fixB[fix_sL1] = (inputSampleL * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sL2]; + fixB[fix_sL2] = (inputSampleL * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleL = temp; // fixed biquad filtering ultrasonics + temp = (inputSampleR * fixB[fix_a0]) + fixB[fix_sR1]; + fixB[fix_sR1] = (inputSampleR * fixB[fix_a1]) - (temp * fixB[fix_b1]) + fixB[fix_sR2]; + fixB[fix_sR2] = (inputSampleR * fixB[fix_a2]) - (temp * fixB[fix_b2]); + inputSampleR = temp; // fixed biquad filtering ultrasonics + + if (wet < 1.0) { + inputSampleL = (inputSampleL * wet) + (drySampleL * (1.0 - wet)); + inputSampleR = (inputSampleR * wet) + (drySampleR * (1.0 - wet)); + } + + // begin 32 bit stereo floating point dither + int expon; + frexpf((float)inputSampleL, &expon); + fpdL ^= fpdL << 13; + fpdL ^= fpdL >> 17; + fpdL ^= fpdL << 5; + inputSampleL += ((double(fpdL) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + frexpf((float)inputSampleR, &expon); + fpdR ^= fpdR << 13; + fpdR ^= fpdR >> 17; + fpdR ^= fpdR << 5; + inputSampleR += ((double(fpdR) - uint32_t(0x7fffffff)) * 5.5e-36l * pow(2, expon + 62)); + // end 32 bit stereo floating point dither + + *out1 = inputSampleL; + *out2 = inputSampleR; + + in1++; + in2++; + out1++; + out2++; + } + } + +private: + double samplerate; + enum { + biq_freq, + biq_reso, + biq_a0, + biq_a1, + biq_a2, + biq_b1, + biq_b2, + biq_aA0, + biq_aA1, + biq_aA2, + biq_bA1, + biq_bA2, + biq_aB0, + biq_aB1, + biq_aB2, + biq_bB1, + biq_bB2, + biq_sL1, + biq_sL2, + biq_sR1, + biq_sR2, + biq_total + }; // coefficient interpolating biquad filter, stereo + std::array biquad; + + double powFactorA; + double powFactorB; + double inTrimA; + double inTrimB; + double outTrimA; + double outTrimB; + + enum { + fix_freq, + fix_reso, + fix_a0, + fix_a1, + fix_a2, + fix_b1, + fix_b2, + fix_sL1, + fix_sL2, + fix_sR1, + fix_sR2, + fix_total + }; // fixed frequency biquad filter for ultrasonics, stereo + std::array fixA; + std::array fixB; + + uint32_t fpdL; + uint32_t fpdR; + // default stuff + + float A; + float B; + float C; + float D; + float E; + float F; // parameters. Always 0-1, and we scale/alter them elsewhere. +}; +} \ No newline at end of file diff --git a/filter/ysvf.h b/filter/ysvf.h new file mode 100644 index 0000000..78110fe --- /dev/null +++ b/filter/ysvf.h @@ -0,0 +1,112 @@ +#pragma once +#include "ylowpass.h" +#include "yhighpass.h" +#include "ybandpass.h" +#include "ynotch.h" + +namespace trnr::lib::filter { + +enum filter_types { + lowpass = 0, + highpass, + bandpass, + notch +}; + +class ysvf { +public: + ysvf(double _samplerate) + : lowpass { _samplerate } + , highpass { _samplerate } + , bandpass { _samplerate } + , notch { _samplerate } + {} + + void set_samplerate(double _samplerate) { + lowpass.set_samplerate(_samplerate); + highpass.set_samplerate(_samplerate); + bandpass.set_samplerate(_samplerate); + notch.set_samplerate(_samplerate); + } + + void set_filter_type(filter_types type) { + filter_type = type; + } + + void set_drive(float value) { + lowpass.set_drive(value); + highpass.set_drive(value); + bandpass.set_drive(value); + notch.set_drive(value); + } + + void set_frequency(float value) { + lowpass.set_frequency(value); + highpass.set_frequency(value); + bandpass.set_frequency(value); + notch.set_frequency(value); + } + + void set_resonance(float value) { + lowpass.set_resonance(value); + highpass.set_resonance(value); + bandpass.set_resonance(value); + notch.set_resonance(value); + } + + void set_edge(float value) { + lowpass.set_edge(value); + highpass.set_edge(value); + bandpass.set_edge(value); + notch.set_edge(value); + } + + void set_output(float value) { + lowpass.set_output(value); + highpass.set_output(value); + bandpass.set_output(value); + notch.set_output(value); + } + + void set_mix(float value) { + lowpass.set_mix(value); + highpass.set_mix(value); + bandpass.set_mix(value); + notch.set_mix(value); + } + + void process_block(double** inputs, double** outputs, int block_size) { + + switch (filter_type) { + case filter_types::lowpass: + lowpass.processblock(inputs, outputs, block_size); + break; + case filter_types::highpass: + highpass.processblock(inputs, outputs, block_size); + break; + case filter_types::bandpass: + bandpass.processblock(inputs, outputs, block_size); + break; + case filter_types::notch: + notch.processblock(inputs, outputs, block_size); + break; + } + } + +private: + filter_types filter_type; + ylowpass lowpass; + yhighpass highpass; + ybandpass bandpass; + ynotch notch; + + double clamp(double& value, double min, double max) { + if (value < min) { + value = min; + } else if (value > max) { + value = max; + } + return value; + } +}; +} \ No newline at end of file diff --git a/synth/tx_envelope.h b/synth/tx_envelope.h new file mode 100644 index 0000000..fb160aa --- /dev/null +++ b/synth/tx_envelope.h @@ -0,0 +1,280 @@ +#pragma once +#include + +namespace trnr::lib::synth { + +enum env_state { + idle = 0, + attack1, + attack2, + hold, + decay1, + decay2, + sustain, + release1, + release2 +}; + +class tx_envelope { +public: + float attack1_rate; + float attack1_level; + float attack2_rate; + float hold_rate; + float decay1_rate; + float decay1_level; + float decay2_rate; + float sustain_level; + float release1_rate; + float release1_level; + float release2_rate; + + tx_envelope(double _samplerate) + : samplerate { _samplerate } + , attack1_rate { 0 } + , attack1_level { 0 } + , attack2_rate { 0 } + , hold_rate { 0 } + , decay1_rate { 0 } + , decay1_level { 0 } + , decay2_rate { 0 } + , sustain_level { 0 } + , release1_rate { 0 } + , release1_level { 0 } + , release2_rate { 0 } + , level { 0.f } + , phase { 0 } + , state { idle } + , start_level { 0.f } + , h1 { 0. } + , h2 { 0. } + , h3 { 0. } + { + } + + float process_sample(bool gate, bool trigger) { + + int attack_mid_x1 = ms_to_samples(attack1_rate); + int attack_mid_x2 = ms_to_samples(attack2_rate); + int hold_samp = ms_to_samples(hold_rate); + int decay_mid_x1 = ms_to_samples(decay1_rate); + int decay_mid_x2 = ms_to_samples(decay2_rate); + int release_mid_x1 = ms_to_samples(release1_rate); + int release_mid_x2 = ms_to_samples(release2_rate); + + // if note on is triggered, transition to attack phase + if (trigger) { + start_level = level; + phase = 0; + state = attack1; + } + // attack 1st half + if (state == attack1) { + // while in attack phase + if (phase < attack_mid_x1) { + level = lerp(0, start_level, attack_mid_x1, attack1_level, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > attack_mid_x1) { + phase = attack_mid_x1; + } + // if attack phase is done, transition to decay phase + if (phase == attack_mid_x1) { + state = attack2; + phase = 0; + } + } + // attack 2nd half + if (state == attack2) { + // while in attack phase + if (phase < attack_mid_x2) { + level = lerp(0, attack1_level, attack_mid_x2, 1, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > attack_mid_x2) { + phase = attack_mid_x2; + } + // if attack phase is done, transition to decay phase + if (phase == attack_mid_x2) { + state = hold; + phase = 0; + } + } + // hold + if (state == hold) { + if (phase < hold_samp) { + level = 1.0; + phase += 1; + } + if (phase > hold_samp) { + phase = hold_samp; + } + if (phase == hold_samp) { + state = decay1; + phase = 0; + } + } + // decay 1st half + if (state == decay1) { + // while in decay phase + if (phase < decay_mid_x1) { + level = lerp(0, 1, decay_mid_x1, decay1_level, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > decay_mid_x1) { + phase = decay_mid_x1; + } + // if decay phase is done, transition to sustain phase + if (phase == decay_mid_x1) { + state = decay2; + phase = 0; + } + } + // decay 2nd half + if (state == decay2) { + // while in decay phase + if (phase < decay_mid_x2) { + level = lerp(0, decay1_level, decay_mid_x2, sustain_level, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > decay_mid_x2) { + phase = decay_mid_x2; + } + // if decay phase is done, transition to sustain phase + if (phase == decay_mid_x2) { + state = sustain; + phase = 0; + level = sustain_level; + } + } + // while sustain phase: if note off is triggered, transition to release phase + if (state == sustain && !gate) { + state = release1; + level = sustain_level; + } + // release 1st half + if (state == release1) { + // while in release phase + if (phase < release_mid_x1) { + level = lerp(0, sustain_level, release_mid_x1, release1_level, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > release_mid_x1) { + phase = release_mid_x1; + } + // transition to 2nd release half + if (phase == release_mid_x1) { + phase = 0; + state = release2; + } + } + // release 2nd half + if (state == release2) { + // while in release phase + if (phase < release_mid_x2) { + level = lerp(0, release1_level, release_mid_x2, 0, phase); + phase += 1; + } + // reset phase if parameter was changed + if (phase > release_mid_x2) { + phase = release_mid_x2; + } + // reset + if (phase == release_mid_x2) { + phase = 0; + state = idle; + level = 0; + } + } + + return smooth(level); + } + + bool is_busy() { return state != 0; } + + void set_samplerate(double sampleRate) { + this->samplerate = sampleRate; + } + + // returns the x/y coordinates of the envelope points as a list for graphical representation. + std::array calc_coordinates() { + + float a_x = 0; + float a_y = 0; + + float b_x = attack1_rate; + float b_y = attack1_level; + + float c_x = b_x + attack2_rate; + float c_y = 1; + + float d_x = c_x + hold_rate; + float d_y = 1; + + float e_x = d_x + decay1_rate; + float e_y = decay1_level; + + float f_x = e_x + decay2_rate; + float f_y = sustain_level; + + float g_x = f_x + 125; + float g_y = sustain_level; + + float h_x = g_x + release1_rate; + float h_y = release1_level; + + float i_x = h_x + release2_rate; + float i_y = 0; + + float total = i_x; + + return { + a_x, + a_y, + b_x / total, + b_y, + c_x / total, + c_y, + d_x / total, + d_y, + e_x / total, + e_y, + f_x / total, + f_y, + g_x / total, + g_y, + h_x / total, + h_y, + i_x / total, + i_y + }; + } + +private: + double samplerate; + int phase; + float level; + env_state state; + float start_level; + float h1; + float h2; + float h3; + + float lerp(float x1, float y1, float x2, float y2, float x) { return y1 + (((x - x1) * (y2 - y1)) / (x2 - x1)); } + + float smooth(float sample) { + h3 = h2; + h2 = h1; + h1 = sample; + + return (h1 + h2 + h3) / 3.f; + } + + float ms_to_samples(float ms) { return ms * samplerate / 1000.f; } +}; +} \ No newline at end of file diff --git a/synth/tx_operator.h b/synth/tx_operator.h new file mode 100644 index 0000000..0070a7d --- /dev/null +++ b/synth/tx_operator.h @@ -0,0 +1,39 @@ +#pragma once +#include "tx_sineosc.h" +#include "tx_envelope.h" + +namespace trnr::lib::synth { +class tx_operator { +public: + tx_operator(double samplerate) + : ratio { 1 } + , amplitude { 1.0f } + , envelope(samplerate) + , oscillator(samplerate) + { + } + + tx_envelope envelope; + tx_sineosc oscillator; + float ratio; + float amplitude; + + float process_sample(const bool& gate, const bool& trigger, const float& frequency, const float& velocity, const float& pm = 0) { + + float env = envelope.process_sample(gate, trigger); + + // drifts and sounds better! + if (envelope.is_busy()) { + double osc = oscillator.process_sample(trigger, frequency, pm); + return osc * env * velocity; + } else { + return 0.; + } + } + + void set_samplerate(double samplerate) { + this->envelope.set_samplerate(samplerate); + this->oscillator.set_samplerate(samplerate); + } +}; +} diff --git a/synth/tx_sineosc.h b/synth/tx_sineosc.h new file mode 100644 index 0000000..3f3f640 --- /dev/null +++ b/synth/tx_sineosc.h @@ -0,0 +1,95 @@ +#pragma once +#include + +namespace trnr::lib::synth { + +class tx_sineosc { +public: + bool phase_reset; + + tx_sineosc(double _samplerate) + : samplerate { _samplerate } + , phase_resolution { 16.f } + , phase { 0. } + , history { 0. } + , phase_reset { false } + { + } + + void set_phase_resolution(float res) { + phase_resolution = powf(2, res); + } + + float process_sample(bool trigger, float frequency, float phase_modulation = 0.f) { + if (trigger && phase_reset) { + phase = 0.0; + } + + float lookup_phase = phase + phase_modulation; + wrap(lookup_phase); + phase += frequency / samplerate; + wrap(phase); + + redux(lookup_phase); + + float output = sine(lookup_phase * 4096.); + + filter(output); + return output; + } + + void set_samplerate(double _samplerate) { + this->samplerate = _samplerate; + } + +private: + double samplerate; + float phase_resolution; + float phase; + float history; + + float sine(float x) { + // x is scaled 0<=x<4096 + const float a = -0.40319426317E-08; + const float b = 0.21683205691E+03; + const float c = 0.28463350538E-04; + const float d = -0.30774648337E-02; + float y; + + bool negate = false; + if (x > 2048) { + negate = true; + x -= 2048; + } + if (x > 1024) + x = 2048 - x; + y = (a + x) / (b + c * x * x) + d * x; + if (negate) + return (float)(-y); + else + return (float)y; + } + + float wrap(float& phase) { + while (phase < 0.) + phase += 1.; + + while (phase >= 1.) + phase -= 1.; + + return phase; + } + + float filter(float& value) { + value = 0.5 * (value + history); + history = value; + return value; + } + + float redux(float& value) + { + value = static_cast(value * phase_resolution) / phase_resolution; + return value; + } +}; +} \ No newline at end of file diff --git a/synth/tx_voice.h b/synth/tx_voice.h new file mode 100644 index 0000000..2f06038 --- /dev/null +++ b/synth/tx_voice.h @@ -0,0 +1,166 @@ +#pragma once +#include "tx_sineosc.h" +#include "tx_envelope.h" +#include "tx_operator.h" + +namespace trnr::lib::synth { + +class tx_voice { +public: + tx_voice(double samplerate) + : algorithm { 0 } + , pitch_env_amt { 0.f } + , feedback_amt { 0.f } + , pitch_env(samplerate) + , feedback_osc(samplerate) + , op1(samplerate) + , op2(samplerate) + , op3(samplerate) + , bit_resolution(12.f) + { + } + + bool gate = false; + bool trigger = false; + float frequency = 100.f; + float velocity = 1.f; + + int algorithm; + float pitch_env_amt; + float feedback_amt; + float bit_resolution; + tx_sineosc feedback_osc; + tx_envelope pitch_env; + tx_operator op1; + tx_operator op2; + tx_operator op3; + + float process_sample() { + float pitch_env_signal = pitch_env.process_sample(gate, trigger) * pitch_env_amt; + float pitched_freq = frequency + pitch_env_signal; + + float output = 0.f; + + // mix operator signals according to selected algorithm + switch (algorithm) { + case 0: + output = calc_algo1(pitched_freq); + break; + case 1: + output = calc_algo2(pitched_freq); + break; + case 2: + output = calc_algo3(pitched_freq); + break; + case 3: + output = calc_algo4(pitched_freq); + break; + default: + output = calc_algo1(pitched_freq); + break; + } + + // reset trigger + trigger = false; + + return redux(output, bit_resolution); + } + + bool is_busy() { return gate || op1.envelope.is_busy() || op2.envelope.is_busy() || op3.envelope.is_busy(); } + + void set_samplerate(double samplerate) { + pitch_env.set_samplerate(samplerate); + feedback_osc.set_samplerate(samplerate); + op1.set_samplerate(samplerate); + op2.set_samplerate(samplerate); + op3.set_samplerate(samplerate); + } + + void set_phase_reset(bool phase_reset) { + op1.oscillator.phase_reset = phase_reset; + op2.oscillator.phase_reset = phase_reset; + op3.oscillator.phase_reset = phase_reset; + feedback_osc.phase_reset = phase_reset; + } + +private: + const float MOD_INDEX_COEFF = 4.f; + + float calc_algo1(const float frequency) { + float fb_freq = frequency * op3.ratio; + float fb_mod_index = (feedback_amt * MOD_INDEX_COEFF); + float fb_signal = feedback_osc.process_sample(trigger, fb_freq) * fb_mod_index; + + float op3_Freq = frequency * op3.ratio; + float op3_mod_index = (op3.amplitude * MOD_INDEX_COEFF); + float op3_signal = op3.process_sample(gate, trigger, op3_Freq, velocity, fb_signal) * op3_mod_index; + + float op2_freq = frequency * op2.ratio; + float op2_mod_index = (op2.amplitude * MOD_INDEX_COEFF); + float op2_signal = op2.process_sample(gate, trigger, op2_freq, velocity, op3_signal) * op2_mod_index; + + float op1_freq = frequency * op1.ratio; + return op1.process_sample(gate, trigger, op1_freq, velocity, op2_signal) * op1.amplitude; + } + + float calc_algo2(const float frequency) { + float fb_freq = frequency * op3.ratio; + float fb_mod_index = (feedback_amt * MOD_INDEX_COEFF); + float fb_signal = feedback_osc.process_sample(trigger, fb_freq) * fb_mod_index; + + float op3_freq = frequency * op3.ratio; + float op3_signal = op3.process_sample(gate, trigger, op3_freq, velocity, fb_signal) * op3.amplitude; + + float op2_freq = frequency * op2.ratio; + float op2_mod_index = (op2.amplitude * MOD_INDEX_COEFF); + float op2_signal = op2.process_sample(gate, trigger, op2_freq, velocity) * op2_mod_index; + + float op1_freq = frequency * op1.ratio; + float op1_signal = op1.process_sample(gate, trigger, op1_freq, velocity, op2_signal) * op1.amplitude; + + return op1_signal + op3_signal; + } + + float calc_algo3(const float frequency) { + float fb_freq = frequency * op3.ratio; + float fb_mod_index = (feedback_amt * MOD_INDEX_COEFF); + float fb_signal = feedback_osc.process_sample(trigger, fb_freq) * fb_mod_index; + + float op3_freq = frequency * op3.ratio; + float op3_signal = op3.process_sample(gate, trigger, op3_freq, velocity, fb_signal) * op3.amplitude; + + float op2_freq = frequency * op2.ratio; + float op2_signal = op2.process_sample(gate, trigger, op2_freq, velocity) * op2.amplitude; + + float op1_freq = frequency * op1.ratio; + float op1_signal = op1.process_sample(gate, trigger, op1_freq, velocity) * op1.amplitude; + + return op1_signal + op2_signal + op3_signal; + } + + float calc_algo4(const float frequency) { + float fb_freq = frequency * op3.ratio; + float fb_mod_index = (feedback_amt * MOD_INDEX_COEFF); + float fb_signal = feedback_osc.process_sample(trigger, fb_freq) * fb_mod_index; + + float op3_freq = frequency * op3.ratio; + float op3_mod_index = (op3.amplitude * MOD_INDEX_COEFF); + float op3_signal = op3.process_sample(gate, trigger, op3_freq, velocity, fb_signal) * op3_mod_index; + + float op2_freq = frequency * op2.ratio; + float op2_mod_index = (op2.amplitude * MOD_INDEX_COEFF); + float op2_signal = op2.process_sample(gate, trigger, op2_freq, velocity) * op2_mod_index; + + float op1_freq = frequency * op1.ratio; + return op1.process_sample(gate, trigger, op1_freq, velocity, op2_signal + op3_signal) * op1.amplitude; + } + + float redux(float& value, float resolution) + { + float res = powf(2, resolution); + value = roundf(value * res) / res; + + return value; + } +}; +} \ No newline at end of file diff --git a/util/audio_math.h b/util/audio_math.h new file mode 100644 index 0000000..13efef5 --- /dev/null +++ b/util/audio_math.h @@ -0,0 +1,11 @@ +#include + +namespace trnr::lib::util { + static inline double lin_2_db(double lin) { + return 20 * log(lin); + } + + static inline double db_2_lin(double db) { + return pow(10, db/20); + } +} \ No newline at end of file