FFmpeg
aap_template.c
Go to the documentation of this file.
1 /*
2  * This file is part of FFmpeg.
3  *
4  * FFmpeg is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU Lesser General Public
6  * License as published by the Free Software Foundation; either
7  * version 2.1 of the License, or (at your option) any later version.
8  *
9  * FFmpeg is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12  * Lesser General Public License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with FFmpeg; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 #undef ZERO
20 #undef ONE
21 #undef ftype
22 #undef SAMPLE_FORMAT
23 #if DEPTH == 32
24 #define SAMPLE_FORMAT float
25 #define ftype float
26 #define ONE 1.f
27 #define ZERO 0.f
28 #else
29 #define SAMPLE_FORMAT double
30 #define ftype double
31 #define ONE 1.0
32 #define ZERO 0.0
33 #endif
34 
35 #define fn3(a,b) a##_##b
36 #define fn2(a,b) fn3(a,b)
37 #define fn(a) fn2(a, SAMPLE_FORMAT)
38 
40  ftype *coeffs, ftype *tmp, int *offset)
41 {
42  const int order = s->order;
43  ftype output;
44 
45  delay[*offset] = sample;
46 
47  memcpy(tmp, coeffs + order - *offset, order * sizeof(ftype));
48 #if DEPTH == 32
49  output = s->fdsp->scalarproduct_float(delay, tmp, s->kernel_size);
50 #else
51  output = s->fdsp->scalarproduct_double(delay, tmp, s->kernel_size);
52 #endif
53 
54  if (--(*offset) < 0)
55  *offset = order - 1;
56 
57  return output;
58 }
59 
60 static int fn(lup_decompose)(ftype **MA, const int N, const ftype tol, int *P)
61 {
62  for (int i = 0; i <= N; i++)
63  P[i] = i;
64 
65  for (int i = 0; i < N; i++) {
66  ftype maxA = ZERO;
67  int imax = i;
68 
69  for (int k = i; k < N; k++) {
70  ftype absA = fabs(MA[k][i]);
71  if (absA > maxA) {
72  maxA = absA;
73  imax = k;
74  }
75  }
76 
77  if (maxA < tol)
78  return 0;
79 
80  if (imax != i) {
81  FFSWAP(int, P[i], P[imax]);
82  FFSWAP(ftype *, MA[i], MA[imax]);
83  P[N]++;
84  }
85 
86  for (int j = i + 1; j < N; j++) {
87  MA[j][i] /= MA[i][i];
88 
89  for (int k = i + 1; k < N; k++)
90  MA[j][k] -= MA[j][i] * MA[i][k];
91  }
92  }
93 
94  return 1;
95 }
96 
97 static void fn(lup_invert)(ftype *const *MA, const int *P, const int N, ftype **IA)
98 {
99  for (int j = 0; j < N; j++) {
100  for (int i = 0; i < N; i++) {
101  IA[i][j] = P[i] == j ? ONE : ZERO;
102 
103  for (int k = 0; k < i; k++)
104  IA[i][j] -= MA[i][k] * IA[k][j];
105  }
106 
107  for (int i = N - 1; i >= 0; i--) {
108  for (int k = i + 1; k < N; k++)
109  IA[i][j] -= MA[i][k] * IA[k][j];
110 
111  IA[i][j] /= MA[i][i];
112  }
113  }
114 }
115 
116 static ftype fn(process_sample)(AudioAPContext *s, ftype input, ftype desired, int ch)
117 {
118  ftype *dcoeffs = (ftype *)s->dcoeffs->extended_data[ch];
119  ftype *coeffs = (ftype *)s->coeffs->extended_data[ch];
120  ftype *delay = (ftype *)s->delay->extended_data[ch];
121  ftype **itmpmp = (ftype **)&s->itmpmp[s->projection * ch];
122  ftype **tmpmp = (ftype **)&s->tmpmp[s->projection * ch];
123  ftype *tmpm = (ftype *)s->tmpm->extended_data[ch];
124  ftype *tmp = (ftype *)s->tmp->extended_data[ch];
125  ftype *e = (ftype *)s->e->extended_data[ch];
126  ftype *x = (ftype *)s->x->extended_data[ch];
127  ftype *w = (ftype *)s->w->extended_data[ch];
128  int *p = (int *)s->p->extended_data[ch];
129  int *offset = (int *)s->offset->extended_data[ch];
130  const int projection = s->projection;
131  const ftype delta = s->delta;
132  const int order = s->order;
133  const int length = projection + order;
134  const ftype mu = s->mu;
135  const ftype tol = 0.00001f;
136  ftype output;
137 
138  x[offset[2] + length] = x[offset[2]] = input;
139  delay[offset[0] + order] = input;
140 
141  output = fn(fir_sample)(s, input, delay, coeffs, tmp, offset);
142  e[offset[1]] = e[offset[1] + projection] = desired - output;
143 
144  for (int i = 0; i < projection; i++) {
145  const int iprojection = i * projection;
146 
147  for (int j = i; j < projection; j++) {
148  ftype sum = ZERO;
149  for (int k = 0; k < order; k++)
150  sum += x[offset[2] + i + k] * x[offset[2] + j + k];
151  tmpm[iprojection + j] = sum;
152  if (i != j)
153  tmpm[j * projection + i] = sum;
154  }
155 
156  tmpm[iprojection + i] += delta;
157  }
158 
159  fn(lup_decompose)(tmpmp, projection, tol, p);
160  fn(lup_invert)(tmpmp, p, projection, itmpmp);
161 
162  for (int i = 0; i < projection; i++) {
163  ftype sum = ZERO;
164  for (int j = 0; j < projection; j++)
165  sum += itmpmp[i][j] * e[j + offset[1]];
166  w[i] = sum;
167  }
168 
169  for (int i = 0; i < order; i++) {
170  ftype sum = ZERO;
171  for (int j = 0; j < projection; j++)
172  sum += x[offset[2] + i + j] * w[j];
173  dcoeffs[i] = sum;
174  }
175 
176  for (int i = 0; i < order; i++)
177  coeffs[i] = coeffs[i + order] = coeffs[i] + mu * dcoeffs[i];
178 
179  if (--offset[1] < 0)
180  offset[1] = projection - 1;
181 
182  if (--offset[2] < 0)
183  offset[2] = length - 1;
184 
185  switch (s->output_mode) {
186  case IN_MODE: output = input; break;
187  case DESIRED_MODE: output = desired; break;
188  case OUT_MODE: output = desired - output; break;
189  case NOISE_MODE: output = input - output; break;
190  case ERROR_MODE: break;
191  }
192  return output;
193 }
194 
195 static int fn(filter_channels)(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
196 {
197  AudioAPContext *s = ctx->priv;
198  AVFrame *out = arg;
199  const int start = (out->ch_layout.nb_channels * jobnr) / nb_jobs;
200  const int end = (out->ch_layout.nb_channels * (jobnr+1)) / nb_jobs;
201 
202  for (int c = start; c < end; c++) {
203  const ftype *input = (const ftype *)s->frame[0]->extended_data[c];
204  const ftype *desired = (const ftype *)s->frame[1]->extended_data[c];
205  ftype *output = (ftype *)out->extended_data[c];
206 
207  for (int n = 0; n < out->nb_samples; n++) {
208  output[n] = fn(process_sample)(s, input[n], desired[n], c);
209  if (ctx->is_disabled)
210  output[n] = input[n];
211  }
212  }
213 
214  return 0;
215 }
out
FILE * out
Definition: movenc.c:55
output
filter_frame For filters that do not use the this method is called when a frame is pushed to the filter s input It can be called at any time except in a reentrant way If the input frame is enough to produce output
Definition: filter_design.txt:225
AVFrame
This structure describes decoded (raw) audio or video data.
Definition: frame.h:374
tmp
static uint8_t tmp[11]
Definition: aes_ctr.c:28
w
uint8_t w
Definition: llviddspenc.c:38
IA
#define IA(x)
Definition: cast5.c:30
ftype
#define ftype
Definition: aap_template.c:30
fir_sample
static ftype fn() fir_sample(AudioAPContext *s, ftype sample, ftype *delay, ftype *coeffs, ftype *tmp, int *offset)
Definition: aap_template.c:39
ERROR_MODE
@ ERROR_MODE
Definition: af_aap.c:38
s
#define s(width, name)
Definition: cbs_vp9.c:198
process_sample
static ftype fn() process_sample(AudioAPContext *s, ftype input, ftype desired, int ch)
Definition: aap_template.c:116
ctx
AVFormatContext * ctx
Definition: movenc.c:49
DESIRED_MODE
@ DESIRED_MODE
Definition: af_aap.c:35
arg
const char * arg
Definition: jacosubdec.c:67
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
lup_decompose
static int fn() lup_decompose(ftype **MA, const int N, const ftype tol, int *P)
Definition: aap_template.c:60
c
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
lup_invert
static void fn() lup_invert(ftype *const *MA, const int *P, const int N, ftype **IA)
Definition: aap_template.c:97
filter_channels
static int fn() filter_channels(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
Definition: aap_template.c:195
P
#define P
for
for(k=2;k<=8;++k)
Definition: h264pred_template.c:425
sample
#define sample
Definition: flacdsp_template.c:44
fn
#define fn(a)
Definition: aap_template.c:37
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
N
#define N
Definition: af_mcompand.c:54
input
and forward the test the status of outputs and forward it to the corresponding return FFERROR_NOT_READY If the filters stores internally one or a few frame for some input
Definition: filter_design.txt:172
ZERO
#define ZERO
Definition: aap_template.c:32
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
delta
float delta
Definition: vorbis_enc_data.h:430
AudioAPContext
Definition: af_aap.c:42
NOISE_MODE
@ NOISE_MODE
Definition: af_aap.c:37
FFSWAP
#define FFSWAP(type, a, b)
Definition: macros.h:52
IN_MODE
@ IN_MODE
Definition: af_aap.c:34
OUT_MODE
@ OUT_MODE
Definition: af_aap.c:36
AVFilterContext
An instance of a filter.
Definition: avfilter.h:407
ONE
#define ONE
Definition: aap_template.c:31