|
| 1 | +// Sallen-Key low-pass filter benchmark |
| 2 | +// Process 2 million samples with per-sample coefficient recomputation |
| 3 | +// Based on Andrew Simper's nodal filter equations (SkfLinearTrapOptimised2.pdf) |
| 4 | + |
| 5 | +#include <math.h> |
| 6 | +#include <stdio.h> |
| 7 | +#include <time.h> |
| 8 | + |
| 9 | +int main(void) { |
| 10 | + int n = 2000000; |
| 11 | + double inv_sr = 1.0 / 44100.0; |
| 12 | + double f_min = log(50.0); |
| 13 | + double f_max = log(16000.0); |
| 14 | + double f_range = f_max - f_min; |
| 15 | + |
| 16 | + double ic1eq = 0.0; |
| 17 | + double ic2eq = 0.0; |
| 18 | + double last_y = 0.0; |
| 19 | + |
| 20 | + double phase = 0.0; |
| 21 | + double freq = 440.0 / 44100.0; |
| 22 | + double two_pi = 2.0 * M_PI; |
| 23 | + |
| 24 | + double fc_val = 0.5; |
| 25 | + double res_val = 0.3; |
| 26 | + |
| 27 | + struct timespec t0, t1; |
| 28 | + clock_gettime(CLOCK_MONOTONIC, &t0); |
| 29 | + |
| 30 | + for (int i = 0; i < n; i++) { |
| 31 | + double x = sin(phase * two_pi); |
| 32 | + phase += freq; |
| 33 | + if (phase > 1.0) phase -= 1.0; |
| 34 | + |
| 35 | + double fc_norm = fc_val + 0.001 * x; |
| 36 | + if (fc_norm < 0.0) fc_norm = 0.0; |
| 37 | + if (fc_norm > 1.0) fc_norm = 1.0; |
| 38 | + double res_c = res_val; |
| 39 | + if (res_c < 0.0) res_c = 0.0; |
| 40 | + if (res_c > 1.0) res_c = 1.0; |
| 41 | + |
| 42 | + double cutoff = exp(fc_norm * f_range + f_min); |
| 43 | + double g = tan(M_PI * cutoff * inv_sr); |
| 44 | + double k = 2.0 * res_c; |
| 45 | + double g1 = 1.0 + g; |
| 46 | + double a0 = 1.0 / (g1 * g1 - g * k); |
| 47 | + double a1 = k * a0; |
| 48 | + double a2 = g1 * a0; |
| 49 | + double a3 = g * a2; |
| 50 | + double a4 = 1.0 / g1; |
| 51 | + double a5 = g * a4; |
| 52 | + |
| 53 | + double v1 = a1 * ic2eq + a2 * ic1eq + a3 * x; |
| 54 | + double v2 = a4 * ic2eq + a5 * v1; |
| 55 | + |
| 56 | + ic1eq = 2.0 * (v1 - k * v2) - ic1eq; |
| 57 | + ic2eq = 2.0 * v2 - ic2eq; |
| 58 | + |
| 59 | + last_y = v2; |
| 60 | + } |
| 61 | + |
| 62 | + clock_gettime(CLOCK_MONOTONIC, &t1); |
| 63 | + double elapsed = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) * 1e-9; |
| 64 | + fprintf(stderr, "exec: %.3fs\n", elapsed); |
| 65 | + printf("%d samples, last_y=%f\n", n, last_y); |
| 66 | + return 0; |
| 67 | +} |
0 commit comments