+static void get_response(int channel, int format, double w,
+ const double *b, const double *a,
+ int nb_b, int nb_a, double *r, double *i)
+{
+ double realz, realp;
+ double imagz, imagp;
+ double real, imag;
+ double div;
+
+ if (format == 0) {
+ realz = 0., realp = 0.;
+ imagz = 0., imagp = 0.;
+ for (int x = 0; x < nb_a; x++) {
+ realz += cos(-x * w) * a[x];
+ imagz += sin(-x * w) * a[x];
+ }
+
+ for (int x = 0; x < nb_b; x++) {
+ realp += cos(-x * w) * b[x];
+ imagp += sin(-x * w) * b[x];
+ }
+
+ div = realp * realp + imagp * imagp;
+ real = (realz * realp + imagz * imagp) / div;
+ imag = (imagz * realp - imagp * realz) / div;
+ } else {
+ real = 1;
+ imag = 0;
+ for (int x = 0; x < nb_a; x++) {
+ double ore, oim, re, im;
+
+ re = cos(w) - a[2 * x];
+ im = sin(w) - a[2 * x + 1];
+
+ ore = real;
+ oim = imag;
+
+ real = ore * re - oim * im;
+ imag = ore * im + oim * re;
+ }
+
+ for (int x = 0; x < nb_b; x++) {
+ double ore, oim, re, im;
+
+ re = cos(w) - b[2 * x];
+ im = sin(w) - b[2 * x + 1];
+
+ ore = real;
+ oim = imag;
+ div = re * re + im * im;
+
+ real = (ore * re + oim * im) / div;
+ imag = (oim * re - ore * im) / div;
+ }
+ }
+
+ *r = real;
+ *i = imag;
+}
+