git.sesse.net Git - ffmpeg/blob - libavfilter/af_afftdn.c

   1 /*
   2  * Copyright (c) 2018 The FFmpeg Project
   3  *
   4  * This file is part of FFmpeg.
   5  *
   6  * FFmpeg is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * FFmpeg is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FFmpeg; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19  */
  20
  21 #include <float.h>
  22
  23 #include "libavutil/audio_fifo.h"
  24 #include "libavutil/avstring.h"
  25 #include "libavutil/channel_layout.h"
  26 #include "libavutil/opt.h"
  27 #include "libavcodec/avfft.h"
  28 #include "avfilter.h"
  29 #include "audio.h"
  30 #include "formats.h"
  31
  32 #define C       (M_LN10 * 0.1)
  33 #define RATIO    0.98
  34 #define RRATIO  (1.0 - RATIO)
  35
  36 enum OutModes {
  37     IN_MODE,
  38     OUT_MODE,
  39     NOISE_MODE,
  40     NB_MODES
  41 };
  42
  43 enum NoiseType {
  44     WHITE_NOISE,
  45     VINYL_NOISE,
  46     SHELLAC_NOISE,
  47     CUSTOM_NOISE,
  48     NB_NOISE
  49 };
  50
  51 typedef struct DeNoiseChannel {
  52     int         band_noise[15];
  53     double      noise_band_auto_var[15];
  54     double      noise_band_sample[15];
  55     double     *amt;
  56     double     *band_amt;
  57     double     *band_excit;
  58     double     *gain;
  59     double     *prior;
  60     double     *prior_band_excit;
  61     double     *clean_data;
  62     double     *noisy_data;
  63     double     *out_samples;
  64     double     *spread_function;
  65     double     *abs_var;
  66     double     *rel_var;
  67     double     *min_abs_var;
  68     FFTComplex *fft_data;
  69     FFTContext *fft, *ifft;
  70
  71     double      noise_band_norm[15];
  72     double      noise_band_avr[15];
  73     double      noise_band_avi[15];
  74     double      noise_band_var[15];
  75
  76     double      sfm_threshold;
  77     double      sfm_alpha;
  78     double      sfm_results[3];
  79     int         sfm_fail_flags[512];
  80     int         sfm_fail_total;
  81 } DeNoiseChannel;
  82
  83 typedef struct AudioFFTDeNoiseContext {
  84     const AVClass *class;
  85
  86     float   noise_reduction;
  87     float   noise_floor;
  88     int     noise_type;
  89     char   *band_noise_str;
  90     float   residual_floor;
  91     int     track_noise;
  92     int     track_residual;
  93     int     output_mode;
  94
  95     float   last_residual_floor;
  96     float   last_noise_floor;
  97     float   last_noise_reduction;
  98     float   last_noise_balance;
  99     int64_t block_count;
 100
 101     int64_t pts;
 102     int     channels;
 103     int     sample_noise;
 104     int     sample_noise_start;
 105     int     sample_noise_end;
 106     float   sample_rate;
 107     int     buffer_length;
 108     int     fft_length;
 109     int     fft_length2;
 110     int     bin_count;
 111     int     window_length;
 112     int     sample_advance;
 113     int     number_of_bands;
 114
 115     int     band_centre[15];
 116
 117     int    *bin2band;
 118     double *window;
 119     double *band_alpha;
 120     double *band_beta;
 121
 122     DeNoiseChannel *dnch;
 123
 124     double  max_gain;
 125     double  max_var;
 126     double  gain_scale;
 127     double  window_weight;
 128     double  floor;
 129     double  sample_floor;
 130     double  auto_floor;
 131
 132     int     noise_band_edge[17];
 133     int     noise_band_count;
 134     double  matrix_a[25];
 135     double  vector_b[5];
 136     double  matrix_b[75];
 137     double  matrix_c[75];
 138
 139     AVAudioFifo *fifo;
 140 } AudioFFTDeNoiseContext;
 141
 142 #define OFFSET(x) offsetof(AudioFFTDeNoiseContext, x)
 143 #define A AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
 144
 145 static const AVOption afftdn_options[] = {
 146     { "nr", "set the noise reduction",    OFFSET(noise_reduction), AV_OPT_TYPE_FLOAT,  {.dbl = 12},          .01, 97, A },
 147     { "nf", "set the noise floor",        OFFSET(noise_floor),     AV_OPT_TYPE_FLOAT,  {.dbl =-50},          -80,-20, A },
 148     { "nt", "set the noise type",         OFFSET(noise_type),      AV_OPT_TYPE_INT,    {.i64 = WHITE_NOISE}, WHITE_NOISE, NB_NOISE-1, A, "type" },
 149     {  "w", "white noise",                0,                       AV_OPT_TYPE_CONST,  {.i64 = WHITE_NOISE},   0,  0, A, "type" },
 150     {  "v", "vinyl noise",                0,                       AV_OPT_TYPE_CONST,  {.i64 = VINYL_NOISE},   0,  0, A, "type" },
 151     {  "s", "shellac noise",              0,                       AV_OPT_TYPE_CONST,  {.i64 = SHELLAC_NOISE}, 0,  0, A, "type" },
 152     {  "c", "custom noise",               0,                       AV_OPT_TYPE_CONST,  {.i64 = CUSTOM_NOISE},  0,  0, A, "type" },
 153     { "bn", "set the custom bands noise", OFFSET(band_noise_str),  AV_OPT_TYPE_STRING, {.str = 0},             0,  0, A },
 154     { "rf", "set the residual floor",     OFFSET(residual_floor),  AV_OPT_TYPE_FLOAT,  {.dbl =-38},          -80,-20, A },
 155     { "tn", "track noise",                OFFSET(track_noise),     AV_OPT_TYPE_BOOL,   {.i64 =  0},            0,  1, A },
 156     { "tr", "track residual",             OFFSET(track_residual),  AV_OPT_TYPE_BOOL,   {.i64 =  0},            0,  1, A },
 157     { "om", "set output mode",            OFFSET(output_mode),     AV_OPT_TYPE_INT,    {.i64 = OUT_MODE},      0,  NB_MODES-1, A, "mode" },
 158     {  "i", "input",                      0,                       AV_OPT_TYPE_CONST,  {.i64 = IN_MODE},       0,  0, A, "mode" },
 159     {  "o", "output",                     0,                       AV_OPT_TYPE_CONST,  {.i64 = OUT_MODE},      0,  0, A, "mode" },
 160     {  "n", "noise",                      0,                       AV_OPT_TYPE_CONST,  {.i64 = NOISE_MODE},    0,  0, A, "mode" },
 161     { NULL }
 162 };
 163
 164 AVFILTER_DEFINE_CLASS(afftdn);
 165
 166 static int get_band_noise(AudioFFTDeNoiseContext *s,
 167                           int band, double a,
 168                           double b, double c)
 169 {
 170     double d1, d2, d3;
 171
 172     d1 = a / s->band_centre[band];
 173     d1 = 10.0 * log(1.0 + d1 * d1) / M_LN10;
 174     d2 = b / s->band_centre[band];
 175     d2 = 10.0 * log(1.0 + d2 * d2) / M_LN10;
 176     d3 = s->band_centre[band] / c;
 177     d3 = 10.0 * log(1.0 + d3 * d3) / M_LN10;
 178
 179     return lrint(-d1 + d2 - d3);
 180 }
 181
 182 static void factor(double *array, int size)
 183 {
 184     for (int i = 0; i < size - 1; i++) {
 185         for (int j = i + 1; j < size; j++) {
 186             double d = array[j + i * size] / array[i + i * size];
 187
 188             array[j + i * size] = d;
 189             for (int k = i + 1; k < size; k++) {
 190                 array[j + k * size] -= d * array[i + k * size];
 191             }
 192         }
 193     }
 194 }
 195
 196 static void solve(double *matrix, double *vector, int size)
 197 {
 198     for (int i = 0; i < size - 1; i++) {
 199         for (int j = i + 1; j < size; j++) {
 200             double d = matrix[j + i * size];
 201             vector[j] -= d * vector[i];
 202         }
 203     }
 204
 205     vector[size - 1] /= matrix[size * size - 1];
 206
 207     for (int i = size - 2; i >= 0; i--) {
 208         double d = vector[i];
 209         for (int j = i + 1; j < size; j++)
 210             d -= matrix[i + j * size] * vector[j];
 211         vector[i] = d / matrix[i + i * size];
 212     }
 213 }
 214
 215 static int process_get_band_noise(AudioFFTDeNoiseContext *s,
 216                                   DeNoiseChannel *dnch,
 217                                   int band)
 218 {
 219     double product, sum, f;
 220     int i = 0;
 221
 222     if (band < 15)
 223         return dnch->band_noise[band];
 224
 225     for (int j = 0; j < 5; j++) {
 226         sum = 0.0;
 227         for (int k = 0; k < 15; k++)
 228             sum += s->matrix_b[i++] * dnch->band_noise[k];
 229         s->vector_b[j] = sum;
 230     }
 231
 232     solve(s->matrix_a, s->vector_b, 5);
 233     f = (0.5 * s->sample_rate) / s->band_centre[14];
 234     f = 15.0 + log(f / 1.5) / log(1.5);
 235     sum = 0.0;
 236     product = 1.0;
 237     for (int j = 0; j < 5; j++) {
 238         sum += product * s->vector_b[j];
 239         product *= f;
 240     }
 241
 242     return lrint(sum);
 243 }
 244
 245 static void calculate_sfm(AudioFFTDeNoiseContext *s,
 246                           DeNoiseChannel *dnch,
 247                           int start, int end)
 248 {
 249     double d1 = 0.0, d2 = 1.0;
 250     int i = 0, j = 0;
 251
 252     for (int k = start; k < end; k++) {
 253         if (dnch->noisy_data[k] > s->sample_floor) {
 254             j++;
 255             d1 += dnch->noisy_data[k];
 256             d2 *= dnch->noisy_data[k];
 257             if (d2 > 1.0E100) {
 258                 d2 *= 1.0E-100;
 259                 i++;
 260             } else if (d2 < 1.0E-100) {
 261                 d2 *= 1.0E100;
 262                 i--;
 263             }
 264         }
 265     }
 266     if (j > 1) {
 267         d1 /= j;
 268         dnch->sfm_results[0] = d1;
 269         d2 = log(d2) + 230.2585 * i;
 270         d2 /= j;
 271         d1 = log(d1);
 272         dnch->sfm_results[1] = d1;
 273         dnch->sfm_results[2] = d1 - d2;
 274     } else {
 275         dnch->sfm_results[0] = s->auto_floor;
 276         dnch->sfm_results[1] = dnch->sfm_threshold;
 277         dnch->sfm_results[2] = dnch->sfm_threshold;
 278     }
 279 }
 280
 281 static double limit_gain(double a, double b)
 282 {
 283     if (a > 1.0)
 284         return (b * a - 1.0) / (b + a - 2.0);
 285     if (a < 1.0)
 286         return (b * a - 2.0 * a + 1.0) / (b - a);
 287     return 1.0;
 288 }
 289
 290 static void process_frame(AudioFFTDeNoiseContext *s, DeNoiseChannel *dnch,
 291                           FFTComplex *fft_data,
 292                           double *prior, double *prior_band_excit, int track_noise)
 293 {
 294     double d1, d2, d3, gain;
 295     int n, i1;
 296
 297     d1 = fft_data[0].re * fft_data[0].re;
 298     dnch->noisy_data[0] = d1;
 299     d2 = d1 / dnch->abs_var[0];
 300     d3 = RATIO * prior[0] + RRATIO * fmax(d2 - 1.0, 0.0);
 301     gain = d3 / (1.0 + d3);
 302     gain *= (gain + M_PI_4 / fmax(d2, 1.0E-6));
 303     prior[0] = (d2 * gain);
 304     dnch->clean_data[0] = (d1 * gain);
 305     gain = sqrt(gain);
 306     dnch->gain[0] = gain;
 307     n = 0;
 308     for (int i = 1; i < s->fft_length2; i++) {
 309         d1 = fft_data[i].re * fft_data[i].re + fft_data[i].im * fft_data[i].im;
 310         if (d1 > s->sample_floor)
 311             n = i;
 312
 313         dnch->noisy_data[i] = d1;
 314         d2 = d1 / dnch->abs_var[i];
 315         d3 = RATIO * prior[i] + RRATIO * fmax(d2 - 1.0, 0.0);
 316         gain = d3 / (1.0 + d3);
 317         gain *= (gain + M_PI_4 / fmax(d2, 1.0E-6));
 318         prior[i] = d2 * gain;
 319         dnch->clean_data[i] = d1 * gain;
 320         gain = sqrt(gain);
 321         dnch->gain[i] = gain;
 322     }
 323     d1 = fft_data[0].im * fft_data[0].im;
 324     if (d1 > s->sample_floor)
 325         n = s->fft_length2;
 326
 327     dnch->noisy_data[s->fft_length2] = d1;
 328     d2 = d1 / dnch->abs_var[s->fft_length2];
 329     d3 = RATIO * prior[s->fft_length2] + RRATIO * fmax(d2 - 1.0, 0.0);
 330     gain = d3 / (1.0 + d3);
 331     gain *= gain + M_PI_4 / fmax(d2, 1.0E-6);
 332     prior[s->fft_length2] = d2 * gain;
 333     dnch->clean_data[s->fft_length2] = d1 * gain;
 334     gain = sqrt(gain);
 335     dnch->gain[s->fft_length2] = gain;
 336     if (n > s->fft_length2 - 2) {
 337         n = s->bin_count;
 338         i1 = s->noise_band_count;
 339     } else {
 340         i1 = 0;
 341         for (int i = 0; i <= s->noise_band_count; i++) {
 342             if (n > 1.1 * s->noise_band_edge[i]) {
 343                 i1 = i;
 344             }
 345         }
 346     }
 347
 348     if (track_noise && (i1 > s->noise_band_count / 2)) {
 349         int j = FFMIN(n, s->noise_band_edge[i1]);
 350         int m = 3, k;
 351
 352         for (k = i1 - 1; k >= 0; k--) {
 353             int i = s->noise_band_edge[k];
 354             calculate_sfm(s, dnch, i, j);
 355             dnch->noise_band_sample[k] = dnch->sfm_results[0];
 356             if (dnch->sfm_results[2] + 0.013 * m * fmax(0.0, dnch->sfm_results[1] - 20.53) >= dnch->sfm_threshold) {
 357                 break;
 358             }
 359             j = i;
 360             m++;
 361         }
 362
 363         if (k < i1 - 1) {
 364             double sum = 0.0, min, max;
 365             int i;
 366
 367             for (i = i1 - 1; i > k; i--) {
 368                 min = log(dnch->noise_band_sample[i] / dnch->noise_band_auto_var[i]);
 369                 sum += min;
 370             }
 371
 372             i = i1 - k - 1;
 373             if (i < 5) {
 374                 min = 3.0E-4 * i * i;
 375             } else {
 376                 min = 3.0E-4 * (8 * i - 16);
 377             }
 378             if (i < 3) {
 379                 max = 2.0E-4 * i * i;
 380             } else {
 381                 max = 2.0E-4 * (4 * i - 4);
 382             }
 383
 384             if (s->track_residual) {
 385                 if (s->last_noise_floor > s->last_residual_floor + 9) {
 386                     min *= 0.5;
 387                     max *= 0.75;
 388                 } else if (s->last_noise_floor > s->last_residual_floor + 6) {
 389                     min *= 0.4;
 390                     max *= 1.0;
 391                 } else if (s->last_noise_floor > s->last_residual_floor + 4) {
 392                     min *= 0.3;
 393                     max *= 1.3;
 394                 } else if (s->last_noise_floor > s->last_residual_floor + 2) {
 395                     min *= 0.2;
 396                     max *= 1.6;
 397                 } else if (s->last_noise_floor > s->last_residual_floor) {
 398                     min *= 0.1;
 399                     max *= 2.0;
 400                 } else {
 401                     min = 0.0;
 402                     max *= 2.5;
 403                 }
 404             }
 405
 406             sum = av_clipd(sum, -min, max);
 407             sum = exp(sum);
 408             for (int i = 0; i < 15; i++)
 409                 dnch->noise_band_auto_var[i] *= sum;
 410         } else if (dnch->sfm_results[2] >= dnch->sfm_threshold) {
 411             dnch->sfm_fail_flags[s->block_count & 0x1FF] = 1;
 412             dnch->sfm_fail_total += 1;
 413         }
 414     }
 415
 416     for (int i = 0; i < s->number_of_bands; i++) {
 417         dnch->band_excit[i] = 0.0;
 418         dnch->band_amt[i] = 0.0;
 419     }
 420
 421     for (int i = 0; i < s->bin_count; i++) {
 422         dnch->band_excit[s->bin2band[i]] += dnch->clean_data[i];
 423     }
 424
 425     for (int i = 0; i < s->number_of_bands; i++) {
 426         dnch->band_excit[i] = fmax(dnch->band_excit[i],
 427                                 s->band_alpha[i] * dnch->band_excit[i] +
 428                                 s->band_beta[i] * prior_band_excit[i]);
 429         prior_band_excit[i] = dnch->band_excit[i];
 430     }
 431
 432     for (int j = 0, i = 0; j < s->number_of_bands; j++) {
 433         for (int k = 0; k < s->number_of_bands; k++) {
 434             dnch->band_amt[j] += dnch->spread_function[i++] * dnch->band_excit[k];
 435         }
 436     }
 437
 438     for (int i = 0; i < s->bin_count; i++)
 439         dnch->amt[i] = dnch->band_amt[s->bin2band[i]];
 440
 441     if (dnch->amt[0] > dnch->abs_var[0]) {
 442         dnch->gain[0] = 1.0;
 443     } else if (dnch->amt[0] > dnch->min_abs_var[0]) {
 444         double limit = sqrt(dnch->abs_var[0] / dnch->amt[0]);
 445         dnch->gain[0] = limit_gain(dnch->gain[0], limit);
 446     } else {
 447         dnch->gain[0] = limit_gain(dnch->gain[0], s->max_gain);
 448     }
 449     if (dnch->amt[s->fft_length2] > dnch->abs_var[s->fft_length2]) {
 450         dnch->gain[s->fft_length2] = 1.0;
 451     } else if (dnch->amt[s->fft_length2] > dnch->min_abs_var[s->fft_length2]) {
 452         double limit = sqrt(dnch->abs_var[s->fft_length2] / dnch->amt[s->fft_length2]);
 453         dnch->gain[s->fft_length2] = limit_gain(dnch->gain[s->fft_length2], limit);
 454     } else {
 455         dnch->gain[s->fft_length2] = limit_gain(dnch->gain[s->fft_length2], s->max_gain);
 456     }
 457
 458     for (int i = 1; i < s->fft_length2; i++) {
 459         if (dnch->amt[i] > dnch->abs_var[i]) {
 460             dnch->gain[i] = 1.0;
 461         } else if (dnch->amt[i] > dnch->min_abs_var[i]) {
 462             double limit = sqrt(dnch->abs_var[i] / dnch->amt[i]);
 463             dnch->gain[i] = limit_gain(dnch->gain[i], limit);
 464         } else {
 465             dnch->gain[i] = limit_gain(dnch->gain[i], s->max_gain);
 466         }
 467     }
 468
 469     gain = dnch->gain[0];
 470     dnch->clean_data[0] = (gain * gain * dnch->noisy_data[0]);
 471     fft_data[0].re *= gain;
 472     gain = dnch->gain[s->fft_length2];
 473     dnch->clean_data[s->fft_length2] = (gain * gain * dnch->noisy_data[s->fft_length2]);
 474     fft_data[0].im *= gain;
 475     for (int i = 1; i < s->fft_length2; i++) {
 476         gain = dnch->gain[i];
 477         dnch->clean_data[i] = (gain * gain * dnch->noisy_data[i]);
 478         fft_data[i].re *= gain;
 479         fft_data[i].im *= gain;
 480     }
 481 }
 482
 483 static double freq2bark(double x)
 484 {
 485     double d = x / 7500.0;
 486
 487     return 13.0 * atan(7.6E-4 * x) + 3.5 * atan(d * d);
 488 }
 489
 490 static int get_band_centre(AudioFFTDeNoiseContext *s, int band)
 491 {
 492     if (band == -1)
 493         return lrint(s->band_centre[0] / 1.5);
 494
 495     return s->band_centre[band];
 496 }
 497
 498 static int get_band_edge(AudioFFTDeNoiseContext *s, int band)
 499 {
 500     int i;
 501
 502     if (band == 15) {
 503         i = lrint(s->band_centre[14] * 1.224745);
 504     } else {
 505         i = lrint(s->band_centre[band] / 1.224745);
 506     }
 507
 508     return FFMIN(i, s->sample_rate / 2);
 509 }
 510
 511 static void set_band_parameters(AudioFFTDeNoiseContext *s,
 512                                 DeNoiseChannel *dnch)
 513 {
 514     double band_noise, d2, d3, d4, d5;
 515     int i = 0, j = 0, k = 0;
 516
 517     d5 = 0.0;
 518     band_noise = process_get_band_noise(s, dnch, 0);
 519     for (int m = j; m <= s->fft_length2; m++) {
 520         if (m == j) {
 521             i = j;
 522             d5 = band_noise;
 523             if (k == 15) {
 524                 j = s->bin_count;
 525             } else {
 526                 j = s->fft_length * get_band_centre(s, k) / s->sample_rate;
 527             }
 528             d2 = j - i;
 529             band_noise = process_get_band_noise(s, dnch, k);
 530             k++;
 531         }
 532         d3 = (j - m) / d2;
 533         d4 = (m - i) / d2;
 534         dnch->rel_var[m] = exp((d5 * d3 + band_noise * d4) * C);
 535     }
 536     dnch->rel_var[s->fft_length2] = exp(band_noise * C);
 537
 538     for (i = 0; i < 15; i++)
 539         dnch->noise_band_auto_var[i] = s->max_var * exp((process_get_band_noise(s, dnch, i) - 2.0) * C);
 540
 541     for (i = 0; i <= s->fft_length2; i++) {
 542         dnch->abs_var[i] = fmax(s->max_var * dnch->rel_var[i], 1.0);
 543         dnch->min_abs_var[i] = s->gain_scale * dnch->abs_var[i];
 544     }
 545 }
 546
 547 static void read_custom_noise(AudioFFTDeNoiseContext *s, int ch)
 548 {
 549     DeNoiseChannel *dnch = &s->dnch[ch];
 550     char *p, *arg, *saveptr = NULL;
 551     int i, ret, band_noise[15] = { 0 };
 552
 553     if (!s->band_noise_str)
 554         return;
 555
 556     p = av_strdup(s->band_noise_str);
 557     if (!p)
 558         return;
 559
 560     for (i = 0; i < 15; i++) {
 561         if (!(arg = av_strtok(p, "| ", &saveptr)))
 562             break;
 563
 564         p = NULL;
 565
 566         ret = sscanf(arg, "%d", &band_noise[i]);
 567         if (ret != 1) {
 568             av_log(s, AV_LOG_ERROR, "Custom band noise must be integer.\n");
 569             break;
 570         }
 571
 572         band_noise[i] = av_clip(band_noise[i], -24, 24);
 573     }
 574
 575     av_free(p);
 576     memcpy(dnch->band_noise, band_noise, sizeof(band_noise));
 577 }
 578
 579 static void set_parameters(AudioFFTDeNoiseContext *s)
 580 {
 581     if (s->last_noise_floor != s->noise_floor)
 582         s->last_noise_floor = s->noise_floor;
 583
 584     if (s->track_residual)
 585         s->last_noise_floor = fmaxf(s->last_noise_floor, s->residual_floor);
 586
 587     s->max_var = s->floor * exp((100.0 + s->last_noise_floor) * C);
 588
 589     if (s->track_residual) {
 590         s->last_residual_floor = s->residual_floor;
 591         s->last_noise_reduction = fmax(s->last_noise_floor - s->last_residual_floor, 0);
 592         s->max_gain = exp(s->last_noise_reduction * (0.5 * C));
 593     } else if (s->noise_reduction != s->last_noise_reduction) {
 594         s->last_noise_reduction = s->noise_reduction;
 595         s->last_residual_floor = av_clipf(s->last_noise_floor - s->last_noise_reduction, -80, -20);
 596         s->max_gain = exp(s->last_noise_reduction * (0.5 * C));
 597     }
 598
 599     s->gain_scale = 1.0 / (s->max_gain * s->max_gain);
 600
 601     for (int ch = 0; ch < s->channels; ch++) {
 602         DeNoiseChannel *dnch = &s->dnch[ch];
 603
 604         set_band_parameters(s, dnch);
 605     }
 606 }
 607
 608 static int config_input(AVFilterLink *inlink)
 609 {
 610     AVFilterContext *ctx = inlink->dst;
 611     AudioFFTDeNoiseContext *s = ctx->priv;
 612     double wscale, sar, sum, sdiv;
 613     int i, j, k, m, n;
 614
 615     s->dnch = av_calloc(inlink->channels, sizeof(*s->dnch));
 616     if (!s->dnch)
 617         return AVERROR(ENOMEM);
 618
 619     s->pts = AV_NOPTS_VALUE;
 620     s->channels = inlink->channels;
 621     s->sample_rate = inlink->sample_rate;
 622     s->sample_advance = s->sample_rate / 80;
 623     s->window_length = 3 * s->sample_advance;
 624     s->fft_length2 = 1 << (32 - ff_clz(s->window_length));
 625     s->fft_length = s->fft_length2 * 2;
 626     s->buffer_length = s->fft_length * 2;
 627     s->bin_count = s->fft_length2 + 1;
 628
 629     s->band_centre[0] = 80;
 630     for (i = 1; i < 15; i++) {
 631         s->band_centre[i] = lrint(1.5 * s->band_centre[i - 1] + 5.0);
 632         if (s->band_centre[i] < 1000) {
 633             s->band_centre[i] = 10 * (s->band_centre[i] / 10);
 634         } else if (s->band_centre[i] < 5000) {
 635             s->band_centre[i] = 50 * ((s->band_centre[i] + 20) / 50);
 636         } else if (s->band_centre[i] < 15000) {
 637             s->band_centre[i] = 100 * ((s->band_centre[i] + 45) / 100);
 638         } else {
 639             s->band_centre[i] = 1000 * ((s->band_centre[i] + 495) / 1000);
 640         }
 641     }
 642
 643     for (j = 0; j < 5; j++) {
 644         for (k = 0; k < 5; k++) {
 645             s->matrix_a[j + k * 5] = 0.0;
 646             for (m = 0; m < 15; m++)
 647                 s->matrix_a[j + k * 5] += pow(m, j + k);
 648         }
 649     }
 650
 651     factor(s->matrix_a, 5);
 652
 653     i = 0;
 654     for (j = 0; j < 5; j++)
 655         for (k = 0; k < 15; k++)
 656             s->matrix_b[i++] = pow(k, j);
 657
 658     i = 0;
 659     for (j = 0; j < 15; j++)
 660         for (k = 0; k < 5; k++)
 661             s->matrix_c[i++] = pow(j, k);
 662
 663     s->window = av_calloc(s->window_length, sizeof(*s->window));
 664     s->bin2band = av_calloc(s->bin_count, sizeof(*s->bin2band));
 665     if (!s->window || !s->bin2band)
 666         return AVERROR(ENOMEM);
 667
 668     sdiv = s->sample_rate / 17640.0;
 669     for (i = 0; i <= s->fft_length2; i++)
 670         s->bin2band[i] = lrint(sdiv * freq2bark((0.5 * i * s->sample_rate) / s->fft_length2));
 671
 672     s->number_of_bands = s->bin2band[s->fft_length2] + 1;
 673
 674     s->band_alpha = av_calloc(s->number_of_bands, sizeof(*s->band_alpha));
 675     s->band_beta = av_calloc(s->number_of_bands, sizeof(*s->band_beta));
 676     if (!s->band_alpha || !s->band_beta)
 677         return AVERROR(ENOMEM);
 678
 679     for (int ch = 0; ch < inlink->channels; ch++) {
 680         DeNoiseChannel *dnch = &s->dnch[ch];
 681
 682         switch (s->noise_type) {
 683         case WHITE_NOISE:
 684             for (i = 0; i < 15; i++)
 685                 dnch->band_noise[i] = 0;
 686             break;
 687         case VINYL_NOISE:
 688             for (i = 0; i < 15; i++)
 689                 dnch->band_noise[i] = get_band_noise(s, i, 50.0, 500.5, 2125.0) + FFMAX(i - 7, 0);
 690             break;
 691         case SHELLAC_NOISE:
 692             for (i = 0; i < 15; i++)
 693                 dnch->band_noise[i] = get_band_noise(s, i, 1.0, 500.0, 1.0E10) + FFMAX(i - 12, -5);
 694         case CUSTOM_NOISE:
 695             read_custom_noise(s, ch);
 696             break;
 697         default:
 698             return AVERROR_BUG;
 699         }
 700
 701
 702         dnch->sfm_threshold = 0.8;
 703         dnch->sfm_alpha = 0.05;
 704         for (i = 0; i < 512; i++)
 705             dnch->sfm_fail_flags[i] = 0;
 706
 707         dnch->sfm_fail_total = 0;
 708         j = FFMAX((int)(10.0 * (1.3 - dnch->sfm_threshold)), 1);
 709
 710         for (i = 0; i < 512; i += j) {
 711             dnch->sfm_fail_flags[i] = 1;
 712             dnch->sfm_fail_total += 1;
 713         }
 714
 715         dnch->amt = av_calloc(s->bin_count, sizeof(*dnch->amt));
 716         dnch->band_amt = av_calloc(s->number_of_bands, sizeof(*dnch->band_amt));
 717         dnch->band_excit = av_calloc(s->number_of_bands, sizeof(*dnch->band_excit));
 718         dnch->gain = av_calloc(s->bin_count, sizeof(*dnch->gain));
 719         dnch->prior = av_calloc(s->bin_count, sizeof(*dnch->prior));
 720         dnch->prior_band_excit = av_calloc(s->number_of_bands, sizeof(*dnch->prior_band_excit));
 721         dnch->clean_data = av_calloc(s->bin_count, sizeof(*dnch->clean_data));
 722         dnch->noisy_data = av_calloc(s->bin_count, sizeof(*dnch->noisy_data));
 723         dnch->out_samples = av_calloc(s->buffer_length, sizeof(*dnch->out_samples));
 724         dnch->abs_var = av_calloc(s->bin_count, sizeof(*dnch->abs_var));
 725         dnch->rel_var = av_calloc(s->bin_count, sizeof(*dnch->rel_var));
 726         dnch->min_abs_var = av_calloc(s->bin_count, sizeof(*dnch->min_abs_var));
 727         dnch->fft_data = av_calloc(s->fft_length2 + 1, sizeof(*dnch->fft_data));
 728         dnch->fft  = av_fft_init(av_log2(s->fft_length2), 0);
 729         dnch->ifft = av_fft_init(av_log2(s->fft_length2), 1);
 730         dnch->spread_function = av_calloc(s->number_of_bands * s->number_of_bands,
 731                                           sizeof(*dnch->spread_function));
 732
 733         if (!dnch->amt ||
 734             !dnch->band_amt ||
 735             !dnch->band_excit ||
 736             !dnch->gain ||
 737             !dnch->prior ||
 738             !dnch->prior_band_excit ||
 739             !dnch->clean_data ||
 740             !dnch->noisy_data ||
 741             !dnch->out_samples ||
 742             !dnch->fft_data ||
 743             !dnch->abs_var ||
 744             !dnch->rel_var ||
 745             !dnch->min_abs_var ||
 746             !dnch->spread_function ||
 747             !dnch->fft ||
 748             !dnch->ifft)
 749             return AVERROR(ENOMEM);
 750     }
 751
 752     for (int ch = 0; ch < inlink->channels; ch++) {
 753         DeNoiseChannel *dnch = &s->dnch[ch];
 754         double *prior_band_excit = dnch->prior_band_excit;
 755         double *prior = dnch->prior;
 756         double min, max;
 757         double p1, p2;
 758
 759         p1 = pow(0.1, 2.5 / sdiv);
 760         p2 = pow(0.1, 1.0 / sdiv);
 761         j = 0;
 762         for (m = 0; m < s->number_of_bands; m++) {
 763             for (n = 0; n < s->number_of_bands; n++) {
 764                 if (n < m) {
 765                     dnch->spread_function[j++] = pow(p2, m - n);
 766                 } else if (n > m) {
 767                     dnch->spread_function[j++] = pow(p1, n - m);
 768                 } else {
 769                     dnch->spread_function[j++] = 1.0;
 770                 }
 771             }
 772         }
 773
 774         for (m = 0; m < s->number_of_bands; m++) {
 775             dnch->band_excit[m] = 0.0;
 776             prior_band_excit[m] = 0.0;
 777         }
 778
 779         for (m = 0; m <= s->fft_length2; m++)
 780             dnch->band_excit[s->bin2band[m]] += 1.0;
 781
 782         j = 0;
 783         for (m = 0; m < s->number_of_bands; m++) {
 784             for (n = 0; n < s->number_of_bands; n++)
 785                 prior_band_excit[m] += dnch->spread_function[j++] * dnch->band_excit[n];
 786         }
 787
 788         min = pow(0.1, 2.5);
 789         max = pow(0.1, 1.0);
 790         for (int i = 0; i < s->number_of_bands; i++) {
 791             if (i < lrint(12.0 * sdiv)) {
 792                 dnch->band_excit[i] = pow(0.1, 1.45 + 0.1 * i / sdiv);
 793             } else {
 794                 dnch->band_excit[i] = pow(0.1, 2.5 - 0.2 * (i / sdiv - 14.0));
 795             }
 796             dnch->band_excit[i] = av_clipd(dnch->band_excit[i], min, max);
 797         }
 798
 799         for (int i = 0; i <= s->fft_length2; i++)
 800             prior[i] = RRATIO;
 801         for (int i = 0; i < s->buffer_length; i++)
 802             dnch->out_samples[i] = 0;
 803
 804         j = 0;
 805         for (int i = 0; i < s->number_of_bands; i++)
 806             for (int k = 0; k < s->number_of_bands; k++)
 807                 dnch->spread_function[j++] *= dnch->band_excit[i] / prior_band_excit[i];
 808     }
 809
 810     j = 0;
 811     sar = s->sample_advance / s->sample_rate;
 812     for (int i = 0; i <= s->fft_length2; i++) {
 813         if ((i == s->fft_length2) || (s->bin2band[i] > j)) {
 814             double d6 = (i - 1) * s->sample_rate / s->fft_length;
 815             double d7 = fmin(0.008 + 2.2 / d6, 0.03);
 816             s->band_alpha[j] = exp(-sar / d7);
 817             s->band_beta[j] = 1.0 - s->band_alpha[j];
 818             j = s->bin2band[i];
 819         }
 820     }
 821
 822     wscale = sqrt(16.0 / (9.0 * s->fft_length));
 823     sum = 0.0;
 824     for (int i = 0; i < s->window_length; i++) {
 825         double d10 = sin(i * M_PI / s->window_length);
 826         d10 *= wscale * d10;
 827         s->window[i] = d10;
 828         sum += d10 * d10;
 829     }
 830
 831     s->window_weight = 0.5 * sum;
 832     s->floor = (1LL << 48) * exp(-23.025558369790467) * s->window_weight;
 833     s->sample_floor = s->floor * exp(4.144600506562284);
 834     s->auto_floor = s->floor * exp(6.907667510937141);
 835
 836     set_parameters(s);
 837
 838     s->noise_band_edge[0] = FFMIN(s->fft_length2, s->fft_length * get_band_edge(s, 0) / s->sample_rate);
 839     i = 0;
 840     for (int j = 1; j < 16; j++) {
 841         s->noise_band_edge[j] = FFMIN(s->fft_length2, s->fft_length * get_band_edge(s, j) / s->sample_rate);
 842         if (s->noise_band_edge[j] > lrint(1.1 * s->noise_band_edge[j - 1]))
 843             i++;
 844         s->noise_band_edge[16] = i;
 845     }
 846     s->noise_band_count = s->noise_band_edge[16];
 847
 848     s->fifo = av_audio_fifo_alloc(inlink->format, inlink->channels, s->fft_length);
 849     if (!s->fifo)
 850         return AVERROR(ENOMEM);
 851
 852     return 0;
 853 }
 854
 855 static void preprocess(FFTComplex *in, int len)
 856 {
 857     double d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
 858     int n, i, k;
 859
 860     d5 = 2.0 * M_PI / len;
 861     d8 = sin(0.5 * d5);
 862     d8 = -2.0 * d8 * d8;
 863     d7 = sin(d5);
 864     d9 = 1.0 + d8;
 865     d6 = d7;
 866     n = len / 2;
 867
 868     for (i = 1; i < len / 4; i++) {
 869         k = n - i;
 870         d2 = 0.5 * (in[i].re + in[k].re);
 871         d1 = 0.5 * (in[i].im - in[k].im);
 872         d4 = 0.5 * (in[i].im + in[k].im);
 873         d3 = 0.5 * (in[k].re - in[i].re);
 874         in[i].re = d2 + d9 * d4 + d6 * d3;
 875         in[i].im = d1 + d9 * d3 - d6 * d4;
 876         in[k].re = d2 - d9 * d4 - d6 * d3;
 877         in[k].im = -d1 + d9 * d3 - d6 * d4;
 878         d10 = d9;
 879         d9 += d9 * d8 - d6 * d7;
 880         d6 += d6 * d8 + d10 * d7;
 881     }
 882
 883     d2 = in[0].re;
 884     in[0].re = d2 + in[0].im;
 885     in[0].im = d2 - in[0].im;
 886 }
 887
 888 static void postprocess(FFTComplex *in, int len)
 889 {
 890     double d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
 891     int n, i, k;
 892
 893     d5 = 2.0 * M_PI / len;
 894     d8 = sin(0.5 * d5);
 895     d8 = -2.0 * d8 * d8;
 896     d7 = sin(d5);
 897     d9 = 1.0 + d8;
 898     d6 = d7;
 899     n = len / 2;
 900     for (i = 1; i < len / 4; i++) {
 901         k = n - i;
 902         d2 = 0.5 * (in[i].re + in[k].re);
 903         d1 = 0.5 * (in[i].im - in[k].im);
 904         d4 = 0.5 * (in[i].re - in[k].re);
 905         d3 = 0.5 * (in[i].im + in[k].im);
 906         in[i].re = d2 - d9 * d3 - d6 * d4;
 907         in[i].im = d1 + d9 * d4 - d6 * d3;
 908         in[k].re = d2 + d9 * d3 + d6 * d4;
 909         in[k].im = -d1 + d9 * d4 - d6 * d3;
 910         d10 = d9;
 911         d9 += d9 * d8 - d6 * d7;
 912         d6 += d6 * d8 + d10 * d7;
 913     }
 914     d2 = in[0].re;
 915     in[0].re = 0.5 * (d2 + in[0].im);
 916     in[0].im = 0.5 * (d2 - in[0].im);
 917 }
 918
 919 static void init_sample_noise(DeNoiseChannel *dnch)
 920 {
 921     for (int i = 0; i < 15; i++) {
 922         dnch->noise_band_norm[i] = 0.0;
 923         dnch->noise_band_avr[i] = 0.0;
 924         dnch->noise_band_avi[i] = 0.0;
 925         dnch->noise_band_var[i] = 0.0;
 926     }
 927 }
 928
 929 static void sample_noise_block(AudioFFTDeNoiseContext *s,
 930                                DeNoiseChannel *dnch,
 931                                AVFrame *in, int ch)
 932 {
 933     float *src = (float *)in->extended_data[ch];
 934     double mag2, var = 0.0, avr = 0.0, avi = 0.0;
 935     int edge, j, k, n, edgemax;
 936
 937     for (int i = 0; i < s->window_length; i++) {
 938         dnch->fft_data[i].re = s->window[i] * src[i] * (1LL << 24);
 939         dnch->fft_data[i].im = 0.0;
 940     }
 941
 942     for (int i = s->window_length; i < s->fft_length2; i++) {
 943         dnch->fft_data[i].re = 0.0;
 944         dnch->fft_data[i].im = 0.0;
 945     }
 946
 947     av_fft_permute(dnch->fft, dnch->fft_data);
 948     av_fft_calc(dnch->fft, dnch->fft_data);
 949
 950     preprocess(dnch->fft_data, s->fft_length);
 951
 952     edge = s->noise_band_edge[0];
 953     j = edge;
 954     k = 0;
 955     n = j;
 956     edgemax = fmin(s->fft_length2, s->noise_band_edge[15]);
 957     dnch->fft_data[s->fft_length2].re = dnch->fft_data[0].im;
 958     dnch->fft_data[0].im = 0.0;
 959     dnch->fft_data[s->fft_length2].im = 0.0;
 960
 961     for (int i = j; i <= edgemax; i++) {
 962         if ((i == j) && (i < edgemax)) {
 963             if (j > edge) {
 964                 dnch->noise_band_norm[k - 1] += j - edge;
 965                 dnch->noise_band_avr[k - 1] += avr;
 966                 dnch->noise_band_avi[k - 1] += avi;
 967                 dnch->noise_band_var[k - 1] += var;
 968             }
 969             k++;
 970             edge = j;
 971             j = s->noise_band_edge[k];
 972             if (k == 15) {
 973                 j++;
 974             }
 975             var = 0.0;
 976             avr = 0.0;
 977             avi = 0.0;
 978         }
 979         avr += dnch->fft_data[n].re;
 980         avi += dnch->fft_data[n].im;
 981         mag2 = dnch->fft_data[n].re * dnch->fft_data[n].re +
 982                dnch->fft_data[n].im * dnch->fft_data[n].im;
 983
 984         mag2 = fmax(mag2, s->sample_floor);
 985
 986         dnch->noisy_data[i] = mag2;
 987         var += mag2;
 988         n++;
 989     }
 990
 991     dnch->noise_band_norm[k - 1] += j - edge;
 992     dnch->noise_band_avr[k - 1] += avr;
 993     dnch->noise_band_avi[k - 1] += avi;
 994     dnch->noise_band_var[k - 1] += var;
 995 }
 996
 997 static void finish_sample_noise(AudioFFTDeNoiseContext *s,
 998                                 DeNoiseChannel *dnch,
 999                                 double *sample_noise)
1000 {
1001     for (int i = 0; i < s->noise_band_count; i++) {
1002         dnch->noise_band_avr[i] /= dnch->noise_band_norm[i];
1003         dnch->noise_band_avi[i] /= dnch->noise_band_norm[i];
1004         dnch->noise_band_var[i] /= dnch->noise_band_norm[i];
1005         dnch->noise_band_var[i] -= dnch->noise_band_avr[i] * dnch->noise_band_avr[i] +
1006                                    dnch->noise_band_avi[i] * dnch->noise_band_avi[i];
1007         dnch->noise_band_auto_var[i] = dnch->noise_band_var[i];
1008         sample_noise[i] = (1.0 / C) * log(dnch->noise_band_var[i] / s->floor) - 100.0;
1009     }
1010     if (s->noise_band_count < 15) {
1011         for (int i = s->noise_band_count; i < 15; i++)
1012             sample_noise[i] = sample_noise[i - 1];
1013     }
1014 }
1015
1016 static void set_noise_profile(AudioFFTDeNoiseContext *s,
1017                               DeNoiseChannel *dnch,
1018                               double *sample_noise,
1019                               int new_profile)
1020 {
1021     int new_band_noise[15];
1022     double temp[15];
1023     double sum = 0.0, d1;
1024     float new_noise_floor;
1025     int i, n;
1026
1027     for (int m = 0; m < 15; m++)
1028         temp[m] = sample_noise[m];
1029
1030     if (new_profile) {
1031         i = 0;
1032         for (int m = 0; m < 5; m++) {
1033             sum = 0.0;
1034             for (n = 0; n < 15; n++)
1035                 sum += s->matrix_b[i++] * temp[n];
1036             s->vector_b[m] = sum;
1037         }
1038         solve(s->matrix_a, s->vector_b, 5);
1039         i = 0;
1040         for (int m = 0; m < 15; m++) {
1041             sum = 0.0;
1042             for (n = 0; n < 5; n++)
1043                 sum += s->matrix_c[i++] * s->vector_b[n];
1044             temp[m] = sum;
1045         }
1046     }
1047
1048     sum = 0.0;
1049     for (int m = 0; m < 15; m++)
1050         sum += temp[m];
1051
1052     d1 = (int)(sum / 15.0 - 0.5);
1053     if (!new_profile)
1054         i = lrint(temp[7] - d1);
1055
1056     for (d1 -= dnch->band_noise[7] - i; d1 > -20.0; d1 -= 1.0)
1057         ;
1058
1059     for (int m = 0; m < 15; m++)
1060         temp[m] -= d1;
1061
1062     new_noise_floor = d1 + 2.5;
1063
1064     if (new_profile) {
1065         av_log(s, AV_LOG_INFO, "bn=");
1066         for (int m = 0; m < 15; m++) {
1067             new_band_noise[m] = lrint(temp[m]);
1068             new_band_noise[m] = av_clip(new_band_noise[m], -24, 24);
1069             av_log(s, AV_LOG_INFO, "%d ", new_band_noise[m]);
1070         }
1071         av_log(s, AV_LOG_INFO, "\n");
1072         memcpy(dnch->band_noise, new_band_noise, sizeof(new_band_noise));
1073     }
1074
1075     if (s->track_noise)
1076         s->noise_floor = new_noise_floor;
1077 }
1078
1079 typedef struct ThreadData {
1080     AVFrame *in;
1081 } ThreadData;
1082
1083 static int filter_channel(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
1084 {
1085     AudioFFTDeNoiseContext *s = ctx->priv;
1086     ThreadData *td = arg;
1087     AVFrame *in = td->in;
1088     const int start = (in->channels * jobnr) / nb_jobs;
1089     const int end = (in->channels * (jobnr+1)) / nb_jobs;
1090
1091     for (int ch = start; ch < end; ch++) {
1092         DeNoiseChannel *dnch = &s->dnch[ch];
1093         const float *src = (const float *)in->extended_data[ch];
1094         double *dst = dnch->out_samples;
1095
1096         if (s->track_noise) {
1097             int i = s->block_count & 0x1FF;
1098
1099             if (dnch->sfm_fail_flags[i])
1100                 dnch->sfm_fail_total--;
1101             dnch->sfm_fail_flags[i] = 0;
1102             dnch->sfm_threshold *= 1.0 - dnch->sfm_alpha;
1103             dnch->sfm_threshold += dnch->sfm_alpha * (0.5 + (1.0 / 640) * dnch->sfm_fail_total);
1104         }
1105
1106         for (int m = 0; m < s->window_length; m++) {
1107             dnch->fft_data[m].re = s->window[m] * src[m] * (1LL << 24);
1108             dnch->fft_data[m].im = 0;
1109         }
1110
1111         for (int m = s->window_length; m < s->fft_length2; m++) {
1112             dnch->fft_data[m].re = 0;
1113             dnch->fft_data[m].im = 0;
1114         }
1115
1116         av_fft_permute(dnch->fft, dnch->fft_data);
1117         av_fft_calc(dnch->fft, dnch->fft_data);
1118
1119         preprocess(dnch->fft_data, s->fft_length);
1120         process_frame(s, dnch, dnch->fft_data,
1121                       dnch->prior,
1122                       dnch->prior_band_excit,
1123                       s->track_noise);
1124         postprocess(dnch->fft_data, s->fft_length);
1125
1126         av_fft_permute(dnch->ifft, dnch->fft_data);
1127         av_fft_calc(dnch->ifft, dnch->fft_data);
1128
1129         for (int m = 0; m < s->window_length; m++)
1130             dst[m] += s->window[m] * dnch->fft_data[m].re / (1LL << 24);
1131     }
1132
1133     return 0;
1134 }
1135
1136 static void get_auto_noise_levels(AudioFFTDeNoiseContext *s,
1137                                   DeNoiseChannel *dnch,
1138                                   double *levels)
1139 {
1140     if (s->noise_band_count > 0) {
1141         for (int i = 0; i < s->noise_band_count; i++) {
1142             levels[i] = (1.0 / C) * log(dnch->noise_band_auto_var[i] / s->floor) - 100.0;
1143         }
1144         if (s->noise_band_count < 15) {
1145             for (int i = s->noise_band_count; i < 15; i++)
1146                 levels[i] = levels[i - 1];
1147         }
1148     } else {
1149         for (int i = 0; i < 15; i++) {
1150             levels[i] = -100.0;
1151         }
1152     }
1153 }
1154
1155 static int filter_frame(AVFilterLink *inlink, AVFrame *frame)
1156 {
1157     AVFilterContext *ctx = inlink->dst;
1158     AVFilterLink *outlink = ctx->outputs[0];
1159     AudioFFTDeNoiseContext *s = ctx->priv;
1160     AVFrame *out = NULL, *in = NULL;
1161     ThreadData td;
1162     int ret = 0;
1163
1164     if (s->pts == AV_NOPTS_VALUE)
1165         s->pts = frame->pts;
1166
1167     ret = av_audio_fifo_write(s->fifo, (void **)frame->extended_data, frame->nb_samples);
1168     av_frame_free(&frame);
1169     if (ret < 0)
1170         return ret;
1171
1172     while (av_audio_fifo_size(s->fifo) >= s->window_length) {
1173         if (!in) {
1174             in = ff_get_audio_buffer(outlink, s->window_length);
1175             if (!in)
1176                 return AVERROR(ENOMEM);
1177         }
1178
1179         ret = av_audio_fifo_peek(s->fifo, (void **)in->extended_data, s->window_length);
1180         if (ret < 0)
1181             break;
1182
1183         if (s->track_noise) {
1184             for (int ch = 0; ch < inlink->channels; ch++) {
1185                 DeNoiseChannel *dnch = &s->dnch[ch];
1186                 double levels[15];
1187
1188                 get_auto_noise_levels(s, dnch, levels);
1189                 set_noise_profile(s, dnch, levels, 0);
1190             }
1191
1192             if (s->noise_floor != s->last_noise_floor)
1193                 set_parameters(s);
1194         }
1195
1196         if (s->sample_noise_start) {
1197             for (int ch = 0; ch < inlink->channels; ch++) {
1198                 DeNoiseChannel *dnch = &s->dnch[ch];
1199
1200                 init_sample_noise(dnch);
1201             }
1202             s->sample_noise_start = 0;
1203             s->sample_noise = 1;
1204         }
1205
1206         if (s->sample_noise) {
1207             for (int ch = 0; ch < inlink->channels; ch++) {
1208                 DeNoiseChannel *dnch = &s->dnch[ch];
1209
1210                 sample_noise_block(s, dnch, in, ch);
1211             }
1212         }
1213
1214         if (s->sample_noise_end) {
1215             for (int ch = 0; ch < inlink->channels; ch++) {
1216                 DeNoiseChannel *dnch = &s->dnch[ch];
1217                 double sample_noise[15];
1218
1219                 finish_sample_noise(s, dnch, sample_noise);
1220                 set_noise_profile(s, dnch, sample_noise, 1);
1221                 set_band_parameters(s, dnch);
1222             }
1223             s->sample_noise = 0;
1224             s->sample_noise_end = 0;
1225         }
1226
1227         s->block_count++;
1228         td.in = in;
1229         ctx->internal->execute(ctx, filter_channel, &td, NULL,
1230                                FFMIN(outlink->channels, ff_filter_get_nb_threads(ctx)));
1231
1232         out = ff_get_audio_buffer(outlink, s->sample_advance);
1233         if (!out) {
1234             ret = AVERROR(ENOMEM);
1235             break;
1236         }
1237
1238         for (int ch = 0; ch < inlink->channels; ch++) {
1239             DeNoiseChannel *dnch = &s->dnch[ch];
1240             double *src = dnch->out_samples;
1241             float *orig = (float *)in->extended_data[ch];
1242             float *dst = (float *)out->extended_data[ch];
1243
1244             switch (s->output_mode) {
1245             case IN_MODE:
1246                 for (int m = 0; m < s->sample_advance; m++)
1247                     dst[m] = orig[m];
1248                 break;
1249             case OUT_MODE:
1250                 for (int m = 0; m < s->sample_advance; m++)
1251                     dst[m] = src[m];
1252                 break;
1253             case NOISE_MODE:
1254                 for (int m = 0; m < s->sample_advance; m++)
1255                     dst[m] = orig[m] - src[m];
1256                 break;
1257             default:
1258                 return AVERROR_BUG;
1259             }
1260             memmove(src, src + s->sample_advance, (s->window_length - s->sample_advance) * sizeof(*src));
1261             memset(src + (s->window_length - s->sample_advance), 0, s->sample_advance * sizeof(*src));
1262         }
1263
1264         av_audio_fifo_drain(s->fifo, s->sample_advance);
1265
1266         out->pts = s->pts;
1267         ret = ff_filter_frame(outlink, out);
1268         if (ret < 0)
1269             break;
1270         s->pts += s->sample_advance;
1271     }
1272     av_frame_free(&in);
1273
1274     return ret;
1275 }
1276
1277 static av_cold void uninit(AVFilterContext *ctx)
1278 {
1279     AudioFFTDeNoiseContext *s = ctx->priv;
1280
1281     av_freep(&s->window);
1282     av_freep(&s->bin2band);
1283     av_freep(&s->band_alpha);
1284     av_freep(&s->band_beta);
1285
1286     if (s->dnch) {
1287         for (int ch = 0; ch < s->channels; ch++) {
1288             DeNoiseChannel *dnch = &s->dnch[ch];
1289             av_freep(&dnch->amt);
1290             av_freep(&dnch->band_amt);
1291             av_freep(&dnch->band_excit);
1292             av_freep(&dnch->gain);
1293             av_freep(&dnch->prior);
1294             av_freep(&dnch->prior_band_excit);
1295             av_freep(&dnch->clean_data);
1296             av_freep(&dnch->noisy_data);
1297             av_freep(&dnch->out_samples);
1298             av_freep(&dnch->spread_function);
1299             av_freep(&dnch->abs_var);
1300             av_freep(&dnch->rel_var);
1301             av_freep(&dnch->min_abs_var);
1302             av_freep(&dnch->fft_data);
1303             av_fft_end(dnch->fft);
1304             dnch->fft = NULL;
1305             av_fft_end(dnch->ifft);
1306             dnch->ifft = NULL;
1307         }
1308         av_freep(&s->dnch);
1309     }
1310
1311     av_audio_fifo_free(s->fifo);
1312 }
1313
1314 static int query_formats(AVFilterContext *ctx)
1315 {
1316     AVFilterFormats *formats = NULL;
1317     AVFilterChannelLayouts *layouts = NULL;
1318     static const enum AVSampleFormat sample_fmts[] = {
1319         AV_SAMPLE_FMT_FLTP,
1320         AV_SAMPLE_FMT_NONE
1321     };
1322     int ret;
1323
1324     formats = ff_make_format_list(sample_fmts);
1325     if (!formats)
1326         return AVERROR(ENOMEM);
1327     ret = ff_set_common_formats(ctx, formats);
1328     if (ret < 0)
1329         return ret;
1330
1331     layouts = ff_all_channel_counts();
1332     if (!layouts)
1333         return AVERROR(ENOMEM);
1334
1335     ret = ff_set_common_channel_layouts(ctx, layouts);
1336     if (ret < 0)
1337         return ret;
1338
1339     formats = ff_all_samplerates();
1340     return ff_set_common_samplerates(ctx, formats);
1341 }
1342
1343 static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
1344                            char *res, int res_len, int flags)
1345 {
1346     AudioFFTDeNoiseContext *s = ctx->priv;
1347     int need_reset = 0;
1348
1349     if (!strcmp(cmd, "sample_noise") ||
1350         !strcmp(cmd, "sn")) {
1351         if (!strcmp(args, "start")) {
1352             s->sample_noise_start = 1;
1353             s->sample_noise_end = 0;
1354         } else if (!strcmp(args, "end")) {
1355             s->sample_noise_start = 0;
1356             s->sample_noise_end = 1;
1357         }
1358     } else if (!strcmp(cmd, "nr") ||
1359                !strcmp(cmd, "noise_reduction")) {
1360         float nr;
1361
1362         if (sscanf(args, "%f", &nr) == 1) {
1363             s->noise_reduction = av_clipf(nr, 0.01, 97);
1364             need_reset = 1;
1365         }
1366     } else if (!strcmp(cmd, "nf") ||
1367                !strcmp(cmd, "noise_floor")) {
1368         float nf;
1369
1370         if (sscanf(args, "%f", &nf) == 1) {
1371             s->noise_floor = av_clipf(nf, -80, -20);
1372             need_reset = 1;
1373         }
1374     } else if (!strcmp(cmd, "output_mode") ||
1375                !strcmp(cmd, "om")) {
1376         if (!strcmp(args, "i")) {
1377             s->output_mode = IN_MODE;
1378         } else if (!strcmp(args, "o")) {
1379             s->output_mode = OUT_MODE;
1380         } else if (!strcmp(args, "n")) {
1381             s->output_mode = NOISE_MODE;
1382         }
1383     }
1384
1385     if (need_reset)
1386         set_parameters(s);
1387
1388     return 0;
1389 }
1390
1391 static const AVFilterPad inputs[] = {
1392     {
1393         .name         = "default",
1394         .type         = AVMEDIA_TYPE_AUDIO,
1395         .filter_frame = filter_frame,
1396         .config_props = config_input,
1397     },
1398     { NULL }
1399 };
1400
1401 static const AVFilterPad outputs[] = {
1402     {
1403         .name = "default",
1404         .type = AVMEDIA_TYPE_AUDIO,
1405     },
1406     { NULL }
1407 };
1408
1409 AVFilter ff_af_afftdn = {
1410     .name            = "afftdn",
1411     .description     = NULL_IF_CONFIG_SMALL("Denoise audio samples using FFT."),
1412     .query_formats   = query_formats,
1413     .priv_size       = sizeof(AudioFFTDeNoiseContext),
1414     .priv_class      = &afftdn_class,
1415     .uninit          = uninit,
1416     .inputs          = inputs,
1417     .outputs         = outputs,
1418     .process_command = process_command,
1419     .flags           = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
1420                        AVFILTER_FLAG_SLICE_THREADS,
1421 };