FFmpeg
afir_template.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2017 Paul B Mahol
3  *
4  * This file is part of FFmpeg.
5  *
6  * FFmpeg is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2.1 of the License, or (at your option) any later version.
10  *
11  * FFmpeg is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with FFmpeg; if not, write to the Free Software
18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19  */
20 
21 #include "avfilter.h"
22 #include "formats.h"
23 #include "internal.h"
24 #include "audio.h"
25 
26 #undef ctype
27 #undef ftype
28 #undef SQRT
29 #undef SAMPLE_FORMAT
30 #if DEPTH == 32
31 #define SAMPLE_FORMAT float
32 #define SQRT sqrtf
33 #define ctype AVComplexFloat
34 #define ftype float
35 #else
36 #define SAMPLE_FORMAT double
37 #define SQRT sqrt
38 #define ctype AVComplexDouble
39 #define ftype double
40 #endif
41 
42 #define fn3(a,b) a##_##b
43 #define fn2(a,b) fn3(a,b)
44 #define fn(a) fn2(a, SAMPLE_FORMAT)
45 
47 {
48  AudioFIRContext *s = ctx->priv;
49  ftype *mag, *phase, *delay, min = FLT_MAX, max = FLT_MIN;
50  ftype min_delay = FLT_MAX, max_delay = FLT_MIN;
51  int prev_ymag = -1, prev_yphase = -1, prev_ydelay = -1;
52  char text[32];
53  int channel, i, x;
54 
55  memset(out->data[0], 0, s->h * out->linesize[0]);
56 
57  phase = av_malloc_array(s->w, sizeof(*phase));
58  mag = av_malloc_array(s->w, sizeof(*mag));
59  delay = av_malloc_array(s->w, sizeof(*delay));
60  if (!mag || !phase || !delay)
61  goto end;
62 
63  channel = av_clip(s->ir_channel, 0, s->ir[s->selir]->ch_layout.nb_channels - 1);
64  for (i = 0; i < s->w; i++) {
65  const ftype *src = (const ftype *)s->ir[s->selir]->extended_data[channel];
66  double w = i * M_PI / (s->w - 1);
67  double div, real_num = 0., imag_num = 0., real = 0., imag = 0.;
68 
69  for (x = 0; x < s->nb_taps; x++) {
70  real += cos(-x * w) * src[x];
71  imag += sin(-x * w) * src[x];
72  real_num += cos(-x * w) * src[x] * x;
73  imag_num += sin(-x * w) * src[x] * x;
74  }
75 
76  mag[i] = hypot(real, imag);
77  phase[i] = atan2(imag, real);
78  div = real * real + imag * imag;
79  delay[i] = (real_num * real + imag_num * imag) / div;
80  min = fminf(min, mag[i]);
81  max = fmaxf(max, mag[i]);
82  min_delay = fminf(min_delay, delay[i]);
83  max_delay = fmaxf(max_delay, delay[i]);
84  }
85 
86  for (i = 0; i < s->w; i++) {
87  int ymag = mag[i] / max * (s->h - 1);
88  int ydelay = (delay[i] - min_delay) / (max_delay - min_delay) * (s->h - 1);
89  int yphase = (0.5 * (1. + phase[i] / M_PI)) * (s->h - 1);
90 
91  ymag = s->h - 1 - av_clip(ymag, 0, s->h - 1);
92  yphase = s->h - 1 - av_clip(yphase, 0, s->h - 1);
93  ydelay = s->h - 1 - av_clip(ydelay, 0, s->h - 1);
94 
95  if (prev_ymag < 0)
96  prev_ymag = ymag;
97  if (prev_yphase < 0)
98  prev_yphase = yphase;
99  if (prev_ydelay < 0)
100  prev_ydelay = ydelay;
101 
102  draw_line(out, i, ymag, FFMAX(i - 1, 0), prev_ymag, 0xFFFF00FF);
103  draw_line(out, i, yphase, FFMAX(i - 1, 0), prev_yphase, 0xFF00FF00);
104  draw_line(out, i, ydelay, FFMAX(i - 1, 0), prev_ydelay, 0xFF00FFFF);
105 
106  prev_ymag = ymag;
107  prev_yphase = yphase;
108  prev_ydelay = ydelay;
109  }
110 
111  if (s->w > 400 && s->h > 100) {
112  drawtext(out, 2, 2, "Max Magnitude:", 0xDDDDDDDD);
113  snprintf(text, sizeof(text), "%.2f", max);
114  drawtext(out, 15 * 8 + 2, 2, text, 0xDDDDDDDD);
115 
116  drawtext(out, 2, 12, "Min Magnitude:", 0xDDDDDDDD);
117  snprintf(text, sizeof(text), "%.2f", min);
118  drawtext(out, 15 * 8 + 2, 12, text, 0xDDDDDDDD);
119 
120  drawtext(out, 2, 22, "Max Delay:", 0xDDDDDDDD);
121  snprintf(text, sizeof(text), "%.2f", max_delay);
122  drawtext(out, 11 * 8 + 2, 22, text, 0xDDDDDDDD);
123 
124  drawtext(out, 2, 32, "Min Delay:", 0xDDDDDDDD);
125  snprintf(text, sizeof(text), "%.2f", min_delay);
126  drawtext(out, 11 * 8 + 2, 32, text, 0xDDDDDDDD);
127  }
128 
129 end:
130  av_free(delay);
131  av_free(phase);
132  av_free(mag);
133 }
134 
136 {
137  for (int ch = 0; ch < ctx->inputs[1 + s->selir]->ch_layout.nb_channels; ch++) {
138  ftype *time = (ftype *)s->ir[s->selir]->extended_data[!s->one2many * ch];
139  int toffset = 0;
140 
141  for (int i = FFMAX(1, s->length * s->nb_taps); i < s->nb_taps; i++)
142  time[i] = 0;
143 
144  av_log(ctx, AV_LOG_DEBUG, "channel: %d\n", ch);
145 
146  for (int segment = 0; segment < s->nb_segments; segment++) {
147  AudioFIRSegment *seg = &s->seg[segment];
148  ftype *blockin = (ftype *)seg->blockin->extended_data[ch];
149  ftype *blockout = (ftype *)seg->blockout->extended_data[ch];
150  ctype *coeff = (ctype *)seg->coeff->extended_data[ch];
151 
152  av_log(ctx, AV_LOG_DEBUG, "segment: %d\n", segment);
153 
154  for (int i = 0; i < seg->nb_partitions; i++) {
155  const int coffset = i * seg->coeff_size;
156  const int remaining = s->nb_taps - toffset;
157  const int size = remaining >= seg->part_size ? seg->part_size : remaining;
158 
159  if (size < 8) {
160  for (int n = 0; n < size; n++)
161  coeff[coffset + n].re = time[toffset + n];
162 
163  toffset += size;
164  continue;
165  }
166 
167  memset(blockin, 0, sizeof(*blockin) * seg->fft_length);
168  memcpy(blockin, time + toffset, size * sizeof(*blockin));
169 
170  seg->tx_fn(seg->tx[0], blockout, blockin, sizeof(ftype));
171 
172  for (int n = 0; n < seg->part_size + 1; n++) {
173  coeff[coffset + n].re = blockout[2 * n];
174  coeff[coffset + n].im = blockout[2 * n + 1];
175  }
176 
177  toffset += size;
178  }
179 
180  av_log(ctx, AV_LOG_DEBUG, "nb_partitions: %d\n", seg->nb_partitions);
181  av_log(ctx, AV_LOG_DEBUG, "partition size: %d\n", seg->part_size);
182  av_log(ctx, AV_LOG_DEBUG, "block size: %d\n", seg->block_size);
183  av_log(ctx, AV_LOG_DEBUG, "fft_length: %d\n", seg->fft_length);
184  av_log(ctx, AV_LOG_DEBUG, "coeff_size: %d\n", seg->coeff_size);
185  av_log(ctx, AV_LOG_DEBUG, "input_size: %d\n", seg->input_size);
186  av_log(ctx, AV_LOG_DEBUG, "input_offset: %d\n", seg->input_offset);
187  }
188  }
189 }
190 
191 static int fn(get_power)(AVFilterContext *ctx, AudioFIRContext *s, int cur_nb_taps)
192 {
193  ftype power = 0;
194  int ch;
195 
196  switch (s->gtype) {
197  case -1:
198  /* nothing to do */
199  break;
200  case 0:
201  for (ch = 0; ch < ctx->inputs[1 + s->selir]->ch_layout.nb_channels; ch++) {
202  ftype *time = (ftype *)s->ir[s->selir]->extended_data[!s->one2many * ch];
203 
204  for (int i = 0; i < cur_nb_taps; i++)
205  power += FFABS(time[i]);
206  }
207  s->gain = ctx->inputs[1 + s->selir]->ch_layout.nb_channels / power;
208  break;
209  case 1:
210  for (ch = 0; ch < ctx->inputs[1 + s->selir]->ch_layout.nb_channels; ch++) {
211  ftype *time = (ftype *)s->ir[s->selir]->extended_data[!s->one2many * ch];
212 
213  for (int i = 0; i < cur_nb_taps; i++)
214  power += time[i];
215  }
216  s->gain = ctx->inputs[1 + s->selir]->ch_layout.nb_channels / power;
217  break;
218  case 2:
219  for (ch = 0; ch < ctx->inputs[1 + s->selir]->ch_layout.nb_channels; ch++) {
220  ftype *time = (ftype *)s->ir[s->selir]->extended_data[!s->one2many * ch];
221 
222  for (int i = 0; i < cur_nb_taps; i++)
223  power += time[i] * time[i];
224  }
225  s->gain = SQRT(ch / power);
226  break;
227  default:
228  return AVERROR_BUG;
229  }
230 
231  s->gain = FFMIN(s->gain * s->ir_gain, 1.);
232 
233  av_log(ctx, AV_LOG_DEBUG, "power %f, gain %f\n", power, s->gain);
234 
235  for (int ch = 0; ch < ctx->inputs[1 + s->selir]->ch_layout.nb_channels; ch++) {
236  ftype *time = (ftype *)s->ir[s->selir]->extended_data[!s->one2many * ch];
237 
238 #if DEPTH == 32
239  s->fdsp->vector_fmul_scalar(time, time, s->gain, FFALIGN(cur_nb_taps, 4));
240 #else
241  s->fdsp->vector_dmul_scalar(time, time, s->gain, FFALIGN(cur_nb_taps, 8));
242 #endif
243  }
244 
245  return 0;
246 }
247 
248 static void fn(direct)(const ftype *in, const ctype *ir, int len, ftype *out)
249 {
250  for (int n = 0; n < len; n++)
251  for (int m = 0; m <= n; m++)
252  out[n] += ir[m].re * in[n - m];
253 }
254 
255 static void fn(fir_fadd)(AudioFIRContext *s, ftype *dst, const ftype *src, int nb_samples)
256 {
257  if ((nb_samples & 15) == 0 && nb_samples >= 16) {
258 #if DEPTH == 32
259  s->fdsp->vector_fmac_scalar(dst, src, 1.f, nb_samples);
260 #else
261  s->fdsp->vector_dmac_scalar(dst, src, 1.0, nb_samples);
262 #endif
263  } else {
264  for (int n = 0; n < nb_samples; n++)
265  dst[n] += src[n];
266  }
267 }
268 
269 static int fn(fir_quantum)(AVFilterContext *ctx, AVFrame *out, int ch, int offset)
270 {
271  AudioFIRContext *s = ctx->priv;
272  const ftype *in = (const ftype *)s->in->extended_data[ch] + offset;
273  ftype *blockin, *blockout, *buf, *ptr = (ftype *)out->extended_data[ch] + offset;
274  const int nb_samples = FFMIN(s->min_part_size, out->nb_samples - offset);
275  int n, i, j;
276 
277  for (int segment = 0; segment < s->nb_segments; segment++) {
278  AudioFIRSegment *seg = &s->seg[segment];
279  ftype *src = (ftype *)seg->input->extended_data[ch];
280  ftype *dst = (ftype *)seg->output->extended_data[ch];
281  ftype *sumin = (ftype *)seg->sumin->extended_data[ch];
282  ftype *sumout = (ftype *)seg->sumout->extended_data[ch];
283 
284  if (s->min_part_size >= 8) {
285 #if DEPTH == 32
286  s->fdsp->vector_fmul_scalar(src + seg->input_offset, in, s->dry_gain, FFALIGN(nb_samples, 4));
287 #else
288  s->fdsp->vector_dmul_scalar(src + seg->input_offset, in, s->dry_gain, FFALIGN(nb_samples, 8));
289 #endif
290  emms_c();
291  } else {
292  for (n = 0; n < nb_samples; n++)
293  src[seg->input_offset + n] = in[n] * s->dry_gain;
294  }
295 
296  seg->output_offset[ch] += s->min_part_size;
297  if (seg->output_offset[ch] == seg->part_size) {
298  seg->output_offset[ch] = 0;
299  } else {
300  memmove(src, src + s->min_part_size, (seg->input_size - s->min_part_size) * sizeof(*src));
301 
302  dst += seg->output_offset[ch];
303  fn(fir_fadd)(s, ptr, dst, nb_samples);
304  continue;
305  }
306 
307  if (seg->part_size < 8) {
308  memset(dst, 0, sizeof(*dst) * seg->part_size * seg->nb_partitions);
309 
310  j = seg->part_index[ch];
311 
312  for (i = 0; i < seg->nb_partitions; i++) {
313  const int coffset = j * seg->coeff_size;
314  const ctype *coeff = (const ctype *)seg->coeff->extended_data[ch * !s->one2many] + coffset;
315 
316  fn(direct)(src, coeff, nb_samples, dst);
317 
318  if (j == 0)
319  j = seg->nb_partitions;
320  j--;
321  }
322 
323  seg->part_index[ch] = (seg->part_index[ch] + 1) % seg->nb_partitions;
324 
325  memmove(src, src + s->min_part_size, (seg->input_size - s->min_part_size) * sizeof(*src));
326 
327  for (n = 0; n < nb_samples; n++) {
328  ptr[n] += dst[n];
329  }
330  continue;
331  }
332 
333  memset(sumin, 0, sizeof(*sumin) * seg->fft_length);
334  blockin = (ftype *)seg->blockin->extended_data[ch] + seg->part_index[ch] * seg->block_size;
335  blockout = (ftype *)seg->blockout->extended_data[ch] + seg->part_index[ch] * seg->block_size;
336  memset(blockin + seg->part_size, 0, sizeof(*blockin) * (seg->fft_length - seg->part_size));
337 
338  memcpy(blockin, src, sizeof(*src) * seg->part_size);
339 
340  seg->tx_fn(seg->tx[ch], blockout, blockin, sizeof(ftype));
341 
342  j = seg->part_index[ch];
343 
344  for (i = 0; i < seg->nb_partitions; i++) {
345  const int coffset = j * seg->coeff_size;
346  const ftype *blockout = (const ftype *)seg->blockout->extended_data[ch] + i * seg->block_size;
347  const ctype *coeff = (const ctype *)seg->coeff->extended_data[ch * !s->one2many] + coffset;
348 
349 #if DEPTH == 32
350  s->afirdsp.fcmul_add(sumin, blockout, (const ftype *)coeff, seg->part_size);
351 #else
352  s->afirdsp.dcmul_add(sumin, blockout, (const ftype *)coeff, seg->part_size);
353 #endif
354 
355  if (j == 0)
356  j = seg->nb_partitions;
357  j--;
358  }
359 
360  seg->itx_fn(seg->itx[ch], sumout, sumin, sizeof(ctype));
361 
362  buf = (ftype *)seg->buffer->extended_data[ch];
363  fn(fir_fadd)(s, buf, sumout, seg->part_size);
364 
365  memcpy(dst, buf, seg->part_size * sizeof(*dst));
366 
367  buf = (ftype *)seg->buffer->extended_data[ch];
368  memcpy(buf, sumout + seg->part_size, seg->part_size * sizeof(*buf));
369 
370  seg->part_index[ch] = (seg->part_index[ch] + 1) % seg->nb_partitions;
371 
372  memmove(src, src + s->min_part_size, (seg->input_size - s->min_part_size) * sizeof(*src));
373 
374  fn(fir_fadd)(s, ptr, dst, nb_samples);
375  }
376 
377  if (s->min_part_size >= 8) {
378 #if DEPTH == 32
379  s->fdsp->vector_fmul_scalar(ptr, ptr, s->wet_gain, FFALIGN(nb_samples, 4));
380 #else
381  s->fdsp->vector_dmul_scalar(ptr, ptr, s->wet_gain, FFALIGN(nb_samples, 8));
382 #endif
383  emms_c();
384  } else {
385  for (n = 0; n < nb_samples; n++)
386  ptr[n] *= s->wet_gain;
387  }
388 
389  return 0;
390 }
391 
392 
av_clip
#define av_clip
Definition: common.h:95
draw_response
static void fn() draw_response(AVFilterContext *ctx, AVFrame *out)
Definition: afir_template.c:46
AudioFIRSegment::block_size
int block_size
Definition: af_afir.h:34
out
FILE * out
Definition: movenc.c:54
ctype
#define ctype
Definition: afir_template.c:38
SQRT
#define SQRT
Definition: afir_template.c:37
AVFrame
This structure describes decoded (raw) audio or video data.
Definition: frame.h:325
AudioFIRSegment::buffer
AVFrame * buffer
Definition: af_afir.h:47
w
uint8_t w
Definition: llviddspenc.c:38
AudioFIRSegment::input_offset
int input_offset
Definition: af_afir.h:38
AudioFIRSegment::tx_fn
av_tx_fn tx_fn
Definition: af_afir.h:53
max
#define max(a, b)
Definition: cuda_runtime.h:33
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
AudioFIRSegment::blockin
AVFrame * blockin
Definition: af_afir.h:45
AudioFIRSegment::part_size
int part_size
Definition: af_afir.h:33
AudioFIRSegment::input_size
int input_size
Definition: af_afir.h:37
formats.h
AudioFIRSegment::coeff
AVFrame * coeff
Definition: af_afir.h:48
fn
#define fn(a)
Definition: afir_template.c:44
AudioFIRSegment::blockout
AVFrame * blockout
Definition: af_afir.h:46
AudioFIRSegment
Definition: af_afir.h:31
ftype
#define ftype
Definition: afir_template.c:39
AudioFIRSegment::tx
AVTXContext ** tx
Definition: af_afir.h:52
s
#define s(width, name)
Definition: cbs_vp9.c:256
fminf
float fminf(float, float)
direct
static void fn() direct(const ftype *in, const ctype *ir, int len, ftype *out)
Definition: afir_template.c:248
AV_LOG_DEBUG
#define AV_LOG_DEBUG
Stuff which is only useful for libav* developers.
Definition: log.h:201
ctx
AVFormatContext * ctx
Definition: movenc.c:48
FFABS
#define FFABS(a)
Absolute value, Note, INT_MIN / INT64_MIN result in undefined behavior as they are not representable ...
Definition: common.h:64
draw_line
static void draw_line(AVFrame *out, int x0, int y0, int x1, int y1, uint32_t color)
Definition: af_afir.c:71
AudioFIRSegment::itx_fn
av_tx_fn itx_fn
Definition: af_afir.h:53
AudioFIRSegment::output
AVFrame * output
Definition: af_afir.h:50
f
f
Definition: af_crystalizer.c:122
for
for(k=2;k<=8;++k)
Definition: h264pred_template.c:425
fmaxf
float fmaxf(float, float)
hypot
static av_const double hypot(double x, double y)
Definition: libm.h:366
size
int size
Definition: twinvq_data.h:10344
AudioFIRSegment::sumin
AVFrame * sumin
Definition: af_afir.h:43
fir_fadd
static void fn() fir_fadd(AudioFIRContext *s, ftype *dst, const ftype *src, int nb_samples)
Definition: afir_template.c:255
offset
it s the only field you need to keep assuming you have a context There is some magic you don t need to care about around this just let it vf offset
Definition: writing_filters.txt:86
M_PI
#define M_PI
Definition: mathematics.h:52
internal.h
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
AudioFIRSegment::input
AVFrame * input
Definition: af_afir.h:49
AudioFIRSegment::coeff_size
int coeff_size
Definition: af_afir.h:36
convert_channels
static void fn() convert_channels(AVFilterContext *ctx, AudioFIRContext *s)
Definition: afir_template.c:135
AVFrame::extended_data
uint8_t ** extended_data
pointers to the data planes/channels.
Definition: frame.h:386
av_malloc_array
#define av_malloc_array(a, b)
Definition: tableprint_vlc.h:31
AudioFIRSegment::nb_partitions
int nb_partitions
Definition: af_afir.h:32
FFMIN
#define FFMIN(a, b)
Definition: macros.h:49
DEPTH
#define DEPTH
Definition: v210enc.c:43
len
int len
Definition: vorbis_enc_data.h:426
AudioFIRSegment::itx
AVTXContext ** itx
Definition: af_afir.h:52
AudioFIRSegment::fft_length
int fft_length
Definition: af_afir.h:35
power
static float power(float r, float g, float b, float max)
Definition: preserve_color.h:45
AudioFIRSegment::sumout
AVFrame * sumout
Definition: af_afir.h:44
AudioFIRContext
Definition: af_afir.h:56
avfilter.h
get_power
static int fn() get_power(AVFilterContext *ctx, AudioFIRContext *s, int cur_nb_taps)
Definition: afir_template.c:191
segment
Definition: hls.c:75
AVFilterContext
An instance of a filter.
Definition: avfilter.h:415
audio.h
av_free
#define av_free(p)
Definition: tableprint_vlc.h:33
FFALIGN
#define FFALIGN(x, a)
Definition: macros.h:78
src
INIT_CLIP pixel * src
Definition: h264pred_template.c:418
fir_quantum
static int fn() fir_quantum(AVFilterContext *ctx, AVFrame *out, int ch, int offset)
Definition: afir_template.c:269
AVERROR_BUG
#define AVERROR_BUG
Internal bug, also see AVERROR_BUG2.
Definition: error.h:52
coeff
static const double coeff[2][5]
Definition: vf_owdenoise.c:78
av_log
#define av_log(a,...)
Definition: tableprint_vlc.h:27
snprintf
#define snprintf
Definition: snprintf.h:34
AudioFIRSegment::output_offset
int * output_offset
Definition: af_afir.h:40
channel
channel
Definition: ebur128.h:39
drawtext
static void drawtext(AVFrame *pic, int x, int y, const char *txt, uint32_t color)
Definition: af_afir.c:48
re
float re
Definition: fft.c:79
min
float min
Definition: vorbis_enc_data.h:429
AudioFIRSegment::part_index
int * part_index
Definition: af_afir.h:41