]> git.sesse.net Git - ffmpeg/blob - libavfilter/af_afreqshift.c
avfilter/af_acrossover: simplify coefficients
[ffmpeg] / libavfilter / af_afreqshift.c
1 /*
2  * Copyright (c) Paul B Mahol
3  * Copyright (c) Laurent de Soras, 2005
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 #define NB_COEFS 16
30
31 typedef struct AFreqShift {
32     const AVClass *class;
33
34     double shift;
35
36     double c[NB_COEFS];
37
38     int64_t in_samples;
39
40     AVFrame *i1, *o1;
41     AVFrame *i2, *o2;
42
43     void (*filter_channel)(AVFilterContext *ctx,
44                            int nb_samples,
45                            int sample_rate,
46                            const double *src, double *dst,
47                            double *i1, double *o1,
48                            double *i2, double *o2);
49 } AFreqShift;
50
51 static int query_formats(AVFilterContext *ctx)
52 {
53     AVFilterFormats *formats = NULL;
54     AVFilterChannelLayouts *layouts = NULL;
55     static const enum AVSampleFormat sample_fmts[] = {
56         AV_SAMPLE_FMT_DBLP,
57         AV_SAMPLE_FMT_NONE
58     };
59     int ret;
60
61     formats = ff_make_format_list(sample_fmts);
62     if (!formats)
63         return AVERROR(ENOMEM);
64     ret = ff_set_common_formats(ctx, formats);
65     if (ret < 0)
66         return ret;
67
68     layouts = ff_all_channel_counts();
69     if (!layouts)
70         return AVERROR(ENOMEM);
71
72     ret = ff_set_common_channel_layouts(ctx, layouts);
73     if (ret < 0)
74         return ret;
75
76     formats = ff_all_samplerates();
77     return ff_set_common_samplerates(ctx, formats);
78 }
79
80 static void pfilter_channel(AVFilterContext *ctx,
81                             int nb_samples,
82                             int sample_rate,
83                             const double *src, double *dst,
84                             double *i1, double *o1,
85                             double *i2, double *o2)
86 {
87     AFreqShift *s = ctx->priv;
88     double *c = s->c;
89     double shift = s->shift * M_PI;
90     double cos_theta = cos(shift);
91     double sin_theta = sin(shift);
92
93     for (int n = 0; n < nb_samples; n++) {
94         double xn1 = src[n], xn2 = src[n];
95         double I, Q;
96
97         for (int j = 0; j < NB_COEFS / 2; j++) {
98             I = c[j] * (xn1 + o2[j]) - i2[j];
99             i2[j] = i1[j];
100             i1[j] = xn1;
101             o2[j] = o1[j];
102             o1[j] = I;
103             xn1 = I;
104         }
105
106         for (int j = NB_COEFS / 2; j < NB_COEFS; j++) {
107             Q = c[j] * (xn2 + o2[j]) - i2[j];
108             i2[j] = i1[j];
109             i1[j] = xn2;
110             o2[j] = o1[j];
111             o1[j] = Q;
112             xn2 = Q;
113         }
114         Q = o2[NB_COEFS - 1];
115
116         dst[n] = I * cos_theta - Q * sin_theta;
117     }
118 }
119
120 static void ffilter_channel(AVFilterContext *ctx,
121                             int nb_samples,
122                             int sample_rate,
123                             const double *src, double *dst,
124                             double *i1, double *o1,
125                             double *i2, double *o2)
126 {
127     AFreqShift *s = ctx->priv;
128     double *c = s->c;
129     double ts = 1. / sample_rate;
130     double shift = s->shift;
131     int64_t N = s->in_samples;
132
133     for (int n = 0; n < nb_samples; n++) {
134         double xn1 = src[n], xn2 = src[n];
135         double I, Q, theta;
136
137         for (int j = 0; j < NB_COEFS / 2; j++) {
138             I = c[j] * (xn1 + o2[j]) - i2[j];
139             i2[j] = i1[j];
140             i1[j] = xn1;
141             o2[j] = o1[j];
142             o1[j] = I;
143             xn1 = I;
144         }
145
146         for (int j = NB_COEFS / 2; j < NB_COEFS; j++) {
147             Q = c[j] * (xn2 + o2[j]) - i2[j];
148             i2[j] = i1[j];
149             i1[j] = xn2;
150             o2[j] = o1[j];
151             o1[j] = Q;
152             xn2 = Q;
153         }
154         Q = o2[NB_COEFS - 1];
155
156         theta = 2. * M_PI * fmod(shift * (N + n) * ts, 1.);
157         dst[n] = I * cos(theta) - Q * sin(theta);
158     }
159 }
160
161 static void compute_transition_param(double *K, double *Q, double transition)
162 {
163     double kksqrt, e, e2, e4, k, q;
164
165     k  = tan((1. - transition * 2.) * M_PI / 4.);
166     k *= k;
167     kksqrt = pow(1 - k * k, 0.25);
168     e = 0.5 * (1. - kksqrt) / (1. + kksqrt);
169     e2 = e * e;
170     e4 = e2 * e2;
171     q = e * (1. + e4 * (2. + e4 * (15. + 150. * e4)));
172
173     *Q = q;
174     *K = k;
175 }
176
177 static double ipowp(double x, int64_t n)
178 {
179     double z = 1.;
180
181     while (n != 0) {
182         if (n & 1)
183             z *= x;
184         n >>= 1;
185         x *= x;
186     }
187
188     return z;
189 }
190
191 static double compute_acc_num(double q, int order, int c)
192 {
193     int64_t i = 0;
194     int j = 1;
195     double acc = 0.;
196     double q_ii1;
197
198     do {
199         q_ii1  = ipowp(q, i * (i + 1));
200         q_ii1 *= sin((i * 2 + 1) * c * M_PI / order) * j;
201         acc   += q_ii1;
202
203         j = -j;
204         i++;
205     } while (fabs(q_ii1) > 1e-100);
206
207     return acc;
208 }
209
210 static double compute_acc_den(double q, int order, int c)
211 {
212     int64_t i = 1;
213     int j = -1;
214     double acc = 0.;
215     double q_i2;
216
217     do {
218         q_i2  = ipowp(q, i * i);
219         q_i2 *= cos(i * 2 * c * M_PI / order) * j;
220         acc  += q_i2;
221
222         j = -j;
223         i++;
224     } while (fabs(q_i2) > 1e-100);
225
226     return acc;
227 }
228
229 static double compute_coef(int index, double k, double q, int order)
230 {
231     const int    c    = index + 1;
232     const double num  = compute_acc_num(q, order, c) * pow(q, 0.25);
233     const double den  = compute_acc_den(q, order, c) + 0.5;
234     const double ww   = num / den;
235     const double wwsq = ww * ww;
236
237     const double x    = sqrt((1 - wwsq * k) * (1 - wwsq / k)) / (1 + wwsq);
238     const double coef = (1 - x) / (1 + x);
239
240     return coef;
241 }
242
243 static void compute_coefs(double *coef_arr, int nbr_coefs, double transition)
244 {
245     const int order = nbr_coefs * 2 + 1;
246     double k, q;
247
248     compute_transition_param(&k, &q, transition);
249
250     for (int n = 0; n < nbr_coefs; n++)
251         coef_arr[(n / 2) + (n & 1) * nbr_coefs / 2] = compute_coef(n, k, q, order);
252 }
253
254 static int config_input(AVFilterLink *inlink)
255 {
256     AVFilterContext *ctx = inlink->dst;
257     AFreqShift *s = ctx->priv;
258
259     compute_coefs(s->c, NB_COEFS, 2. * 20. / inlink->sample_rate);
260
261     s->i1 = ff_get_audio_buffer(inlink, NB_COEFS);
262     s->o1 = ff_get_audio_buffer(inlink, NB_COEFS);
263     s->i2 = ff_get_audio_buffer(inlink, NB_COEFS);
264     s->o2 = ff_get_audio_buffer(inlink, NB_COEFS);
265     if (!s->i1 || !s->o1 || !s->i2 || !s->o2)
266         return AVERROR(ENOMEM);
267
268     if (!strcmp(ctx->filter->name, "afreqshift"))
269         s->filter_channel = ffilter_channel;
270     else
271         s->filter_channel = pfilter_channel;
272
273     return 0;
274 }
275
276 typedef struct ThreadData {
277     AVFrame *in, *out;
278 } ThreadData;
279
280 static int filter_channels(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
281 {
282     AFreqShift *s = ctx->priv;
283     ThreadData *td = arg;
284     AVFrame *out = td->out;
285     AVFrame *in = td->in;
286     const int start = (in->channels * jobnr) / nb_jobs;
287     const int end = (in->channels * (jobnr+1)) / nb_jobs;
288
289     for (int ch = start; ch < end; ch++) {
290         s->filter_channel(ctx, in->nb_samples,
291                           in->sample_rate,
292                           (const double *)in->extended_data[ch],
293                           (double *)out->extended_data[ch],
294                           (double *)s->i1->extended_data[ch],
295                           (double *)s->o1->extended_data[ch],
296                           (double *)s->i2->extended_data[ch],
297                           (double *)s->o2->extended_data[ch]);
298     }
299
300     return 0;
301 }
302
303 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
304 {
305     AVFilterContext *ctx = inlink->dst;
306     AVFilterLink *outlink = ctx->outputs[0];
307     AFreqShift *s = ctx->priv;
308     AVFrame *out;
309     ThreadData td;
310
311     if (av_frame_is_writable(in)) {
312         out = in;
313     } else {
314         out = ff_get_audio_buffer(outlink, in->nb_samples);
315         if (!out) {
316             av_frame_free(&in);
317             return AVERROR(ENOMEM);
318         }
319         av_frame_copy_props(out, in);
320     }
321
322     td.in = in; td.out = out;
323     ctx->internal->execute(ctx, filter_channels, &td, NULL, FFMIN(inlink->channels,
324                                                             ff_filter_get_nb_threads(ctx)));
325
326     s->in_samples += in->nb_samples;
327
328     if (out != in)
329         av_frame_free(&in);
330     return ff_filter_frame(outlink, out);
331 }
332
333 static av_cold void uninit(AVFilterContext *ctx)
334 {
335     AFreqShift *s = ctx->priv;
336
337     av_frame_free(&s->i1);
338     av_frame_free(&s->o1);
339     av_frame_free(&s->i2);
340     av_frame_free(&s->o2);
341 }
342
343 #define OFFSET(x) offsetof(AFreqShift, x)
344 #define FLAGS AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_RUNTIME_PARAM
345
346 static const AVOption afreqshift_options[] = {
347     { "shift", "set frequency shift", OFFSET(shift), AV_OPT_TYPE_DOUBLE, {.dbl=0}, -INT_MAX, INT_MAX, FLAGS },
348     { NULL }
349 };
350
351 AVFILTER_DEFINE_CLASS(afreqshift);
352
353 static const AVFilterPad inputs[] = {
354     {
355         .name         = "default",
356         .type         = AVMEDIA_TYPE_AUDIO,
357         .filter_frame = filter_frame,
358         .config_props = config_input,
359     },
360     { NULL }
361 };
362
363 static const AVFilterPad outputs[] = {
364     {
365         .name = "default",
366         .type = AVMEDIA_TYPE_AUDIO,
367     },
368     { NULL }
369 };
370
371 AVFilter ff_af_afreqshift = {
372     .name            = "afreqshift",
373     .description     = NULL_IF_CONFIG_SMALL("Apply frequency shifting to input audio."),
374     .query_formats   = query_formats,
375     .priv_size       = sizeof(AFreqShift),
376     .priv_class      = &afreqshift_class,
377     .uninit          = uninit,
378     .inputs          = inputs,
379     .outputs         = outputs,
380     .process_command = ff_filter_process_command,
381     .flags           = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
382                        AVFILTER_FLAG_SLICE_THREADS,
383 };
384
385 static const AVOption aphaseshift_options[] = {
386     { "shift", "set phase shift", OFFSET(shift), AV_OPT_TYPE_DOUBLE, {.dbl=0}, -1.0, 1.0, FLAGS },
387     { NULL }
388 };
389
390 AVFILTER_DEFINE_CLASS(aphaseshift);
391
392 AVFilter ff_af_aphaseshift = {
393     .name            = "aphaseshift",
394     .description     = NULL_IF_CONFIG_SMALL("Apply phase shifting to input audio."),
395     .query_formats   = query_formats,
396     .priv_size       = sizeof(AFreqShift),
397     .priv_class      = &aphaseshift_class,
398     .uninit          = uninit,
399     .inputs          = inputs,
400     .outputs         = outputs,
401     .process_command = ff_filter_process_command,
402     .flags           = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
403                        AVFILTER_FLAG_SLICE_THREADS,
404 };