1 /*
2 * Copyright (c) 2005 Boðaç Topaktaþ
3 * Copyright (c) 2020 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/channel_layout.h"
23 #include "libavutil/ffmath.h"
24 #include "libavutil/opt.h"
25 #include "avfilter.h"
26 #include "audio.h"
27 #include "formats.h"
28
29 typedef struct BiquadCoeffs {
30 double a1, a2;
31 double b0, b1, b2;
32 } BiquadCoeffs;
33
34 typedef struct ASuperCutContext {
35 const AVClass *class;
36
37 double cutoff;
38 double level;
39 double qfactor;
40 int order;
41
42 int filter_count;
43 int bypass;
44
45 BiquadCoeffs coeffs[10];
46
47 AVFrame *w;
48
49 int (*filter_channels)(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs);
50 } ASuperCutContext;
51
52 static const enum AVSampleFormat sample_fmts[] = {
53 AV_SAMPLE_FMT_FLTP, AV_SAMPLE_FMT_DBLP, AV_SAMPLE_FMT_NONE
54 };
55
calc_q_factors(int n,double * q)56 static void calc_q_factors(int n, double *q)
57 {
58 for (int i = 0; i < n / 2; i++)
59 q[i] = 1. / (-2. * cos(M_PI * (2. * (i + 1) + n - 1.) / (2. * n)));
60 }
61
get_coeffs(AVFilterContext * ctx)62 static int get_coeffs(AVFilterContext *ctx)
63 {
64 ASuperCutContext *s = ctx->priv;
65 AVFilterLink *inlink = ctx->inputs[0];
66 double w0 = s->cutoff / inlink->sample_rate;
67 double K = tan(M_PI * w0);
68 double q[10];
69
70 s->bypass = w0 >= 0.5;
71 if (s->bypass)
72 return 0;
73
74 if (!strcmp(ctx->filter->name, "asubcut")) {
75 s->filter_count = s->order / 2 + (s->order & 1);
76
77 calc_q_factors(s->order, q);
78
79 if (s->order & 1) {
80 BiquadCoeffs *coeffs = &s->coeffs[0];
81 double omega = 2. * tan(M_PI * w0);
82
83 coeffs->b0 = 2. / (2. + omega);
84 coeffs->b1 = -coeffs->b0;
85 coeffs->b2 = 0.;
86 coeffs->a1 = -(omega - 2.) / (2. + omega);
87 coeffs->a2 = 0.;
88 }
89
90 for (int b = (s->order & 1); b < s->filter_count; b++) {
91 BiquadCoeffs *coeffs = &s->coeffs[b];
92 const int idx = b - (s->order & 1);
93 double norm = 1.0 / (1.0 + K / q[idx] + K * K);
94
95 coeffs->b0 = norm;
96 coeffs->b1 = -2.0 * coeffs->b0;
97 coeffs->b2 = coeffs->b0;
98 coeffs->a1 = -2.0 * (K * K - 1.0) * norm;
99 coeffs->a2 = -(1.0 - K / q[idx] + K * K) * norm;
100 }
101 } else if (!strcmp(ctx->filter->name, "asupercut")) {
102 s->filter_count = s->order / 2 + (s->order & 1);
103
104 calc_q_factors(s->order, q);
105
106 if (s->order & 1) {
107 BiquadCoeffs *coeffs = &s->coeffs[0];
108 double omega = 2. * tan(M_PI * w0);
109
110 coeffs->b0 = omega / (2. + omega);
111 coeffs->b1 = coeffs->b0;
112 coeffs->b2 = 0.;
113 coeffs->a1 = -(omega - 2.) / (2. + omega);
114 coeffs->a2 = 0.;
115 }
116
117 for (int b = (s->order & 1); b < s->filter_count; b++) {
118 BiquadCoeffs *coeffs = &s->coeffs[b];
119 const int idx = b - (s->order & 1);
120 double norm = 1.0 / (1.0 + K / q[idx] + K * K);
121
122 coeffs->b0 = K * K * norm;
123 coeffs->b1 = 2.0 * coeffs->b0;
124 coeffs->b2 = coeffs->b0;
125 coeffs->a1 = -2.0 * (K * K - 1.0) * norm;
126 coeffs->a2 = -(1.0 - K / q[idx] + K * K) * norm;
127 }
128 } else if (!strcmp(ctx->filter->name, "asuperpass")) {
129 double alpha, beta, gamma, theta;
130 double theta_0 = 2. * M_PI * (s->cutoff / inlink->sample_rate);
131 double d_E;
132
133 s->filter_count = s->order / 2;
134 d_E = (2. * tan(theta_0 / (2. * s->qfactor))) / sin(theta_0);
135
136 for (int b = 0; b < s->filter_count; b += 2) {
137 double D = 2. * sin(((b + 1) * M_PI) / (2. * s->filter_count));
138 double A = (1. + pow((d_E / 2.), 2)) / (D * d_E / 2.);
139 double d = sqrt((d_E * D) / (A + sqrt(A * A - 1.)));
140 double B = D * (d_E / 2.) / d;
141 double W = B + sqrt(B * B - 1.);
142
143 for (int j = 0; j < 2; j++) {
144 BiquadCoeffs *coeffs = &s->coeffs[b + j];
145
146 if (j == 1)
147 theta = 2. * atan(tan(theta_0 / 2.) / W);
148 else
149 theta = 2. * atan(W * tan(theta_0 / 2.));
150
151 beta = 0.5 * ((1. - (d / 2.) * sin(theta)) / (1. + (d / 2.) * sin(theta)));
152 gamma = (0.5 + beta) * cos(theta);
153 alpha = 0.5 * (0.5 - beta) * sqrt(1. + pow((W - (1. / W)) / d, 2.));
154
155 coeffs->a1 = 2. * gamma;
156 coeffs->a2 = -2. * beta;
157 coeffs->b0 = 2. * alpha;
158 coeffs->b1 = 0.;
159 coeffs->b2 = -2. * alpha;
160 }
161 }
162 } else if (!strcmp(ctx->filter->name, "asuperstop")) {
163 double alpha, beta, gamma, theta;
164 double theta_0 = 2. * M_PI * (s->cutoff / inlink->sample_rate);
165 double d_E;
166
167 s->filter_count = s->order / 2;
168 d_E = (2. * tan(theta_0 / (2. * s->qfactor))) / sin(theta_0);
169
170 for (int b = 0; b < s->filter_count; b += 2) {
171 double D = 2. * sin(((b + 1) * M_PI) / (2. * s->filter_count));
172 double A = (1. + pow((d_E / 2.), 2)) / (D * d_E / 2.);
173 double d = sqrt((d_E * D) / (A + sqrt(A * A - 1.)));
174 double B = D * (d_E / 2.) / d;
175 double W = B + sqrt(B * B - 1.);
176
177 for (int j = 0; j < 2; j++) {
178 BiquadCoeffs *coeffs = &s->coeffs[b + j];
179
180 if (j == 1)
181 theta = 2. * atan(tan(theta_0 / 2.) / W);
182 else
183 theta = 2. * atan(W * tan(theta_0 / 2.));
184
185 beta = 0.5 * ((1. - (d / 2.) * sin(theta)) / (1. + (d / 2.) * sin(theta)));
186 gamma = (0.5 + beta) * cos(theta);
187 alpha = 0.5 * (0.5 + beta) * ((1. - cos(theta)) / (1. - cos(theta_0)));
188
189 coeffs->a1 = 2. * gamma;
190 coeffs->a2 = -2. * beta;
191 coeffs->b0 = 2. * alpha;
192 coeffs->b1 = -4. * alpha * cos(theta_0);
193 coeffs->b2 = 2. * alpha;
194 }
195 }
196 }
197
198 return 0;
199 }
200
201 typedef struct ThreadData {
202 AVFrame *in, *out;
203 } ThreadData;
204
205 #define FILTER(name, type) \
206 static int filter_channels_## name(AVFilterContext *ctx, void *arg, \
207 int jobnr, int nb_jobs) \
208 { \
209 ASuperCutContext *s = ctx->priv; \
210 ThreadData *td = arg; \
211 AVFrame *out = td->out; \
212 AVFrame *in = td->in; \
213 const int start = (in->ch_layout.nb_channels * jobnr) / nb_jobs; \
214 const int end = (in->ch_layout.nb_channels * (jobnr+1)) / nb_jobs; \
215 const double level = s->level; \
216 \
217 for (int ch = start; ch < end; ch++) { \
218 const type *src = (const type *)in->extended_data[ch]; \
219 type *dst = (type *)out->extended_data[ch]; \
220 \
221 for (int b = 0; b < s->filter_count; b++) { \
222 BiquadCoeffs *coeffs = &s->coeffs[b]; \
223 const type a1 = coeffs->a1; \
224 const type a2 = coeffs->a2; \
225 const type b0 = coeffs->b0; \
226 const type b1 = coeffs->b1; \
227 const type b2 = coeffs->b2; \
228 type *w = ((type *)s->w->extended_data[ch]) + b * 2; \
229 \
230 for (int n = 0; n < in->nb_samples; n++) { \
231 type sin = b ? dst[n] : src[n] * level; \
232 type sout = sin * b0 + w[0]; \
233 \
234 w[0] = b1 * sin + w[1] + a1 * sout; \
235 w[1] = b2 * sin + a2 * sout; \
236 \
237 dst[n] = sout; \
238 } \
239 } \
240 } \
241 \
242 return 0; \
243 }
244
FILTER(fltp,float)245 FILTER(fltp, float)
246 FILTER(dblp, double)
247
248 static int config_input(AVFilterLink *inlink)
249 {
250 AVFilterContext *ctx = inlink->dst;
251 ASuperCutContext *s = ctx->priv;
252
253 switch (inlink->format) {
254 case AV_SAMPLE_FMT_FLTP: s->filter_channels = filter_channels_fltp; break;
255 case AV_SAMPLE_FMT_DBLP: s->filter_channels = filter_channels_dblp; break;
256 }
257
258 s->w = ff_get_audio_buffer(inlink, 2 * 10);
259 if (!s->w)
260 return AVERROR(ENOMEM);
261
262 return get_coeffs(ctx);
263 }
264
filter_frame(AVFilterLink * inlink,AVFrame * in)265 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
266 {
267 AVFilterContext *ctx = inlink->dst;
268 ASuperCutContext *s = ctx->priv;
269 AVFilterLink *outlink = ctx->outputs[0];
270 ThreadData td;
271 AVFrame *out;
272
273 if (s->bypass)
274 return ff_filter_frame(outlink, in);
275
276 if (av_frame_is_writable(in)) {
277 out = in;
278 } else {
279 out = ff_get_audio_buffer(outlink, in->nb_samples);
280 if (!out) {
281 av_frame_free(&in);
282 return AVERROR(ENOMEM);
283 }
284 av_frame_copy_props(out, in);
285 }
286
287 td.in = in; td.out = out;
288 ff_filter_execute(ctx, s->filter_channels, &td, NULL,
289 FFMIN(inlink->ch_layout.nb_channels, ff_filter_get_nb_threads(ctx)));
290
291 if (out != in)
292 av_frame_free(&in);
293 return ff_filter_frame(outlink, out);
294 }
295
process_command(AVFilterContext * ctx,const char * cmd,const char * args,char * res,int res_len,int flags)296 static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
297 char *res, int res_len, int flags)
298 {
299 int ret;
300
301 ret = ff_filter_process_command(ctx, cmd, args, res, res_len, flags);
302 if (ret < 0)
303 return ret;
304
305 return get_coeffs(ctx);
306 }
307
uninit(AVFilterContext * ctx)308 static av_cold void uninit(AVFilterContext *ctx)
309 {
310 ASuperCutContext *s = ctx->priv;
311
312 av_frame_free(&s->w);
313 }
314
315 #define OFFSET(x) offsetof(ASuperCutContext, x)
316 #define FLAGS AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_RUNTIME_PARAM
317
318 static const AVOption asupercut_options[] = {
319 { "cutoff", "set cutoff frequency", OFFSET(cutoff), AV_OPT_TYPE_DOUBLE, {.dbl=20000}, 20000, 192000, FLAGS },
320 { "order", "set filter order", OFFSET(order), AV_OPT_TYPE_INT, {.i64=10}, 3, 20, FLAGS },
321 { "level", "set input level", OFFSET(level), AV_OPT_TYPE_DOUBLE, {.dbl=1.}, 0., 1., FLAGS },
322 { NULL }
323 };
324
325 AVFILTER_DEFINE_CLASS(asupercut);
326
327 static const AVFilterPad inputs[] = {
328 {
329 .name = "default",
330 .type = AVMEDIA_TYPE_AUDIO,
331 .filter_frame = filter_frame,
332 .config_props = config_input,
333 },
334 };
335
336 static const AVFilterPad outputs[] = {
337 {
338 .name = "default",
339 .type = AVMEDIA_TYPE_AUDIO,
340 },
341 };
342
343 const AVFilter ff_af_asupercut = {
344 .name = "asupercut",
345 .description = NULL_IF_CONFIG_SMALL("Cut super frequencies."),
346 .priv_size = sizeof(ASuperCutContext),
347 .priv_class = &asupercut_class,
348 .uninit = uninit,
349 FILTER_INPUTS(inputs),
350 FILTER_OUTPUTS(outputs),
351 FILTER_SAMPLEFMTS_ARRAY(sample_fmts),
352 .process_command = process_command,
353 .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
354 AVFILTER_FLAG_SLICE_THREADS,
355 };
356
357 static const AVOption asubcut_options[] = {
358 { "cutoff", "set cutoff frequency", OFFSET(cutoff), AV_OPT_TYPE_DOUBLE, {.dbl=20}, 2, 200, FLAGS },
359 { "order", "set filter order", OFFSET(order), AV_OPT_TYPE_INT, {.i64=10}, 3, 20, FLAGS },
360 { "level", "set input level", OFFSET(level), AV_OPT_TYPE_DOUBLE, {.dbl=1.}, 0., 1., FLAGS },
361 { NULL }
362 };
363
364 AVFILTER_DEFINE_CLASS(asubcut);
365
366 const AVFilter ff_af_asubcut = {
367 .name = "asubcut",
368 .description = NULL_IF_CONFIG_SMALL("Cut subwoofer frequencies."),
369 .priv_size = sizeof(ASuperCutContext),
370 .priv_class = &asubcut_class,
371 .uninit = uninit,
372 FILTER_INPUTS(inputs),
373 FILTER_OUTPUTS(outputs),
374 FILTER_SAMPLEFMTS_ARRAY(sample_fmts),
375 .process_command = process_command,
376 .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
377 AVFILTER_FLAG_SLICE_THREADS,
378 };
379
380 static const AVOption asuperpass_asuperstop_options[] = {
381 { "centerf","set center frequency", OFFSET(cutoff), AV_OPT_TYPE_DOUBLE, {.dbl=1000}, 2, 999999, FLAGS },
382 { "order", "set filter order", OFFSET(order), AV_OPT_TYPE_INT, {.i64=4}, 4, 20, FLAGS },
383 { "qfactor","set Q-factor", OFFSET(qfactor),AV_OPT_TYPE_DOUBLE, {.dbl=1.},0.01, 100., FLAGS },
384 { "level", "set input level", OFFSET(level), AV_OPT_TYPE_DOUBLE, {.dbl=1.}, 0., 2., FLAGS },
385 { NULL }
386 };
387
388 AVFILTER_DEFINE_CLASS_EXT(asuperpass_asuperstop, "asuperpass/asuperstop",
389 asuperpass_asuperstop_options);
390
391 const AVFilter ff_af_asuperpass = {
392 .name = "asuperpass",
393 .description = NULL_IF_CONFIG_SMALL("Apply high order Butterworth band-pass filter."),
394 .priv_class = &asuperpass_asuperstop_class,
395 .priv_size = sizeof(ASuperCutContext),
396 .uninit = uninit,
397 FILTER_INPUTS(inputs),
398 FILTER_OUTPUTS(outputs),
399 FILTER_SAMPLEFMTS_ARRAY(sample_fmts),
400 .process_command = process_command,
401 .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
402 AVFILTER_FLAG_SLICE_THREADS,
403 };
404
405 const AVFilter ff_af_asuperstop = {
406 .name = "asuperstop",
407 .description = NULL_IF_CONFIG_SMALL("Apply high order Butterworth band-stop filter."),
408 .priv_class = &asuperpass_asuperstop_class,
409 .priv_size = sizeof(ASuperCutContext),
410 .uninit = uninit,
411 FILTER_INPUTS(inputs),
412 FILTER_OUTPUTS(outputs),
413 FILTER_SAMPLEFMTS_ARRAY(sample_fmts),
414 .process_command = process_command,
415 .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
416 AVFILTER_FLAG_SLICE_THREADS,
417 };
418