From patchwork Tue Aug 22 18:38:15 2017 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Muhammad Faiz X-Patchwork-Id: 4798 Delivered-To: ffmpegpatchwork@gmail.com Received: by 10.2.137.29 with SMTP id o29csp1142832jaj; Tue, 22 Aug 2017 11:45:28 -0700 (PDT) X-Received: by 10.223.142.144 with SMTP id q16mr30739wrb.331.1503427528227; Tue, 22 Aug 2017 11:45:28 -0700 (PDT) ARC-Seal: i=1; a=rsa-sha256; t=1503427528; cv=none; d=google.com; s=arc-20160816; b=tEjgA8nFQ/qXIed8cDcF6xjPWJjcyY5CqBk0tX8gAtNJgdGTDKkeLG1ezOCvoexe2t BtFKg+9vFM+RrXZhDsPfugHEB7v4aUiNRhQDOF1JNoijDFdS4N3b/fyMHFdnGH7GYJvb nguS2QB6z8qxnMGGOzY09sfCnjMhO6NlP0y/Jkyuit8xF4yUjBbQGBu9REyUyoAY5xGk LjJJcY9Zo0oi8iaO0pT06bmj4g9sT/zZs7ZRSkETDiYXAMPk4Q/PZqr+j1K7dlTqIGkG I+FIeiNw2W9TSsNfcTY8wxmH22kYM0/7oSBzkV+3DQLdFVLcX7AjJ573PfuO9aMelZXx bRsQ== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=arc-20160816; h=sender:errors-to:content-transfer-encoding:reply-to:list-subscribe :list-help:list-post:list-archive:list-unsubscribe:list-id :precedence:subject:mime-version:message-id:date:to:from :dkim-signature:delivered-to:arc-authentication-results; bh=Bz4LZi5K0txu0xYwjt+BROPUsXEaIWPyHry000Hgy9I=; b=NfHHUof/ZHlbS0qzt8/89KheaLAZ0eRl4EBajJTVbNeIc76qLY7XprzkmlGUaVEhxG 9mYDMRVUAGnO2uqz+H5X/aDs4JDsuEtnMZ6xf1CK1geEmsSnNlgGLY44zY0ksjR9Xo0Q J0sdHmoLN5EiSFZL0i8wFLwPQzrJXnTvXJ62b+1kA0nl1Hqz90XMrEUa57qlROl/FhDM ov6pkXP2bxJfJVx9yZ8m6/Xd7XQt/0U2N9NlMJydyiHrZYRVlWtM5spuXcjv5FMifDmV bCYA/dDUCoPEqlJftxWqSG/MUutA8SjlcHohjD9XWvJ//Ob7W0cpajdIWXGhv7MFQshB R6xA== ARC-Authentication-Results: i=1; mx.google.com; dkim=neutral (body hash did not verify) header.i=@gmail.com header.s=20161025 header.b=bmzTaIaJ; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org; dmarc=fail (p=NONE sp=NONE dis=NONE) header.from=gmail.com Return-Path: Received: from ffbox0-bg.mplayerhq.hu (ffbox0-bg.ffmpeg.org. [79.124.17.100]) by mx.google.com with ESMTP id m199si296568wma.178.2017.08.22.11.45.25; Tue, 22 Aug 2017 11:45:28 -0700 (PDT) Received-SPF: pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) client-ip=79.124.17.100; Authentication-Results: mx.google.com; dkim=neutral (body hash did not verify) header.i=@gmail.com header.s=20161025 header.b=bmzTaIaJ; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org; dmarc=fail (p=NONE sp=NONE dis=NONE) header.from=gmail.com Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id 4656E689BEC; Tue, 22 Aug 2017 21:45:15 +0300 (EEST) X-Original-To: ffmpeg-devel@ffmpeg.org Delivered-To: ffmpeg-devel@ffmpeg.org Received: from mail-pg0-f66.google.com (mail-pg0-f66.google.com [74.125.83.66]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id 67339689824 for ; Tue, 22 Aug 2017 21:45:08 +0300 (EEST) Received: by mail-pg0-f66.google.com with SMTP id m133so8664624pga.2 for ; Tue, 22 Aug 2017 11:45:16 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=from:to:subject:date:message-id:mime-version :content-transfer-encoding; bh=rmZ7s62iDcQrU+SIINXviqi+t47NS6Z6G1YcqB79svA=; b=bmzTaIaJTFt1rAZ6V2MBi95r3x/jMHRiF+1cKVJ6WgXXqW9nvwqCW8qcO03GV1AfdH +mAxWBBYvdvE1lrlMHxMJ0poC4kpQwrQ0VeNOqLdBY+ZRJiI8z1q2/ww67dpU2BH3ppm 6U3dzjqSUDDGpuUfTg0Piw9ey66nMcNTdFmd7MlymU8SxdHMRf55GMTiz5pEjnRhCin2 utom8EuXyM6exU3Bq92ODD7MkYdUTZVSCUwo2DgOUkp+MKVlw6FvDqOZgR4lN1KsOwrG HAbC1WaXzbCeO1/NFpvvlJXE3QySpFPyvHQ7fWFaTL3MhUEh7ScWh2hgIEUo0cFzsTkX u84Q== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:from:to:subject:date:message-id:mime-version :content-transfer-encoding; bh=rmZ7s62iDcQrU+SIINXviqi+t47NS6Z6G1YcqB79svA=; b=Jukc1vWMyXmM6eb15ug3gf0ts5paJ40vgY7z1DXUaquggCi2h6mQgl0vMOSc2VwXSW pt4pDKti1cKMSAJRHDadu3uBDQDiBn4YTaCwaCl7GVvT3caDiRKLc5a2ppRg4ja/qxGJ Rftv090R1pVs7GoBk+OUC0/vPWrxOF7tDa4O2uNco54EqpuEDA+krDFxQpWLfWV3RIbc gE0WexVOkSXcs4SdfAq2Ci5T5yG4i+PJEyEytWLN0gdcb2ZWaF7xSjfT/YNAVfllmS/L 1mjj5k9Cy9WgGu2sOYxLDaQh8Q707jX4/7/pJbmi7xn01bRYwdqw8mCzaWXcjFDllEfT 38Ww== X-Gm-Message-State: AHYfb5glxzH3wA+kbYlyNRkx2BrMndlKNVFuCWcAGQP3mPuDgQVcrFHJ ZbMfwLan/mUukCvqNmM= X-Received: by 10.84.233.139 with SMTP id l11mr116895plk.298.1503427109614; Tue, 22 Aug 2017 11:38:29 -0700 (PDT) Received: from localhost.localdomain ([114.124.205.114]) by smtp.gmail.com with ESMTPSA id d184sm3662697pfc.110.2017.08.22.11.38.27 for (version=TLS1_2 cipher=ECDHE-RSA-AES128-GCM-SHA256 bits=128/128); Tue, 22 Aug 2017 11:38:28 -0700 (PDT) From: Muhammad Faiz To: ffmpeg-devel@ffmpeg.org Date: Wed, 23 Aug 2017 01:38:15 +0700 Message-Id: <20170822183816.4724-1-mfcc64@gmail.com> X-Mailer: git-send-email 2.13.2 MIME-Version: 1.0 Subject: [FFmpeg-devel] [PATCH] avfilter/af_firequalizer: add min_phase option X-BeenThere: ffmpeg-devel@ffmpeg.org X-Mailman-Version: 2.1.20 Precedence: list List-Id: FFmpeg development discussions and patches List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Reply-To: FFmpeg development discussions and patches Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" Signed-off-by: Muhammad Faiz --- doc/filters.texi | 3 + libavfilter/af_firequalizer.c | 147 +++++++++++++++++++++++++++++++++++++++++- 2 files changed, 147 insertions(+), 3 deletions(-) diff --git a/doc/filters.texi b/doc/filters.texi index 3b5a38fc9f..bd88665c89 100644 --- a/doc/filters.texi +++ b/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 diff --git a/libavfilter/af_firequalizer.c b/libavfilter/af_firequalizer.c index 7741057a65..8348e7558b 100644 --- a/libavfilter/af_firequalizer.c +++ b/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);