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
27
32
35
40
44
48
50 {
53 const float *coeffs =
s->coeffs;
55 int nb_samples;
56
59
60 nb_samples =
FFMIN(
s->nb_samples,
s->n -
s->pts);
61 if (nb_samples <= 0) {
63 return 0;
64 }
65
68
69 memcpy(
frame->data[0], coeffs +
s->pts, nb_samples *
sizeof(
float));
70
73
75 }
76
80 {
89
93
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;
103
106
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;
116 }
117 }
118
119 for (
i = 0; dc_norm &&
i < num_taps;
i++)
121
123 }
124
126 {
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 / .0005
f) / logf(2.
f);
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 }
149 return .1102f * (att - 8.7f);
151 return .58417f *
powf(att - 20.96
f, .4
f) + .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 * .5
f / Fc): *beta;
158 att = att < 60.f ? (att - 7.95f) / (2.285
f *
M_PI * 2.
f) :
159 ((.0007528358f-1.577737e-05 * *beta) * *beta + 0.6248022
f) * *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;
170 }
171
172 att = att ? att : 120.f;
173
175
176 if (!n) {
177 n = *num_taps;
178 *num_taps =
av_clip(n, 11, 32767);
180 *num_taps = 1 + 2 * (int)((
int)((*num_taps / 2) * Fc + .5
f) / Fc + .5f);
181 }
182
183 return make_lpf(*num_taps |= 1, Fc, *beta, 0.
f, 1.
f, 0);
184 }
185
187 {
188 for (
int i = 0;
i < n;
i++)
190
192 }
193
194 #define SQR(a) ((a) * (a))
195
197 {
199 if (x)
200 return logf(x);
201 return -26;
202 }
203
205 {
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));
218 pi_wraps = &
work[work_len + 2];
219
221
230
231 s->tx_fn(
s->tx,
work,
work,
sizeof(
float));
/* Cepstral: */
232
233 for (
i = 0;
i <= work_len;
i += 2) {
235 float detect = 2 *
M_PI;
236 float delta = angle - prev_angle2;
238
239 prev_angle2 = angle;
241 angle += cum_2pi;
243 delta = angle - prev_angle1;
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
251 }
252
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 */
260 work[
i + work_len / 2] = 0;
261 }
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
269 for (
i = 2;
i < work_len;
i += 2) {
271
274 }
275
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 + .5
f);
i++) {
283 if (
fabs(imp_sum) >
fabs(peak_imp_sum)) {
284 peak_imp_sum = imp_sum;
286 }
287 if (
work[
i] >
work[imp_peak])
/* For debug check only */
289 }
290
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) * .22
f) * *
len + .5f;
301 end = (.997f + (0 - phase1) * .22
f) * *
len + .5f;
302 begin = peak - (begin & ~3);
303 end = peak + 1 + ((end + 3) & ~3);
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
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
323
325 }
326
328 {
331 float Fn =
s->sample_rate * .5f;
333 int i, n, post_peak, longer,
ret;
334
337
338 if (
s->Fc0 >= Fn ||
s->Fc1 >= Fn) {
340 "filter frequency must be less than %d/2.\n",
s->sample_rate);
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
349
350 longer =
s->num_taps[1] >
s->num_taps[0];
351 n =
s->num_taps[longer];
352
354 for (
i = 0;
i <
s->num_taps[!longer];
i++)
355 h[longer][
i + (n -
s->num_taps[!longer]) / 2] +=
h[!longer][
i];
356
359
361 }
362
363 if (
s->phase != 50.f) {
367 } else {
368 post_peak = n >> 1;
369 }
370
377 }
378
379 for (
i = 0;
i < n;
i++)
380 s->coeffs[
i] =
h[longer][
i];
381
385
389 }
390
392 {
394
398 }
399
401 {
405 },
406 };
407
408 #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
409 #define OFFSET(x) offsetof(SincContext, x)
410
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 },
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 },
425 };
426
428
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,
438 };