FFmpeg coverage

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_common_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(SincContext *s, float **h, int *len, int *post_len, float phase)
205			{
206		✗	float *pi_wraps, *work, phase1 = (phase > 50.f ? 100.f - phase : phase) / 50.f;
207		✗	int i, work_len, begin, end, imp_peak = 0, peak = 0, ret;
208		✗	float imp_sum = 0, peak_imp_sum = 0, scale = 1.f;
209		✗	float prev_angle2 = 0, cum_2pi = 0, prev_angle1 = 0, cum_1pi = 0;
210
211		✗	for (i = *len, work_len = 2 * 2 * 8; i > 1; work_len <<= 1, i >>= 1);
212
213			/* The first part is for work (+2 for (UN)PACK), the latter for pi_wraps. */
214		✗	work = av_calloc((work_len + 2) + (work_len / 2 + 1), sizeof(float));
215		✗	if (!work)
216		✗	return AVERROR(ENOMEM);
217		✗	pi_wraps = &work[work_len + 2];
218
219		✗	memcpy(work, *h, *len * sizeof(*work));
220
221		✗	av_tx_uninit(&s->tx);
222		✗	av_tx_uninit(&s->itx);
223		✗	ret = av_tx_init(&s->tx, &s->tx_fn, AV_TX_FLOAT_RDFT, 0, work_len, &scale, AV_TX_INPLACE);
224		✗	if (ret < 0)
225		✗	goto fail;
226		✗	ret = av_tx_init(&s->itx, &s->itx_fn, AV_TX_FLOAT_RDFT, 1, work_len, &scale, AV_TX_INPLACE);
227		✗	if (ret < 0)
228		✗	goto fail;
229
230		✗	s->tx_fn(s->tx, work, work, sizeof(float)); /* Cepstral: */
231
232		✗	for (i = 0; i <= work_len; i += 2) {
233		✗	float angle = atan2f(work[i + 1], work[i]);
234		✗	float detect = 2 * M_PI;
235		✗	float delta = angle - prev_angle2;
236		✗	float adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
237
238		✗	prev_angle2 = angle;
239		✗	cum_2pi += adjust;
240		✗	angle += cum_2pi;
241		✗	detect = M_PI;
242		✗	delta = angle - prev_angle1;
243		✗	adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
244		✗	prev_angle1 = angle;
245		✗	cum_1pi += fabsf(adjust); /* fabs for when 2pi and 1pi have combined */
246		✗	pi_wraps[i >> 1] = cum_1pi;
247
248		✗	work[i] = safe_log(sqrtf(SQR(work[i]) + SQR(work[i + 1])));
249		✗	work[i + 1] = 0;
250			}
251
252		✗	s->itx_fn(s->itx, work, work, sizeof(AVComplexFloat));
253
254		✗	for (i = 0; i < work_len; i++)
255		✗	work[i] *= 2.f / work_len;
256
257		✗	for (i = 1; i < work_len / 2; i++) { /* Window to reject acausal components */
258		✗	work[i] *= 2;
259		✗	work[i + work_len / 2] = 0;
260			}
261		✗	s->tx_fn(s->tx, work, work, sizeof(float));
262
263		✗	for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
264		✗	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];
265
266		✗	work[0] = exp(work[0]);
267		✗	work[1] = exp(work[1]);
268		✗	for (i = 2; i < work_len; i += 2) {
269		✗	float x = expf(work[i]);
270
271		✗	work[i ] = x * cosf(work[i + 1]);
272		✗	work[i + 1] = x * sinf(work[i + 1]);
273			}
274
275		✗	s->itx_fn(s->itx, work, work, sizeof(AVComplexFloat));
276		✗	for (i = 0; i < work_len; i++)
277		✗	work[i] *= 2.f / work_len;
278
279			/* Find peak pos. */
280		✗	for (i = 0; i <= (int) (pi_wraps[work_len >> 1] / M_PI + .5f); i++) {
281		✗	imp_sum += work[i];
282		✗	if (fabs(imp_sum) > fabs(peak_imp_sum)) {
283		✗	peak_imp_sum = imp_sum;
284		✗	peak = i;
285			}
286		✗	if (work[i] > work[imp_peak]) /* For debug check only */
287		✗	imp_peak = i;
288			}
289
290		✗	while (peak && fabsf(work[peak - 1]) > fabsf(work[peak]) && (work[peak - 1] * work[peak] > 0)) {
291		✗	peak--;
292			}
293
294		✗	if (!phase1) {
295		✗	begin = 0;
296		✗	} else if (phase1 == 1) {
297		✗	begin = peak - *len / 2;
298			} else {
299		✗	begin = (.997f - (2 - phase1) * .22f) * *len + .5f;
300		✗	end = (.997f + (0 - phase1) * .22f) * *len + .5f;
301		✗	begin = peak - (begin & ~3);
302		✗	end = peak + 1 + ((end + 3) & ~3);
303		✗	*len = end - begin;
304		✗	*h = av_realloc_f(*h, *len, sizeof(**h));
305		✗	if (!*h) {
306		✗	av_free(work);
307		✗	return AVERROR(ENOMEM);
308			}
309			}
310
311		✗	for (i = 0; i < *len; i++) {
312		✗	(*h)[i] = work[(begin + (phase > 50.f ? *len - 1 - i : i) + work_len) & (work_len - 1)];
313			}
314		✗	*post_len = phase > 50 ? peak - begin : begin + *len - (peak + 1);
315
316		✗	av_log(s, AV_LOG_DEBUG, "%d nPI=%g peak-sum@%i=%g (val@%i=%g); len=%i post=%i (%g%%)\n",
317		✗	work_len, pi_wraps[work_len >> 1] / M_PI, peak, peak_imp_sum, imp_peak,
318		✗	work[imp_peak], *len, *post_len, 100.f - 100.f * *post_len / (*len - 1));
319
320		✗	fail:
321		✗	av_free(work);
322
323		✗	return ret;
324			}
325
326		✗	static int config_output(AVFilterLink *outlink)
327			{
328		✗	AVFilterContext *ctx = outlink->src;
329		✗	SincContext *s = ctx->priv;
330		✗	float Fn = s->sample_rate * .5f;
331			float *h[2];
332			int i, n, post_peak, longer;
333
334		✗	outlink->sample_rate = s->sample_rate;
335		✗	s->pts = 0;
336
337		✗	if (s->Fc0 >= Fn || s->Fc1 >= Fn) {
338		✗	av_log(ctx, AV_LOG_ERROR,
339			"filter frequency must be less than %d/2.\n", s->sample_rate);
340		✗	return AVERROR(EINVAL);
341			}
342
343		✗	h[0] = lpf(Fn, s->Fc0, s->tbw0, &s->num_taps[0], s->att, &s->beta, s->round);
344		✗	h[1] = lpf(Fn, s->Fc1, s->tbw1, &s->num_taps[1], s->att, &s->beta, s->round);
345
346		✗	if (h[0])
347		✗	invert(h[0], s->num_taps[0]);
348
349		✗	longer = s->num_taps[1] > s->num_taps[0];
350		✗	n = s->num_taps[longer];
351
352		✗	if (h[0] && h[1]) {
353		✗	for (i = 0; i < s->num_taps[!longer]; i++)
354		✗	h[longer][i + (n - s->num_taps[!longer]) / 2] += h[!longer][i];
355
356		✗	if (s->Fc0 < s->Fc1)
357		✗	invert(h[longer], n);
358
359		✗	av_free(h[!longer]);
360			}
361
362		✗	if (s->phase != 50.f) {
363		✗	int ret = fir_to_phase(s, &h[longer], &n, &post_peak, s->phase);
364		✗	if (ret < 0)
365		✗	return ret;
366			} else {
367		✗	post_peak = n >> 1;
368			}
369
370		✗	s->n = 1 << (av_log2(n) + 1);
371		✗	s->rdft_len = 1 << av_log2(n);
372		✗	s->coeffs = av_calloc(s->n, sizeof(*s->coeffs));
373		✗	if (!s->coeffs)
374		✗	return AVERROR(ENOMEM);
375
376		✗	for (i = 0; i < n; i++)
377		✗	s->coeffs[i] = h[longer][i];
378		✗	av_free(h[longer]);
379
380		✗	av_tx_uninit(&s->tx);
381		✗	av_tx_uninit(&s->itx);
382
383		✗	return 0;
384			}
385
386		✗	static av_cold void uninit(AVFilterContext *ctx)
387			{
388		✗	SincContext *s = ctx->priv;
389
390		✗	av_freep(&s->coeffs);
391		✗	av_tx_uninit(&s->tx);
392		✗	av_tx_uninit(&s->itx);
393		✗	}
394
395			static const AVFilterPad sinc_outputs[] = {
396			{
397			.name = "default",
398			.type = AVMEDIA_TYPE_AUDIO,
399			.config_props = config_output,
400			},
401			};
402
403			#define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
404			#define OFFSET(x) offsetof(SincContext, x)
405
406			static const AVOption sinc_options[] = {
407			{ "sample_rate", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
408			{ "r", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
409			{ "nb_samples", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
410			{ "n", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
411			{ "hp", "set high-pass filter frequency", OFFSET(Fc0), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
412			{ "lp", "set low-pass filter frequency", OFFSET(Fc1), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
413			{ "phase", "set filter phase response", OFFSET(phase), AV_OPT_TYPE_FLOAT, {.dbl=50}, 0, 100, AF },
414			{ "beta", "set kaiser window beta", OFFSET(beta), AV_OPT_TYPE_FLOAT, {.dbl=-1}, -1, 256, AF },
415			{ "att", "set stop-band attenuation", OFFSET(att), AV_OPT_TYPE_FLOAT, {.dbl=120}, 40, 180, AF },
416			{ "round", "enable rounding", OFFSET(round), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, AF },
417			{ "hptaps", "set number of taps for high-pass filter", OFFSET(num_taps[0]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
418			{ "lptaps", "set number of taps for low-pass filter", OFFSET(num_taps[1]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
419			{ NULL }
420			};
421
422			AVFILTER_DEFINE_CLASS(sinc);
423
424			const AVFilter ff_asrc_sinc = {
425			.name = "sinc",
426			.description = NULL_IF_CONFIG_SMALL("Generate a sinc kaiser-windowed low-pass, high-pass, band-pass, or band-reject FIR coefficients."),
427			.priv_size = sizeof(SincContext),
428			.priv_class = &sinc_class,
429			.uninit = uninit,
430			.activate = activate,
431			.inputs = NULL,
432			FILTER_OUTPUTS(sinc_outputs),
433			FILTER_QUERY_FUNC2(query_formats),
434			};
435