I've spent 4+ hours trying to convert to C++ the Matlab code from this excellent work. To speed-up the process I decided to use good old friend ChatGPT and here I am, with headache.
After a 4+ hours I have some code that works at test program, but gives me NaNs wtih actual sample values, altough they are the same. What I need is a function that takes float x, and distortion coef (0..1), and return me the distorted x. I call it in the cycle with other computations so I need exactly such function.
If someone interesting, me with ChatGPT ended up with the following code. Actually, I gave up and will make more simply distortion that I can understand, not the physical modelling, etc.:
Code: Select all
#include <cmath>
#include <cmath>
#include <vector>
#include <iostream>
template <typename T>
T clampTo01(const T& value) {
if (value < 0.0) {
return 0.0;
} else if (value > 1.0) {
return 1.0;
} else {
return value;
}
}
class DistortionProcessor {
public:
float fs;
DistortionProcessor() : fs(44100) {}
//Matlab's filter works in that way?
float filter(const std::vector<double>& b, const std::vector<double>& a, const float x, std::vector<double>& state) {
double y = b[0] * x;
for (size_t i = 1; i < b.size(); i++) {
if (i <= state.size()) {
y += b[i] * x - a[i] * state[i - 1];
}
}
for (size_t i = 1; i < state.size(); i++) {
state[i] = state[i - 1];
}
if (state.size() > 0) {
state[0] = x;
}
if (a[0] != 1.0) {
for (size_t i = 0; i < state.size(); i++) {
state[i] /= a[0];
}
}
return static_cast<float>(y);
}
float bjtfilt(float x, float fs) {
double coeff = M_PI / fs;
// filter's coefs
std::vector<double> B = {1800.0 * pow(coeff, 2) + 603.0 * coeff + 1.0,
3600.0 * pow(coeff, 2) - 2.0,
1800.0 * pow(coeff, 2) - 603.0 * coeff + 1.0};
double amp = pow(10, (36.0 / 20.0));
for (size_t i = 0; i < B.size(); i++) {
B[i] *= amp;
}
std::vector<double> state(B.size() - 1, 0.0);
return filter(B, B, x, state);
}
// Op-Amp Gain Stage
float opampfilt(float x, float DIST) {
double Rt = 100000 * DIST;
double Rb = 100000 * (1 - DIST) + 4700;
double Cz = 0.000001;
double Cc = 0.000000000250;
double c = 2 * fs;
double ab0 = 1 / (Rt * Cc * Rb * Cz);
double a1 = 1 / (Rb * Cz) + 1 / (Rt * Cc);
double b1 = a1 + 1 / (Rb * Cc);
double B0 = ab0 + b1 * c + pow(c, 2);
double B1 = 2 * ab0 - 2 * pow(c, 2);
double B2 = ab0 - b1 * c + pow(c, 2);
double A0 = ab0 + a1 * c + pow(c, 2);
double A1 = B1;
double A2 = ab0 - a1 * c + pow(c, 2);
double y = B0 * x;
if (prev_1_initialized_) y += B1 * prev_1_;
if (prev_2_initialized_) y += B2 * prev_2_;
if (prev_1_initialized_) y -= A1 * prev_1_;
if (prev_2_initialized_) y -= A2 * prev_2_;
prev_2_ = prev_1_;
prev_1_ = y;
prev_1_initialized_ = true;
prev_2_initialized_ = true;
return static_cast<float>(y / A0);
}
// Diode clipper
float diodeclip(float x) {
float n = 2.5;
return static_cast<float>(x / pow(1 + pow(fabs(x), n), 1.0 / n));
}
float processSignal(float input, float DIST)
{
float a = bjtfilt(input, fs);
std::cout << "a:" << a << std::endl;
float b = opampfilt(a, DIST);
std::cout << "b:" << a << std::endl;
// Ограничиваем значения в диапазоне от 0 до 1
float c = clampTo01(b);
std::cout << "c:" << c << std::endl;
float d = diodeclip(c);
return d;
}
private:
float prev_1_ = 0.0;
float prev_2_ = 0.0;
bool prev_1_initialized_ = false;
bool prev_2_initialized_ = false;
};
#include <iostream>
int main() {
DistortionProcessor processor;
processor.fs = 44100;
float x[5] = {0.00897738, 0.00941722, 0.0124338, 0.0139674, 0.0133841};
float DIST = 0.5;
for (size_t i = 0; i < 5; i++) {
float processed = processor.processSignal(x[i], DIST);
std::cout << "Processed signal: " << processed << std::endl;
}
return 0;
}