From 270068b5af281dae9280306b1621d1f2b7920adb Mon Sep 17 00:00:00 2001
From: Paul B Mahol <onemda@gmail.com>
Date: Sun, 14 Apr 2019 19:03:39 +0200
Subject: [PATCH] avfilter/af_acrossover: improve filter output

Makes sum always flat. Also faster.
---
 libavfilter/af_acrossover.c | 58 ++++++++++++-------------------------
 1 file changed, 19 insertions(+), 39 deletions(-)

diff --git a/libavfilter/af_acrossover.c b/libavfilter/af_acrossover.c
index b1a1155c29..31d7568abc 100644
--- a/libavfilter/af_acrossover.c
+++ b/libavfilter/af_acrossover.c
@@ -47,7 +47,6 @@ typedef struct BiquadContext {
 
 typedef struct CrossoverChannel {
     BiquadContext lp[MAX_BANDS][4];
-    BiquadContext hp[MAX_BANDS][4];
 } CrossoverChannel;
 
 typedef struct AudioCrossoverContext {
@@ -131,30 +130,16 @@ static av_cold int init(AVFilterContext *ctx)
     return ret;
 }
 
-static void set_lp(BiquadContext *b, float fc, float q, float sr)
+static void set_lp(BiquadContext *b, double fc, double q, double sr)
 {
-    double omega = (2.0 * M_PI * fc / sr);
+    double omega = 2.0 * M_PI * fc / sr;
     double sn = sin(omega);
     double cs = cos(omega);
-    double alpha = (sn / (2 * q));
-    double inv = (1.0 / (1.0 + alpha));
-
-    b->a2 = b->a0 = (inv * (1.0 - cs) * 0.5);
-    b->a1 = b->a0 + b->a0;
-    b->b1 = -2. * cs * inv;
-    b->b2 = (1. - alpha) * inv;
-}
-
-static void set_hp(BiquadContext *b, float fc, float q, float sr)
-{
-    double omega = 2 * M_PI * fc / sr;
-    double sn = sin(omega);
-    double cs = cos(omega);
-    double alpha = sn / (2 * q);
+    double alpha = sn / (2. * q);
     double inv = 1.0 / (1.0 + alpha);
 
-    b->a0 = inv * (1. + cs) / 2.;
-    b->a1 = -2. * b->a0;
+    b->a0 = (1. - cs) * 0.5 * inv;
+    b->a1 = (1. - cs) * inv;
     b->a2 = b->a0;
     b->b1 = -2. * cs * inv;
     b->b2 = (1. - alpha) * inv;
@@ -189,18 +174,13 @@ static int config_input(AVFilterLink *inlink)
     for (ch = 0; ch < inlink->channels; ch++) {
         for (band = 0; band <= s->nb_splits; band++) {
             set_lp(&s->xover[ch].lp[band][0], s->splits[band], q, sample_rate);
-            set_hp(&s->xover[ch].hp[band][0], s->splits[band], q, sample_rate);
 
             if (s->order > 1) {
                 set_lp(&s->xover[ch].lp[band][1], s->splits[band], 1.34, sample_rate);
-                set_hp(&s->xover[ch].hp[band][1], s->splits[band], 1.34, sample_rate);
                 set_lp(&s->xover[ch].lp[band][2], s->splits[band],    q, sample_rate);
-                set_hp(&s->xover[ch].hp[band][2], s->splits[band],    q, sample_rate);
                 set_lp(&s->xover[ch].lp[band][3], s->splits[band], 1.34, sample_rate);
-                set_hp(&s->xover[ch].hp[band][3], s->splits[band], 1.34, sample_rate);
             } else {
                 set_lp(&s->xover[ch].lp[band][1], s->splits[band], q, sample_rate);
-                set_hp(&s->xover[ch].hp[band][1], s->splits[band], q, sample_rate);
             }
         }
     }
@@ -275,23 +255,23 @@ static int filter_frame(AVFilterLink *inlink, AVFrame *in)
         const double *src = (const double *)in->extended_data[ch];
         CrossoverChannel *xover = &s->xover[ch];
 
-        for (band = 0; band < ctx->nb_outputs; band++) {
-            double *dst = (double *)frames[band]->extended_data[ch];
+        for (i = 0; i < in->nb_samples; i++) {
+            double sample = src[i], lo, hi;
 
-            for (i = 0; i < in->nb_samples; i++) {
-                dst[i] = src[i];
+            for (band = 0; band < ctx->nb_outputs; band++) {
+                double *dst = (double *)frames[band]->extended_data[ch];
 
-                for (f = 0; f < s->filter_count; f++) {
-                    if (band + 1 < ctx->nb_outputs) {
-                        BiquadContext *lp = &xover->lp[band][f];
-                        dst[i] = biquad_process(lp, dst[i]);
-                    }
-
-                    if (band - 1 >= 0) {
-                        BiquadContext *hp = &xover->hp[band - 1][f];
-                        dst[i] = biquad_process(hp, dst[i]);
-                    }
+                lo = sample;
+                for (f = 0; band + 1 < ctx->nb_outputs && f < s->filter_count; f++) {
+                    BiquadContext *lp = &xover->lp[band][f];
+                    lo = biquad_process(lp, lo);
                 }
+
+                hi = sample - lo;
+
+                dst[i] = lo;
+
+                sample = hi;
             }
         }
     }