FFmpeg coverage


Directory: ../../../ffmpeg/
File: src/libavfilter/asrc_sinc.c
Date: 2022-12-06 04:36:11
Exec Total Coverage
Lines: 0 217 0.0%
Functions: 0 12 0.0%
Branches: 0 126 0.0%

Line Branch Exec Source
1 /*
2 * Copyright (c) 2008-2009 Rob Sykes <robs@users.sourceforge.net>
3 * Copyright (c) 2017 Paul B Mahol
4 *
5 * This file is part of FFmpeg.
6 *
7 * FFmpeg is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * FFmpeg is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with FFmpeg; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22 #include "libavutil/avassert.h"
23 #include "libavutil/channel_layout.h"
24 #include "libavutil/opt.h"
25 #include "libavutil/tx.h"
26
27 #include "audio.h"
28 #include "avfilter.h"
29 #include "filters.h"
30 #include "internal.h"
31
32 typedef struct SincContext {
33 const AVClass *class;
34
35 int sample_rate, nb_samples;
36 float att, beta, phase, Fc0, Fc1, tbw0, tbw1;
37 int num_taps[2];
38 int round;
39
40 int n, rdft_len;
41 float *coeffs;
42 int64_t pts;
43
44 AVTXContext *tx, *itx;
45 av_tx_fn tx_fn, itx_fn;
46 } SincContext;
47
48 static int activate(AVFilterContext *ctx)
49 {
50 AVFilterLink *outlink = ctx->outputs[0];
51 SincContext *s = ctx->priv;
52 const float *coeffs = s->coeffs;
53 AVFrame *frame = NULL;
54 int nb_samples;
55
56 if (!ff_outlink_frame_wanted(outlink))
57 return FFERROR_NOT_READY;
58
59 nb_samples = FFMIN(s->nb_samples, s->n - s->pts);
60 if (nb_samples <= 0) {
61 ff_outlink_set_status(outlink, AVERROR_EOF, s->pts);
62 return 0;
63 }
64
65 if (!(frame = ff_get_audio_buffer(outlink, nb_samples)))
66 return AVERROR(ENOMEM);
67
68 memcpy(frame->data[0], coeffs + s->pts, nb_samples * sizeof(float));
69
70 frame->pts = s->pts;
71 s->pts += nb_samples;
72
73 return ff_filter_frame(outlink, frame);
74 }
75
76 static int query_formats(AVFilterContext *ctx)
77 {
78 SincContext *s = ctx->priv;
79 static const AVChannelLayout chlayouts[] = { AV_CHANNEL_LAYOUT_MONO, { 0 } };
80 int sample_rates[] = { s->sample_rate, -1 };
81 static const enum AVSampleFormat sample_fmts[] = { AV_SAMPLE_FMT_FLT,
82 AV_SAMPLE_FMT_NONE };
83 int ret = ff_set_common_formats_from_list(ctx, sample_fmts);
84 if (ret < 0)
85 return ret;
86
87 ret = ff_set_common_channel_layouts_from_list(ctx, chlayouts);
88 if (ret < 0)
89 return ret;
90
91 return ff_set_common_samplerates_from_list(ctx, sample_rates);
92 }
93
94 static float bessel_I_0(float x)
95 {
96 float term = 1, sum = 1, last_sum, x2 = x / 2;
97 int i = 1;
98
99 do {
100 float y = x2 / i++;
101
102 last_sum = sum;
103 sum += term *= y * y;
104 } while (sum != last_sum);
105
106 return sum;
107 }
108
109 static float *make_lpf(int num_taps, float Fc, float beta, float rho,
110 float scale, int dc_norm)
111 {
112 int i, m = num_taps - 1;
113 float *h = av_calloc(num_taps, sizeof(*h)), sum = 0;
114 float mult = scale / bessel_I_0(beta), mult1 = 1.f / (.5f * m + rho);
115
116 if (!h)
117 return NULL;
118
119 av_assert0(Fc >= 0 && Fc <= 1);
120
121 for (i = 0; i <= m / 2; i++) {
122 float z = i - .5f * m, x = z * M_PI, y = z * mult1;
123 h[i] = x ? sinf(Fc * x) / x : Fc;
124 sum += h[i] *= bessel_I_0(beta * sqrtf(1.f - y * y)) * mult;
125 if (m - i != i) {
126 h[m - i] = h[i];
127 sum += h[i];
128 }
129 }
130
131 for (i = 0; dc_norm && i < num_taps; i++)
132 h[i] *= scale / sum;
133
134 return h;
135 }
136
137 static float kaiser_beta(float att, float tr_bw)
138 {
139 if (att >= 60.f) {
140 static const float coefs[][4] = {
141 {-6.784957e-10, 1.02856e-05, 0.1087556, -0.8988365 + .001},
142 {-6.897885e-10, 1.027433e-05, 0.10876, -0.8994658 + .002},
143 {-1.000683e-09, 1.030092e-05, 0.1087677, -0.9007898 + .003},
144 {-3.654474e-10, 1.040631e-05, 0.1087085, -0.8977766 + .006},
145 {8.106988e-09, 6.983091e-06, 0.1091387, -0.9172048 + .015},
146 {9.519571e-09, 7.272678e-06, 0.1090068, -0.9140768 + .025},
147 {-5.626821e-09, 1.342186e-05, 0.1083999, -0.9065452 + .05},
148 {-9.965946e-08, 5.073548e-05, 0.1040967, -0.7672778 + .085},
149 {1.604808e-07, -5.856462e-05, 0.1185998, -1.34824 + .1},
150 {-1.511964e-07, 6.363034e-05, 0.1064627, -0.9876665 + .18},
151 };
152 float realm = logf(tr_bw / .0005f) / logf(2.f);
153 float const *c0 = coefs[av_clip((int)realm, 0, FF_ARRAY_ELEMS(coefs) - 1)];
154 float const *c1 = coefs[av_clip(1 + (int)realm, 0, FF_ARRAY_ELEMS(coefs) - 1)];
155 float b0 = ((c0[0] * att + c0[1]) * att + c0[2]) * att + c0[3];
156 float b1 = ((c1[0] * att + c1[1]) * att + c1[2]) * att + c1[3];
157
158 return b0 + (b1 - b0) * (realm - (int)realm);
159 }
160 if (att > 50.f)
161 return .1102f * (att - 8.7f);
162 if (att > 20.96f)
163 return .58417f * powf(att - 20.96f, .4f) + .07886f * (att - 20.96f);
164 return 0;
165 }
166
167 static void kaiser_params(float att, float Fc, float tr_bw, float *beta, int *num_taps)
168 {
169 *beta = *beta < 0.f ? kaiser_beta(att, tr_bw * .5f / Fc): *beta;
170 att = att < 60.f ? (att - 7.95f) / (2.285f * M_PI * 2.f) :
171 ((.0007528358f-1.577737e-05 * *beta) * *beta + 0.6248022f) * *beta + .06186902f;
172 *num_taps = !*num_taps ? ceilf(att/tr_bw + 1) : *num_taps;
173 }
174
175 static float *lpf(float Fn, float Fc, float tbw, int *num_taps, float att, float *beta, int round)
176 {
177 int n = *num_taps;
178
179 if ((Fc /= Fn) <= 0.f || Fc >= 1.f) {
180 *num_taps = 0;
181 return NULL;
182 }
183
184 att = att ? att : 120.f;
185
186 kaiser_params(att, Fc, (tbw ? tbw / Fn : .05f) * .5f, beta, num_taps);
187
188 if (!n) {
189 n = *num_taps;
190 *num_taps = av_clip(n, 11, 32767);
191 if (round)
192 *num_taps = 1 + 2 * (int)((int)((*num_taps / 2) * Fc + .5f) / Fc + .5f);
193 }
194
195 return make_lpf(*num_taps |= 1, Fc, *beta, 0.f, 1.f, 0);
196 }
197
198 static void invert(float *h, int n)
199 {
200 for (int i = 0; i < n; i++)
201 h[i] = -h[i];
202
203 h[(n - 1) / 2] += 1;
204 }
205
206 #define SQR(a) ((a) * (a))
207
208 static float safe_log(float x)
209 {
210 av_assert0(x >= 0);
211 if (x)
212 return logf(x);
213 return -26;
214 }
215
216 static int fir_to_phase(SincContext *s, float **h, int *len, int *post_len, float phase)
217 {
218 float *pi_wraps, *work, phase1 = (phase > 50.f ? 100.f - phase : phase) / 50.f;
219 int i, work_len, begin, end, imp_peak = 0, peak = 0;
220 float imp_sum = 0, peak_imp_sum = 0, scale = 1.f;
221 float prev_angle2 = 0, cum_2pi = 0, prev_angle1 = 0, cum_1pi = 0;
222
223 for (i = *len, work_len = 2 * 2 * 8; i > 1; work_len <<= 1, i >>= 1);
224
225 /* The first part is for work (+2 for (UN)PACK), the latter for pi_wraps. */
226 work = av_calloc((work_len + 2) + (work_len / 2 + 1), sizeof(float));
227 if (!work)
228 return AVERROR(ENOMEM);
229 pi_wraps = &work[work_len + 2];
230
231 memcpy(work, *h, *len * sizeof(*work));
232
233 av_tx_uninit(&s->tx);
234 av_tx_uninit(&s->itx);
235 av_tx_init(&s->tx, &s->tx_fn, AV_TX_FLOAT_RDFT, 0, work_len, &scale, AV_TX_INPLACE);
236 av_tx_init(&s->itx, &s->itx_fn, AV_TX_FLOAT_RDFT, 1, work_len, &scale, AV_TX_INPLACE);
237 if (!s->tx || !s->itx) {
238 av_free(work);
239 return AVERROR(ENOMEM);
240 }
241
242 s->tx_fn(s->tx, work, work, sizeof(float)); /* Cepstral: */
243
244 for (i = 0; i <= work_len; i += 2) {
245 float angle = atan2f(work[i + 1], work[i]);
246 float detect = 2 * M_PI;
247 float delta = angle - prev_angle2;
248 float adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
249
250 prev_angle2 = angle;
251 cum_2pi += adjust;
252 angle += cum_2pi;
253 detect = M_PI;
254 delta = angle - prev_angle1;
255 adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
256 prev_angle1 = angle;
257 cum_1pi += fabsf(adjust); /* fabs for when 2pi and 1pi have combined */
258 pi_wraps[i >> 1] = cum_1pi;
259
260 work[i] = safe_log(sqrtf(SQR(work[i]) + SQR(work[i + 1])));
261 work[i + 1] = 0;
262 }
263
264 s->itx_fn(s->itx, work, work, sizeof(AVComplexFloat));
265
266 for (i = 0; i < work_len; i++)
267 work[i] *= 2.f / work_len;
268
269 for (i = 1; i < work_len / 2; i++) { /* Window to reject acausal components */
270 work[i] *= 2;
271 work[i + work_len / 2] = 0;
272 }
273 s->tx_fn(s->tx, work, work, sizeof(float));
274
275 for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
276 work[i + 1] = phase1 * i / work_len * pi_wraps[work_len >> 1] + (1 - phase1) * (work[i + 1] + pi_wraps[i >> 1]) - pi_wraps[i >> 1];
277
278 work[0] = exp(work[0]);
279 work[1] = exp(work[1]);
280 for (i = 2; i < work_len; i += 2) {
281 float x = expf(work[i]);
282
283 work[i ] = x * cosf(work[i + 1]);
284 work[i + 1] = x * sinf(work[i + 1]);
285 }
286
287 s->itx_fn(s->itx, work, work, sizeof(AVComplexFloat));
288 for (i = 0; i < work_len; i++)
289 work[i] *= 2.f / work_len;
290
291 /* Find peak pos. */
292 for (i = 0; i <= (int) (pi_wraps[work_len >> 1] / M_PI + .5f); i++) {
293 imp_sum += work[i];
294 if (fabs(imp_sum) > fabs(peak_imp_sum)) {
295 peak_imp_sum = imp_sum;
296 peak = i;
297 }
298 if (work[i] > work[imp_peak]) /* For debug check only */
299 imp_peak = i;
300 }
301
302 while (peak && fabsf(work[peak - 1]) > fabsf(work[peak]) && (work[peak - 1] * work[peak] > 0)) {
303 peak--;
304 }
305
306 if (!phase1) {
307 begin = 0;
308 } else if (phase1 == 1) {
309 begin = peak - *len / 2;
310 } else {
311 begin = (.997f - (2 - phase1) * .22f) * *len + .5f;
312 end = (.997f + (0 - phase1) * .22f) * *len + .5f;
313 begin = peak - (begin & ~3);
314 end = peak + 1 + ((end + 3) & ~3);
315 *len = end - begin;
316 *h = av_realloc_f(*h, *len, sizeof(**h));
317 if (!*h) {
318 av_free(work);
319 return AVERROR(ENOMEM);
320 }
321 }
322
323 for (i = 0; i < *len; i++) {
324 (*h)[i] = work[(begin + (phase > 50.f ? *len - 1 - i : i) + work_len) & (work_len - 1)];
325 }
326 *post_len = phase > 50 ? peak - begin : begin + *len - (peak + 1);
327
328 av_log(s, AV_LOG_DEBUG, "%d nPI=%g peak-sum@%i=%g (val@%i=%g); len=%i post=%i (%g%%)\n",
329 work_len, pi_wraps[work_len >> 1] / M_PI, peak, peak_imp_sum, imp_peak,
330 work[imp_peak], *len, *post_len, 100.f - 100.f * *post_len / (*len - 1));
331
332 av_free(work);
333
334 return 0;
335 }
336
337 static int config_output(AVFilterLink *outlink)
338 {
339 AVFilterContext *ctx = outlink->src;
340 SincContext *s = ctx->priv;
341 float Fn = s->sample_rate * .5f;
342 float *h[2];
343 int i, n, post_peak, longer;
344
345 outlink->sample_rate = s->sample_rate;
346 s->pts = 0;
347
348 if (s->Fc0 >= Fn || s->Fc1 >= Fn) {
349 av_log(ctx, AV_LOG_ERROR,
350 "filter frequency must be less than %d/2.\n", s->sample_rate);
351 return AVERROR(EINVAL);
352 }
353
354 h[0] = lpf(Fn, s->Fc0, s->tbw0, &s->num_taps[0], s->att, &s->beta, s->round);
355 h[1] = lpf(Fn, s->Fc1, s->tbw1, &s->num_taps[1], s->att, &s->beta, s->round);
356
357 if (h[0])
358 invert(h[0], s->num_taps[0]);
359
360 longer = s->num_taps[1] > s->num_taps[0];
361 n = s->num_taps[longer];
362
363 if (h[0] && h[1]) {
364 for (i = 0; i < s->num_taps[!longer]; i++)
365 h[longer][i + (n - s->num_taps[!longer]) / 2] += h[!longer][i];
366
367 if (s->Fc0 < s->Fc1)
368 invert(h[longer], n);
369
370 av_free(h[!longer]);
371 }
372
373 if (s->phase != 50.f) {
374 int ret = fir_to_phase(s, &h[longer], &n, &post_peak, s->phase);
375 if (ret < 0)
376 return ret;
377 } else {
378 post_peak = n >> 1;
379 }
380
381 s->n = 1 << (av_log2(n) + 1);
382 s->rdft_len = 1 << av_log2(n);
383 s->coeffs = av_calloc(s->n, sizeof(*s->coeffs));
384 if (!s->coeffs)
385 return AVERROR(ENOMEM);
386
387 for (i = 0; i < n; i++)
388 s->coeffs[i] = h[longer][i];
389 av_free(h[longer]);
390
391 av_tx_uninit(&s->tx);
392 av_tx_uninit(&s->itx);
393
394 return 0;
395 }
396
397 static av_cold void uninit(AVFilterContext *ctx)
398 {
399 SincContext *s = ctx->priv;
400
401 av_freep(&s->coeffs);
402 av_tx_uninit(&s->tx);
403 av_tx_uninit(&s->itx);
404 }
405
406 static const AVFilterPad sinc_outputs[] = {
407 {
408 .name = "default",
409 .type = AVMEDIA_TYPE_AUDIO,
410 .config_props = config_output,
411 },
412 };
413
414 #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
415 #define OFFSET(x) offsetof(SincContext, x)
416
417 static const AVOption sinc_options[] = {
418 { "sample_rate", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
419 { "r", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
420 { "nb_samples", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
421 { "n", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
422 { "hp", "set high-pass filter frequency", OFFSET(Fc0), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
423 { "lp", "set low-pass filter frequency", OFFSET(Fc1), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
424 { "phase", "set filter phase response", OFFSET(phase), AV_OPT_TYPE_FLOAT, {.dbl=50}, 0, 100, AF },
425 { "beta", "set kaiser window beta", OFFSET(beta), AV_OPT_TYPE_FLOAT, {.dbl=-1}, -1, 256, AF },
426 { "att", "set stop-band attenuation", OFFSET(att), AV_OPT_TYPE_FLOAT, {.dbl=120}, 40, 180, AF },
427 { "round", "enable rounding", OFFSET(round), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, AF },
428 { "hptaps", "set number of taps for high-pass filter", OFFSET(num_taps[0]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
429 { "lptaps", "set number of taps for low-pass filter", OFFSET(num_taps[1]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
430 { NULL }
431 };
432
433 AVFILTER_DEFINE_CLASS(sinc);
434
435 const AVFilter ff_asrc_sinc = {
436 .name = "sinc",
437 .description = NULL_IF_CONFIG_SMALL("Generate a sinc kaiser-windowed low-pass, high-pass, band-pass, or band-reject FIR coefficients."),
438 .priv_size = sizeof(SincContext),
439 .priv_class = &sinc_class,
440 .uninit = uninit,
441 .activate = activate,
442 .inputs = NULL,
443 FILTER_OUTPUTS(sinc_outputs),
444 FILTER_QUERY_FUNC(query_formats),
445 };
446