+static double distance(double x0, double x1, double y0, double y1)
+{
+ return hypot(x0 - x1, y0 - y1);
+}
+
+static void get_response(int channel, int format, double w,
+ const double *b, const double *a,
+ int nb_b, int nb_a, double *magnitude, double *phase)
+{
+ 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;
+
+ *magnitude = hypot(real, imag);
+ *phase = atan2(imag, real);
+ } else {
+ double p = 1., z = 1.;
+ double acc = 0.;
+
+ for (int x = 0; x < nb_a; x++) {
+ z *= distance(cos(w), a[2 * x], sin(w), a[2 * x + 1]);
+ acc += atan2(sin(w) - a[2 * x + 1], cos(w) - a[2 * x]);
+ }
+
+ for (int x = 0; x < nb_b; x++) {
+ p *= distance(cos(w), b[2 * x], sin(w), b[2 * x + 1]);
+ acc -= atan2(sin(w) - b[2 * x + 1], cos(w) - b[2 * x]);
+ }
+
+ *magnitude = z / p;
+ *phase = acc;
+ }
+}
+
+static void draw_response(AVFilterContext *ctx, AVFrame *out, int sample_rate)