FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavfilter/af_arnndn.c
Date: 2026-09-25 23:13:28
Exec Total Coverage
Lines: 0 761 0.0%
Functions: 0 41 0.0%
Branches: 0 698 0.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2018 Gregor Richards
3 * Copyright (c) 2017 Mozilla
4 * Copyright (c) 2005-2009 Xiph.Org Foundation
5 * Copyright (c) 2007-2008 CSIRO
6 * Copyright (c) 2008-2011 Octasic Inc.
7 * Copyright (c) Jean-Marc Valin
8 * Copyright (c) 2019 Paul B Mahol
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 *
14 * - Redistributions of source code must retain the above copyright
15 * notice, this list of conditions and the following disclaimer.
16 *
17 * - Redistributions in binary form must reproduce the above copyright
18 * notice, this list of conditions and the following disclaimer in the
19 * documentation and/or other materials provided with the distribution.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
25 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #include "libavutil/avassert.h"
35 #include "libavutil/file_open.h"
36 #include "libavutil/float_dsp.h"
37 #include "libavutil/mem.h"
38 #include "libavutil/mem_internal.h"
39 #include "libavutil/opt.h"
40 #include "libavutil/tx.h"
41 #include "avfilter.h"
42 #include "audio.h"
43 #include "filters.h"
44 #include "formats.h"
45
46 #define FRAME_SIZE_SHIFT 2
47 #define FRAME_SIZE (120<<FRAME_SIZE_SHIFT)
48 #define WINDOW_SIZE (2*FRAME_SIZE)
49 #define FREQ_SIZE (FRAME_SIZE + 1)
50
51 #define PITCH_MIN_PERIOD 60
52 #define PITCH_MAX_PERIOD 768
53 #define PITCH_FRAME_SIZE 960
54 #define PITCH_BUF_SIZE (PITCH_MAX_PERIOD+PITCH_FRAME_SIZE)
55
56 #define SQUARE(x) ((x)*(x))
57
58 #define NB_BANDS 22
59
60 #define CEPS_MEM 8
61 #define NB_DELTA_CEPS 6
62
63 #define NB_FEATURES (NB_BANDS+3*NB_DELTA_CEPS+2)
64
65 #define WEIGHTS_SCALE (1.f/256)
66
67 #define MAX_NEURONS 128
68
69 #define ACTIVATION_TANH 0
70 #define ACTIVATION_SIGMOID 1
71 #define ACTIVATION_RELU 2
72
73 #define Q15ONE 1.0f
74
75 typedef struct DenseLayer {
76 const float *bias;
77 const float *input_weights;
78 int nb_inputs;
79 int nb_neurons;
80 int activation;
81 } DenseLayer;
82
83 typedef struct GRULayer {
84 const float *bias;
85 const float *input_weights;
86 const float *recurrent_weights;
87 int nb_inputs;
88 int nb_neurons;
89 int activation;
90 } GRULayer;
91
92 typedef struct RNNModel {
93 int input_dense_size;
94 const DenseLayer *input_dense;
95
96 int vad_gru_size;
97 const GRULayer *vad_gru;
98
99 int noise_gru_size;
100 const GRULayer *noise_gru;
101
102 int denoise_gru_size;
103 const GRULayer *denoise_gru;
104
105 int denoise_output_size;
106 const DenseLayer *denoise_output;
107
108 int vad_output_size;
109 const DenseLayer *vad_output;
110 } RNNModel;
111
112 typedef struct RNNState {
113 float *vad_gru_state;
114 float *noise_gru_state;
115 float *denoise_gru_state;
116 RNNModel *model;
117 } RNNState;
118
119 typedef struct DenoiseState {
120 float analysis_mem[FRAME_SIZE];
121 float cepstral_mem[CEPS_MEM][NB_BANDS];
122 int memid;
123 DECLARE_ALIGNED(32, float, synthesis_mem)[FRAME_SIZE];
124 float pitch_buf[PITCH_BUF_SIZE];
125 float pitch_enh_buf[PITCH_BUF_SIZE];
126 float last_gain;
127 int last_period;
128 float mem_hp_x[2];
129 float lastg[NB_BANDS];
130 float history[FRAME_SIZE];
131 RNNState rnn[2];
132 AVTXContext *tx, *txi;
133 av_tx_fn tx_fn, txi_fn;
134 } DenoiseState;
135
136 typedef struct AudioRNNContext {
137 const AVClass *class;
138
139 char *model_name;
140 float mix;
141
142 int channels;
143 DenoiseState *st;
144
145 DECLARE_ALIGNED(32, float, window)[WINDOW_SIZE];
146 DECLARE_ALIGNED(32, float, dct_table)[FFALIGN(NB_BANDS, 4)][FFALIGN(NB_BANDS, 4)];
147
148 RNNModel *model[2];
149
150 AVFloatDSPContext *fdsp;
151 } AudioRNNContext;
152
153 #define F_ACTIVATION_TANH 0
154 #define F_ACTIVATION_SIGMOID 1
155 #define F_ACTIVATION_RELU 2
156
157 ✗ static void rnnoise_model_free(RNNModel *model)
158 {
159 #define FREE_MAYBE(ptr) do { if (ptr) free(ptr); } while (0)
160 #define FREE_DENSE(name) do { \
161 if (model->name) { \
162 av_free((void *) model->name->input_weights); \
163 av_free((void *) model->name->bias); \
164 av_free((void *) model->name); \
165 } \
166 } while (0)
167 #define FREE_GRU(name) do { \
168 if (model->name) { \
169 av_free((void *) model->name->input_weights); \
170 av_free((void *) model->name->recurrent_weights); \
171 av_free((void *) model->name->bias); \
172 av_free((void *) model->name); \
173 } \
174 } while (0)
175
176 ✗ if (!model)
177 ✗ return;
178 ✗ FREE_DENSE(input_dense);
179 ✗ FREE_GRU(vad_gru);
180 ✗ FREE_GRU(noise_gru);
181 ✗ FREE_GRU(denoise_gru);
182 ✗ FREE_DENSE(denoise_output);
183 ✗ FREE_DENSE(vad_output);
184 ✗ av_free(model);
185 }
186
187 ✗ static int rnnoise_model_from_file(FILE *f, RNNModel **rnn)
188 {
189 ✗ RNNModel *ret = NULL;
190 DenseLayer *input_dense;
191 GRULayer *vad_gru;
192 GRULayer *noise_gru;
193 GRULayer *denoise_gru;
194 DenseLayer *denoise_output;
195 DenseLayer *vad_output;
196 int in;
197
198 ✗ if (fscanf(f, "rnnoise-nu model file version %d\n", &in) != 1 || in != 1)
199 ✗ return AVERROR_INVALIDDATA;
200
201 ✗ ret = av_calloc(1, sizeof(RNNModel));
202 ✗ if (!ret)
203 ✗ return AVERROR(ENOMEM);
204
205 #define ALLOC_LAYER(type, name) \
206 name = av_calloc(1, sizeof(type)); \
207 if (!name) { \
208 rnnoise_model_free(ret); \
209 return AVERROR(ENOMEM); \
210 } \
211 ret->name = name
212
213 ✗ ALLOC_LAYER(DenseLayer, input_dense);
214 ✗ ALLOC_LAYER(GRULayer, vad_gru);
215 ✗ ALLOC_LAYER(GRULayer, noise_gru);
216 ✗ ALLOC_LAYER(GRULayer, denoise_gru);
217 ✗ ALLOC_LAYER(DenseLayer, denoise_output);
218 ✗ ALLOC_LAYER(DenseLayer, vad_output);
219
220 #define INPUT_VAL(name) do { \
221 if (fscanf(f, "%d", &in) != 1 || in < 0 || in > 128) { \
222 rnnoise_model_free(ret); \
223 return AVERROR(EINVAL); \
224 } \
225 name = in; \
226 } while (0)
227
228 #define INPUT_ACTIVATION(name) do { \
229 int activation; \
230 INPUT_VAL(activation); \
231 switch (activation) { \
232 case F_ACTIVATION_SIGMOID: \
233 name = ACTIVATION_SIGMOID; \
234 break; \
235 case F_ACTIVATION_RELU: \
236 name = ACTIVATION_RELU; \
237 break; \
238 default: \
239 name = ACTIVATION_TANH; \
240 } \
241 } while (0)
242
243 #define INPUT_ARRAY(name, len) do { \
244 float *values = av_calloc((len), sizeof(float)); \
245 if (!values) { \
246 rnnoise_model_free(ret); \
247 return AVERROR(ENOMEM); \
248 } \
249 name = values; \
250 for (int i = 0; i < (len); i++) { \
251 if (fscanf(f, "%d", &in) != 1) { \
252 rnnoise_model_free(ret); \
253 return AVERROR(EINVAL); \
254 } \
255 values[i] = in; \
256 } \
257 } while (0)
258
259 #define INPUT_ARRAY3(name, len0, len1, len2) do { \
260 float *values = av_calloc(FFALIGN((len0), 4) * FFALIGN((len1), 4) * (len2), sizeof(float)); \
261 if (!values) { \
262 rnnoise_model_free(ret); \
263 return AVERROR(ENOMEM); \
264 } \
265 name = values; \
266 for (int k = 0; k < (len0); k++) { \
267 for (int i = 0; i < (len2); i++) { \
268 for (int j = 0; j < (len1); j++) { \
269 if (fscanf(f, "%d", &in) != 1) { \
270 rnnoise_model_free(ret); \
271 return AVERROR(EINVAL); \
272 } \
273 values[j * (len2) * FFALIGN((len0), 4) + i * FFALIGN((len0), 4) + k] = in; \
274 } \
275 } \
276 } \
277 } while (0)
278
279 #define NEW_LINE() do { \
280 int c; \
281 while ((c = fgetc(f)) != EOF) { \
282 if (c == '\n') \
283 break; \
284 } \
285 } while (0)
286
287 #define INPUT_DENSE(name) do { \
288 INPUT_VAL(name->nb_inputs); \
289 INPUT_VAL(name->nb_neurons); \
290 ret->name ## _size = name->nb_neurons; \
291 INPUT_ACTIVATION(name->activation); \
292 NEW_LINE(); \
293 INPUT_ARRAY(name->input_weights, name->nb_inputs * name->nb_neurons); \
294 NEW_LINE(); \
295 INPUT_ARRAY(name->bias, name->nb_neurons); \
296 NEW_LINE(); \
297 } while (0)
298
299 #define INPUT_GRU(name) do { \
300 INPUT_VAL(name->nb_inputs); \
301 INPUT_VAL(name->nb_neurons); \
302 ret->name ## _size = name->nb_neurons; \
303 INPUT_ACTIVATION(name->activation); \
304 NEW_LINE(); \
305 INPUT_ARRAY3(name->input_weights, name->nb_inputs, name->nb_neurons, 3); \
306 NEW_LINE(); \
307 INPUT_ARRAY3(name->recurrent_weights, name->nb_neurons, name->nb_neurons, 3); \
308 NEW_LINE(); \
309 INPUT_ARRAY(name->bias, name->nb_neurons * 3); \
310 NEW_LINE(); \
311 } while (0)
312
313 ✗ INPUT_DENSE(input_dense);
314 ✗ INPUT_GRU(vad_gru);
315 ✗ INPUT_GRU(noise_gru);
316 ✗ INPUT_GRU(denoise_gru);
317 ✗ INPUT_DENSE(denoise_output);
318 ✗ INPUT_DENSE(vad_output);
319
320 ✗ if (vad_output->nb_neurons != 1) {
321 ✗ rnnoise_model_free(ret);
322 ✗ return AVERROR(EINVAL);
323 }
324
325 ✗ *rnn = ret;
326
327 ✗ return 0;
328 }
329
330 ✗ static int query_formats(const AVFilterContext *ctx,
331 AVFilterFormatsConfig **cfg_in,
332 AVFilterFormatsConfig **cfg_out)
333 {
334 static const enum AVSampleFormat sample_fmts[] = {
335 AV_SAMPLE_FMT_FLTP,
336 AV_SAMPLE_FMT_NONE
337 };
338 ✗ int ret, sample_rates[] = { 48000, -1 };
339
340 ✗ ret = ff_set_sample_formats_from_list2(ctx, cfg_in, cfg_out, sample_fmts);
341 ✗ if (ret < 0)
342 ✗ return ret;
343
344 ✗ return ff_set_common_samplerates_from_list2(ctx, cfg_in, cfg_out, sample_rates);
345 }
346
347 ✗ static int config_input(AVFilterLink *inlink)
348 {
349 ✗ AVFilterContext *ctx = inlink->dst;
350 ✗ AudioRNNContext *s = ctx->priv;
351 ✗ int ret = 0;
352
353 ✗ s->channels = inlink->ch_layout.nb_channels;
354
355 ✗ if (!s->st)
356 ✗ s->st = av_calloc(s->channels, sizeof(DenoiseState));
357 ✗ if (!s->st)
358 ✗ return AVERROR(ENOMEM);
359
360 ✗ for (int i = 0; i < s->channels; i++) {
361 ✗ DenoiseState *st = &s->st[i];
362
363 ✗ st->rnn[0].model = s->model[0];
364 ✗ st->rnn[0].vad_gru_state = av_calloc(FFALIGN(s->model[0]->vad_gru_size, 16), sizeof(float));
365 ✗ st->rnn[0].noise_gru_state = av_calloc(FFALIGN(s->model[0]->noise_gru_size, 16), sizeof(float));
366 ✗ st->rnn[0].denoise_gru_state = av_calloc(FFALIGN(s->model[0]->denoise_gru_size, 16), sizeof(float));
367 ✗ if (!st->rnn[0].vad_gru_state ||
368 ✗ !st->rnn[0].noise_gru_state ||
369 ✗ !st->rnn[0].denoise_gru_state)
370 ✗ return AVERROR(ENOMEM);
371 }
372
373 ✗ for (int i = 0; i < s->channels; i++) {
374 ✗ DenoiseState *st = &s->st[i];
375 ✗ float scale = 1.f;
376
377 ✗ if (!st->tx)
378 ✗ ret = av_tx_init(&st->tx, &st->tx_fn, AV_TX_FLOAT_FFT, 0, WINDOW_SIZE, &scale, 0);
379 ✗ if (ret < 0)
380 ✗ return ret;
381
382 ✗ if (!st->txi)
383 ✗ ret = av_tx_init(&st->txi, &st->txi_fn, AV_TX_FLOAT_FFT, 1, WINDOW_SIZE, &scale, 0);
384 ✗ if (ret < 0)
385 ✗ return ret;
386 }
387
388 ✗ return ret;
389 }
390
391 ✗ static void biquad(float *y, float mem[2], const float *x,
392 const float *b, const float *a, int N)
393 {
394 ✗ for (int i = 0; i < N; i++) {
395 float xi, yi;
396
397 ✗ xi = x[i];
398 ✗ yi = x[i] + mem[0];
399 ✗ mem[0] = mem[1] + (b[0]*xi - a[0]*yi);
400 ✗ mem[1] = (b[1]*xi - a[1]*yi);
401 ✗ y[i] = yi;
402 }
403 ✗ }
404
405 #define RNN_MOVE(dst, src, n) (memmove((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) ))
406 #define RNN_CLEAR(dst, n) (memset((dst), 0, (n)*sizeof(*(dst))))
407 #define RNN_COPY(dst, src, n) (memcpy((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) ))
408
409 ✗ static void forward_transform(DenoiseState *st, AVComplexFloat *out, const float *in)
410 {
411 AVComplexFloat x[WINDOW_SIZE];
412 AVComplexFloat y[WINDOW_SIZE];
413
414 ✗ for (int i = 0; i < WINDOW_SIZE; i++) {
415 ✗ x[i].re = in[i];
416 ✗ x[i].im = 0;
417 }
418
419 ✗ st->tx_fn(st->tx, y, x, sizeof(AVComplexFloat));
420
421 ✗ RNN_COPY(out, y, FREQ_SIZE);
422 ✗ }
423
424 ✗ static void inverse_transform(DenoiseState *st, float *out, const AVComplexFloat *in)
425 {
426 AVComplexFloat x[WINDOW_SIZE];
427 AVComplexFloat y[WINDOW_SIZE];
428
429 ✗ RNN_COPY(x, in, FREQ_SIZE);
430
431 ✗ for (int i = FREQ_SIZE; i < WINDOW_SIZE; i++) {
432 ✗ x[i].re = x[WINDOW_SIZE - i].re;
433 ✗ x[i].im = -x[WINDOW_SIZE - i].im;
434 }
435
436 ✗ st->txi_fn(st->txi, y, x, sizeof(AVComplexFloat));
437
438 ✗ for (int i = 0; i < WINDOW_SIZE; i++)
439 ✗ out[i] = y[i].re / WINDOW_SIZE;
440 ✗ }
441
442 static const uint8_t eband5ms[] = {
443 /*0 200 400 600 800 1k 1.2 1.4 1.6 2k 2.4 2.8 3.2 4k 4.8 5.6 6.8 8k 9.6 12k 15.6 20k*/
444 0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28, 34, 40, 48, 60, 78, 100
445 };
446
447 ✗ static void compute_band_energy(float *bandE, const AVComplexFloat *X)
448 {
449 ✗ float sum[NB_BANDS] = {0};
450
451 ✗ for (int i = 0; i < NB_BANDS - 1; i++) {
452 int band_size;
453
454 ✗ band_size = (eband5ms[i + 1] - eband5ms[i]) << FRAME_SIZE_SHIFT;
455 ✗ for (int j = 0; j < band_size; j++) {
456 ✗ float tmp, frac = (float)j / band_size;
457
458 ✗ tmp = SQUARE(X[(eband5ms[i] << FRAME_SIZE_SHIFT) + j].re);
459 ✗ tmp += SQUARE(X[(eband5ms[i] << FRAME_SIZE_SHIFT) + j].im);
460 ✗ sum[i] += (1.f - frac) * tmp;
461 ✗ sum[i + 1] += frac * tmp;
462 }
463 }
464
465 ✗ sum[0] *= 2;
466 ✗ sum[NB_BANDS - 1] *= 2;
467
468 ✗ for (int i = 0; i < NB_BANDS; i++)
469 ✗ bandE[i] = sum[i];
470 ✗ }
471
472 ✗ static void compute_band_corr(float *bandE, const AVComplexFloat *X, const AVComplexFloat *P)
473 {
474 ✗ float sum[NB_BANDS] = { 0 };
475
476 ✗ for (int i = 0; i < NB_BANDS - 1; i++) {
477 int band_size;
478
479 ✗ band_size = (eband5ms[i + 1] - eband5ms[i]) << FRAME_SIZE_SHIFT;
480 ✗ for (int j = 0; j < band_size; j++) {
481 ✗ float tmp, frac = (float)j / band_size;
482
483 ✗ tmp = X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].re * P[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].re;
484 ✗ tmp += X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].im * P[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].im;
485 ✗ sum[i] += (1 - frac) * tmp;
486 ✗ sum[i + 1] += frac * tmp;
487 }
488 }
489
490 ✗ sum[0] *= 2;
491 ✗ sum[NB_BANDS-1] *= 2;
492
493 ✗ for (int i = 0; i < NB_BANDS; i++)
494 ✗ bandE[i] = sum[i];
495 ✗ }
496
497 ✗ static void frame_analysis(AudioRNNContext *s, DenoiseState *st, AVComplexFloat *X, float *Ex, const float *in)
498 {
499 ✗ LOCAL_ALIGNED_32(float, x, [WINDOW_SIZE]);
500
501 ✗ RNN_COPY(x, st->analysis_mem, FRAME_SIZE);
502 ✗ RNN_COPY(x + FRAME_SIZE, in, FRAME_SIZE);
503 ✗ RNN_COPY(st->analysis_mem, in, FRAME_SIZE);
504 ✗ s->fdsp->vector_fmul(x, x, s->window, WINDOW_SIZE);
505 ✗ forward_transform(st, X, x);
506 ✗ compute_band_energy(Ex, X);
507 ✗ }
508
509 ✗ static void frame_synthesis(AudioRNNContext *s, DenoiseState *st, float *out, const AVComplexFloat *y)
510 {
511 ✗ LOCAL_ALIGNED_32(float, x, [WINDOW_SIZE]);
512 ✗ const float *src = st->history;
513 ✗ const float mix = s->mix;
514 ✗ const float imix = 1.f - FFMAX(mix, 0.f);
515
516 ✗ inverse_transform(st, x, y);
517 ✗ s->fdsp->vector_fmul(x, x, s->window, WINDOW_SIZE);
518 ✗ s->fdsp->vector_fmac_scalar(x, st->synthesis_mem, 1.f, FRAME_SIZE);
519 ✗ RNN_COPY(out, x, FRAME_SIZE);
520 ✗ RNN_COPY(st->synthesis_mem, &x[FRAME_SIZE], FRAME_SIZE);
521
522 ✗ for (int n = 0; n < FRAME_SIZE; n++)
523 ✗ out[n] = out[n] * mix + src[n] * imix;
524 ✗ }
525
526 ✗ static inline void xcorr_kernel(const float *x, const float *y, float sum[4], int len)
527 {
528 ✗ float y_0, y_1, y_2, y_3 = 0;
529 int j;
530
531 ✗ y_0 = *y++;
532 ✗ y_1 = *y++;
533 ✗ y_2 = *y++;
534
535 ✗ for (j = 0; j < len - 3; j += 4) {
536 float tmp;
537
538 ✗ tmp = *x++;
539 ✗ y_3 = *y++;
540 ✗ sum[0] += tmp * y_0;
541 ✗ sum[1] += tmp * y_1;
542 ✗ sum[2] += tmp * y_2;
543 ✗ sum[3] += tmp * y_3;
544 ✗ tmp = *x++;
545 ✗ y_0 = *y++;
546 ✗ sum[0] += tmp * y_1;
547 ✗ sum[1] += tmp * y_2;
548 ✗ sum[2] += tmp * y_3;
549 ✗ sum[3] += tmp * y_0;
550 ✗ tmp = *x++;
551 ✗ y_1 = *y++;
552 ✗ sum[0] += tmp * y_2;
553 ✗ sum[1] += tmp * y_3;
554 ✗ sum[2] += tmp * y_0;
555 ✗ sum[3] += tmp * y_1;
556 ✗ tmp = *x++;
557 ✗ y_2 = *y++;
558 ✗ sum[0] += tmp * y_3;
559 ✗ sum[1] += tmp * y_0;
560 ✗ sum[2] += tmp * y_1;
561 ✗ sum[3] += tmp * y_2;
562 }
563
564 ✗ if (j++ < len) {
565 ✗ float tmp = *x++;
566
567 ✗ y_3 = *y++;
568 ✗ sum[0] += tmp * y_0;
569 ✗ sum[1] += tmp * y_1;
570 ✗ sum[2] += tmp * y_2;
571 ✗ sum[3] += tmp * y_3;
572 }
573
574 ✗ if (j++ < len) {
575 ✗ float tmp=*x++;
576
577 ✗ y_0 = *y++;
578 ✗ sum[0] += tmp * y_1;
579 ✗ sum[1] += tmp * y_2;
580 ✗ sum[2] += tmp * y_3;
581 ✗ sum[3] += tmp * y_0;
582 }
583
584 ✗ if (j < len) {
585 ✗ float tmp=*x++;
586
587 ✗ y_1 = *y++;
588 ✗ sum[0] += tmp * y_2;
589 ✗ sum[1] += tmp * y_3;
590 ✗ sum[2] += tmp * y_0;
591 ✗ sum[3] += tmp * y_1;
592 }
593 ✗ }
594
595 ✗ static inline float celt_inner_prod(const float *x,
596 const float *y, int N)
597 {
598 ✗ float xy = 0.f;
599
600 ✗ for (int i = 0; i < N; i++)
601 ✗ xy += x[i] * y[i];
602
603 ✗ return xy;
604 }
605
606 ✗ static void celt_pitch_xcorr(const float *x, const float *y,
607 float *xcorr, int len, int max_pitch)
608 {
609 int i;
610
611 ✗ for (i = 0; i < max_pitch - 3; i += 4) {
612 ✗ float sum[4] = { 0, 0, 0, 0};
613
614 ✗ xcorr_kernel(x, y + i, sum, len);
615
616 ✗ xcorr[i] = sum[0];
617 ✗ xcorr[i + 1] = sum[1];
618 ✗ xcorr[i + 2] = sum[2];
619 ✗ xcorr[i + 3] = sum[3];
620 }
621 /* In case max_pitch isn't a multiple of 4, do non-unrolled version. */
622 ✗ for (; i < max_pitch; i++) {
623 ✗ xcorr[i] = celt_inner_prod(x, y + i, len);
624 }
625 ✗ }
626
627 ✗ static int celt_autocorr(const float *x, /* in: [0...n-1] samples x */
628 float *ac, /* out: [0...lag-1] ac values */
629 const float *window,
630 int overlap,
631 int lag,
632 int n)
633 {
634 ✗ int fastN = n - lag;
635 int shift;
636 const float *xptr;
637 float xx[PITCH_BUF_SIZE>>1];
638
639 ✗ if (overlap == 0) {
640 ✗ xptr = x;
641 } else {
642 ✗ for (int i = 0; i < n; i++)
643 ✗ xx[i] = x[i];
644 ✗ for (int i = 0; i < overlap; i++) {
645 ✗ xx[i] = x[i] * window[i];
646 ✗ xx[n-i-1] = x[n-i-1] * window[i];
647 }
648 ✗ xptr = xx;
649 }
650
651 ✗ shift = 0;
652 ✗ celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1);
653
654 ✗ for (int k = 0; k <= lag; k++) {
655 ✗ float d = 0.f;
656
657 ✗ for (int i = k + fastN; i < n; i++)
658 ✗ d += xptr[i] * xptr[i-k];
659 ✗ ac[k] += d;
660 }
661
662 ✗ return shift;
663 }
664
665 ✗ static void celt_lpc(float *lpc, /* out: [0...p-1] LPC coefficients */
666 const float *ac, /* in: [0...p] autocorrelation values */
667 int p)
668 {
669 ✗ float r, error = ac[0];
670
671 ✗ RNN_CLEAR(lpc, p);
672 ✗ if (ac[0] != 0) {
673 ✗ for (int i = 0; i < p; i++) {
674 /* Sum up this iteration's reflection coefficient */
675 ✗ float rr = 0;
676 ✗ for (int j = 0; j < i; j++)
677 ✗ rr += (lpc[j] * ac[i - j]);
678 ✗ rr += ac[i + 1];
679 ✗ r = -rr/error;
680 /* Update LPC coefficients and total error */
681 ✗ lpc[i] = r;
682 ✗ for (int j = 0; j < (i + 1) >> 1; j++) {
683 float tmp1, tmp2;
684 ✗ tmp1 = lpc[j];
685 ✗ tmp2 = lpc[i-1-j];
686 ✗ lpc[j] = tmp1 + (r*tmp2);
687 ✗ lpc[i-1-j] = tmp2 + (r*tmp1);
688 }
689
690 ✗ error = error - (r * r *error);
691 /* Bail out once we get 30 dB gain */
692 ✗ if (error < .001f * ac[0])
693 ✗ break;
694 }
695 }
696 ✗ }
697
698 ✗ static void celt_fir5(const float *x,
699 const float *num,
700 float *y,
701 int N,
702 float *mem)
703 {
704 float num0, num1, num2, num3, num4;
705 float mem0, mem1, mem2, mem3, mem4;
706
707 ✗ num0 = num[0];
708 ✗ num1 = num[1];
709 ✗ num2 = num[2];
710 ✗ num3 = num[3];
711 ✗ num4 = num[4];
712 ✗ mem0 = mem[0];
713 ✗ mem1 = mem[1];
714 ✗ mem2 = mem[2];
715 ✗ mem3 = mem[3];
716 ✗ mem4 = mem[4];
717
718 ✗ for (int i = 0; i < N; i++) {
719 ✗ float sum = x[i];
720
721 ✗ sum += (num0*mem0);
722 ✗ sum += (num1*mem1);
723 ✗ sum += (num2*mem2);
724 ✗ sum += (num3*mem3);
725 ✗ sum += (num4*mem4);
726 ✗ mem4 = mem3;
727 ✗ mem3 = mem2;
728 ✗ mem2 = mem1;
729 ✗ mem1 = mem0;
730 ✗ mem0 = x[i];
731 ✗ y[i] = sum;
732 }
733
734 ✗ mem[0] = mem0;
735 ✗ mem[1] = mem1;
736 ✗ mem[2] = mem2;
737 ✗ mem[3] = mem3;
738 ✗ mem[4] = mem4;
739 ✗ }
740
741 ✗ static void pitch_downsample(float *x[], float *x_lp,
742 int len, int C)
743 {
744 float ac[5];
745 ✗ float tmp=Q15ONE;
746 ✗ float lpc[4], mem[5]={0,0,0,0,0};
747 float lpc2[5];
748 ✗ float c1 = .8f;
749
750 ✗ for (int i = 1; i < len >> 1; i++)
751 ✗ x_lp[i] = .5f * (.5f * (x[0][(2*i-1)]+x[0][(2*i+1)])+x[0][2*i]);
752 ✗ x_lp[0] = .5f * (.5f * (x[0][1])+x[0][0]);
753 ✗ if (C==2) {
754 ✗ for (int i = 1; i < len >> 1; i++)
755 ✗ x_lp[i] += (.5f * (.5f * (x[1][(2*i-1)]+x[1][(2*i+1)])+x[1][2*i]));
756 ✗ x_lp[0] += .5f * (.5f * (x[1][1])+x[1][0]);
757 }
758
759 ✗ celt_autocorr(x_lp, ac, NULL, 0, 4, len>>1);
760
761 /* Noise floor -40 dB */
762 ✗ ac[0] *= 1.0001f;
763 /* Lag windowing */
764 ✗ for (int i = 1; i <= 4; i++) {
765 /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
766 ✗ ac[i] -= ac[i]*(.008f*i)*(.008f*i);
767 }
768
769 ✗ celt_lpc(lpc, ac, 4);
770 ✗ for (int i = 0; i < 4; i++) {
771 ✗ tmp = .9f * tmp;
772 ✗ lpc[i] = (lpc[i] * tmp);
773 }
774 /* Add a zero */
775 ✗ lpc2[0] = lpc[0] + .8f;
776 ✗ lpc2[1] = lpc[1] + (c1 * lpc[0]);
777 ✗ lpc2[2] = lpc[2] + (c1 * lpc[1]);
778 ✗ lpc2[3] = lpc[3] + (c1 * lpc[2]);
779 ✗ lpc2[4] = (c1 * lpc[3]);
780 ✗ celt_fir5(x_lp, lpc2, x_lp, len>>1, mem);
781 ✗ }
782
783 ✗ static inline void dual_inner_prod(const float *x, const float *y01, const float *y02,
784 int N, float *xy1, float *xy2)
785 {
786 ✗ float xy01 = 0, xy02 = 0;
787
788 ✗ for (int i = 0; i < N; i++) {
789 ✗ xy01 += (x[i] * y01[i]);
790 ✗ xy02 += (x[i] * y02[i]);
791 }
792
793 ✗ *xy1 = xy01;
794 ✗ *xy2 = xy02;
795 ✗ }
796
797 ✗ static float compute_pitch_gain(float xy, float xx, float yy)
798 {
799 ✗ return xy / sqrtf(1.f + xx * yy);
800 }
801
802 static const uint8_t second_check[16] = {0, 0, 3, 2, 3, 2, 5, 2, 3, 2, 3, 2, 5, 2, 3, 2};
803 ✗ static float remove_doubling(float *x, int maxperiod, int minperiod, int N,
804 int *T0_, int prev_period, float prev_gain)
805 {
806 int k, i, T, T0;
807 float g, g0;
808 float pg;
809 float xy,xx,yy,xy2;
810 float xcorr[3];
811 float best_xy, best_yy;
812 int offset;
813 int minperiod0;
814 float yy_lookup[PITCH_MAX_PERIOD+1];
815
816 ✗ minperiod0 = minperiod;
817 ✗ maxperiod /= 2;
818 ✗ minperiod /= 2;
819 ✗ *T0_ /= 2;
820 ✗ prev_period /= 2;
821 ✗ N /= 2;
822 ✗ x += maxperiod;
823 ✗ if (*T0_>=maxperiod)
824 ✗ *T0_=maxperiod-1;
825
826 ✗ T = T0 = *T0_;
827 ✗ dual_inner_prod(x, x, x-T0, N, &xx, &xy);
828 ✗ yy_lookup[0] = xx;
829 ✗ yy=xx;
830 ✗ for (i = 1; i <= maxperiod; i++) {
831 ✗ yy = yy+(x[-i] * x[-i])-(x[N-i] * x[N-i]);
832 ✗ yy_lookup[i] = FFMAX(0, yy);
833 }
834 ✗ yy = yy_lookup[T0];
835 ✗ best_xy = xy;
836 ✗ best_yy = yy;
837 ✗ g = g0 = compute_pitch_gain(xy, xx, yy);
838 /* Look for any pitch at T/k */
839 ✗ for (k = 2; k <= 15; k++) {
840 int T1, T1b;
841 float g1;
842 ✗ float cont=0;
843 float thresh;
844 ✗ T1 = (2*T0+k)/(2*k);
845 ✗ if (T1 < minperiod)
846 ✗ break;
847 /* Look for another strong correlation at T1b */
848 ✗ if (k==2)
849 {
850 ✗ if (T1+T0>maxperiod)
851 ✗ T1b = T0;
852 else
853 ✗ T1b = T0+T1;
854 } else
855 {
856 ✗ T1b = (2*second_check[k]*T0+k)/(2*k);
857 }
858 ✗ dual_inner_prod(x, &x[-T1], &x[-T1b], N, &xy, &xy2);
859 ✗ xy = .5f * (xy + xy2);
860 ✗ yy = .5f * (yy_lookup[T1] + yy_lookup[T1b]);
861 ✗ g1 = compute_pitch_gain(xy, xx, yy);
862 ✗ if (FFABS(T1-prev_period)<=1)
863 ✗ cont = prev_gain;
864 ✗ else if (FFABS(T1-prev_period)<=2 && 5 * k * k < T0)
865 ✗ cont = prev_gain * .5f;
866 else
867 ✗ cont = 0;
868 ✗ thresh = FFMAX(.3f, (.7f * g0) - cont);
869 /* Bias against very high pitch (very short period) to avoid false-positives
870 due to short-term correlation */
871 ✗ if (T1<3*minperiod)
872 ✗ thresh = FFMAX(.4f, (.85f * g0) - cont);
873 ✗ else if (T1<2*minperiod)
874 ✗ thresh = FFMAX(.5f, (.9f * g0) - cont);
875 ✗ if (g1 > thresh)
876 {
877 ✗ best_xy = xy;
878 ✗ best_yy = yy;
879 ✗ T = T1;
880 ✗ g = g1;
881 }
882 }
883 ✗ best_xy = FFMAX(0, best_xy);
884 ✗ if (best_yy <= best_xy)
885 ✗ pg = Q15ONE;
886 else
887 ✗ pg = best_xy/(best_yy + 1);
888
889 ✗ for (k = 0; k < 3; k++)
890 ✗ xcorr[k] = celt_inner_prod(x, x-(T+k-1), N);
891 ✗ if ((xcorr[2]-xcorr[0]) > .7f * (xcorr[1]-xcorr[0]))
892 ✗ offset = 1;
893 ✗ else if ((xcorr[0]-xcorr[2]) > (.7f * (xcorr[1] - xcorr[2])))
894 ✗ offset = -1;
895 else
896 ✗ offset = 0;
897 ✗ if (pg > g)
898 ✗ pg = g;
899 ✗ *T0_ = 2*T+offset;
900
901 ✗ if (*T0_<minperiod0)
902 ✗ *T0_=minperiod0;
903 ✗ return pg;
904 }
905
906 ✗ static void find_best_pitch(float *xcorr, float *y, int len,
907 int max_pitch, int *best_pitch)
908 {
909 float best_num[2];
910 float best_den[2];
911 ✗ float Syy = 1.f;
912
913 ✗ best_num[0] = -1;
914 ✗ best_num[1] = -1;
915 ✗ best_den[0] = 0;
916 ✗ best_den[1] = 0;
917 ✗ best_pitch[0] = 0;
918 ✗ best_pitch[1] = 1;
919
920 ✗ for (int j = 0; j < len; j++)
921 ✗ Syy += y[j] * y[j];
922
923 ✗ for (int i = 0; i < max_pitch; i++) {
924 ✗ if (xcorr[i]>0) {
925 float num;
926 float xcorr16;
927
928 ✗ xcorr16 = xcorr[i];
929 /* Considering the range of xcorr16, this should avoid both underflows
930 and overflows (inf) when squaring xcorr16 */
931 ✗ xcorr16 *= 1e-12f;
932 ✗ num = xcorr16 * xcorr16;
933 ✗ if ((num * best_den[1]) > (best_num[1] * Syy)) {
934 ✗ if ((num * best_den[0]) > (best_num[0] * Syy)) {
935 ✗ best_num[1] = best_num[0];
936 ✗ best_den[1] = best_den[0];
937 ✗ best_pitch[1] = best_pitch[0];
938 ✗ best_num[0] = num;
939 ✗ best_den[0] = Syy;
940 ✗ best_pitch[0] = i;
941 } else {
942 ✗ best_num[1] = num;
943 ✗ best_den[1] = Syy;
944 ✗ best_pitch[1] = i;
945 }
946 }
947 }
948 ✗ Syy += y[i+len]*y[i+len] - y[i] * y[i];
949 ✗ Syy = FFMAX(1, Syy);
950 }
951 ✗ }
952
953 ✗ static void pitch_search(const float *x_lp, float *y,
954 int len, int max_pitch, int *pitch)
955 {
956 int lag;
957 ✗ int best_pitch[2]={0,0};
958 int offset;
959
960 float x_lp4[WINDOW_SIZE];
961 float y_lp4[WINDOW_SIZE];
962 float xcorr[WINDOW_SIZE];
963
964 ✗ lag = len+max_pitch;
965
966 /* Downsample by 2 again */
967 ✗ for (int j = 0; j < len >> 2; j++)
968 ✗ x_lp4[j] = x_lp[2*j];
969 ✗ for (int j = 0; j < lag >> 2; j++)
970 ✗ y_lp4[j] = y[2*j];
971
972 /* Coarse search with 4x decimation */
973
974 ✗ celt_pitch_xcorr(x_lp4, y_lp4, xcorr, len>>2, max_pitch>>2);
975
976 ✗ find_best_pitch(xcorr, y_lp4, len>>2, max_pitch>>2, best_pitch);
977
978 /* Finer search with 2x decimation */
979 ✗ for (int i = 0; i < max_pitch >> 1; i++) {
980 float sum;
981 ✗ xcorr[i] = 0;
982 ✗ if (FFABS(i-2*best_pitch[0])>2 && FFABS(i-2*best_pitch[1])>2)
983 ✗ continue;
984 ✗ sum = celt_inner_prod(x_lp, y+i, len>>1);
985 ✗ xcorr[i] = FFMAX(-1, sum);
986 }
987
988 ✗ find_best_pitch(xcorr, y, len>>1, max_pitch>>1, best_pitch);
989
990 /* Refine by pseudo-interpolation */
991 ✗ if (best_pitch[0] > 0 && best_pitch[0] < (max_pitch >> 1) - 1) {
992 float a, b, c;
993
994 ✗ a = xcorr[best_pitch[0] - 1];
995 ✗ b = xcorr[best_pitch[0]];
996 ✗ c = xcorr[best_pitch[0] + 1];
997 ✗ if (c - a > .7f * (b - a))
998 ✗ offset = 1;
999 ✗ else if (a - c > .7f * (b-c))
1000 ✗ offset = -1;
1001 else
1002 ✗ offset = 0;
1003 } else {
1004 ✗ offset = 0;
1005 }
1006
1007 ✗ *pitch = 2 * best_pitch[0] - offset;
1008 ✗ }
1009
1010 ✗ static void dct(AudioRNNContext *s, float *out, const float *in)
1011 {
1012 ✗ for (int i = 0; i < NB_BANDS; i++) {
1013 float sum;
1014
1015 ✗ sum = s->fdsp->scalarproduct_float(in, s->dct_table[i], FFALIGN(NB_BANDS, 4));
1016 ✗ out[i] = sum * sqrtf(2.f / 22);
1017 }
1018 ✗ }
1019
1020 ✗ static int compute_frame_features(AudioRNNContext *s, DenoiseState *st, AVComplexFloat *X, AVComplexFloat *P,
1021 float *Ex, float *Ep, float *Exp, float *features, const float *in)
1022 {
1023 ✗ float E = 0;
1024 float *ceps_0, *ceps_1, *ceps_2;
1025 ✗ float spec_variability = 0;
1026 ✗ LOCAL_ALIGNED_32(float, Ly, [FFALIGN(NB_BANDS, 4)]);
1027 ✗ LOCAL_ALIGNED_32(float, p, [WINDOW_SIZE]);
1028 float pitch_buf[PITCH_BUF_SIZE>>1];
1029 int pitch_index;
1030 float gain;
1031 float *(pre[1]);
1032 float tmp[NB_BANDS];
1033 float follow, logMax;
1034
1035 ✗ frame_analysis(s, st, X, Ex, in);
1036 ✗ RNN_MOVE(st->pitch_buf, &st->pitch_buf[FRAME_SIZE], PITCH_BUF_SIZE-FRAME_SIZE);
1037 ✗ RNN_COPY(&st->pitch_buf[PITCH_BUF_SIZE-FRAME_SIZE], in, FRAME_SIZE);
1038 ✗ pre[0] = &st->pitch_buf[0];
1039 ✗ pitch_downsample(pre, pitch_buf, PITCH_BUF_SIZE, 1);
1040 ✗ pitch_search(pitch_buf+(PITCH_MAX_PERIOD>>1), pitch_buf, PITCH_FRAME_SIZE,
1041 PITCH_MAX_PERIOD-3*PITCH_MIN_PERIOD, &pitch_index);
1042 ✗ pitch_index = PITCH_MAX_PERIOD-pitch_index;
1043
1044 ✗ gain = remove_doubling(pitch_buf, PITCH_MAX_PERIOD, PITCH_MIN_PERIOD,
1045 PITCH_FRAME_SIZE, &pitch_index, st->last_period, st->last_gain);
1046 ✗ st->last_period = pitch_index;
1047 ✗ st->last_gain = gain;
1048
1049 ✗ for (int i = 0; i < WINDOW_SIZE; i++)
1050 ✗ p[i] = st->pitch_buf[PITCH_BUF_SIZE-WINDOW_SIZE-pitch_index+i];
1051
1052 ✗ s->fdsp->vector_fmul(p, p, s->window, WINDOW_SIZE);
1053 ✗ forward_transform(st, P, p);
1054 ✗ compute_band_energy(Ep, P);
1055 ✗ compute_band_corr(Exp, X, P);
1056
1057 ✗ for (int i = 0; i < NB_BANDS; i++)
1058 ✗ Exp[i] = Exp[i] / sqrtf(.001f+Ex[i]*Ep[i]);
1059
1060 ✗ memset(Exp + NB_BANDS, 0, (FFALIGN(NB_BANDS, 4) - NB_BANDS) * sizeof(*Exp));
1061 ✗ dct(s, tmp, Exp);
1062
1063 ✗ for (int i = 0; i < NB_DELTA_CEPS; i++)
1064 ✗ features[NB_BANDS+2*NB_DELTA_CEPS+i] = tmp[i];
1065
1066 ✗ features[NB_BANDS+2*NB_DELTA_CEPS] -= 1.3;
1067 ✗ features[NB_BANDS+2*NB_DELTA_CEPS+1] -= 0.9;
1068 ✗ features[NB_BANDS+3*NB_DELTA_CEPS] = .01*(pitch_index-300);
1069 ✗ logMax = -2;
1070 ✗ follow = -2;
1071
1072 ✗ for (int i = 0; i < NB_BANDS; i++) {
1073 ✗ Ly[i] = log10f(1e-2f + Ex[i]);
1074 ✗ Ly[i] = FFMAX(logMax-7, FFMAX(follow-1.5, Ly[i]));
1075 ✗ logMax = FFMAX(logMax, Ly[i]);
1076 ✗ follow = FFMAX(follow-1.5, Ly[i]);
1077 ✗ E += Ex[i];
1078 }
1079
1080 ✗ if (E < 0.04f) {
1081 /* If there's no audio, avoid messing up the state. */
1082 ✗ RNN_CLEAR(features, NB_FEATURES);
1083 ✗ return 1;
1084 }
1085
1086 ✗ memset(Ly + NB_BANDS, 0, (FFALIGN(NB_BANDS, 4) - NB_BANDS) * sizeof(*Ly));
1087 ✗ dct(s, features, Ly);
1088 ✗ features[0] -= 12;
1089 ✗ features[1] -= 4;
1090 ✗ ceps_0 = st->cepstral_mem[st->memid];
1091 ✗ ceps_1 = (st->memid < 1) ? st->cepstral_mem[CEPS_MEM+st->memid-1] : st->cepstral_mem[st->memid-1];
1092 ✗ ceps_2 = (st->memid < 2) ? st->cepstral_mem[CEPS_MEM+st->memid-2] : st->cepstral_mem[st->memid-2];
1093
1094 ✗ for (int i = 0; i < NB_BANDS; i++)
1095 ✗ ceps_0[i] = features[i];
1096
1097 ✗ st->memid++;
1098 ✗ for (int i = 0; i < NB_DELTA_CEPS; i++) {
1099 ✗ features[i] = ceps_0[i] + ceps_1[i] + ceps_2[i];
1100 ✗ features[NB_BANDS+i] = ceps_0[i] - ceps_2[i];
1101 ✗ features[NB_BANDS+NB_DELTA_CEPS+i] = ceps_0[i] - 2*ceps_1[i] + ceps_2[i];
1102 }
1103 /* Spectral variability features. */
1104 ✗ if (st->memid == CEPS_MEM)
1105 ✗ st->memid = 0;
1106
1107 ✗ for (int i = 0; i < CEPS_MEM; i++) {
1108 ✗ float mindist = 1e15f;
1109 ✗ for (int j = 0; j < CEPS_MEM; j++) {
1110 ✗ float dist = 0.f;
1111 ✗ for (int k = 0; k < NB_BANDS; k++) {
1112 float tmp;
1113
1114 ✗ tmp = st->cepstral_mem[i][k] - st->cepstral_mem[j][k];
1115 ✗ dist += tmp*tmp;
1116 }
1117
1118 ✗ if (j != i)
1119 ✗ mindist = FFMIN(mindist, dist);
1120 }
1121
1122 ✗ spec_variability += mindist;
1123 }
1124
1125 ✗ features[NB_BANDS+3*NB_DELTA_CEPS+1] = spec_variability/CEPS_MEM-2.1;
1126
1127 ✗ return 0;
1128 }
1129
1130 ✗ static void interp_band_gain(float *g, const float *bandE)
1131 {
1132 ✗ memset(g, 0, sizeof(*g) * FREQ_SIZE);
1133
1134 ✗ for (int i = 0; i < NB_BANDS - 1; i++) {
1135 ✗ const int band_size = (eband5ms[i + 1] - eband5ms[i]) << FRAME_SIZE_SHIFT;
1136
1137 ✗ for (int j = 0; j < band_size; j++) {
1138 ✗ float frac = (float)j / band_size;
1139
1140 ✗ g[(eband5ms[i] << FRAME_SIZE_SHIFT) + j] = (1.f - frac) * bandE[i] + frac * bandE[i + 1];
1141 }
1142 }
1143 ✗ }
1144
1145 ✗ static void pitch_filter(AVComplexFloat *X, const AVComplexFloat *P, const float *Ex, const float *Ep,
1146 const float *Exp, const float *g)
1147 {
1148 float newE[NB_BANDS];
1149 float r[NB_BANDS];
1150 float norm[NB_BANDS];
1151 ✗ float rf[FREQ_SIZE] = {0};
1152 ✗ float normf[FREQ_SIZE]={0};
1153
1154 ✗ for (int i = 0; i < NB_BANDS; i++) {
1155 ✗ if (Exp[i]>g[i]) r[i] = 1;
1156 ✗ else r[i] = SQUARE(Exp[i])*(1-SQUARE(g[i]))/(.001 + SQUARE(g[i])*(1-SQUARE(Exp[i])));
1157 ✗ r[i] = sqrtf(av_clipf(r[i], 0, 1));
1158 ✗ r[i] *= sqrtf(Ex[i]/(1e-8+Ep[i]));
1159 }
1160 ✗ interp_band_gain(rf, r);
1161 ✗ for (int i = 0; i < FREQ_SIZE; i++) {
1162 ✗ X[i].re += rf[i]*P[i].re;
1163 ✗ X[i].im += rf[i]*P[i].im;
1164 }
1165 ✗ compute_band_energy(newE, X);
1166 ✗ for (int i = 0; i < NB_BANDS; i++) {
1167 ✗ norm[i] = sqrtf(Ex[i] / (1e-8+newE[i]));
1168 }
1169 ✗ interp_band_gain(normf, norm);
1170 ✗ for (int i = 0; i < FREQ_SIZE; i++) {
1171 ✗ X[i].re *= normf[i];
1172 ✗ X[i].im *= normf[i];
1173 }
1174 ✗ }
1175
1176 static const float tansig_table[201] = {
1177 0.000000f, 0.039979f, 0.079830f, 0.119427f, 0.158649f,
1178 0.197375f, 0.235496f, 0.272905f, 0.309507f, 0.345214f,
1179 0.379949f, 0.413644f, 0.446244f, 0.477700f, 0.507977f,
1180 0.537050f, 0.564900f, 0.591519f, 0.616909f, 0.641077f,
1181 0.664037f, 0.685809f, 0.706419f, 0.725897f, 0.744277f,
1182 0.761594f, 0.777888f, 0.793199f, 0.807569f, 0.821040f,
1183 0.833655f, 0.845456f, 0.856485f, 0.866784f, 0.876393f,
1184 0.885352f, 0.893698f, 0.901468f, 0.908698f, 0.915420f,
1185 0.921669f, 0.927473f, 0.932862f, 0.937863f, 0.942503f,
1186 0.946806f, 0.950795f, 0.954492f, 0.957917f, 0.961090f,
1187 0.964028f, 0.966747f, 0.969265f, 0.971594f, 0.973749f,
1188 0.975743f, 0.977587f, 0.979293f, 0.980869f, 0.982327f,
1189 0.983675f, 0.984921f, 0.986072f, 0.987136f, 0.988119f,
1190 0.989027f, 0.989867f, 0.990642f, 0.991359f, 0.992020f,
1191 0.992631f, 0.993196f, 0.993718f, 0.994199f, 0.994644f,
1192 0.995055f, 0.995434f, 0.995784f, 0.996108f, 0.996407f,
1193 0.996682f, 0.996937f, 0.997172f, 0.997389f, 0.997590f,
1194 0.997775f, 0.997946f, 0.998104f, 0.998249f, 0.998384f,
1195 0.998508f, 0.998623f, 0.998728f, 0.998826f, 0.998916f,
1196 0.999000f, 0.999076f, 0.999147f, 0.999213f, 0.999273f,
1197 0.999329f, 0.999381f, 0.999428f, 0.999472f, 0.999513f,
1198 0.999550f, 0.999585f, 0.999617f, 0.999646f, 0.999673f,
1199 0.999699f, 0.999722f, 0.999743f, 0.999763f, 0.999781f,
1200 0.999798f, 0.999813f, 0.999828f, 0.999841f, 0.999853f,
1201 0.999865f, 0.999875f, 0.999885f, 0.999893f, 0.999902f,
1202 0.999909f, 0.999916f, 0.999923f, 0.999929f, 0.999934f,
1203 0.999939f, 0.999944f, 0.999948f, 0.999952f, 0.999956f,
1204 0.999959f, 0.999962f, 0.999965f, 0.999968f, 0.999970f,
1205 0.999973f, 0.999975f, 0.999977f, 0.999978f, 0.999980f,
1206 0.999982f, 0.999983f, 0.999984f, 0.999986f, 0.999987f,
1207 0.999988f, 0.999989f, 0.999990f, 0.999990f, 0.999991f,
1208 0.999992f, 0.999992f, 0.999993f, 0.999994f, 0.999994f,
1209 0.999994f, 0.999995f, 0.999995f, 0.999996f, 0.999996f,
1210 0.999996f, 0.999997f, 0.999997f, 0.999997f, 0.999997f,
1211 0.999997f, 0.999998f, 0.999998f, 0.999998f, 0.999998f,
1212 0.999998f, 0.999998f, 0.999999f, 0.999999f, 0.999999f,
1213 0.999999f, 0.999999f, 0.999999f, 0.999999f, 0.999999f,
1214 0.999999f, 0.999999f, 0.999999f, 0.999999f, 0.999999f,
1215 1.000000f, 1.000000f, 1.000000f, 1.000000f, 1.000000f,
1216 1.000000f, 1.000000f, 1.000000f, 1.000000f, 1.000000f,
1217 1.000000f,
1218 };
1219
1220 ✗ static inline float tansig_approx(float x)
1221 {
1222 float y, dy;
1223 ✗ float sign=1;
1224 int i;
1225
1226 /* Tests are reversed to catch NaNs */
1227 ✗ if (!(x<8))
1228 ✗ return 1;
1229 ✗ if (!(x>-8))
1230 ✗ return -1;
1231 /* Another check in case of -ffast-math */
1232
1233 ✗ if (isnan(x))
1234 ✗ return 0;
1235
1236 ✗ if (x < 0) {
1237 ✗ x=-x;
1238 ✗ sign=-1;
1239 }
1240 ✗ i = (int)floor(.5f+25*x);
1241 ✗ x -= .04f*i;
1242 ✗ y = tansig_table[i];
1243 ✗ dy = 1-y*y;
1244 ✗ y = y + x*dy*(1 - y*x);
1245 ✗ return sign*y;
1246 }
1247
1248 ✗ static inline float sigmoid_approx(float x)
1249 {
1250 ✗ return .5f + .5f*tansig_approx(.5f*x);
1251 }
1252
1253 ✗ static void compute_dense(const DenseLayer *layer, float *output, const float *input)
1254 {
1255 ✗ const int N = layer->nb_neurons, M = layer->nb_inputs, stride = N;
1256
1257 ✗ for (int i = 0; i < N; i++) {
1258 /* Compute update gate. */
1259 ✗ float sum = layer->bias[i];
1260
1261 ✗ for (int j = 0; j < M; j++)
1262 ✗ sum += layer->input_weights[j * stride + i] * input[j];
1263
1264 ✗ output[i] = WEIGHTS_SCALE * sum;
1265 }
1266
1267 ✗ if (layer->activation == ACTIVATION_SIGMOID) {
1268 ✗ for (int i = 0; i < N; i++)
1269 ✗ output[i] = sigmoid_approx(output[i]);
1270 ✗ } else if (layer->activation == ACTIVATION_TANH) {
1271 ✗ for (int i = 0; i < N; i++)
1272 ✗ output[i] = tansig_approx(output[i]);
1273 ✗ } else if (layer->activation == ACTIVATION_RELU) {
1274 ✗ for (int i = 0; i < N; i++)
1275 ✗ output[i] = FFMAX(0, output[i]);
1276 } else {
1277 ✗ av_assert0(0);
1278 }
1279 ✗ }
1280
1281 ✗ static void compute_gru(AudioRNNContext *s, const GRULayer *gru, float *state, const float *input)
1282 {
1283 ✗ LOCAL_ALIGNED_32(float, z, [MAX_NEURONS]);
1284 ✗ LOCAL_ALIGNED_32(float, r, [MAX_NEURONS]);
1285 ✗ LOCAL_ALIGNED_32(float, h, [MAX_NEURONS]);
1286 ✗ const int M = gru->nb_inputs;
1287 ✗ const int N = gru->nb_neurons;
1288 ✗ const int AN = FFALIGN(N, 4);
1289 ✗ const int AM = FFALIGN(M, 4);
1290 ✗ const int stride = 3 * AN, istride = 3 * AM;
1291
1292 ✗ for (int i = 0; i < N; i++) {
1293 /* Compute update gate. */
1294 ✗ float sum = gru->bias[i];
1295
1296 ✗ sum += s->fdsp->scalarproduct_float(gru->input_weights + i * istride, input, AM);
1297 ✗ sum += s->fdsp->scalarproduct_float(gru->recurrent_weights + i * stride, state, AN);
1298 ✗ z[i] = sigmoid_approx(WEIGHTS_SCALE * sum);
1299 }
1300
1301 ✗ for (int i = 0; i < N; i++) {
1302 /* Compute reset gate. */
1303 ✗ float sum = gru->bias[N + i];
1304
1305 ✗ sum += s->fdsp->scalarproduct_float(gru->input_weights + AM + i * istride, input, AM);
1306 ✗ sum += s->fdsp->scalarproduct_float(gru->recurrent_weights + AN + i * stride, state, AN);
1307 ✗ r[i] = sigmoid_approx(WEIGHTS_SCALE * sum);
1308 }
1309
1310 ✗ for (int i = 0; i < N; i++) {
1311 /* Compute output. */
1312 ✗ float sum = gru->bias[2 * N + i];
1313
1314 ✗ sum += s->fdsp->scalarproduct_float(gru->input_weights + 2 * AM + i * istride, input, AM);
1315 ✗ for (int j = 0; j < N; j++)
1316 ✗ sum += gru->recurrent_weights[2 * AN + i * stride + j] * state[j] * r[j];
1317
1318 ✗ if (gru->activation == ACTIVATION_SIGMOID)
1319 ✗ sum = sigmoid_approx(WEIGHTS_SCALE * sum);
1320 ✗ else if (gru->activation == ACTIVATION_TANH)
1321 ✗ sum = tansig_approx(WEIGHTS_SCALE * sum);
1322 ✗ else if (gru->activation == ACTIVATION_RELU)
1323 ✗ sum = FFMAX(0, WEIGHTS_SCALE * sum);
1324 else
1325 ✗ av_assert0(0);
1326 ✗ h[i] = z[i] * state[i] + (1.f - z[i]) * sum;
1327 }
1328
1329 ✗ RNN_COPY(state, h, N);
1330 ✗ }
1331
1332 #define INPUT_SIZE 42
1333
1334 ✗ static void compute_rnn(AudioRNNContext *s, RNNState *rnn, float *gains, float *vad, const float *input)
1335 {
1336 ✗ LOCAL_ALIGNED_32(float, dense_out, [MAX_NEURONS]);
1337 ✗ LOCAL_ALIGNED_32(float, noise_input, [MAX_NEURONS * 3]);
1338 ✗ LOCAL_ALIGNED_32(float, denoise_input, [MAX_NEURONS * 3]);
1339
1340 ✗ compute_dense(rnn->model->input_dense, dense_out, input);
1341 ✗ compute_gru(s, rnn->model->vad_gru, rnn->vad_gru_state, dense_out);
1342 ✗ compute_dense(rnn->model->vad_output, vad, rnn->vad_gru_state);
1343
1344 ✗ memcpy(noise_input, dense_out, rnn->model->input_dense_size * sizeof(float));
1345 ✗ memcpy(noise_input + rnn->model->input_dense_size,
1346 ✗ rnn->vad_gru_state, rnn->model->vad_gru_size * sizeof(float));
1347 ✗ memcpy(noise_input + rnn->model->input_dense_size + rnn->model->vad_gru_size,
1348 input, INPUT_SIZE * sizeof(float));
1349
1350 ✗ compute_gru(s, rnn->model->noise_gru, rnn->noise_gru_state, noise_input);
1351
1352 ✗ memcpy(denoise_input, rnn->vad_gru_state, rnn->model->vad_gru_size * sizeof(float));
1353 ✗ memcpy(denoise_input + rnn->model->vad_gru_size,
1354 ✗ rnn->noise_gru_state, rnn->model->noise_gru_size * sizeof(float));
1355 ✗ memcpy(denoise_input + rnn->model->vad_gru_size + rnn->model->noise_gru_size,
1356 input, INPUT_SIZE * sizeof(float));
1357
1358 ✗ compute_gru(s, rnn->model->denoise_gru, rnn->denoise_gru_state, denoise_input);
1359 ✗ compute_dense(rnn->model->denoise_output, gains, rnn->denoise_gru_state);
1360 ✗ }
1361
1362 ✗ static float rnnoise_channel(AudioRNNContext *s, DenoiseState *st, float *out, const float *in,
1363 int disabled)
1364 {
1365 AVComplexFloat X[FREQ_SIZE];
1366 AVComplexFloat P[WINDOW_SIZE];
1367 float x[FRAME_SIZE];
1368 float Ex[NB_BANDS], Ep[NB_BANDS];
1369 ✗ LOCAL_ALIGNED_32(float, Exp, [FFALIGN(NB_BANDS, 4)]);
1370 float features[NB_FEATURES];
1371 float g[NB_BANDS];
1372 float gf[FREQ_SIZE];
1373 ✗ float vad_prob = 0;
1374 ✗ float *history = st->history;
1375 static const float a_hp[2] = {-1.99599, 0.99600};
1376 static const float b_hp[2] = {-2, 1};
1377 int silence;
1378
1379 ✗ biquad(x, st->mem_hp_x, in, b_hp, a_hp, FRAME_SIZE);
1380 ✗ silence = compute_frame_features(s, st, X, P, Ex, Ep, Exp, features, x);
1381
1382 ✗ if (!silence && !disabled) {
1383 ✗ compute_rnn(s, &st->rnn[0], g, &vad_prob, features);
1384 ✗ pitch_filter(X, P, Ex, Ep, Exp, g);
1385 ✗ for (int i = 0; i < NB_BANDS; i++) {
1386 ✗ float alpha = .6f;
1387
1388 ✗ g[i] = FFMAX(g[i], alpha * st->lastg[i]);
1389 ✗ st->lastg[i] = g[i];
1390 }
1391
1392 ✗ interp_band_gain(gf, g);
1393
1394 ✗ for (int i = 0; i < FREQ_SIZE; i++) {
1395 ✗ X[i].re *= gf[i];
1396 ✗ X[i].im *= gf[i];
1397 }
1398 }
1399
1400 ✗ frame_synthesis(s, st, out, X);
1401 ✗ memcpy(history, in, FRAME_SIZE * sizeof(*history));
1402
1403 ✗ return vad_prob;
1404 }
1405
1406 typedef struct ThreadData {
1407 AVFrame *in, *out;
1408 } ThreadData;
1409
1410 ✗ static int rnnoise_channels(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
1411 {
1412 ✗ AudioRNNContext *s = ctx->priv;
1413 ✗ ThreadData *td = arg;
1414 ✗ AVFrame *in = td->in;
1415 ✗ AVFrame *out = td->out;
1416 ✗ const int start = ff_slice_pos(out->ch_layout.nb_channels, jobnr, nb_jobs);
1417 ✗ const int end = ff_slice_pos(out->ch_layout.nb_channels, jobnr + 1, nb_jobs);
1418
1419 ✗ for (int ch = start; ch < end; ch++) {
1420 ✗ rnnoise_channel(s, &s->st[ch],
1421 ✗ (float *)out->extended_data[ch],
1422 ✗ (const float *)in->extended_data[ch],
1423 ctx->is_disabled);
1424 }
1425
1426 ✗ return 0;
1427 }
1428
1429 ✗ static int filter_frame(AVFilterLink *inlink, AVFrame *in)
1430 {
1431 ✗ AVFilterContext *ctx = inlink->dst;
1432 ✗ AVFilterLink *outlink = ctx->outputs[0];
1433 ✗ AVFrame *out = NULL;
1434 ThreadData td;
1435
1436 ✗ out = ff_get_audio_buffer(outlink, FRAME_SIZE);
1437 ✗ if (!out) {
1438 ✗ av_frame_free(&in);
1439 ✗ return AVERROR(ENOMEM);
1440 }
1441 ✗ av_frame_copy_props(out, in);
1442
1443 ✗ td.in = in; td.out = out;
1444 ✗ ff_filter_execute(ctx, rnnoise_channels, &td, NULL,
1445 ✗ FFMIN(outlink->ch_layout.nb_channels, ff_filter_get_nb_threads(ctx)));
1446
1447 ✗ av_frame_free(&in);
1448 ✗ return ff_filter_frame(outlink, out);
1449 }
1450
1451 ✗ static int activate(AVFilterContext *ctx)
1452 {
1453 ✗ AVFilterLink *inlink = ctx->inputs[0];
1454 ✗ AVFilterLink *outlink = ctx->outputs[0];
1455 ✗ AVFrame *in = NULL;
1456 int ret;
1457
1458 ✗ FF_FILTER_FORWARD_STATUS_BACK(outlink, inlink);
1459
1460 ✗ ret = ff_inlink_consume_samples(inlink, FRAME_SIZE, FRAME_SIZE, &in);
1461 ✗ if (ret < 0)
1462 ✗ return ret;
1463
1464 ✗ if (ret > 0)
1465 ✗ return filter_frame(inlink, in);
1466
1467 ✗ FF_FILTER_FORWARD_STATUS(inlink, outlink);
1468 ✗ FF_FILTER_FORWARD_WANTED(outlink, inlink);
1469
1470 ✗ return FFERROR_NOT_READY;
1471 }
1472
1473 ✗ static int open_model(AVFilterContext *ctx, RNNModel **model)
1474 {
1475 ✗ AudioRNNContext *s = ctx->priv;
1476 int ret;
1477 FILE *f;
1478
1479 ✗ if (!s->model_name)
1480 ✗ return AVERROR(EINVAL);
1481 ✗ f = avpriv_fopen_utf8(s->model_name, "r");
1482 ✗ if (!f) {
1483 ✗ av_log(ctx, AV_LOG_ERROR, "Failed to open model file: %s\n", s->model_name);
1484 ✗ return AVERROR(EINVAL);
1485 }
1486
1487 ✗ ret = rnnoise_model_from_file(f, model);
1488 ✗ fclose(f);
1489 ✗ if (!*model || ret < 0)
1490 ✗ return ret;
1491
1492 ✗ return 0;
1493 }
1494
1495 ✗ static av_cold int init(AVFilterContext *ctx)
1496 {
1497 ✗ AudioRNNContext *s = ctx->priv;
1498 int ret;
1499
1500 ✗ s->fdsp = avpriv_float_dsp_alloc(0);
1501 ✗ if (!s->fdsp)
1502 ✗ return AVERROR(ENOMEM);
1503
1504 ✗ ret = open_model(ctx, &s->model[0]);
1505 ✗ if (ret < 0)
1506 ✗ return ret;
1507
1508 ✗ for (int i = 0; i < FRAME_SIZE; i++) {
1509 ✗ s->window[i] = sin(.5*M_PI*sin(.5*M_PI*(i+.5)/FRAME_SIZE) * sin(.5*M_PI*(i+.5)/FRAME_SIZE));
1510 ✗ s->window[WINDOW_SIZE - 1 - i] = s->window[i];
1511 }
1512
1513 ✗ for (int i = 0; i < NB_BANDS; i++) {
1514 ✗ for (int j = 0; j < NB_BANDS; j++) {
1515 ✗ s->dct_table[j][i] = cosf((i + .5f) * j * M_PI / NB_BANDS);
1516 ✗ if (j == 0)
1517 ✗ s->dct_table[j][i] *= sqrtf(.5);
1518 }
1519 }
1520
1521 ✗ return 0;
1522 }
1523
1524 ✗ static void free_model(AVFilterContext *ctx, int n)
1525 {
1526 ✗ AudioRNNContext *s = ctx->priv;
1527
1528 ✗ rnnoise_model_free(s->model[n]);
1529 ✗ s->model[n] = NULL;
1530
1531 ✗ for (int ch = 0; ch < s->channels && s->st; ch++) {
1532 ✗ av_freep(&s->st[ch].rnn[n].vad_gru_state);
1533 ✗ av_freep(&s->st[ch].rnn[n].noise_gru_state);
1534 ✗ av_freep(&s->st[ch].rnn[n].denoise_gru_state);
1535 }
1536 ✗ }
1537
1538 ✗ static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
1539 char *res, int res_len, int flags)
1540 {
1541 ✗ AudioRNNContext *s = ctx->priv;
1542 int ret;
1543
1544 ✗ ret = ff_filter_process_command(ctx, cmd, args, res, res_len, flags);
1545 ✗ if (ret < 0)
1546 ✗ return ret;
1547
1548 ✗ ret = open_model(ctx, &s->model[1]);
1549 ✗ if (ret < 0)
1550 ✗ return ret;
1551
1552 ✗ FFSWAP(RNNModel *, s->model[0], s->model[1]);
1553 ✗ for (int ch = 0; ch < s->channels; ch++)
1554 ✗ FFSWAP(RNNState, s->st[ch].rnn[0], s->st[ch].rnn[1]);
1555
1556 ✗ ret = config_input(ctx->inputs[0]);
1557 ✗ if (ret < 0) {
1558 ✗ for (int ch = 0; ch < s->channels; ch++)
1559 ✗ FFSWAP(RNNState, s->st[ch].rnn[0], s->st[ch].rnn[1]);
1560 ✗ FFSWAP(RNNModel *, s->model[0], s->model[1]);
1561 ✗ return ret;
1562 }
1563
1564 ✗ free_model(ctx, 1);
1565 ✗ return 0;
1566 }
1567
1568 ✗ static av_cold void uninit(AVFilterContext *ctx)
1569 {
1570 ✗ AudioRNNContext *s = ctx->priv;
1571
1572 ✗ av_freep(&s->fdsp);
1573 ✗ free_model(ctx, 0);
1574 ✗ for (int ch = 0; ch < s->channels && s->st; ch++) {
1575 ✗ av_tx_uninit(&s->st[ch].tx);
1576 ✗ av_tx_uninit(&s->st[ch].txi);
1577 }
1578 ✗ av_freep(&s->st);
1579 ✗ }
1580
1581 static const AVFilterPad inputs[] = {
1582 {
1583 .name = "default",
1584 .type = AVMEDIA_TYPE_AUDIO,
1585 .config_props = config_input,
1586 },
1587 };
1588
1589 #define OFFSET(x) offsetof(AudioRNNContext, x)
1590 #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_RUNTIME_PARAM
1591
1592 static const AVOption arnndn_options[] = {
1593 { "model", "set model name", OFFSET(model_name), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, AF },
1594 { "m", "set model name", OFFSET(model_name), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, AF },
1595 { "mix", "set output vs input mix", OFFSET(mix), AV_OPT_TYPE_FLOAT, {.dbl=1.0},-1, 1, AF },
1596 { NULL }
1597 };
1598
1599 AVFILTER_DEFINE_CLASS(arnndn);
1600
1601 const FFFilter ff_af_arnndn = {
1602 .p.name = "arnndn",
1603 .p.description = NULL_IF_CONFIG_SMALL("Reduce noise from speech using Recurrent Neural Networks."),
1604 .p.priv_class = &arnndn_class,
1605 .p.flags = AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL |
1606 AVFILTER_FLAG_SLICE_THREADS,
1607 .priv_size = sizeof(AudioRNNContext),
1608 .activate = activate,
1609 .init = init,
1610 .uninit = uninit,
1611 FILTER_INPUTS(inputs),
1612 FILTER_OUTPUTS(ff_audio_default_filterpad),
1613 FILTER_QUERY_FUNC2(query_formats),
1614 .process_command = process_command,
1615 };
1616