avfilter/af_silenceremove: prefer outlink instead of inlink
[ffmpeg.git] / 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             break;
695         case CUSTOM_NOISE:
696             read_custom_noise(s, ch);
697             break;
698         default:
699             return AVERROR_BUG;
700         }
701
702
703         dnch->sfm_threshold = 0.8;
704         dnch->sfm_alpha = 0.05;
705         for (i = 0; i < 512; i++)
706             dnch->sfm_fail_flags[i] = 0;
707
708         dnch->sfm_fail_total = 0;
709         j = FFMAX((int)(10.0 * (1.3 - dnch->sfm_threshold)), 1);
710
711         for (i = 0; i < 512; i += j) {
712             dnch->sfm_fail_flags[i] = 1;
713             dnch->sfm_fail_total += 1;
714         }
715
716         dnch->amt = av_calloc(s->bin_count, sizeof(*dnch->amt));
717         dnch->band_amt = av_calloc(s->number_of_bands, sizeof(*dnch->band_amt));
718         dnch->band_excit = av_calloc(s->number_of_bands, sizeof(*dnch->band_excit));
719         dnch->gain = av_calloc(s->bin_count, sizeof(*dnch->gain));
720         dnch->prior = av_calloc(s->bin_count, sizeof(*dnch->prior));
721         dnch->prior_band_excit = av_calloc(s->number_of_bands, sizeof(*dnch->prior_band_excit));
722         dnch->clean_data = av_calloc(s->bin_count, sizeof(*dnch->clean_data));
723         dnch->noisy_data = av_calloc(s->bin_count, sizeof(*dnch->noisy_data));
724         dnch->out_samples = av_calloc(s->buffer_length, sizeof(*dnch->out_samples));
725         dnch->abs_var = av_calloc(s->bin_count, sizeof(*dnch->abs_var));
726         dnch->rel_var = av_calloc(s->bin_count, sizeof(*dnch->rel_var));
727         dnch->min_abs_var = av_calloc(s->bin_count, sizeof(*dnch->min_abs_var));
728         dnch->fft_data = av_calloc(s->fft_length2 + 1, sizeof(*dnch->fft_data));
729         dnch->fft  = av_fft_init(av_log2(s->fft_length2), 0);
730         dnch->ifft = av_fft_init(av_log2(s->fft_length2), 1);
731         dnch->spread_function = av_calloc(s->number_of_bands * s->number_of_bands,
732                                           sizeof(*dnch->spread_function));
733
734         if (!dnch->amt ||
735             !dnch->band_amt ||
736             !dnch->band_excit ||
737             !dnch->gain ||
738             !dnch->prior ||
739             !dnch->prior_band_excit ||
740             !dnch->clean_data ||
741             !dnch->noisy_data ||
742             !dnch->out_samples ||
743             !dnch->fft_data ||
744             !dnch->abs_var ||
745             !dnch->rel_var ||
746             !dnch->min_abs_var ||
747             !dnch->spread_function ||
748             !dnch->fft ||
749             !dnch->ifft)
750             return AVERROR(ENOMEM);
751     }
752
753     for (int ch = 0; ch < inlink->channels; ch++) {
754         DeNoiseChannel *dnch = &s->dnch[ch];
755         double *prior_band_excit = dnch->prior_band_excit;
756         double *prior = dnch->prior;
757         double min, max;
758         double p1, p2;
759
760         p1 = pow(0.1, 2.5 / sdiv);
761         p2 = pow(0.1, 1.0 / sdiv);
762         j = 0;
763         for (m = 0; m < s->number_of_bands; m++) {
764             for (n = 0; n < s->number_of_bands; n++) {
765                 if (n < m) {
766                     dnch->spread_function[j++] = pow(p2, m - n);
767                 } else if (n > m) {
768                     dnch->spread_function[j++] = pow(p1, n - m);
769                 } else {
770                     dnch->spread_function[j++] = 1.0;
771                 }
772             }
773         }
774
775         for (m = 0; m < s->number_of_bands; m++) {
776             dnch->band_excit[m] = 0.0;
777             prior_band_excit[m] = 0.0;
778         }
779
780         for (m = 0; m <= s->fft_length2; m++)
781             dnch->band_excit[s->bin2band[m]] += 1.0;
782
783         j = 0;
784         for (m = 0; m < s->number_of_bands; m++) {
785             for (n = 0; n < s->number_of_bands; n++)
786                 prior_band_excit[m] += dnch->spread_function[j++] * dnch->band_excit[n];
787         }
788
789         min = pow(0.1, 2.5);
790         max = pow(0.1, 1.0);
791         for (int i = 0; i < s->number_of_bands; i++) {
792             if (i < lrint(12.0 * sdiv)) {
793                 dnch->band_excit[i] = pow(0.1, 1.45 + 0.1 * i / sdiv);
794             } else {
795                 dnch->band_excit[i] = pow(0.1, 2.5 - 0.2 * (i / sdiv - 14.0));
796             }
797             dnch->band_excit[i] = av_clipd(dnch->band_excit[i], min, max);
798         }
799
800         for (int i = 0; i <= s->fft_length2; i++)
801             prior[i] = RRATIO;
802         for (int i = 0; i < s->buffer_length; i++)
803             dnch->out_samples[i] = 0;
804
805         j = 0;
806         for (int i = 0; i < s->number_of_bands; i++)
807             for (int k = 0; k < s->number_of_bands; k++)
808                 dnch->spread_function[j++] *= dnch->band_excit[i] / prior_band_excit[i];
809     }
810
811     j = 0;
812     sar = s->sample_advance / s->sample_rate;
813     for (int i = 0; i <= s->fft_length2; i++) {
814         if ((i == s->fft_length2) || (s->bin2band[i] > j)) {
815             double d6 = (i - 1) * s->sample_rate / s->fft_length;
816             double d7 = fmin(0.008 + 2.2 / d6, 0.03);
817             s->band_alpha[j] = exp(-sar / d7);
818             s->band_beta[j] = 1.0 - s->band_alpha[j];
819             j = s->bin2band[i];
820         }
821     }
822
823     wscale = sqrt(16.0 / (9.0 * s->fft_length));
824     sum = 0.0;
825     for (int i = 0; i < s->window_length; i++) {
826         double d10 = sin(i * M_PI / s->window_length);
827         d10 *= wscale * d10;
828         s->window[i] = d10;
829         sum += d10 * d10;
830     }
831
832     s->window_weight = 0.5 * sum;
833     s->floor = (1LL << 48) * exp(-23.025558369790467) * s->window_weight;
834     s->sample_floor = s->floor * exp(4.144600506562284);
835     s->auto_floor = s->floor * exp(6.907667510937141);
836
837     set_parameters(s);
838
839     s->noise_band_edge[0] = FFMIN(s->fft_length2, s->fft_length * get_band_edge(s, 0) / s->sample_rate);
840     i = 0;
841     for (int j = 1; j < 16; j++) {
842         s->noise_band_edge[j] = FFMIN(s->fft_length2, s->fft_length * get_band_edge(s, j) / s->sample_rate);
843         if (s->noise_band_edge[j] > lrint(1.1 * s->noise_band_edge[j - 1]))
844             i++;
845         s->noise_band_edge[16] = i;
846     }
847     s->noise_band_count = s->noise_band_edge[16];
848
849     s->fifo = av_audio_fifo_alloc(inlink->format, inlink->channels, s->fft_length);
850     if (!s->fifo)
851         return AVERROR(ENOMEM);
852
853     return 0;
854 }
855
856 static void preprocess(FFTComplex *in, int len)
857 {
858     double d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
859     int n, i, k;
860
861     d5 = 2.0 * M_PI / len;
862     d8 = sin(0.5 * d5);
863     d8 = -2.0 * d8 * d8;
864     d7 = sin(d5);
865     d9 = 1.0 + d8;
866     d6 = d7;
867     n = len / 2;
868
869     for (i = 1; i < len / 4; i++) {
870         k = n - i;
871         d2 = 0.5 * (in[i].re + in[k].re);
872         d1 = 0.5 * (in[i].im - in[k].im);
873         d4 = 0.5 * (in[i].im + in[k].im);
874         d3 = 0.5 * (in[k].re - in[i].re);
875         in[i].re = d2 + d9 * d4 + d6 * d3;
876         in[i].im = d1 + d9 * d3 - d6 * d4;
877         in[k].re = d2 - d9 * d4 - d6 * d3;
878         in[k].im = -d1 + d9 * d3 - d6 * d4;
879         d10 = d9;
880         d9 += d9 * d8 - d6 * d7;
881         d6 += d6 * d8 + d10 * d7;
882     }
883
884     d2 = in[0].re;
885     in[0].re = d2 + in[0].im;
886     in[0].im = d2 - in[0].im;
887 }
888
889 static void postprocess(FFTComplex *in, int len)
890 {
891     double d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
892     int n, i, k;
893
894     d5 = 2.0 * M_PI / len;
895     d8 = sin(0.5 * d5);
896     d8 = -2.0 * d8 * d8;
897     d7 = sin(d5);
898     d9 = 1.0 + d8;
899     d6 = d7;
900     n = len / 2;
901     for (i = 1; i < len / 4; i++) {
902         k = n - i;
903         d2 = 0.5 * (in[i].re + in[k].re);
904         d1 = 0.5 * (in[i].im - in[k].im);
905         d4 = 0.5 * (in[i].re - in[k].re);
906         d3 = 0.5 * (in[i].im + in[k].im);
907         in[i].re = d2 - d9 * d3 - d6 * d4;
908         in[i].im = d1 + d9 * d4 - d6 * d3;
909         in[k].re = d2 + d9 * d3 + d6 * d4;
910         in[k].im = -d1 + d9 * d4 - d6 * d3;
911         d10 = d9;
912         d9 += d9 * d8 - d6 * d7;
913         d6 += d6 * d8 + d10 * d7;
914     }
915     d2 = in[0].re;
916     in[0].re = 0.5 * (d2 + in[0].im);
917     in[0].im = 0.5 * (d2 - in[0].im);
918 }
919
920 static void init_sample_noise(DeNoiseChannel *dnch)
921 {
922     for (int i = 0; i < 15; i++) {
923         dnch->noise_band_norm[i] = 0.0;
924         dnch->noise_band_avr[i] = 0.0;
925         dnch->noise_band_avi[i] = 0.0;
926         dnch->noise_band_var[i] = 0.0;
927     }
928 }
929
930 static void sample_noise_block(AudioFFTDeNoiseContext *s,
931                                DeNoiseChannel *dnch,
932                                AVFrame *in, int ch)
933 {
934     float *src = (float *)in->extended_data[ch];
935     double mag2, var = 0.0, avr = 0.0, avi = 0.0;
936     int edge, j, k, n, edgemax;
937
938     for (int i = 0; i < s->window_length; i++) {
939         dnch->fft_data[i].re = s->window[i] * src[i] * (1LL << 24);
940         dnch->fft_data[i].im = 0.0;
941     }
942
943     for (int i = s->window_length; i < s->fft_length2; i++) {
944         dnch->fft_data[i].re = 0.0;
945         dnch->fft_data[i].im = 0.0;
946     }
947
948     av_fft_permute(dnch->fft, dnch->fft_data);
949     av_fft_calc(dnch->fft, dnch->fft_data);
950
951     preprocess(dnch->fft_data, s->fft_length);
952
953     edge = s->noise_band_edge[0];
954     j = edge;
955     k = 0;
956     n = j;
957     edgemax = fmin(s->fft_length2, s->noise_band_edge[15]);
958     dnch->fft_data[s->fft_length2].re = dnch->fft_data[0].im;
959     dnch->fft_data[0].im = 0.0;
960     dnch->fft_data[s->fft_length2].im = 0.0;
961
962     for (int i = j; i <= edgemax; i++) {
963         if ((i == j) && (i < edgemax)) {
964             if (j > edge) {
965                 dnch->noise_band_norm[k - 1] += j - edge;
966                 dnch->noise_band_avr[k - 1] += avr;
967                 dnch->noise_band_avi[k - 1] += avi;
968                 dnch->noise_band_var[k - 1] += var;
969             }
970             k++;
971             edge = j;
972             j = s->noise_band_edge[k];
973             if (k == 15) {
974                 j++;
975             }
976             var = 0.0;
977             avr = 0.0;
978             avi = 0.0;
979         }
980         avr += dnch->fft_data[n].re;
981         avi += dnch->fft_data[n].im;
982         mag2 = dnch->fft_data[n].re * dnch->fft_data[n].re +
983                dnch->fft_data[n].im * dnch->fft_data[n].im;
984
985         mag2 = fmax(mag2, s->sample_floor);
986
987         dnch->noisy_data[i] = mag2;
988         var += mag2;
989         n++;
990     }
991
992     dnch->noise_band_norm[k - 1] += j - edge;
993     dnch->noise_band_avr[k - 1] += avr;
994     dnch->noise_band_avi[k - 1] += avi;
995     dnch->noise_band_var[k - 1] += var;
996 }
997
998 static void finish_sample_noise(AudioFFTDeNoiseContext *s,
999                                 DeNoiseChannel *dnch,
1000                                 double *sample_noise)
1001 {
1002     for (int i = 0; i < s->noise_band_count; i++) {
1003         dnch->noise_band_avr[i] /= dnch->noise_band_norm[i];
1004         dnch->noise_band_avi[i] /= dnch->noise_band_norm[i];
1005         dnch->noise_band_var[i] /= dnch->noise_band_norm[i];
1006         dnch->noise_band_var[i] -= dnch->noise_band_avr[i] * dnch->noise_band_avr[i] +
1007                                    dnch->noise_band_avi[i] * dnch->noise_band_avi[i];
1008         dnch->noise_band_auto_var[i] = dnch->noise_band_var[i];
1009         sample_noise[i] = (1.0 / C) * log(dnch->noise_band_var[i] / s->floor) - 100.0;
1010     }
1011     if (s->noise_band_count < 15) {
1012         for (int i = s->noise_band_count; i < 15; i++)
1013             sample_noise[i] = sample_noise[i - 1];
1014     }
1015 }
1016
1017 static void set_noise_profile(AudioFFTDeNoiseContext *s,
1018                               DeNoiseChannel *dnch,
1019                               double *sample_noise,
1020                               int new_profile)
1021 {
1022     int new_band_noise[15];
1023     double temp[15];
1024     double sum = 0.0, d1;
1025     float new_noise_floor;
1026     int i, n;
1027
1028     for (int m = 0; m < 15; m++)
1029         temp[m] = sample_noise[m];
1030
1031     if (new_profile) {
1032         i = 0;
1033         for (int m = 0; m < 5; m++) {
1034             sum = 0.0;
1035             for (n = 0; n < 15; n++)
1036                 sum += s->matrix_b[i++] * temp[n];
1037             s->vector_b[m] = sum;
1038         }
1039         solve(s->matrix_a, s->vector_b, 5);
1040         i = 0;
1041         for (int m = 0; m < 15; m++) {
1042             sum = 0.0;
1043             for (n = 0; n < 5; n++)
1044                 sum += s->matrix_c[i++] * s->vector_b[n];
1045             temp[m] = sum;
1046         }
1047     }
1048
1049     sum = 0.0;
1050     for (int m = 0; m < 15; m++)
1051         sum += temp[m];
1052
1053     d1 = (int)(sum / 15.0 - 0.5);
1054     if (!new_profile)
1055         i = lrint(temp[7] - d1);
1056
1057     for (d1 -= dnch->band_noise[7] - i; d1 > -20.0; d1 -= 1.0)
1058         ;
1059
1060     for (int m = 0; m < 15; m++)
1061         temp[m] -= d1;
1062
1063     new_noise_floor = d1 + 2.5;
1064
1065     if (new_profile) {
1066         av_log(s, AV_LOG_INFO, "bn=");
1067         for (int m = 0; m < 15; m++) {
1068             new_band_noise[m] = lrint(temp[m]);
1069             new_band_noise[m] = av_clip(new_band_noise[m], -24, 24);
1070             av_log(s, AV_LOG_INFO, "%d ", new_band_noise[m]);
1071         }
1072         av_log(s, AV_LOG_INFO, "\n");
1073         memcpy(dnch->band_noise, new_band_noise, sizeof(new_band_noise));
1074     }
1075
1076     if (s->track_noise)
1077         s->noise_floor = new_noise_floor;
1078 }
1079
1080 typedef struct ThreadData {
1081     AVFrame *in;
1082 } ThreadData;
1083
1084 static int filter_channel(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
1085 {
1086     AudioFFTDeNoiseContext *s = ctx->priv;
1087     ThreadData *td = arg;
1088     AVFrame *in = td->in;
1089     const int start = (in->channels * jobnr) / nb_jobs;
1090     const int end = (in->channels * (jobnr+1)) / nb_jobs;
1091
1092     for (int ch = start; ch < end; ch++) {
1093         DeNoiseChannel *dnch = &s->dnch[ch];
1094         const float *src = (const float *)in->extended_data[ch];
1095         double *dst = dnch->out_samples;
1096
1097         if (s->track_noise) {
1098             int i = s->block_count & 0x1FF;
1099
1100             if (dnch->sfm_fail_flags[i])
1101                 dnch->sfm_fail_total--;
1102             dnch->sfm_fail_flags[i] = 0;
1103             dnch->sfm_threshold *= 1.0 - dnch->sfm_alpha;
1104             dnch->sfm_threshold += dnch->sfm_alpha * (0.5 + (1.0 / 640) * dnch->sfm_fail_total);
1105         }
1106
1107         for (int m = 0; m < s->window_length; m++) {
1108             dnch->fft_data[m].re = s->window[m] * src[m] * (1LL << 24);
1109             dnch->fft_data[m].im = 0;
1110         }
1111
1112         for (int m = s->window_length; m < s->fft_length2; m++) {
1113             dnch->fft_data[m].re = 0;
1114             dnch->fft_data[m].im = 0;
1115         }
1116
1117         av_fft_permute(dnch->fft, dnch->fft_data);
1118         av_fft_calc(dnch->fft, dnch->fft_data);
1119
1120         preprocess(dnch->fft_data, s->fft_length);
1121         process_frame(s, dnch, dnch->fft_data,
1122                       dnch->prior,
1123                       dnch->prior_band_excit,
1124                       s->track_noise);
1125         postprocess(dnch->fft_data, s->fft_length);
1126
1127         av_fft_permute(dnch->ifft, dnch->fft_data);
1128         av_fft_calc(dnch->ifft, dnch->fft_data);
1129
1130         for (int m = 0; m < s->window_length; m++)
1131             dst[m] += s->window[m] * dnch->fft_data[m].re / (1LL << 24);
1132     }
1133
1134     return 0;
1135 }
1136
1137 static void get_auto_noise_levels(AudioFFTDeNoiseContext *s,
1138                                   DeNoiseChannel *dnch,
1139                                   double *levels)
1140 {
1141     if (s->noise_band_count > 0) {
1142         for (int i = 0; i < s->noise_band_count; i++) {
1143             levels[i] = (1.0 / C) * log(dnch->noise_band_auto_var[i] / s->floor) - 100.0;
1144         }
1145         if (s->noise_band_count < 15) {
1146             for (int i = s->noise_band_count; i < 15; i++)
1147                 levels[i] = levels[i - 1];
1148         }
1149     } else {
1150         for (int i = 0; i < 15; i++) {
1151             levels[i] = -100.0;
1152         }
1153     }
1154 }
1155
1156 static int filter_frame(AVFilterLink *inlink, AVFrame *frame)
1157 {
1158     AVFilterContext *ctx = inlink->dst;
1159     AVFilterLink *outlink = ctx->outputs[0];
1160     AudioFFTDeNoiseContext *s = ctx->priv;
1161     AVFrame *out = NULL, *in = NULL;
1162     ThreadData td;
1163     int ret = 0;
1164
1165     if (s->pts == AV_NOPTS_VALUE)
1166         s->pts = frame->pts;
1167
1168     ret = av_audio_fifo_write(s->fifo, (void **)frame->extended_data, frame->nb_samples);
1169     av_frame_free(&frame);
1170     if (ret < 0)
1171         return ret;
1172
1173     while (av_audio_fifo_size(s->fifo) >= s->window_length) {
1174         if (!in) {
1175             in = ff_get_audio_buffer(outlink, s->window_length);
1176             if (!in)
1177                 return AVERROR(ENOMEM);
1178         }
1179
1180         ret = av_audio_fifo_peek(s->fifo, (void **)in->extended_data, s->window_length);
1181         if (ret < 0)
1182             break;
1183
1184         if (s->track_noise) {
1185             for (int ch = 0; ch < inlink->channels; ch++) {
1186                 DeNoiseChannel *dnch = &s->dnch[ch];
1187                 double levels[15];
1188
1189                 get_auto_noise_levels(s, dnch, levels);
1190                 set_noise_profile(s, dnch, levels, 0);
1191             }
1192
1193             if (s->noise_floor != s->last_noise_floor)
1194                 set_parameters(s);
1195         }
1196
1197         if (s->sample_noise_start) {
1198             for (int ch = 0; ch < inlink->channels; ch++) {
1199                 DeNoiseChannel *dnch = &s->dnch[ch];
1200
1201                 init_sample_noise(dnch);
1202             }
1203             s->sample_noise_start = 0;
1204             s->sample_noise = 1;
1205         }
1206
1207         if (s->sample_noise) {
1208             for (int ch = 0; ch < inlink->channels; ch++) {
1209                 DeNoiseChannel *dnch = &s->dnch[ch];
1210
1211                 sample_noise_block(s, dnch, in, ch);
1212             }
1213         }
1214
1215         if (s->sample_noise_end) {
1216             for (int ch = 0; ch < inlink->channels; ch++) {
1217                 DeNoiseChannel *dnch = &s->dnch[ch];
1218                 double sample_noise[15];
1219
1220                 finish_sample_noise(s, dnch, sample_noise);
1221                 set_noise_profile(s, dnch, sample_noise, 1);
1222                 set_band_parameters(s, dnch);
1223             }
1224             s->sample_noise = 0;
1225             s->sample_noise_end = 0;
1226         }
1227
1228         s->block_count++;
1229         td.in = in;
1230         ctx->internal->execute(ctx, filter_channel, &td, NULL,
1231                                FFMIN(outlink->channels, ff_filter_get_nb_threads(ctx)));
1232
1233         out = ff_get_audio_buffer(outlink, s->sample_advance);
1234         if (!out) {
1235             ret = AVERROR(ENOMEM);
1236             break;
1237         }
1238
1239         for (int ch = 0; ch < inlink->channels; ch++) {
1240             DeNoiseChannel *dnch = &s->dnch[ch];
1241             double *src = dnch->out_samples;
1242             float *orig = (float *)in->extended_data[ch];
1243             float *dst = (float *)out->extended_data[ch];
1244
1245             switch (s->output_mode) {
1246             case IN_MODE:
1247                 for (int m = 0; m < s->sample_advance; m++)
1248                     dst[m] = orig[m];
1249                 break;
1250             case OUT_MODE:
1251                 for (int m = 0; m < s->sample_advance; m++)
1252                     dst[m] = src[m];
1253                 break;
1254             case NOISE_MODE:
1255                 for (int m = 0; m < s->sample_advance; m++)
1256                     dst[m] = orig[m] - src[m];
1257                 break;
1258             default:
1259                 return AVERROR_BUG;
1260             }
1261             memmove(src, src + s->sample_advance, (s->window_length - s->sample_advance) * sizeof(*src));
1262             memset(src + (s->window_length - s->sample_advance), 0, s->sample_advance * sizeof(*src));
1263         }
1264
1265         av_audio_fifo_drain(s->fifo, s->sample_advance);
1266
1267         out->pts = s->pts;
1268         ret = ff_filter_frame(outlink, out);
1269         if (ret < 0)
1270             break;
1271         s->pts += s->sample_advance;
1272     }
1273     av_frame_free(&in);
1274
1275     return ret;
1276 }
1277
1278 static av_cold void uninit(AVFilterContext *ctx)
1279 {
1280     AudioFFTDeNoiseContext *s = ctx->priv;
1281
1282     av_freep(&s->window);
1283     av_freep(&s->bin2band);
1284     av_freep(&s->band_alpha);
1285     av_freep(&s->band_beta);
1286
1287     if (s->dnch) {
1288         for (int ch = 0; ch < s->channels; ch++) {
1289             DeNoiseChannel *dnch = &s->dnch[ch];
1290             av_freep(&dnch->amt);
1291             av_freep(&dnch->band_amt);
1292             av_freep(&dnch->band_excit);
1293             av_freep(&dnch->gain);
1294             av_freep(&dnch->prior);
1295             av_freep(&dnch->prior_band_excit);
1296             av_freep(&dnch->clean_data);
1297             av_freep(&dnch->noisy_data);
1298             av_freep(&dnch->out_samples);
1299             av_freep(&dnch->spread_function);
1300             av_freep(&dnch->abs_var);
1301             av_freep(&dnch->rel_var);
1302             av_freep(&dnch->min_abs_var);
1303             av_freep(&dnch->fft_data);
1304             av_fft_end(dnch->fft);
1305             dnch->fft = NULL;
1306             av_fft_end(dnch->ifft);
1307             dnch->ifft = NULL;
1308         }
1309         av_freep(&s->dnch);
1310     }
1311
1312     av_audio_fifo_free(s->fifo);
1313 }
1314
1315 static int query_formats(AVFilterContext *ctx)
1316 {
1317     AVFilterFormats *formats = NULL;
1318     AVFilterChannelLayouts *layouts = NULL;
1319     static const enum AVSampleFormat sample_fmts[] = {
1320         AV_SAMPLE_FMT_FLTP,
1321         AV_SAMPLE_FMT_NONE
1322     };
1323     int ret;
1324
1325     formats = ff_make_format_list(sample_fmts);
1326     if (!formats)
1327         return AVERROR(ENOMEM);
1328     ret = ff_set_common_formats(ctx, formats);
1329     if (ret < 0)
1330         return ret;
1331
1332     layouts = ff_all_channel_counts();
1333     if (!layouts)
1334         return AVERROR(ENOMEM);
1335
1336     ret = ff_set_common_channel_layouts(ctx, layouts);
1337     if (ret < 0)
1338         return ret;
1339
1340     formats = ff_all_samplerates();
1341     return ff_set_common_samplerates(ctx, formats);
1342 }
1343
1344 static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
1345                            char *res, int res_len, int flags)
1346 {
1347     AudioFFTDeNoiseContext *s = ctx->priv;
1348     int need_reset = 0;
1349
1350     if (!strcmp(cmd, "sample_noise") ||
1351         !strcmp(cmd, "sn")) {
1352         if (!strcmp(args, "start")) {
1353             s->sample_noise_start = 1;
1354             s->sample_noise_end = 0;
1355         } else if (!strcmp(args, "end")) {
1356             s->sample_noise_start = 0;
1357             s->sample_noise_end = 1;
1358         }
1359     } else if (!strcmp(cmd, "nr") ||
1360                !strcmp(cmd, "noise_reduction")) {
1361         float nr;
1362
1363         if (sscanf(args, "%f", &nr) == 1) {
1364             s->noise_reduction = av_clipf(nr, 0.01, 97);
1365             need_reset = 1;
1366         }
1367     } else if (!strcmp(cmd, "nf") ||
1368                !strcmp(cmd, "noise_floor")) {
1369         float nf;
1370
1371         if (sscanf(args, "%f", &nf) == 1) {
1372             s->noise_floor = av_clipf(nf, -80, -20);
1373             need_reset = 1;
1374         }
1375     } else if (!strcmp(cmd, "output_mode") ||
1376                !strcmp(cmd, "om")) {
1377         if (!strcmp(args, "i")) {
1378             s->output_mode = IN_MODE;
1379         } else if (!strcmp(args, "o")) {
1380             s->output_mode = OUT_MODE;
1381         } else if (!strcmp(args, "n")) {
1382             s->output_mode = NOISE_MODE;
1383         }
1384     }
1385
1386     if (need_reset)
1387         set_parameters(s);
1388
1389     return 0;
1390 }
1391
1392 static const AVFilterPad inputs[] = {
1393     {
1394         .name         = "default",
1395         .type         = AVMEDIA_TYPE_AUDIO,
1396         .filter_frame = filter_frame,
1397         .config_props = config_input,
1398     },
1399     { NULL }
1400 };
1401
1402 static const AVFilterPad outputs[] = {
1403     {
1404         .name = "default",
1405         .type = AVMEDIA_TYPE_AUDIO,
1406     },
1407     { NULL }
1408 };
1409
1410 AVFilter ff_af_afftdn = {
1411     .name            = "afftdn",
1412     .description     = NULL_IF_CONFIG_SMALL("Denoise audio samples using FFT."),
1413     .query_formats   = query_formats,
1414     .priv_size       = sizeof(AudioFFTDeNoiseContext),
1415     .priv_class      = &afftdn_class,
1416     .uninit          = uninit,
1417     .inputs          = inputs,
1418     .outputs         = outputs,
1419     .process_command = process_command,
1420     .flags           = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
1421                        AVFILTER_FLAG_SLICE_THREADS,
1422 };