avfilter/af_firequalizer: add min_phase option

Signed-off-by: Muhammad Faiz <mfcc64@gmail.com>

doc/filters.texi

@@ -2667,6 +2667,9 @@ Default is linlog.
 @item fft2
 Enable 2-channel convolution using complex FFT. This improves speed significantly.
 Default is disabled.
+
+@item min_phase
+Enable minimum phase impulse response. Default is disabled.
 @end table
 
 @subsection Examples

libavfilter/af_firequalizer.c

@@ -70,13 +70,17 @@ typedef struct FIREqualizerContext {
     RDFTContext   *rdft;
     RDFTContext   *irdft;
     FFTContext    *fft_ctx;
+    RDFTContext   *cepstrum_rdft;
+    RDFTContext   *cepstrum_irdft;
     int           analysis_rdft_len;
     int           rdft_len;
+    int           cepstrum_len;
 
     float         *analysis_buf;
     float         *dump_buf;
     float         *kernel_tmp_buf;
     float         *kernel_buf;
+    float         *cepstrum_buf;
     float         *conv_buf;
     OverlapIndex  *conv_idx;
     int           fir_len;
@@ -99,6 +103,7 @@ typedef struct FIREqualizerContext {
     char          *dumpfile;
     int           dumpscale;
     int           fft2;
+    int           min_phase;
 
     int           nb_gain_entry;
     int           gain_entry_err;
@@ -135,6 +140,7 @@ static const AVOption firequalizer_options[] = {
     { "dumpfile", "set dump file", OFFSET(dumpfile), AV_OPT_TYPE_STRING, { .str = NULL }, 0, 0, FLAGS },
     { "dumpscale", "set dump scale", OFFSET(dumpscale), AV_OPT_TYPE_INT, { .i64 = SCALE_LINLOG }, 0, NB_SCALE-1, FLAGS, "scale" },
     { "fft2", "set 2-channels fft", OFFSET(fft2), AV_OPT_TYPE_BOOL, { .i64 = 0 }, 0, 1, FLAGS },
+    { "min_phase", "set minimum phase mode", OFFSET(min_phase), AV_OPT_TYPE_BOOL, { .i64 = 0 }, 0, 1, FLAGS },
     { NULL }
 };
 
@@ -147,13 +153,18 @@ static void common_uninit(FIREqualizerContext *s)
     av_rdft_end(s->rdft);
     av_rdft_end(s->irdft);
     av_fft_end(s->fft_ctx);
+    av_rdft_end(s->cepstrum_rdft);
+    av_rdft_end(s->cepstrum_irdft);
     s->analysis_rdft = s->analysis_irdft = s->rdft = s->irdft = NULL;
     s->fft_ctx = NULL;
+    s->cepstrum_rdft = NULL;
+    s->cepstrum_irdft = NULL;
 
     av_freep(&s->analysis_buf);
     av_freep(&s->dump_buf);
     av_freep(&s->kernel_tmp_buf);
     av_freep(&s->kernel_buf);
+    av_freep(&s->cepstrum_buf);
     av_freep(&s->conv_buf);
     av_freep(&s->conv_idx);
 }
@@ -235,6 +246,46 @@ static void fast_convolute(FIREqualizerContext *av_restrict s, const float *av_r
     }
 }
 
+static void fast_convolute_nonlinear(FIREqualizerContext *av_restrict s, const float *av_restrict kernel_buf,
+                                     float *av_restrict conv_buf, OverlapIndex *av_restrict idx,
+                                     float *av_restrict data, int nsamples)
+{
+    if (nsamples <= s->nsamples_max) {
+        float *buf = conv_buf + idx->buf_idx * s->rdft_len;
+        float *obuf = conv_buf + !idx->buf_idx * s->rdft_len + idx->overlap_idx;
+        int k;
+
+        memcpy(buf, data, nsamples * sizeof(*data));
+        memset(buf + nsamples, 0, (s->rdft_len - nsamples) * sizeof(*data));
+        av_rdft_calc(s->rdft, buf);
+
+        buf[0] *= kernel_buf[0];
+        buf[1] *= kernel_buf[1];
+        for (k = 2; k < s->rdft_len; k += 2) {
+            float re, im;
+            re = buf[k] * kernel_buf[k] - buf[k+1] * kernel_buf[k+1];
+            im = buf[k] * kernel_buf[k+1] + buf[k+1] * kernel_buf[k];
+            buf[k] = re;
+            buf[k+1] = im;
+        }
+
+        av_rdft_calc(s->irdft, buf);
+        for (k = 0; k < s->rdft_len - idx->overlap_idx; k++)
+            buf[k] += obuf[k];
+        memcpy(data, buf, nsamples * sizeof(*data));
+        idx->buf_idx = !idx->buf_idx;
+        idx->overlap_idx = nsamples;
+    } else {
+        while (nsamples > s->nsamples_max * 2) {
+            fast_convolute_nonlinear(s, kernel_buf, conv_buf, idx, data, s->nsamples_max);
+            data += s->nsamples_max;
+            nsamples -= s->nsamples_max;
+        }
+        fast_convolute_nonlinear(s, kernel_buf, conv_buf, idx, data, nsamples/2);
+        fast_convolute_nonlinear(s, kernel_buf, conv_buf, idx, data + nsamples/2, nsamples - nsamples/2);
+    }
+}
+
 static void fast_convolute2(FIREqualizerContext *av_restrict s, const float *av_restrict kernel_buf, FFTComplex *av_restrict conv_buf,
                             OverlapIndex *av_restrict idx, float *av_restrict data0, float *av_restrict data1, int nsamples)
 {
@@ -310,22 +361,32 @@ static void dump_fir(AVFilterContext *ctx, FILE *fp, int ch)
     double delay = s->zero_phase ? 0.0 : (double) center / rate;
     double vx, ya, yb;
 
+    if (!s->min_phase) {
     s->analysis_buf[0] *= s->rdft_len/2;
     for (x = 1; x <= center; x++) {
         s->analysis_buf[x] *= s->rdft_len/2;
         s->analysis_buf[s->analysis_rdft_len - x] *= s->rdft_len/2;
     }
+    } else {
+        for (x = 0; x < s->fir_len; x++)
+            s->analysis_buf[x] *= s->rdft_len/2;
+    }
 
     if (ch)
         fprintf(fp, "\n\n");
 
     fprintf(fp, "# time[%d] (time amplitude)\n", ch);
 
+    if (!s->min_phase) {
     for (x = center; x > 0; x--)
         fprintf(fp, "%15.10f %15.10f\n", delay - (double) x / rate, (double) s->analysis_buf[s->analysis_rdft_len - x]);
 
     for (x = 0; x <= center; x++)
         fprintf(fp, "%15.10f %15.10f\n", delay + (double)x / rate , (double) s->analysis_buf[x]);
+    } else {
+        for (x = 0; x < s->fir_len; x++)
+            fprintf(fp, "%15.10f %15.10f\n", (double)x / rate, (double) s->analysis_buf[x]);
+    }
 
     av_rdft_calc(s->analysis_rdft, s->analysis_buf);
 
@@ -337,7 +398,9 @@ static void dump_fir(AVFilterContext *ctx, FILE *fp, int ch)
         if (xlog)
             vx = log2(0.05*vx);
         ya = s->dump_buf[i];
-        yb = s->analysis_buf[i];
+        yb = s->min_phase && (i > 1) ? hypotf(s->analysis_buf[i], s->analysis_buf[i+1]) : s->analysis_buf[i];
+        if (s->min_phase)
+            yb = fabs(yb);
         if (ylog) {
             ya = 20.0 * log10(fabs(ya));
             yb = 20.0 * log10(fabs(yb));
@@ -487,6 +550,53 @@ enum VarOffset {
     VAR_NB
 };
 
+static void generate_min_phase_kernel(FIREqualizerContext *s, float *rdft_buf)
+{
+    int k, cepstrum_len = s->cepstrum_len, rdft_len = s->rdft_len;
+    double norm = 2.0 / cepstrum_len;
+
+    memset(s->cepstrum_buf, 0, cepstrum_len * sizeof(*s->cepstrum_buf));
+    memcpy(s->cepstrum_buf, rdft_buf, rdft_len/2 * sizeof(*rdft_buf));
+    memcpy(s->cepstrum_buf + cepstrum_len - rdft_len/2, rdft_buf + rdft_len/2, rdft_len/2  * sizeof(*rdft_buf));
+
+    av_rdft_calc(s->cepstrum_rdft, s->cepstrum_buf);
+
+    s->cepstrum_buf[0] = log(FFMAX(s->cepstrum_buf[0], 1e-10));
+    s->cepstrum_buf[1] = log(FFMAX(s->cepstrum_buf[1], 1e-10));
+
+    for (k = 2; k < cepstrum_len; k += 2) {
+        s->cepstrum_buf[k] = log(FFMAX(s->cepstrum_buf[k], 1e-10));
+        s->cepstrum_buf[k+1] = 0;
+    }
+
+    av_rdft_calc(s->cepstrum_irdft, s->cepstrum_buf);
+
+    memset(s->cepstrum_buf + cepstrum_len/2 + 1, 0, (cepstrum_len/2 - 1) * sizeof(*s->cepstrum_buf));
+    for (k = 1; k < cepstrum_len/2; k++)
+        s->cepstrum_buf[k] *= 2;
+
+    av_rdft_calc(s->cepstrum_rdft, s->cepstrum_buf);
+
+    s->cepstrum_buf[0] = exp(s->cepstrum_buf[0] * norm) * norm;
+    s->cepstrum_buf[1] = exp(s->cepstrum_buf[1] * norm) * norm;
+    for (k = 2; k < cepstrum_len; k += 2) {
+        double mag = exp(s->cepstrum_buf[k] * norm) * norm;
+        double ph = s->cepstrum_buf[k+1] * norm;
+        s->cepstrum_buf[k] = mag * cos(ph);
+        s->cepstrum_buf[k+1] = mag * sin(ph);
+    }
+
+    av_rdft_calc(s->cepstrum_irdft, s->cepstrum_buf);
+    memset(rdft_buf, 0, s->rdft_len * sizeof(*rdft_buf));
+    memcpy(rdft_buf, s->cepstrum_buf, s->fir_len * sizeof(*rdft_buf));
+
+    if (s->dumpfile) {
+        memset(s->analysis_buf, 0, s->analysis_rdft_len * sizeof(*s->analysis_buf));
+        memcpy(s->analysis_buf, s->cepstrum_buf, s->fir_len * sizeof(*s->analysis_buf));
+    }
+
+}
+
 static int generate_kernel(AVFilterContext *ctx, const char *gain, const char *gain_entry)
 {
     FIREqualizerContext *s = ctx->priv;
@@ -549,7 +659,7 @@ static int generate_kernel(AVFilterContext *ctx, const char *gain, const char *g
             if (xlog)
                 vars[VAR_F] = log2(0.05 * vars[VAR_F]);
             result = av_expr_eval(gain_expr, vars, ctx);
-            s->analysis_buf[2*k] = ylog ? pow(10.0, 0.05 * result) : result;
+            s->analysis_buf[2*k] = ylog ? pow(10.0, 0.05 * result) : s->min_phase ? fabs(result) : result;
             s->analysis_buf[2*k+1] = 0.0;
         }
 
@@ -604,6 +714,8 @@ static int generate_kernel(AVFilterContext *ctx, const char *gain, const char *g
         memset(s->analysis_buf + center + 1, 0, (s->analysis_rdft_len - s->fir_len) * sizeof(*s->analysis_buf));
         memcpy(rdft_buf, s->analysis_buf, s->rdft_len/2 * sizeof(*s->analysis_buf));
         memcpy(rdft_buf + s->rdft_len/2, s->analysis_buf + s->analysis_rdft_len - s->rdft_len/2, s->rdft_len/2 * sizeof(*s->analysis_buf));
+        if (s->min_phase)
+            generate_min_phase_kernel(s, rdft_buf);
         av_rdft_calc(s->rdft, rdft_buf);
 
         for (k = 0; k < s->rdft_len; k++) {
@@ -616,10 +728,12 @@ static int generate_kernel(AVFilterContext *ctx, const char *gain, const char *g
             }
         }
 
+        if (!s->min_phase) {
         rdft_buf[s->rdft_len-1] = rdft_buf[1];
         for (k = 0; k < s->rdft_len/2; k++)
             rdft_buf[k] = rdft_buf[2*k];
         rdft_buf[s->rdft_len/2] = rdft_buf[s->rdft_len-1];
+        }
 
         if (dump_fp)
             dump_fir(ctx, dump_fp, ch);
@@ -670,6 +784,25 @@ static int config_input(AVFilterLink *inlink)
     if (s->fft2 && !s->multi && inlink->channels > 1 && !(s->fft_ctx = av_fft_init(rdft_bits, 0)))
         return AVERROR(ENOMEM);
 
+    if (s->min_phase) {
+        int cepstrum_bits = rdft_bits + 2;
+        if (cepstrum_bits > RDFT_BITS_MAX) {
+            av_log(ctx, AV_LOG_ERROR, "too large delay, please decrease it.\n");
+            return AVERROR(EINVAL);
+        }
+
+        cepstrum_bits = FFMIN(RDFT_BITS_MAX, cepstrum_bits + 1);
+        s->cepstrum_rdft = av_rdft_init(cepstrum_bits, DFT_R2C);
+        s->cepstrum_irdft = av_rdft_init(cepstrum_bits, IDFT_C2R);
+        if (!s->cepstrum_rdft || !s->cepstrum_irdft)
+            return AVERROR(ENOMEM);
+
+        s->cepstrum_len = 1 << cepstrum_bits;
+        s->cepstrum_buf = av_malloc_array(s->cepstrum_len, sizeof(*s->cepstrum_buf));
+        if (!s->cepstrum_buf)
+            return AVERROR(ENOMEM);
+    }
+
     for ( ; rdft_bits <= RDFT_BITS_MAX; rdft_bits++) {
         s->analysis_rdft_len = 1 << rdft_bits;
         if (inlink->sample_rate <= s->accuracy * s->analysis_rdft_len)
@@ -712,6 +845,7 @@ static int filter_frame(AVFilterLink *inlink, AVFrame *frame)
     FIREqualizerContext *s = ctx->priv;
     int ch;
 
+    if (!s->min_phase) {
     for (ch = 0; ch + 1 < inlink->channels && s->fft_ctx; ch += 2) {
         fast_convolute2(s, s->kernel_buf, (FFTComplex *)(s->conv_buf + 2 * ch * s->rdft_len),
                         s->conv_idx + ch, (float *) frame->extended_data[ch],
@@ -723,11 +857,18 @@ static int filter_frame(AVFilterLink *inlink, AVFrame *frame)
                        s->conv_buf + 2 * ch * s->rdft_len, s->conv_idx + ch,
                        (float *) frame->extended_data[ch], frame->nb_samples);
     }
+    } else {
+        for (ch = 0; ch < inlink->channels; ch++) {
+            fast_convolute_nonlinear(s, s->kernel_buf + (s->multi ? ch * s->rdft_len : 0),
+                                     s->conv_buf + 2 * ch * s->rdft_len, s->conv_idx + ch,
+                                     (float *) frame->extended_data[ch], frame->nb_samples);
+        }
+    }
 
     s->next_pts = AV_NOPTS_VALUE;
     if (frame->pts != AV_NOPTS_VALUE) {
         s->next_pts = frame->pts + av_rescale_q(frame->nb_samples, av_make_q(1, inlink->sample_rate), inlink->time_base);
-        if (s->zero_phase)
+        if (s->zero_phase && !s->min_phase)
             frame->pts -= av_rescale_q(s->fir_len/2, av_make_q(1, inlink->sample_rate), inlink->time_base);
     }
     s->frame_nsamples_max = FFMAX(s->frame_nsamples_max, frame->nb_samples);