• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (c) 2020 Paul B Mahol
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/avassert.h"
24 #include "libavutil/avstring.h"
25 #include "libavutil/opt.h"
26 #include "avfilter.h"
27 #include "audio.h"
28 #include "filters.h"
29 #include "formats.h"
30 
31 enum WaveletTypes {
32     SYM2,
33     SYM4,
34     RBIOR68,
35     DEB10,
36     SYM10,
37     COIF5,
38     BL3,
39     NB_WAVELET_TYPES,
40 };
41 
42 /*
43  * All wavelets coefficients are taken from: http://wavelets.pybytes.com/
44  */
45 
46 static const double bl3_lp[42] = {
47     0.000146098, -0.000232304, -0.000285414, 0.000462093, 0.000559952,
48     -0.000927187, -0.001103748, 0.00188212, 0.002186714, -0.003882426,
49     -0.00435384, 0.008201477, 0.008685294, -0.017982291, -0.017176331,
50     0.042068328, 0.032080869, -0.110036987, -0.050201753, 0.433923147,
51     0.766130398, 0.433923147, -0.050201753, -0.110036987, 0.032080869,
52     0.042068328, -0.017176331, -0.017982291, 0.008685294, 0.008201477,
53     -0.00435384, -0.003882426, 0.002186714, 0.00188212, -0.001103748,
54     -0.000927187, 0.000559952, 0.000462093, -0.000285414, -0.000232304,
55     0.000146098, 0.0,
56 };
57 
58 static const double bl3_hp[42] = {
59     0.0, 0.000146098, 0.000232304, -0.000285414, -0.000462093, 0.000559952,
60     0.000927187, -0.001103748, -0.00188212, 0.002186714, 0.003882426,
61     -0.00435384, -0.008201477, 0.008685294, 0.017982291, -0.017176331,
62     -0.042068328, 0.032080869, 0.110036987, -0.050201753, -0.433923147,
63     0.766130398, -0.433923147, -0.050201753, 0.110036987, 0.032080869,
64     -0.042068328, -0.017176331, 0.017982291, 0.008685294, -0.008201477,
65     -0.00435384, 0.003882426, 0.002186714, -0.00188212, -0.001103748,
66     0.000927187, 0.000559952, -0.000462093, -0.000285414, 0.000232304,
67     0.000146098,
68 };
69 
70 static const double bl3_ilp[42] = {
71     0.0, 0.000146098, -0.000232304, -0.000285414, 0.000462093, 0.000559952,
72     -0.000927187, -0.001103748, 0.00188212, 0.002186714, -0.003882426,
73     -0.00435384, 0.008201477, 0.008685294, -0.017982291, -0.017176331,
74     0.042068328, 0.032080869, -0.110036987, -0.050201753, 0.433923147,
75     0.766130398, 0.433923147, -0.050201753, -0.110036987, 0.032080869,
76     0.042068328, -0.017176331, -0.017982291, 0.008685294, 0.008201477,
77     -0.00435384, -0.003882426, 0.002186714, 0.00188212, -0.001103748,
78     -0.000927187, 0.000559952, 0.000462093, -0.000285414, -0.000232304,
79     0.000146098,
80 };
81 
82 static const double bl3_ihp[42] = {
83     0.000146098, 0.000232304, -0.000285414, -0.000462093, 0.000559952,
84     0.000927187, -0.001103748, -0.00188212, 0.002186714, 0.003882426,
85     -0.00435384, -0.008201477, 0.008685294, 0.017982291, -0.017176331,
86     -0.042068328, 0.032080869, 0.110036987, -0.050201753, -0.433923147,
87     0.766130398, -0.433923147, -0.050201753, 0.110036987, 0.032080869,
88     -0.042068328, -0.017176331, 0.017982291, 0.008685294, -0.008201477,
89     -0.00435384, 0.003882426, 0.002186714, -0.00188212, -0.001103748,
90     0.000927187, 0.000559952, -0.000462093, -0.000285414, 0.000232304,
91     0.000146098,
92 };
93 
94 static const double sym10_lp[20] = {
95     0.0007701598091144901, 9.563267072289475e-05,
96     -0.008641299277022422, -0.0014653825813050513,
97     0.0459272392310922, 0.011609893903711381,
98     -0.15949427888491757, -0.07088053578324385,
99     0.47169066693843925, 0.7695100370211071,
100     0.38382676106708546, -0.03553674047381755,
101     -0.0319900568824278, 0.04999497207737669,
102     0.005764912033581909, -0.02035493981231129,
103     -0.0008043589320165449, 0.004593173585311828,
104     5.7036083618494284e-05, -0.0004593294210046588,
105 };
106 
107 static const double sym10_hp[20] = {
108     0.0004593294210046588, 5.7036083618494284e-05,
109     -0.004593173585311828, -0.0008043589320165449,
110     0.02035493981231129, 0.005764912033581909,
111     -0.04999497207737669, -0.0319900568824278,
112     0.03553674047381755, 0.38382676106708546,
113     -0.7695100370211071, 0.47169066693843925,
114     0.07088053578324385, -0.15949427888491757,
115     -0.011609893903711381, 0.0459272392310922,
116     0.0014653825813050513, -0.008641299277022422,
117     -9.563267072289475e-05, 0.0007701598091144901,
118 };
119 
120 static const double sym10_ilp[20] = {
121     -0.0004593294210046588, 5.7036083618494284e-05,
122     0.004593173585311828, -0.0008043589320165449,
123     -0.02035493981231129, 0.005764912033581909,
124     0.04999497207737669, -0.0319900568824278,
125     -0.03553674047381755, 0.38382676106708546,
126     0.7695100370211071, 0.47169066693843925,
127     -0.07088053578324385, -0.15949427888491757,
128     0.011609893903711381, 0.0459272392310922,
129     -0.0014653825813050513, -0.008641299277022422,
130     9.563267072289475e-05, 0.0007701598091144901,
131 };
132 
133 static const double sym10_ihp[20] = {
134     0.0007701598091144901, -9.563267072289475e-05,
135     -0.008641299277022422, 0.0014653825813050513,
136     0.0459272392310922, -0.011609893903711381,
137     -0.15949427888491757, 0.07088053578324385,
138     0.47169066693843925, -0.7695100370211071,
139     0.38382676106708546, 0.03553674047381755,
140     -0.0319900568824278, -0.04999497207737669,
141     0.005764912033581909, 0.02035493981231129,
142     -0.0008043589320165449, -0.004593173585311828,
143     5.7036083618494284e-05, 0.0004593294210046588,
144 };
145 
146 static const double rbior68_lp[18] = {
147     0.0, 0.0, 0.0, 0.0,
148     0.014426282505624435, 0.014467504896790148,
149     -0.07872200106262882, -0.04036797903033992,
150     0.41784910915027457, 0.7589077294536541,
151     0.41784910915027457, -0.04036797903033992,
152     -0.07872200106262882, 0.014467504896790148,
153     0.014426282505624435, 0.0, 0.0, 0.0,
154 };
155 
156 static const double rbior68_hp[18] = {
157     -0.0019088317364812906, -0.0019142861290887667,
158     0.016990639867602342, 0.01193456527972926,
159     -0.04973290349094079, -0.07726317316720414,
160     0.09405920349573646, 0.4207962846098268,
161     -0.8259229974584023, 0.4207962846098268,
162     0.09405920349573646, -0.07726317316720414,
163     -0.04973290349094079, 0.01193456527972926,
164     0.016990639867602342, -0.0019142861290887667,
165     -0.0019088317364812906, 0.0,
166 };
167 
168 static const double rbior68_ilp[18] = {
169     0.0019088317364812906, -0.0019142861290887667,
170     -0.016990639867602342, 0.01193456527972926,
171     0.04973290349094079, -0.07726317316720414,
172     -0.09405920349573646, 0.4207962846098268,
173     0.8259229974584023, 0.4207962846098268,
174     -0.09405920349573646, -0.07726317316720414,
175     0.04973290349094079, 0.01193456527972926,
176     -0.016990639867602342, -0.0019142861290887667,
177     0.0019088317364812906, 0.0,
178 };
179 
180 static const double rbior68_ihp[18] = {
181     0.0, 0.0, 0.0, 0.0,
182     0.014426282505624435, -0.014467504896790148,
183     -0.07872200106262882, 0.04036797903033992,
184     0.41784910915027457, -0.7589077294536541,
185     0.41784910915027457, 0.04036797903033992,
186     -0.07872200106262882, -0.014467504896790148,
187     0.014426282505624435, 0.0, 0.0, 0.0,
188 };
189 
190 static const double coif5_lp[30] = {
191     -9.517657273819165e-08, -1.6744288576823017e-07,
192     2.0637618513646814e-06, 3.7346551751414047e-06,
193     -2.1315026809955787e-05, -4.134043227251251e-05,
194     0.00014054114970203437, 0.00030225958181306315,
195     -0.0006381313430451114, -0.0016628637020130838,
196     0.0024333732126576722, 0.006764185448053083,
197     -0.009164231162481846, -0.01976177894257264,
198     0.03268357426711183, 0.0412892087501817,
199     -0.10557420870333893, -0.06203596396290357,
200     0.4379916261718371, 0.7742896036529562,
201     0.4215662066908515, -0.05204316317624377,
202     -0.09192001055969624, 0.02816802897093635,
203     0.023408156785839195, -0.010131117519849788,
204     -0.004159358781386048, 0.0021782363581090178,
205     0.00035858968789573785, -0.00021208083980379827,
206 };
207 
208 static const double coif5_hp[30] = {
209     0.00021208083980379827, 0.00035858968789573785,
210     -0.0021782363581090178, -0.004159358781386048,
211     0.010131117519849788, 0.023408156785839195,
212     -0.02816802897093635, -0.09192001055969624,
213     0.05204316317624377, 0.4215662066908515,
214     -0.7742896036529562, 0.4379916261718371,
215     0.06203596396290357, -0.10557420870333893,
216     -0.0412892087501817, 0.03268357426711183,
217     0.01976177894257264, -0.009164231162481846,
218     -0.006764185448053083, 0.0024333732126576722,
219     0.0016628637020130838, -0.0006381313430451114,
220     -0.00030225958181306315, 0.00014054114970203437,
221     4.134043227251251e-05, -2.1315026809955787e-05,
222     -3.7346551751414047e-06, 2.0637618513646814e-06,
223     1.6744288576823017e-07, -9.517657273819165e-08,
224 };
225 
226 static const double coif5_ilp[30] = {
227     -0.00021208083980379827, 0.00035858968789573785,
228     0.0021782363581090178, -0.004159358781386048,
229     -0.010131117519849788, 0.023408156785839195,
230     0.02816802897093635, -0.09192001055969624,
231     -0.05204316317624377, 0.4215662066908515,
232     0.7742896036529562, 0.4379916261718371,
233     -0.06203596396290357, -0.10557420870333893,
234     0.0412892087501817, 0.03268357426711183,
235     -0.01976177894257264, -0.009164231162481846,
236     0.006764185448053083, 0.0024333732126576722,
237     -0.0016628637020130838, -0.0006381313430451114,
238     0.00030225958181306315, 0.00014054114970203437,
239     -4.134043227251251e-05, -2.1315026809955787e-05,
240     3.7346551751414047e-06, 2.0637618513646814e-06,
241     -1.6744288576823017e-07, -9.517657273819165e-08,
242 };
243 
244 static const double coif5_ihp[30] = {
245     -9.517657273819165e-08, 1.6744288576823017e-07,
246     2.0637618513646814e-06, -3.7346551751414047e-06,
247     -2.1315026809955787e-05, 4.134043227251251e-05,
248     0.00014054114970203437, -0.00030225958181306315,
249     -0.0006381313430451114, 0.0016628637020130838,
250     0.0024333732126576722, -0.006764185448053083,
251     -0.009164231162481846, 0.01976177894257264,
252     0.03268357426711183, -0.0412892087501817,
253     -0.10557420870333893, 0.06203596396290357,
254     0.4379916261718371, -0.7742896036529562,
255     0.4215662066908515, 0.05204316317624377,
256     -0.09192001055969624, -0.02816802897093635,
257     0.023408156785839195, 0.010131117519849788,
258     -0.004159358781386048, -0.0021782363581090178,
259     0.00035858968789573785, 0.00021208083980379827,
260 };
261 
262 static const double deb10_lp[20] = {
263     -1.326420300235487e-05, 9.358867000108985e-05,
264     -0.0001164668549943862, -0.0006858566950046825,
265     0.00199240529499085, 0.0013953517469940798,
266     -0.010733175482979604, 0.0036065535669883944,
267     0.03321267405893324, -0.02945753682194567,
268     -0.07139414716586077, 0.09305736460380659,
269     0.12736934033574265, -0.19594627437659665,
270     -0.24984642432648865, 0.2811723436604265,
271     0.6884590394525921, 0.5272011889309198,
272     0.18817680007762133, 0.026670057900950818,
273 };
274 
275 static const double deb10_hp[20] = {
276     -0.026670057900950818, 0.18817680007762133,
277     -0.5272011889309198, 0.6884590394525921,
278     -0.2811723436604265, -0.24984642432648865,
279     0.19594627437659665, 0.12736934033574265,
280     -0.09305736460380659, -0.07139414716586077,
281     0.02945753682194567, 0.03321267405893324,
282     -0.0036065535669883944, -0.010733175482979604,
283     -0.0013953517469940798, 0.00199240529499085,
284     0.0006858566950046825, -0.0001164668549943862,
285     -9.358867000108985e-05, -1.326420300235487e-05,
286 };
287 
288 static const double deb10_ilp[20] = {
289     0.026670057900950818, 0.18817680007762133,
290     0.5272011889309198, 0.6884590394525921,
291     0.2811723436604265, -0.24984642432648865,
292     -0.19594627437659665, 0.12736934033574265,
293     0.09305736460380659, -0.07139414716586077,
294     -0.02945753682194567, 0.03321267405893324,
295     0.0036065535669883944, -0.010733175482979604,
296     0.0013953517469940798, 0.00199240529499085,
297     -0.0006858566950046825, -0.0001164668549943862,
298     9.358867000108985e-05, -1.326420300235487e-05,
299 };
300 
301 static const double deb10_ihp[20] = {
302     -1.326420300235487e-05, -9.358867000108985e-05,
303     -0.0001164668549943862, 0.0006858566950046825,
304     0.00199240529499085, -0.0013953517469940798,
305     -0.010733175482979604, -0.0036065535669883944,
306     0.03321267405893324, 0.02945753682194567,
307     -0.07139414716586077, -0.09305736460380659,
308     0.12736934033574265, 0.19594627437659665,
309     -0.24984642432648865, -0.2811723436604265,
310     0.6884590394525921, -0.5272011889309198,
311     0.18817680007762133, -0.026670057900950818,
312 };
313 
314 static const double sym4_lp[8] = {
315     -0.07576571478927333,
316     -0.02963552764599851,
317     0.49761866763201545,
318     0.8037387518059161,
319     0.29785779560527736,
320     -0.09921954357684722,
321     -0.012603967262037833,
322     0.0322231006040427,
323 };
324 
325 static const double sym4_hp[8] = {
326     -0.0322231006040427,
327     -0.012603967262037833,
328     0.09921954357684722,
329     0.29785779560527736,
330     -0.8037387518059161,
331     0.49761866763201545,
332     0.02963552764599851,
333     -0.07576571478927333,
334 };
335 
336 static const double sym4_ilp[8] = {
337     0.0322231006040427,
338     -0.012603967262037833,
339     -0.09921954357684722,
340     0.29785779560527736,
341     0.8037387518059161,
342     0.49761866763201545,
343     -0.02963552764599851,
344     -0.07576571478927333,
345 };
346 
347 static const double sym4_ihp[8] = {
348     -0.07576571478927333,
349     0.02963552764599851,
350     0.49761866763201545,
351     -0.8037387518059161,
352     0.29785779560527736,
353     0.09921954357684722,
354     -0.012603967262037833,
355     -0.0322231006040427,
356 };
357 
358 static const double sym2_lp[4] = {
359     -0.12940952255092145, 0.22414386804185735,
360     0.836516303737469, 0.48296291314469025,
361 };
362 
363 static const double sym2_hp[4] = {
364     -0.48296291314469025, 0.836516303737469,
365     -0.22414386804185735, -0.12940952255092145,
366 };
367 
368 static const double sym2_ilp[4] = {
369     0.48296291314469025, 0.836516303737469,
370     0.22414386804185735, -0.12940952255092145,
371 };
372 
373 static const double sym2_ihp[4] = {
374     -0.12940952255092145, -0.22414386804185735,
375     0.836516303737469, -0.48296291314469025,
376 };
377 
378 #define MAX_LEVELS 13
379 
380 typedef struct ChannelParams {
381     int *output_length;
382     int *filter_length;
383     double **output_coefs;
384     double **subbands_to_free;
385     double **filter_coefs;
386 
387     int tempa_length;
388     int tempa_len_max;
389     int temp_in_length;
390     int temp_in_max_length;
391     int buffer_length;
392     int min_left_ext;
393     int max_left_ext;
394 
395     double *tempa;
396     double *tempd;
397     double *temp_in;
398     double *buffer;
399     double *buffer2;
400     double *prev;
401     double *overlap;
402 } ChannelParams;
403 
404 typedef struct AudioFWTDNContext {
405     const AVClass *class;
406 
407     double sigma;
408     double percent;
409     double softness;
410 
411     uint64_t sn;
412     int64_t eof_pts;
413     int eof;
414 
415     int wavelet_type;
416     int channels;
417     int nb_samples;
418     int levels;
419     int wavelet_length;
420     int need_profile;
421     int got_profile;
422     int adaptive;
423 
424     int delay;
425     int drop_samples;
426     int padd_samples;
427     int overlap_length;
428     int prev_length;
429     ChannelParams *cp;
430 
431     const double *lp, *hp;
432     const double *ilp, *ihp;
433 
434     AVFrame *stddev, *absmean, *filter;
435     AVFrame *new_stddev, *new_absmean;
436 
437     int (*filter_channel)(AVFilterContext *ctx, void *arg, int ch, int nb_jobs);
438 } AudioFWTDNContext;
439 
440 #define OFFSET(x) offsetof(AudioFWTDNContext, x)
441 #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
442 #define AFR AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_RUNTIME_PARAM
443 
444 static const AVOption afwtdn_options[] = {
445     { "sigma", "set noise sigma", OFFSET(sigma), AV_OPT_TYPE_DOUBLE, {.dbl=0}, 0, 1, AFR },
446     { "levels", "set number of wavelet levels", OFFSET(levels), AV_OPT_TYPE_INT, {.i64=10}, 1, MAX_LEVELS-1, AF },
447     { "wavet", "set wavelet type", OFFSET(wavelet_type), AV_OPT_TYPE_INT, {.i64=SYM10}, 0, NB_WAVELET_TYPES - 1, AF, "wavet" },
448     { "sym2", "sym2", 0, AV_OPT_TYPE_CONST, {.i64=SYM2}, 0, 0, AF, "wavet" },
449     { "sym4", "sym4", 0, AV_OPT_TYPE_CONST, {.i64=SYM4}, 0, 0, AF, "wavet" },
450     { "rbior68", "rbior68", 0, AV_OPT_TYPE_CONST, {.i64=RBIOR68}, 0, 0, AF, "wavet" },
451     { "deb10", "deb10", 0, AV_OPT_TYPE_CONST, {.i64=DEB10}, 0, 0, AF, "wavet" },
452     { "sym10", "sym10", 0, AV_OPT_TYPE_CONST, {.i64=SYM10}, 0, 0, AF, "wavet" },
453     { "coif5", "coif5", 0, AV_OPT_TYPE_CONST, {.i64=COIF5}, 0, 0, AF, "wavet" },
454     { "bl3", "bl3", 0, AV_OPT_TYPE_CONST, {.i64=BL3}, 0, 0, AF, "wavet" },
455     { "percent", "set percent of full denoising", OFFSET(percent),AV_OPT_TYPE_DOUBLE, {.dbl=85}, 0, 100, AFR },
456     { "profile", "profile noise", OFFSET(need_profile), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, AFR },
457     { "adaptive", "adaptive profiling of noise", OFFSET(adaptive), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, AFR },
458     { "samples", "set frame size in number of samples", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=8192}, 512, 65536, AF },
459     { "softness", "set thresholding softness", OFFSET(softness), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 10, AFR },
460     { NULL }
461 };
462 
463 AVFILTER_DEFINE_CLASS(afwtdn);
464 
465 #define pow2(x) (1U << (x))
466 #define mod_pow2(x, power_of_two) ((x) & ((power_of_two) - 1))
467 
conv_down(double * in,int in_length,double * low,double * high,int out_length,const double * lp,const double * hp,int wavelet_length,int skip,double * buffer,int buffer_length)468 static void conv_down(double *in, int in_length, double *low, double *high,
469                       int out_length, const double *lp, const double *hp,
470                       int wavelet_length, int skip,
471                       double *buffer, int buffer_length)
472 {
473     double thigh = 0.0, tlow = 0.0;
474     int buff_idx = 1 + skip;
475 
476     memcpy(buffer, in, buff_idx * sizeof(*buffer));
477     memset(buffer + buff_idx, 0, (buffer_length - buff_idx) * sizeof(*buffer));
478 
479     for (int i = 0; i < out_length - 1; i++) {
480         double thigh = 0.0, tlow = 0.0;
481 
482         for (int j = 0; j < wavelet_length; j++) {
483             const int idx = mod_pow2(-j + buff_idx - 1, buffer_length);
484             const double btemp = buffer[idx];
485 
486             thigh += btemp * hp[j];
487             tlow += btemp * lp[j];
488         }
489 
490         high[i] = thigh;
491         low[i] = tlow;
492         buffer[buff_idx++] = in[2 * i + 1 + skip];
493         buffer[buff_idx++] = in[2 * i + 2 + skip];
494         buff_idx = mod_pow2(buff_idx, buffer_length);
495     }
496 
497     for (int i = 0; i < wavelet_length; i++) {
498         const int idx = mod_pow2(-i + buff_idx - 1, buffer_length);
499         const double btemp = buffer[idx];
500 
501         thigh += btemp * hp[i];
502         tlow += btemp * lp[i];
503     }
504 
505     high[out_length - 1] = thigh;
506     low[out_length - 1] = tlow;
507 }
508 
left_ext(int wavelet_length,int levels,uint64_t sn)509 static int left_ext(int wavelet_length, int levels, uint64_t sn)
510 {
511     if (!sn)
512         return 0;
513     return (pow2(levels) - 1) * (wavelet_length - 2) + mod_pow2(sn, pow2(levels));
514 }
515 
nb_coefs(int length,int level,uint64_t sn)516 static int nb_coefs(int length, int level, uint64_t sn)
517 {
518     const int pow2_level = pow2(level);
519 
520     return (sn + length) / pow2_level - sn / pow2_level;
521 }
522 
reallocate_inputs(double ** out,int * out_length,int in_length,int levels,int ch,uint64_t sn)523 static int reallocate_inputs(double **out, int *out_length,
524                              int in_length, int levels, int ch, uint64_t sn)
525 {
526     const int temp_length = nb_coefs(in_length, levels, sn);
527 
528     for (int level = 0; level < levels; level++) {
529         const int temp_length = nb_coefs(in_length, level + 1, sn);
530 
531         if (temp_length > out_length[level]) {
532             av_freep(&out[level]);
533             out_length[level] = 0;
534 
535             out[level] = av_calloc(temp_length + 1, sizeof(**out));
536             if (!out[level])
537                 return AVERROR(ENOMEM);
538             out_length[level] = temp_length + 1;
539         }
540 
541         memset(out[level] + temp_length, 0,
542                (out_length[level] - temp_length) * sizeof(**out));
543         out_length[level] = temp_length;
544     }
545 
546     if (temp_length > out_length[levels]) {
547         av_freep(&out[levels]);
548         out_length[levels] = 0;
549 
550         out[levels] = av_calloc(temp_length + 1, sizeof(**out));
551         if (!out[levels])
552             return AVERROR(ENOMEM);
553         out_length[levels] = temp_length + 1;
554     }
555 
556     memset(out[levels] + temp_length, 0,
557            (out_length[levels] - temp_length) * sizeof(**out));
558     out_length[levels] = temp_length;
559 
560     return 0;
561 }
562 
max_left_zeros_inverse(int levels,int level,int wavelet_length)563 static int max_left_zeros_inverse(int levels, int level, int wavelet_length)
564 {
565     return (pow2(levels - level) - 1) * (wavelet_length - 1);
566 }
567 
reallocate_outputs(AudioFWTDNContext * s,double ** out,int * out_length,int in_length,int levels,int ch,uint64_t sn)568 static int reallocate_outputs(AudioFWTDNContext *s,
569                               double **out, int *out_length,
570                               int in_length, int levels, int ch, uint64_t sn)
571 {
572     ChannelParams *cp = &s->cp[ch];
573     int temp_length = 0;
574     int add = 0;
575 
576     for (int level = 0; level < levels; level++) {
577         temp_length = nb_coefs(in_length, level + 1, sn);
578         if (temp_length > out_length[level]) {
579             av_freep(&cp->subbands_to_free[level]);
580             out_length[level] = 0;
581 
582             add = max_left_zeros_inverse(levels, level + 1, s->wavelet_length);
583             cp->subbands_to_free[level] = av_calloc(add + temp_length + 1, sizeof(**out));
584             if (!cp->subbands_to_free[level])
585                 return AVERROR(ENOMEM);
586             out_length[level] = add + temp_length + 1;
587             out[level] = cp->subbands_to_free[level] + add;
588         }
589 
590         memset(out[level] + temp_length, 0,
591                FFMAX(out_length[level] - temp_length - add, 0) * sizeof(**out));
592         out_length[level] = temp_length;
593     }
594 
595     temp_length = nb_coefs(in_length, levels, sn);
596     if (temp_length > out_length[levels]) {
597         av_freep(&cp->subbands_to_free[levels]);
598         out_length[levels] = 0;
599 
600         cp->subbands_to_free[levels] = av_calloc(temp_length + 1, sizeof(**out));
601         if (!cp->subbands_to_free[levels])
602             return AVERROR(ENOMEM);
603         out_length[levels] = temp_length + 1;
604         out[levels] = cp->subbands_to_free[levels];
605     }
606 
607     memset(out[levels] + temp_length, 0,
608            (out_length[levels] - temp_length) * sizeof(**out));
609     out_length[levels] = temp_length;
610 
611     return 0;
612 }
613 
discard_left_ext(int wavelet_length,int levels,int level,uint64_t sn)614 static int discard_left_ext(int wavelet_length, int levels, int level, uint64_t sn)
615 {
616     if (levels == level || sn == 0)
617         return 0;
618     return (pow2(levels - level) - 1) * (wavelet_length - 2) + mod_pow2(sn, pow2(levels)) / pow2(level);
619 }
620 
forward(AudioFWTDNContext * s,const double * in,int in_length,double ** out,int * out_length,int ch,uint64_t sn)621 static int forward(AudioFWTDNContext *s,
622                    const double *in, int in_length,
623                    double **out, int *out_length, int ch, uint64_t sn)
624 {
625     ChannelParams *cp = &s->cp[ch];
626     int levels = s->levels;
627     int skip = sn ? s->wavelet_length - 1 : 1;
628     int leftext, ret;
629 
630     ret = reallocate_inputs(out, out_length, in_length, levels, ch, sn);
631     if (ret < 0)
632         return ret;
633     ret = reallocate_outputs(s, cp->filter_coefs, cp->filter_length,
634                              in_length, levels, ch, sn);
635     if (ret < 0)
636         return ret;
637 
638     leftext = left_ext(s->wavelet_length, levels, sn);
639 
640     if (cp->temp_in_max_length < in_length + cp->max_left_ext + skip) {
641         av_freep(&cp->temp_in);
642         cp->temp_in_max_length = in_length + cp->max_left_ext + skip;
643         cp->temp_in = av_calloc(cp->temp_in_max_length, sizeof(*cp->temp_in));
644         if (!cp->temp_in) {
645             cp->temp_in_max_length = 0;
646             return AVERROR(ENOMEM);
647         }
648     }
649 
650     memset(cp->temp_in, 0, cp->temp_in_max_length * sizeof(*cp->temp_in));
651     cp->temp_in_length = in_length + leftext;
652 
653     if (leftext)
654         memcpy(cp->temp_in, cp->prev + s->prev_length - leftext, leftext * sizeof(*cp->temp_in));
655     memcpy(cp->temp_in + leftext, in, in_length * sizeof(*in));
656 
657     if (levels == 1) {
658         conv_down(cp->temp_in, cp->temp_in_length, out[1], out[0], out_length[1],
659                  s->lp, s->hp, s->wavelet_length, skip,
660                  cp->buffer, cp->buffer_length);
661     } else {
662         int discard = discard_left_ext(s->wavelet_length, levels, 1, sn);
663         int tempa_length_prev;
664 
665         if (cp->tempa_len_max < (in_length + cp->max_left_ext + s->wavelet_length - 1) / 2) {
666             av_freep(&cp->tempa);
667             av_freep(&cp->tempd);
668             cp->tempa_len_max = (in_length + cp->max_left_ext + s->wavelet_length - 1) / 2;
669             cp->tempa = av_calloc(cp->tempa_len_max, sizeof(*cp->tempa));
670             cp->tempd = av_calloc(cp->tempa_len_max, sizeof(*cp->tempd));
671             if (!cp->tempa || !cp->tempd) {
672                 cp->tempa_len_max = 0;
673                 return AVERROR(ENOMEM);
674             }
675         }
676 
677         memset(cp->tempa, 0, cp->tempa_len_max * sizeof(*cp->tempa));
678         memset(cp->tempd, 0, cp->tempa_len_max * sizeof(*cp->tempd));
679 
680         cp->tempa_length = out_length[0] + discard;
681         conv_down(cp->temp_in, cp->temp_in_length,
682                  cp->tempa, cp->tempd, cp->tempa_length,
683                  s->lp, s->hp, s->wavelet_length, skip,
684                  cp->buffer, cp->buffer_length);
685         memcpy(out[0], cp->tempd + discard, out_length[0] * sizeof(**out));
686         tempa_length_prev = cp->tempa_length;
687 
688         for (int level = 1; level < levels - 1; level++) {
689             if (out_length[level] == 0)
690                 return 0;
691             discard = discard_left_ext(s->wavelet_length, levels, level + 1, sn);
692             cp->tempa_length = out_length[level] + discard;
693             conv_down(cp->tempa, tempa_length_prev,
694                      cp->tempa, cp->tempd, cp->tempa_length,
695                      s->lp, s->hp, s->wavelet_length, skip,
696                      cp->buffer, cp->buffer_length);
697             memcpy(out[level], cp->tempd + discard, out_length[level] * sizeof(**out));
698             tempa_length_prev = cp->tempa_length;
699         }
700 
701         if (out_length[levels] == 0)
702             return 0;
703         conv_down(cp->tempa, cp->tempa_length, out[levels], out[levels - 1], out_length[levels],
704                  s->lp, s->hp, s->wavelet_length, skip,
705                  cp->buffer, cp->buffer_length);
706     }
707 
708     if (s->prev_length < in_length) {
709         memcpy(cp->prev, in + in_length - cp->max_left_ext, cp->max_left_ext * sizeof(*cp->prev));
710     } else {
711         memmove(cp->prev, cp->prev + in_length, (s->prev_length - in_length) * sizeof(*cp->prev));
712         memcpy(cp->prev + s->prev_length - in_length, in, in_length * sizeof(*cp->prev));
713     }
714 
715     return 0;
716 }
717 
conv_up(double * low,double * high,int in_length,double * out,int out_length,const double * lp,const double * hp,int filter_length,double * buffer,double * buffer2,int buffer_length)718 static void conv_up(double *low, double *high, int in_length, double *out, int out_length,
719                     const double *lp, const double *hp, int filter_length,
720                     double *buffer, double *buffer2, int buffer_length)
721 {
722     int shift = 0, buff_idx = 0, in_idx = 0;
723 
724     memset(buffer, 0, buffer_length * sizeof(*buffer));
725     memset(buffer2, 0, buffer_length * sizeof(*buffer2));
726 
727     for (int i = 0; i < out_length; i++) {
728         double sum = 0.0;
729 
730         if ((i & 1) == 0) {
731             if (in_idx < in_length) {
732                 buffer[buff_idx] = low[in_idx];
733                 buffer2[buff_idx] = high[in_idx++];
734             } else {
735                 buffer[buff_idx] = 0;
736                 buffer2[buff_idx] = 0;
737             }
738             buff_idx++;
739             if (buff_idx >= buffer_length)
740                 buff_idx = 0;
741             shift = 0;
742         }
743 
744         for (int j = 0; j < (filter_length - shift + 1) / 2; j++) {
745             const int idx = mod_pow2(-j + buff_idx - 1, buffer_length);
746 
747             sum += buffer[idx] * lp[j * 2 + shift] + buffer2[idx] * hp[j * 2 + shift];
748         }
749         out[i] = sum;
750         shift = 1;
751     }
752 }
753 
append_left_ext(int wavelet_length,int levels,int level,uint64_t sn)754 static int append_left_ext(int wavelet_length, int levels, int level, uint64_t sn)
755 {
756     if (levels == level)
757         return 0;
758 
759     return (pow2(levels - level) - 1) * (wavelet_length - 2) +
760             mod_pow2(sn, pow2(levels)) / pow2(level);
761 }
762 
inverse(AudioFWTDNContext * s,double ** in,int * in_length,double * out,int out_length,int ch,uint64_t sn)763 static int inverse(AudioFWTDNContext *s,
764                    double **in, int *in_length,
765                    double *out, int out_length, int ch, uint64_t sn)
766 {
767     ChannelParams *cp = &s->cp[ch];
768     const int levels = s->levels;
769     int leftext = left_ext(s->wavelet_length, levels, sn);
770     int temp_skip = 0;
771 
772     if (sn == 0)
773         temp_skip = cp->min_left_ext;
774 
775     memset(out, 0, out_length * sizeof(*out));
776 
777     if (cp->temp_in_max_length < out_length + cp->max_left_ext + s->wavelet_length - 1) {
778         av_freep(&cp->temp_in);
779         cp->temp_in_max_length = out_length + cp->max_left_ext + s->wavelet_length - 1;
780         cp->temp_in = av_calloc(cp->temp_in_max_length, sizeof(*cp->temp_in));
781         if (!cp->temp_in) {
782             cp->temp_in_max_length = 0;
783             return AVERROR(ENOMEM);
784         }
785     }
786 
787     memset(cp->temp_in, 0, cp->temp_in_max_length * sizeof(*cp->temp_in));
788     cp->temp_in_length = out_length + cp->max_left_ext;
789 
790     if (levels == 1) {
791         conv_up(in[1], in[0], in_length[1], cp->temp_in, cp->temp_in_length,
792                 s->ilp, s->ihp, s->wavelet_length,
793                 cp->buffer, cp->buffer2, cp->buffer_length);
794         memcpy(out + cp->max_left_ext - leftext, cp->temp_in + temp_skip,
795                FFMAX(0, out_length - (cp->max_left_ext - leftext)) * sizeof(*out));
796     } else {
797         double *hp1, *hp2;
798         int add, add2;
799 
800         if (cp->tempa_len_max < (out_length + cp->max_left_ext + s->wavelet_length - 1) / 2) {
801             av_freep(&cp->tempa);
802             cp->tempa_len_max = (out_length + cp->max_left_ext + s->wavelet_length - 1) / 2;
803             cp->tempa = av_calloc(cp->tempa_len_max, sizeof(*cp->tempa));
804             if (!cp->tempa) {
805                 cp->tempa_len_max = 0;
806                 return AVERROR(ENOMEM);
807             }
808         }
809 
810         memset(cp->tempa, 0, cp->tempa_len_max * sizeof(*cp->tempa));
811 
812         hp1 = levels & 1 ? cp->temp_in : cp->tempa;
813         hp2 = levels & 1 ? cp->tempa : cp->temp_in;
814 
815         add = append_left_ext(s->wavelet_length, levels, levels - 1, sn);
816         conv_up(in[levels], in[levels - 1], in_length[levels], hp1, in_length[levels - 2] + add,
817                s->ilp, s->ihp, s->wavelet_length, cp->buffer, cp->buffer2, cp->buffer_length);
818 
819         for (int level = levels - 1; level > 1; level--) {
820             add2 = append_left_ext(s->wavelet_length, levels, level - 1, sn);
821             add = append_left_ext(s->wavelet_length, levels, level, sn);
822             conv_up(hp1, in[level - 1] - add, in_length[level - 1] + add,
823                     hp2, in_length[level - 2] + add2,
824                     s->ilp, s->ihp, s->wavelet_length,
825                     cp->buffer, cp->buffer2, cp->buffer_length);
826             FFSWAP(double *, hp1, hp2);
827         }
828 
829         add = append_left_ext(s->wavelet_length, levels, 1, sn);
830         conv_up(hp1, in[0] - add, in_length[0] + add, cp->temp_in, cp->temp_in_length,
831                s->ilp, s->ihp, s->wavelet_length,
832                cp->buffer, cp->buffer2, cp->buffer_length);
833     }
834 
835     memset(cp->temp_in, 0, temp_skip * sizeof(*cp->temp_in));
836     if (s->overlap_length <= out_length) {
837         memcpy(out + cp->max_left_ext - leftext, cp->temp_in + temp_skip,
838                FFMAX(0, out_length - (cp->max_left_ext - leftext)) * sizeof(*out));
839         for (int i = 0;i < FFMIN(s->overlap_length, out_length); i++)
840             out[i] += cp->overlap[i];
841 
842         memcpy(cp->overlap, cp->temp_in + out_length - (cp->max_left_ext - leftext),
843                s->overlap_length * sizeof(*cp->overlap));
844     } else {
845         for (int i = 0;i < s->overlap_length - (cp->max_left_ext - leftext); i++)
846             cp->overlap[i + cp->max_left_ext - leftext] += cp->temp_in[i];
847         memcpy(out, cp->overlap, out_length * sizeof(*out));
848         memmove(cp->overlap, cp->overlap + out_length,
849                 (s->overlap_length - out_length) * sizeof(*cp->overlap));
850         memcpy(cp->overlap + s->overlap_length - out_length, cp->temp_in + leftext,
851                out_length * sizeof(*cp->overlap));
852     }
853 
854     return 0;
855 }
856 
next_pow2(int in)857 static int next_pow2(int in)
858 {
859     return 1 << (av_log2(in) + 1);
860 }
861 
denoise_level(double * out,const double * in,const double * filter,double percent,int length)862 static void denoise_level(double *out, const double *in,
863                           const double *filter,
864                           double percent, int length)
865 {
866     const double x = percent * 0.01;
867     const double y = 1.0 - x;
868 
869     for (int i = 0; i < length; i++)
870         out[i] = x * filter[i] + in[i] * y;
871 }
872 
sqr(double in)873 static double sqr(double in)
874 {
875     return in * in;
876 }
877 
measure_mean(const double * in,int length)878 static double measure_mean(const double *in, int length)
879 {
880     double sum = 0.0;
881 
882     for (int i = 0; i < length; i++)
883         sum += in[i];
884 
885     return sum / length;
886 }
887 
measure_absmean(const double * in,int length)888 static double measure_absmean(const double *in, int length)
889 {
890     double sum = 0.0;
891 
892     for (int i = 0; i < length; i++)
893         sum += fabs(in[i]);
894 
895     return sum / length;
896 }
897 
measure_stddev(const double * in,int length,double mean)898 static double measure_stddev(const double *in, int length, double mean)
899 {
900     double sum = 0.;
901 
902     for (int i = 0; i < length; i++) {
903         sum += sqr(in[i] - mean);
904     }
905 
906     return sqrt(sum / length);
907 }
908 
noise_filter(const double stddev,const double * in,double * out,double absmean,double softness,double new_stddev,int length)909 static void noise_filter(const double stddev, const double *in,
910                          double *out, double absmean, double softness,
911                          double new_stddev, int length)
912 {
913     for (int i = 0; i < length; i++) {
914         if (new_stddev <= stddev)
915             out[i] = 0.0;
916         else if (fabs(in[i]) <= absmean)
917             out[i] = 0.0;
918         else
919             out[i] = in[i] - FFSIGN(in[i]) * absmean / exp(3.0 * softness * (fabs(in[i]) - absmean) / absmean);
920     }
921 }
922 
923 typedef struct ThreadData {
924     AVFrame *in, *out;
925 } ThreadData;
926 
filter_channel(AVFilterContext * ctx,void * arg,int ch,int nb_jobs)927 static int filter_channel(AVFilterContext *ctx, void *arg, int ch, int nb_jobs)
928 {
929     AudioFWTDNContext *s = ctx->priv;
930     ThreadData *td = arg;
931     AVFrame *in = td->in;
932     AVFrame *out = td->out;
933     ChannelParams *cp = &s->cp[ch];
934     const double *src = (const double *)(in->extended_data[ch]);
935     double *dst = (double *)out->extended_data[ch];
936     double *absmean = (double *)s->absmean->extended_data[ch];
937     double *new_absmean = (double *)s->new_absmean->extended_data[ch];
938     double *stddev = (double *)s->stddev->extended_data[ch];
939     double *new_stddev = (double *)s->new_stddev->extended_data[ch];
940     double *filter = (double *)s->filter->extended_data[ch];
941     double is_noise = 0.0;
942     int ret;
943 
944     ret = forward(s, src, in->nb_samples, cp->output_coefs, cp->output_length, ch, s->sn);
945     if (ret < 0)
946         return ret;
947 
948     if (!s->got_profile && s->need_profile) {
949         for (int level = 0; level <= s->levels; level++) {
950             const int length = cp->output_length[level];
951             const double scale = sqrt(2.0 * log(length));
952 
953             stddev[level] = measure_stddev(cp->output_coefs[level], length,
954                             measure_mean(cp->output_coefs[level], length)) * scale;
955             absmean[level] = measure_absmean(cp->output_coefs[level], length) * scale;
956         }
957     } else if (!s->got_profile && !s->need_profile && !s->adaptive) {
958         for (int level = 0; level <= s->levels; level++) {
959             const int length = cp->output_length[level];
960             const double scale = sqrt(2.0 * log(length));
961 
962             stddev[level] = 0.5 * s->sigma * scale;
963             absmean[level] = 0.5 * s->sigma * scale;
964         }
965     }
966 
967     for (int level = 0; level <= s->levels; level++) {
968         const int length = cp->output_length[level];
969         double vad;
970 
971         new_stddev[level] = measure_stddev(cp->output_coefs[level], length,
972                             measure_mean(cp->output_coefs[level], length));
973         new_absmean[level] = measure_absmean(cp->output_coefs[level], length);
974         if (new_absmean[level] <= FLT_EPSILON)
975             vad = 1.0;
976         else
977             vad = new_stddev[level] / new_absmean[level];
978         if (level < s->levels)
979             is_noise += sqr(vad - 1.232);
980     }
981 
982     is_noise *= in->sample_rate;
983     is_noise /= s->nb_samples;
984     for (int level = 0; level <= s->levels; level++) {
985         const double percent = ctx->is_disabled ? 0. : s->percent;
986         const int length = cp->output_length[level];
987         const double scale = sqrt(2.0 * log(length));
988 
989         if (is_noise < 0.05 && s->adaptive) {
990             stddev[level] = new_stddev[level] * scale;
991             absmean[level] = new_absmean[level] * scale;
992         }
993 
994         noise_filter(stddev[level], cp->output_coefs[level], filter, absmean[level],
995                      s->softness, new_stddev[level], length);
996         denoise_level(cp->filter_coefs[level], cp->output_coefs[level], filter, percent, length);
997     }
998 
999     ret = inverse(s, cp->filter_coefs, cp->filter_length, dst, out->nb_samples, ch, s->sn);
1000     if (ret < 0)
1001         return ret;
1002 
1003     return 0;
1004 }
1005 
filter_frame(AVFilterLink * inlink,AVFrame * in)1006 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
1007 {
1008     AVFilterContext *ctx = inlink->dst;
1009     AudioFWTDNContext *s = ctx->priv;
1010     AVFilterLink *outlink = ctx->outputs[0];
1011     ThreadData td;
1012     AVFrame *out;
1013     int eof = in == NULL;
1014 
1015     out = ff_get_audio_buffer(outlink, s->nb_samples);
1016     if (!out) {
1017         av_frame_free(&in);
1018         return AVERROR(ENOMEM);
1019     }
1020     if (in) {
1021         av_frame_copy_props(out, in);
1022         s->eof_pts = in->pts + in->nb_samples;
1023     }
1024     if (eof)
1025         out->pts = s->eof_pts - s->padd_samples;
1026 
1027     if (!in || in->nb_samples < s->nb_samples) {
1028         AVFrame *new_in = ff_get_audio_buffer(outlink, s->nb_samples);
1029 
1030         if (!new_in) {
1031             av_frame_free(&in);
1032             av_frame_free(&out);
1033             return AVERROR(ENOMEM);
1034         }
1035         if (in)
1036             av_frame_copy_props(new_in, in);
1037 
1038         s->padd_samples -= s->nb_samples - (in ? in->nb_samples: 0);
1039         if (in)
1040             av_samples_copy(new_in->extended_data, in->extended_data, 0, 0,
1041                             in->nb_samples, in->ch_layout.nb_channels, in->format);
1042         av_frame_free(&in);
1043         in = new_in;
1044     }
1045 
1046     td.in  = in;
1047     td.out = out;
1048     ff_filter_execute(ctx, s->filter_channel, &td, NULL, inlink->ch_layout.nb_channels);
1049     if (s->need_profile)
1050         s->got_profile = 1;
1051 
1052     s->sn += s->nb_samples;
1053 
1054     if (s->drop_samples >= in->nb_samples) {
1055         s->drop_samples -= in->nb_samples;
1056         s->delay += in->nb_samples;
1057         av_frame_free(&in);
1058         av_frame_free(&out);
1059         FF_FILTER_FORWARD_STATUS(inlink, outlink);
1060         FF_FILTER_FORWARD_WANTED(outlink, inlink);
1061         return 0;
1062     } else if (s->drop_samples > 0) {
1063         for (int ch = 0; ch < out->ch_layout.nb_channels; ch++) {
1064             memmove(out->extended_data[ch],
1065                     out->extended_data[ch] + s->drop_samples * sizeof(double),
1066                     (in->nb_samples - s->drop_samples) * sizeof(double));
1067         }
1068 
1069         out->nb_samples = in->nb_samples - s->drop_samples;
1070         out->pts = in->pts - av_rescale_q(s->delay, (AVRational){1, outlink->sample_rate}, outlink->time_base);
1071         s->delay += s->drop_samples;
1072         s->drop_samples = 0;
1073     } else {
1074         if (s->padd_samples < 0 && eof) {
1075             out->nb_samples = FFMAX(0, out->nb_samples + s->padd_samples);
1076             s->padd_samples = 0;
1077         }
1078         if (!eof)
1079             out->pts = in->pts - av_rescale_q(s->delay, (AVRational){1, outlink->sample_rate}, outlink->time_base);
1080     }
1081 
1082     av_frame_free(&in);
1083     return ff_filter_frame(outlink, out);
1084 }
1085 
max_left_ext(int wavelet_length,int levels)1086 static int max_left_ext(int wavelet_length, int levels)
1087 {
1088     return (pow2(levels) - 1) * (wavelet_length - 1);
1089 }
1090 
min_left_ext(int wavelet_length,int levels)1091 static int min_left_ext(int wavelet_length, int levels)
1092 {
1093     return (pow2(levels) - 1) * (wavelet_length - 2);
1094 }
1095 
config_output(AVFilterLink * outlink)1096 static int config_output(AVFilterLink *outlink)
1097 {
1098     AVFilterContext *ctx = outlink->src;
1099     AudioFWTDNContext *s = ctx->priv;
1100 
1101     switch (s->wavelet_type) {
1102     case SYM2:
1103         s->wavelet_length = 4;
1104         s->lp  = sym2_lp;
1105         s->hp  = sym2_hp;
1106         s->ilp = sym2_ilp;
1107         s->ihp = sym2_ihp;
1108         break;
1109     case SYM4:
1110         s->wavelet_length = 8;
1111         s->lp  = sym4_lp;
1112         s->hp  = sym4_hp;
1113         s->ilp = sym4_ilp;
1114         s->ihp = sym4_ihp;
1115         break;
1116     case RBIOR68:
1117         s->wavelet_length = 18;
1118         s->lp  = rbior68_lp;
1119         s->hp  = rbior68_hp;
1120         s->ilp = rbior68_ilp;
1121         s->ihp = rbior68_ihp;
1122         break;
1123     case DEB10:
1124         s->wavelet_length = 20;
1125         s->lp  = deb10_lp;
1126         s->hp  = deb10_hp;
1127         s->ilp = deb10_ilp;
1128         s->ihp = deb10_ihp;
1129         break;
1130     case SYM10:
1131         s->wavelet_length = 20;
1132         s->lp  = sym10_lp;
1133         s->hp  = sym10_hp;
1134         s->ilp = sym10_ilp;
1135         s->ihp = sym10_ihp;
1136         break;
1137     case COIF5:
1138         s->wavelet_length = 30;
1139         s->lp  = coif5_lp;
1140         s->hp  = coif5_hp;
1141         s->ilp = coif5_ilp;
1142         s->ihp = coif5_ihp;
1143         break;
1144     case BL3:
1145         s->wavelet_length = 42;
1146         s->lp  = bl3_lp;
1147         s->hp  = bl3_hp;
1148         s->ilp = bl3_ilp;
1149         s->ihp = bl3_ihp;
1150         break;
1151     default:
1152         av_assert0(0);
1153     }
1154 
1155     s->levels = FFMIN(s->levels, lrint(log(s->nb_samples / (s->wavelet_length - 1.0)) / M_LN2));
1156     av_log(ctx, AV_LOG_VERBOSE, "levels: %d\n", s->levels);
1157     s->filter_channel = filter_channel;
1158 
1159     s->stddev = ff_get_audio_buffer(outlink, MAX_LEVELS);
1160     s->new_stddev = ff_get_audio_buffer(outlink, MAX_LEVELS);
1161     s->filter = ff_get_audio_buffer(outlink, s->nb_samples);
1162     s->absmean = ff_get_audio_buffer(outlink, MAX_LEVELS);
1163     s->new_absmean = ff_get_audio_buffer(outlink, MAX_LEVELS);
1164     if (!s->stddev || !s->absmean || !s->filter ||
1165         !s->new_stddev || !s->new_absmean)
1166         return AVERROR(ENOMEM);
1167 
1168     s->channels = outlink->ch_layout.nb_channels;
1169     s->overlap_length = max_left_ext(s->wavelet_length, s->levels);
1170     s->prev_length = s->overlap_length;
1171     s->drop_samples = s->overlap_length;
1172     s->padd_samples = s->overlap_length;
1173     s->sn = 1;
1174 
1175     s->cp = av_calloc(s->channels, sizeof(*s->cp));
1176     if (!s->cp)
1177         return AVERROR(ENOMEM);
1178 
1179     for (int ch = 0; ch < s->channels; ch++) {
1180         ChannelParams *cp = &s->cp[ch];
1181 
1182         cp->output_coefs = av_calloc(s->levels + 1, sizeof(*cp->output_coefs));
1183         cp->filter_coefs = av_calloc(s->levels + 1, sizeof(*cp->filter_coefs));
1184         cp->output_length = av_calloc(s->levels + 1, sizeof(*cp->output_length));
1185         cp->filter_length = av_calloc(s->levels + 1, sizeof(*cp->filter_length));
1186         cp->buffer_length = next_pow2(s->wavelet_length);
1187         cp->buffer = av_calloc(cp->buffer_length, sizeof(*cp->buffer));
1188         cp->buffer2 = av_calloc(cp->buffer_length, sizeof(*cp->buffer2));
1189         cp->subbands_to_free = av_calloc(s->levels + 1, sizeof(*cp->subbands_to_free));
1190         cp->prev = av_calloc(s->prev_length, sizeof(*cp->prev));
1191         cp->overlap = av_calloc(s->overlap_length, sizeof(*cp->overlap));
1192         cp->max_left_ext = max_left_ext(s->wavelet_length, s->levels);
1193         cp->min_left_ext = min_left_ext(s->wavelet_length, s->levels);
1194         if (!cp->output_coefs || !cp->filter_coefs || !cp->output_length ||
1195             !cp->filter_length || !cp->subbands_to_free || !cp->prev || !cp->overlap ||
1196             !cp->buffer || !cp->buffer2)
1197             return AVERROR(ENOMEM);
1198     }
1199 
1200     return 0;
1201 }
1202 
activate(AVFilterContext * ctx)1203 static int activate(AVFilterContext *ctx)
1204 {
1205     AVFilterLink *inlink = ctx->inputs[0];
1206     AVFilterLink *outlink = ctx->outputs[0];
1207     AudioFWTDNContext *s = ctx->priv;
1208     AVFrame *in = NULL;
1209     int ret, status;
1210     int64_t pts;
1211 
1212     FF_FILTER_FORWARD_STATUS_BACK(outlink, inlink);
1213 
1214     if (!s->eof) {
1215         ret = ff_inlink_consume_samples(inlink, s->nb_samples, s->nb_samples, &in);
1216         if (ret < 0)
1217             return ret;
1218         if (ret > 0)
1219             return filter_frame(inlink, in);
1220     }
1221 
1222     if (ff_inlink_acknowledge_status(inlink, &status, &pts)) {
1223         if (status == AVERROR_EOF)
1224             s->eof = 1;
1225     }
1226 
1227     if (s->eof && s->padd_samples != 0) {
1228         return filter_frame(inlink, NULL);
1229     } else if (s->eof) {
1230         ff_outlink_set_status(outlink, AVERROR_EOF, s->eof_pts);
1231         return 0;
1232     }
1233     FF_FILTER_FORWARD_WANTED(outlink, inlink);
1234 
1235     return FFERROR_NOT_READY;
1236 }
1237 
uninit(AVFilterContext * ctx)1238 static av_cold void uninit(AVFilterContext *ctx)
1239 {
1240     AudioFWTDNContext *s = ctx->priv;
1241 
1242     av_frame_free(&s->filter);
1243     av_frame_free(&s->new_stddev);
1244     av_frame_free(&s->stddev);
1245     av_frame_free(&s->new_absmean);
1246     av_frame_free(&s->absmean);
1247 
1248     for (int ch = 0; s->cp && ch < s->channels; ch++) {
1249         ChannelParams *cp = &s->cp[ch];
1250 
1251         av_freep(&cp->tempa);
1252         av_freep(&cp->tempd);
1253         av_freep(&cp->temp_in);
1254         av_freep(&cp->buffer);
1255         av_freep(&cp->buffer2);
1256         av_freep(&cp->prev);
1257         av_freep(&cp->overlap);
1258 
1259         av_freep(&cp->output_length);
1260         av_freep(&cp->filter_length);
1261 
1262         if (cp->output_coefs) {
1263             for (int level = 0; level <= s->levels; level++)
1264                 av_freep(&cp->output_coefs[level]);
1265         }
1266 
1267         if (cp->subbands_to_free) {
1268             for (int level = 0; level <= s->levels; level++)
1269                 av_freep(&cp->subbands_to_free[level]);
1270         }
1271 
1272         av_freep(&cp->subbands_to_free);
1273         av_freep(&cp->output_coefs);
1274         av_freep(&cp->filter_coefs);
1275     }
1276 
1277     av_freep(&s->cp);
1278 }
1279 
process_command(AVFilterContext * ctx,const char * cmd,const char * args,char * res,int res_len,int flags)1280 static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
1281                            char *res, int res_len, int flags)
1282 {
1283     AudioFWTDNContext *s = ctx->priv;
1284     int ret;
1285 
1286     ret = ff_filter_process_command(ctx, cmd, args, res, res_len, flags);
1287     if (ret < 0)
1288         return ret;
1289 
1290     if (!strcmp(cmd, "profile") && s->need_profile)
1291         s->got_profile = 0;
1292 
1293     return 0;
1294 }
1295 
1296 static const AVFilterPad inputs[] = {
1297     {
1298         .name         = "default",
1299         .type         = AVMEDIA_TYPE_AUDIO,
1300     },
1301 };
1302 
1303 static const AVFilterPad outputs[] = {
1304     {
1305         .name          = "default",
1306         .type          = AVMEDIA_TYPE_AUDIO,
1307         .config_props  = config_output,
1308     },
1309 };
1310 
1311 const AVFilter ff_af_afwtdn = {
1312     .name            = "afwtdn",
1313     .description     = NULL_IF_CONFIG_SMALL("Denoise audio stream using Wavelets."),
1314     .priv_size       = sizeof(AudioFWTDNContext),
1315     .priv_class      = &afwtdn_class,
1316     .activate        = activate,
1317     .uninit          = uninit,
1318     FILTER_INPUTS(inputs),
1319     FILTER_OUTPUTS(outputs),
1320     FILTER_SINGLE_SAMPLEFMT(AV_SAMPLE_FMT_DBLP),
1321     .process_command = process_command,
1322     .flags           = AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL |
1323                        AVFILTER_FLAG_SLICE_THREADS,
1324 };
1325