FFmpeg coverage

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