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 "libavutil/tx.h"
22 #include "avfilter.h"
23 #include "internal.h"
24 #include "audio.h"
25 
26 #undef ctype
27 #undef ftype
28 #undef SQRT
29 #undef HYPOT
30 #undef SAMPLE_FORMAT
31 #undef TX_TYPE
32 #if DEPTH == 32
33 #define SAMPLE_FORMAT float
34 #define SQRT sqrtf
35 #define HYPOT hypotf
36 #define ctype AVComplexFloat
37 #define ftype float
38 #define TX_TYPE AV_TX_FLOAT_RDFT
39 #else
40 #define SAMPLE_FORMAT double
41 #define SQRT sqrt
42 #define HYPOT hypot
43 #define ctype AVComplexDouble
44 #define ftype double
45 #define TX_TYPE AV_TX_DOUBLE_RDFT
46 #endif
47 
48 #define fn3(a,b) a##_##b
49 #define fn2(a,b) fn3(a,b)
50 #define fn(a) fn2(a, SAMPLE_FORMAT)
51 
53 {
54  AudioFIRContext *s = ctx->priv;
55  ftype *mag, *phase, *delay, min = FLT_MAX, max = FLT_MIN;
56  ftype min_delay = FLT_MAX, max_delay = FLT_MIN;
57  int prev_ymag = -1, prev_yphase = -1, prev_ydelay = -1;
58  char text[32];
59  int channel, i, x;
60 
61  for (int y = 0; y < s->h; y++)
62  memset(out->data[0] + y * out->linesize[0], 0, s->w * 4);
63 
64  phase = av_malloc_array(s->w, sizeof(*phase));
65  mag = av_malloc_array(s->w, sizeof(*mag));
66  delay = av_malloc_array(s->w, sizeof(*delay));
67  if (!mag || !phase || !delay)
68  goto end;
69 
70  channel = av_clip(s->ir_channel, 0, s->ir[s->selir]->ch_layout.nb_channels - 1);
71  for (i = 0; i < s->w; i++) {
72  const ftype *src = (const ftype *)s->ir[s->selir]->extended_data[channel];
73  double w = i * M_PI / (s->w - 1);
74  double div, real_num = 0., imag_num = 0., real = 0., imag = 0.;
75 
76  for (x = 0; x < s->nb_taps[s->selir]; x++) {
77  real += cos(-x * w) * src[x];
78  imag += sin(-x * w) * src[x];
79  real_num += cos(-x * w) * src[x] * x;
80  imag_num += sin(-x * w) * src[x] * x;
81  }
82 
83  mag[i] = hypot(real, imag);
84  phase[i] = atan2(imag, real);
85  div = real * real + imag * imag;
86  delay[i] = (real_num * real + imag_num * imag) / div;
87  min = fminf(min, mag[i]);
88  max = fmaxf(max, mag[i]);
89  min_delay = fminf(min_delay, delay[i]);
90  max_delay = fmaxf(max_delay, delay[i]);
91  }
92 
93  for (i = 0; i < s->w; i++) {
94  int ymag = mag[i] / max * (s->h - 1);
95  int ydelay = (delay[i] - min_delay) / (max_delay - min_delay) * (s->h - 1);
96  int yphase = (0.5 * (1. + phase[i] / M_PI)) * (s->h - 1);
97 
98  ymag = s->h - 1 - av_clip(ymag, 0, s->h - 1);
99  yphase = s->h - 1 - av_clip(yphase, 0, s->h - 1);
100  ydelay = s->h - 1 - av_clip(ydelay, 0, s->h - 1);
101 
102  if (prev_ymag < 0)
103  prev_ymag = ymag;
104  if (prev_yphase < 0)
105  prev_yphase = yphase;
106  if (prev_ydelay < 0)
107  prev_ydelay = ydelay;
108 
109  draw_line(out, i, ymag, FFMAX(i - 1, 0), prev_ymag, 0xFFFF00FF);
110  draw_line(out, i, yphase, FFMAX(i - 1, 0), prev_yphase, 0xFF00FF00);
111  draw_line(out, i, ydelay, FFMAX(i - 1, 0), prev_ydelay, 0xFF00FFFF);
112 
113  prev_ymag = ymag;
114  prev_yphase = yphase;
115  prev_ydelay = ydelay;
116  }
117 
118  if (s->w > 400 && s->h > 100) {
119  drawtext(out, 2, 2, "Max Magnitude:", 0xDDDDDDDD);
120  snprintf(text, sizeof(text), "%.2f", max);
121  drawtext(out, 15 * 8 + 2, 2, text, 0xDDDDDDDD);
122 
123  drawtext(out, 2, 12, "Min Magnitude:", 0xDDDDDDDD);
124  snprintf(text, sizeof(text), "%.2f", min);
125  drawtext(out, 15 * 8 + 2, 12, text, 0xDDDDDDDD);
126 
127  drawtext(out, 2, 22, "Max Delay:", 0xDDDDDDDD);
128  snprintf(text, sizeof(text), "%.2f", max_delay);
129  drawtext(out, 11 * 8 + 2, 22, text, 0xDDDDDDDD);
130 
131  drawtext(out, 2, 32, "Min Delay:", 0xDDDDDDDD);
132  snprintf(text, sizeof(text), "%.2f", min_delay);
133  drawtext(out, 11 * 8 + 2, 32, text, 0xDDDDDDDD);
134  }
135 
136 end:
137  av_free(delay);
138  av_free(phase);
139  av_free(mag);
140 }
141 
143  int cur_nb_taps, int ch,
144  ftype *time)
145 {
146  ftype ch_gain = 1;
147 
148  switch (s->gtype) {
149  case -1:
150  ch_gain = 1;
151  break;
152  case 0:
153  {
154  ftype sum = 0;
155 
156  for (int i = 0; i < cur_nb_taps; i++)
157  sum += FFABS(time[i]);
158  ch_gain = 1. / sum;
159  }
160  break;
161  case 1:
162  {
163  ftype sum = 0;
164 
165  for (int i = 0; i < cur_nb_taps; i++)
166  sum += time[i];
167  ch_gain = 1. / sum;
168  }
169  break;
170  case 2:
171  {
172  ftype sum = 0;
173 
174  for (int i = 0; i < cur_nb_taps; i++)
175  sum += time[i] * time[i];
176  ch_gain = 1. / SQRT(sum);
177  }
178  break;
179  case 3:
180  case 4:
181  {
182  ftype *inc, *outc, scale, power;
183  AVTXContext *tx;
184  av_tx_fn tx_fn;
185  int ret, size;
186 
187  size = 1 << av_ceil_log2_c(cur_nb_taps);
188  inc = av_calloc(size + 2, sizeof(SAMPLE_FORMAT));
189  outc = av_calloc(size + 2, sizeof(SAMPLE_FORMAT));
190  if (!inc || !outc) {
191  av_free(outc);
192  av_free(inc);
193  break;
194  }
195 
196  scale = 1.;
197  ret = av_tx_init(&tx, &tx_fn, TX_TYPE, 0, size, &scale, 0);
198  if (ret < 0) {
199  av_free(outc);
200  av_free(inc);
201  break;
202  }
203 
204  {
205  memcpy(inc, time, cur_nb_taps * sizeof(SAMPLE_FORMAT));
206  tx_fn(tx, outc, inc, sizeof(SAMPLE_FORMAT));
207 
208  power = 0;
209  if (s->gtype == 3) {
210  for (int i = 0; i < size / 2 + 1; i++)
211  power = FFMAX(power, HYPOT(outc[i * 2], outc[i * 2 + 1]));
212  } else {
213  ftype sum = 0;
214  for (int i = 0; i < size / 2 + 1; i++)
215  sum += HYPOT(outc[i * 2], outc[i * 2 + 1]);
216  power = SQRT(sum / (size / 2 + 1));
217  }
218 
219  ch_gain = 1. / power;
220  }
221 
222  av_tx_uninit(&tx);
223  av_free(outc);
224  av_free(inc);
225  }
226  break;
227  default:
228  return AVERROR_BUG;
229  }
230 
231  if (ch_gain != 1. || s->ir_gain != 1.) {
232  ftype gain = ch_gain * s->ir_gain;
233 
234  av_log(ctx, AV_LOG_DEBUG, "ch%d gain %f\n", ch, gain);
235 #if DEPTH == 32
236  s->fdsp->vector_fmul_scalar(time, time, gain, FFALIGN(cur_nb_taps, 4));
237 #else
238  s->fdsp->vector_dmul_scalar(time, time, gain, FFALIGN(cur_nb_taps, 8));
239 #endif
240  }
241 
242  return 0;
243 }
244 
246  AudioFIRSegment *seg, int coeff_partition, int selir)
247 {
248  const int coffset = coeff_partition * seg->coeff_size;
249  const int nb_taps = s->nb_taps[selir];
250  ftype *time = (ftype *)s->norm_ir[selir]->extended_data[ch];
251  ftype *tempin = (ftype *)seg->tempin->extended_data[ch];
252  ftype *tempout = (ftype *)seg->tempout->extended_data[ch];
253  ctype *coeff = (ctype *)seg->coeff->extended_data[ch];
254  const int remaining = nb_taps - (seg->input_offset + coeff_partition * seg->part_size);
255  const int size = remaining >= seg->part_size ? seg->part_size : remaining;
256 
257  memset(tempin + size, 0, sizeof(*tempin) * (seg->block_size - size));
258  memcpy(tempin, time + seg->input_offset + coeff_partition * seg->part_size,
259  size * sizeof(*tempin));
260  seg->ctx_fn(seg->ctx[ch], tempout, tempin, sizeof(*tempin));
261  memcpy(coeff + coffset, tempout, seg->coeff_size * sizeof(*coeff));
262 
263  av_log(ctx, AV_LOG_DEBUG, "channel: %d\n", ch);
264  av_log(ctx, AV_LOG_DEBUG, "nb_partitions: %d\n", seg->nb_partitions);
265  av_log(ctx, AV_LOG_DEBUG, "partition size: %d\n", seg->part_size);
266  av_log(ctx, AV_LOG_DEBUG, "block size: %d\n", seg->block_size);
267  av_log(ctx, AV_LOG_DEBUG, "fft_length: %d\n", seg->fft_length);
268  av_log(ctx, AV_LOG_DEBUG, "coeff_size: %d\n", seg->coeff_size);
269  av_log(ctx, AV_LOG_DEBUG, "input_size: %d\n", seg->input_size);
270  av_log(ctx, AV_LOG_DEBUG, "input_offset: %d\n", seg->input_offset);
271 }
272 
273 static void fn(fir_fadd)(AudioFIRContext *s, ftype *dst, const ftype *src, int nb_samples)
274 {
275  if ((nb_samples & 15) == 0 && nb_samples >= 8) {
276 #if DEPTH == 32
277  s->fdsp->vector_fmac_scalar(dst, src, 1.f, nb_samples);
278 #else
279  s->fdsp->vector_dmac_scalar(dst, src, 1.0, nb_samples);
280 #endif
281  } else {
282  for (int n = 0; n < nb_samples; n++)
283  dst[n] += src[n];
284  }
285 }
286 
287 static int fn(fir_quantum)(AVFilterContext *ctx, AVFrame *out, int ch, int ioffset, int offset, int selir)
288 {
289  AudioFIRContext *s = ctx->priv;
290  const ftype *in = (const ftype *)s->in->extended_data[ch] + ioffset;
291  ftype *blockout, *ptr = (ftype *)out->extended_data[ch] + offset;
292  const int min_part_size = s->min_part_size;
293  const int nb_samples = FFMIN(min_part_size, out->nb_samples - offset);
294  const int nb_segments = s->nb_segments[selir];
295  const float dry_gain = s->dry_gain;
296  const float wet_gain = s->wet_gain;
297 
298  for (int segment = 0; segment < nb_segments; segment++) {
299  AudioFIRSegment *seg = &s->seg[selir][segment];
300  ftype *src = (ftype *)seg->input->extended_data[ch];
301  ftype *dst = (ftype *)seg->output->extended_data[ch];
302  ftype *sumin = (ftype *)seg->sumin->extended_data[ch];
303  ftype *sumout = (ftype *)seg->sumout->extended_data[ch];
304  ftype *tempin = (ftype *)seg->tempin->extended_data[ch];
305  ftype *buf = (ftype *)seg->buffer->extended_data[ch];
306  int *output_offset = &seg->output_offset[ch];
307  const int nb_partitions = seg->nb_partitions;
308  const int input_offset = seg->input_offset;
309  const int part_size = seg->part_size;
310  int j;
311 
312  seg->part_index[ch] = seg->part_index[ch] % nb_partitions;
313  if (dry_gain == 1.f) {
314  memcpy(src + input_offset, in, nb_samples * sizeof(*src));
315  } else if (min_part_size >= 8) {
316 #if DEPTH == 32
317  s->fdsp->vector_fmul_scalar(src + input_offset, in, dry_gain, FFALIGN(nb_samples, 4));
318 #else
319  s->fdsp->vector_dmul_scalar(src + input_offset, in, dry_gain, FFALIGN(nb_samples, 8));
320 #endif
321  } else {
322  ftype *src2 = src + input_offset;
323  for (int n = 0; n < nb_samples; n++)
324  src2[n] = in[n] * dry_gain;
325  }
326 
327  output_offset[0] += min_part_size;
328  if (output_offset[0] >= part_size) {
329  output_offset[0] = 0;
330  } else {
331  memmove(src, src + min_part_size, (seg->input_size - min_part_size) * sizeof(*src));
332 
333  dst += output_offset[0];
334  fn(fir_fadd)(s, ptr, dst, nb_samples);
335  continue;
336  }
337 
338  memset(sumin, 0, sizeof(*sumin) * seg->fft_length);
339 
340  blockout = (ftype *)seg->blockout->extended_data[ch] + seg->part_index[ch] * seg->block_size;
341  memset(tempin + part_size, 0, sizeof(*tempin) * (seg->block_size - part_size));
342  memcpy(tempin, src, sizeof(*src) * part_size);
343  seg->tx_fn(seg->tx[ch], blockout, tempin, sizeof(ftype));
344 
345  j = seg->part_index[ch];
346  for (int i = 0; i < nb_partitions; i++) {
347  const int input_partition = j;
348  const int coeff_partition = i;
349  const int coffset = coeff_partition * seg->coeff_size;
350  const ftype *blockout = (const ftype *)seg->blockout->extended_data[ch] + input_partition * seg->block_size;
351  const ctype *coeff = ((const ctype *)seg->coeff->extended_data[ch]) + coffset;
352 
353  if (j == 0)
354  j = nb_partitions;
355  j--;
356 
357 #if DEPTH == 32
358  s->afirdsp.fcmul_add(sumin, blockout, (const ftype *)coeff, part_size);
359 #else
360  s->afirdsp.dcmul_add(sumin, blockout, (const ftype *)coeff, part_size);
361 #endif
362  }
363 
364  seg->itx_fn(seg->itx[ch], sumout, sumin, sizeof(ctype));
365 
366  fn(fir_fadd)(s, buf, sumout, part_size);
367  memcpy(dst, buf, part_size * sizeof(*dst));
368  memcpy(buf, sumout + part_size, part_size * sizeof(*buf));
369 
370  fn(fir_fadd)(s, ptr, dst, nb_samples);
371 
372  if (part_size != min_part_size)
373  memmove(src, src + min_part_size, (seg->input_size - min_part_size) * sizeof(*src));
374 
375  seg->part_index[ch] = (seg->part_index[ch] + 1) % nb_partitions;
376  }
377 
378  if (wet_gain == 1.f)
379  return 0;
380 
381  if (min_part_size >= 8) {
382 #if DEPTH == 32
383  s->fdsp->vector_fmul_scalar(ptr, ptr, wet_gain, FFALIGN(nb_samples, 4));
384 #else
385  s->fdsp->vector_dmul_scalar(ptr, ptr, wet_gain, FFALIGN(nb_samples, 8));
386 #endif
387  } else {
388  for (int n = 0; n < nb_samples; n++)
389  ptr[n] *= wet_gain;
390  }
391 
392  return 0;
393 }
av_clip
#define av_clip
Definition: common.h:96
draw_response
static void fn() draw_response(AVFilterContext *ctx, AVFrame *out)
Definition: afir_template.c:52
convert_channel
static void fn() convert_channel(AVFilterContext *ctx, AudioFIRContext *s, int ch, AudioFIRSegment *seg, int coeff_partition, int selir)
Definition: afir_template.c:245
AudioFIRSegment::block_size
int block_size
Definition: af_afir.h:36
out
FILE * out
Definition: movenc.c:54
ctype
#define ctype
Definition: afir_template.c:43
TX_TYPE
#define TX_TYPE
Definition: afir_template.c:45
AVTXContext
Definition: tx_priv.h:235
SQRT
#define SQRT
Definition: afir_template.c:41
AVFrame
This structure describes decoded (raw) audio or video data.
Definition: frame.h:340
AudioFIRSegment::buffer
AVFrame * buffer
Definition: af_afir.h:50
w
uint8_t w
Definition: llviddspenc.c:38
HYPOT
#define HYPOT
Definition: afir_template.c:42
AudioFIRSegment::input_offset
int input_offset
Definition: af_afir.h:40
AudioFIRSegment::tx_fn
av_tx_fn tx_fn
Definition: af_afir.h:56
max
#define max(a, b)
Definition: cuda_runtime.h:33
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
AudioFIRSegment::part_size
int part_size
Definition: af_afir.h:35
AudioFIRSegment::input_size
int input_size
Definition: af_afir.h:39
av_tx_init
av_cold int av_tx_init(AVTXContext **ctx, av_tx_fn *tx, enum AVTXType type, int inv, int len, const void *scale, uint64_t flags)
Initialize a transform context with the given configuration (i)MDCTs with an odd length are currently...
Definition: tx.c:901
AudioFIRSegment::coeff
AVFrame * coeff
Definition: af_afir.h:51
fn
#define fn(a)
Definition: afir_template.c:50
scale
static av_always_inline float scale(float x, float s)
Definition: vf_v360.c:1389
AudioFIRSegment::blockout
AVFrame * blockout
Definition: af_afir.h:47
AudioFIRSegment
Definition: af_afir.h:33
AudioFIRSegment::tx
AVTXContext ** tx
Definition: af_afir.h:55
ftype
#define ftype
Definition: afir_template.c:44
av_tx_fn
void(* av_tx_fn)(AVTXContext *s, void *out, void *in, ptrdiff_t stride)
Function pointer to a function to perform the transform.
Definition: tx.h:151
s
#define s(width, name)
Definition: cbs_vp9.c:198
fminf
float fminf(float, float)
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:65
draw_line
static void draw_line(AVFrame *out, int x0, int y0, int x1, int y1, uint32_t color)
Definition: af_afir.c:73
AudioFIRSegment::itx_fn
av_tx_fn itx_fn
Definition: af_afir.h:56
AudioFIRSegment::output
AVFrame * output
Definition: af_afir.h:53
f
f
Definition: af_crystalizer.c:121
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:45
fir_fadd
static void fn() fir_fadd(AudioFIRContext *s, ftype *dst, const ftype *src, int nb_samples)
Definition: afir_template.c:273
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
AudioFIRSegment::tempin
AVFrame * tempin
Definition: af_afir.h:48
M_PI
#define M_PI
Definition: mathematics.h:67
av_tx_uninit
av_cold void av_tx_uninit(AVTXContext **ctx)
Frees a context and sets *ctx to NULL, does nothing when *ctx == NULL.
Definition: tx.c:294
internal.h
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:255
AudioFIRSegment::input
AVFrame * input
Definition: af_afir.h:52
get_power
static int fn() get_power(AVFilterContext *ctx, AudioFIRContext *s, int cur_nb_taps, int ch, ftype *time)
Definition: afir_template.c:142
AudioFIRSegment::coeff_size
int coeff_size
Definition: af_afir.h:38
AVFrame::extended_data
uint8_t ** extended_data
pointers to the data planes/channels.
Definition: frame.h:401
src2
const pixel * src2
Definition: h264pred_template.c:422
av_malloc_array
#define av_malloc_array(a, b)
Definition: tableprint_vlc.h:31
AudioFIRSegment::nb_partitions
int nb_partitions
Definition: af_afir.h:34
FFMIN
#define FFMIN(a, b)
Definition: macros.h:49
av_calloc
void * av_calloc(size_t nmemb, size_t size)
Definition: mem.c:262
ret
ret
Definition: filter_design.txt:187
AudioFIRSegment::itx
AVTXContext ** itx
Definition: af_afir.h:55
AudioFIRSegment::fft_length
int fft_length
Definition: af_afir.h:37
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:46
AudioFIRContext
Definition: af_afir.h:59
avfilter.h
segment
Definition: hls.c:76
AVFilterContext
An instance of a filter.
Definition: avfilter.h:397
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
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:79
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:42
channel
channel
Definition: ebur128.h:39
av_ceil_log2_c
static av_always_inline av_const int av_ceil_log2_c(int x)
Compute ceil(log2(x)).
Definition: common.h:418
drawtext
static void drawtext(AVFrame *pic, int x, int y, const char *txt, uint32_t color)
Definition: af_afir.c:50
tx.h
SAMPLE_FORMAT
#define SAMPLE_FORMAT
Definition: afir_template.c:40
fir_quantum
static int fn() fir_quantum(AVFilterContext *ctx, AVFrame *out, int ch, int ioffset, int offset, int selir)
Definition: afir_template.c:287
min
float min
Definition: vorbis_enc_data.h:429
AudioFIRSegment::part_index
int * part_index
Definition: af_afir.h:43