FFmpeg
iirfilter.c
Go to the documentation of this file.
1 /*
2  * IIR filter
3  * Copyright (c) 2008 Konstantin Shishkov
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 /**
23  * @file
24  * different IIR filters implementation
25  */
26 
27 #include <math.h>
28 
29 #include "config.h"
30 
31 #include "libavutil/attributes.h"
32 #include "libavutil/common.h"
33 #include "libavutil/log.h"
34 #include "libavutil/mem.h"
35 
36 #include "iirfilter.h"
37 
38 /**
39  * IIR filter global parameters
40  */
41 typedef struct FFIIRFilterCoeffs {
42  int order;
43  float gain;
44  int *cx;
45  float *cy;
47 
48 /**
49  * IIR filter state
50  */
51 typedef struct FFIIRFilterState {
52  float x[1];
54 
55 /// maximum supported filter order
56 #define MAXORDER 30
57 
58 static av_cold int butterworth_init_coeffs(void *avc,
59  struct FFIIRFilterCoeffs *c,
60  enum IIRFilterMode filt_mode,
61  int order, float cutoff_ratio,
62  float stopband)
63 {
64  int i, j;
65  double wa;
66  double p[MAXORDER + 1][2];
67 
68  if (filt_mode != FF_FILTER_MODE_LOWPASS) {
69  av_log(avc, AV_LOG_ERROR, "Butterworth filter currently only supports "
70  "low-pass filter mode\n");
71  return -1;
72  }
73  if (order & 1) {
74  av_log(avc, AV_LOG_ERROR, "Butterworth filter currently only supports "
75  "even filter orders\n");
76  return -1;
77  }
78 
79  wa = 2 * tan(M_PI * 0.5 * cutoff_ratio);
80 
81  c->cx[0] = 1;
82  for (i = 1; i < (order >> 1) + 1; i++)
83  c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i;
84 
85  p[0][0] = 1.0;
86  p[0][1] = 0.0;
87  for (i = 1; i <= order; i++)
88  p[i][0] = p[i][1] = 0.0;
89  for (i = 0; i < order; i++) {
90  double zp[2];
91  double th = (i + (order >> 1) + 0.5) * M_PI / order;
92  double a_re, a_im, c_re, c_im;
93  zp[0] = cos(th) * wa;
94  zp[1] = sin(th) * wa;
95  a_re = zp[0] + 2.0;
96  c_re = zp[0] - 2.0;
97  a_im =
98  c_im = zp[1];
99  zp[0] = (a_re * c_re + a_im * c_im) / (c_re * c_re + c_im * c_im);
100  zp[1] = (a_im * c_re - a_re * c_im) / (c_re * c_re + c_im * c_im);
101 
102  for (j = order; j >= 1; j--) {
103  a_re = p[j][0];
104  a_im = p[j][1];
105  p[j][0] = a_re * zp[0] - a_im * zp[1] + p[j - 1][0];
106  p[j][1] = a_re * zp[1] + a_im * zp[0] + p[j - 1][1];
107  }
108  a_re = p[0][0] * zp[0] - p[0][1] * zp[1];
109  p[0][1] = p[0][0] * zp[1] + p[0][1] * zp[0];
110  p[0][0] = a_re;
111  }
112  c->gain = p[order][0];
113  for (i = 0; i < order; i++) {
114  c->gain += p[i][0];
115  c->cy[i] = (-p[i][0] * p[order][0] + -p[i][1] * p[order][1]) /
116  (p[order][0] * p[order][0] + p[order][1] * p[order][1]);
117  }
118  c->gain /= 1 << order;
119 
120  return 0;
121 }
122 
123 static av_cold int biquad_init_coeffs(void *avc, struct FFIIRFilterCoeffs *c,
124  enum IIRFilterMode filt_mode, int order,
125  float cutoff_ratio, float stopband)
126 {
127  double cos_w0, sin_w0;
128  double a0, x0, x1;
129 
130  if (filt_mode != FF_FILTER_MODE_HIGHPASS &&
131  filt_mode != FF_FILTER_MODE_LOWPASS) {
132  av_log(avc, AV_LOG_ERROR, "Biquad filter currently only supports "
133  "high-pass and low-pass filter modes\n");
134  return -1;
135  }
136  if (order != 2) {
137  av_log(avc, AV_LOG_ERROR, "Biquad filter must have order of 2\n");
138  return -1;
139  }
140 
141  cos_w0 = cos(M_PI * cutoff_ratio);
142  sin_w0 = sin(M_PI * cutoff_ratio);
143 
144  a0 = 1.0 + (sin_w0 / 2.0);
145 
146  if (filt_mode == FF_FILTER_MODE_HIGHPASS) {
147  c->gain = ((1.0 + cos_w0) / 2.0) / a0;
148  x0 = ((1.0 + cos_w0) / 2.0) / a0;
149  x1 = (-(1.0 + cos_w0)) / a0;
150  } else { // FF_FILTER_MODE_LOWPASS
151  c->gain = ((1.0 - cos_w0) / 2.0) / a0;
152  x0 = ((1.0 - cos_w0) / 2.0) / a0;
153  x1 = (1.0 - cos_w0) / a0;
154  }
155  c->cy[0] = (-1.0 + (sin_w0 / 2.0)) / a0;
156  c->cy[1] = (2.0 * cos_w0) / a0;
157 
158  // divide by gain to make the x coeffs integers.
159  // during filtering, the delay state will include the gain multiplication
160  c->cx[0] = lrintf(x0 / c->gain);
161  c->cx[1] = lrintf(x1 / c->gain);
162 
163  return 0;
164 }
165 
167  enum IIRFilterType filt_type,
168  enum IIRFilterMode filt_mode,
169  int order, float cutoff_ratio,
170  float stopband, float ripple)
171 {
173  int ret = 0;
174 
175  if (order <= 0 || order > MAXORDER || cutoff_ratio >= 1.0)
176  return NULL;
177 
178  if (!(c = av_mallocz(sizeof(*c))) ||
179  !(c->cx = av_malloc (sizeof(c->cx[0]) * ((order >> 1) + 1))) ||
180  !(c->cy = av_malloc (sizeof(c->cy[0]) * order)))
181  goto free;
182  c->order = order;
183 
184  switch (filt_type) {
186  ret = butterworth_init_coeffs(avc, c, filt_mode, order, cutoff_ratio,
187  stopband);
188  break;
190  ret = biquad_init_coeffs(avc, c, filt_mode, order, cutoff_ratio,
191  stopband);
192  break;
193  default:
194  av_log(avc, AV_LOG_ERROR, "filter type is not currently implemented\n");
195  goto free;
196  }
197 
198  if (!ret)
199  return c;
200 free:
202  return NULL;
203 }
204 
206 {
207  FFIIRFilterState *s = av_mallocz(sizeof(FFIIRFilterState) + sizeof(s->x[0]) * (order - 1));
208  return s;
209 }
210 
211 #define CONV_S16(dest, source) dest = av_clip_int16(lrintf(source));
212 
213 #define CONV_FLT(dest, source) dest = source;
214 
215 #define FILTER_BW_O4_1(i0, i1, i2, i3, fmt) \
216  in = *src0 * c->gain + \
217  c->cy[0] * s->x[i0] + \
218  c->cy[1] * s->x[i1] + \
219  c->cy[2] * s->x[i2] + \
220  c->cy[3] * s->x[i3]; \
221  res = (s->x[i0] + in) * 1 + \
222  (s->x[i1] + s->x[i3]) * 4 + \
223  s->x[i2] * 6; \
224  CONV_ ## fmt(*dst0, res) \
225  s->x[i0] = in; \
226  src0 += sstep; \
227  dst0 += dstep;
228 
229 #define FILTER_BW_O4(type, fmt) { \
230  int i; \
231  const type *src0 = src; \
232  type *dst0 = dst; \
233  for (i = 0; i < size; i += 4) { \
234  float in, res; \
235  FILTER_BW_O4_1(0, 1, 2, 3, fmt); \
236  FILTER_BW_O4_1(1, 2, 3, 0, fmt); \
237  FILTER_BW_O4_1(2, 3, 0, 1, fmt); \
238  FILTER_BW_O4_1(3, 0, 1, 2, fmt); \
239  } \
240 }
241 
242 #define FILTER_DIRECT_FORM_II(type, fmt) { \
243  int i; \
244  const type *src0 = src; \
245  type *dst0 = dst; \
246  for (i = 0; i < size; i++) { \
247  int j; \
248  float in, res; \
249  in = *src0 * c->gain; \
250  for (j = 0; j < c->order; j++) \
251  in += c->cy[j] * s->x[j]; \
252  res = s->x[0] + in + s->x[c->order >> 1] * c->cx[c->order >> 1]; \
253  for (j = 1; j < c->order >> 1; j++) \
254  res += (s->x[j] + s->x[c->order - j]) * c->cx[j]; \
255  for (j = 0; j < c->order - 1; j++) \
256  s->x[j] = s->x[j + 1]; \
257  CONV_ ## fmt(*dst0, res) \
258  s->x[c->order - 1] = in; \
259  src0 += sstep; \
260  dst0 += dstep; \
261  } \
262 }
263 
264 #define FILTER_O2(type, fmt) { \
265  int i; \
266  const type *src0 = src; \
267  type *dst0 = dst; \
268  for (i = 0; i < size; i++) { \
269  float in = *src0 * c->gain + \
270  s->x[0] * c->cy[0] + \
271  s->x[1] * c->cy[1]; \
272  CONV_ ## fmt(*dst0, s->x[0] + in + s->x[1] * c->cx[1]) \
273  s->x[0] = s->x[1]; \
274  s->x[1] = in; \
275  src0 += sstep; \
276  dst0 += dstep; \
277  } \
278 }
279 
280 void ff_iir_filter(const struct FFIIRFilterCoeffs *c,
281  struct FFIIRFilterState *s, int size,
282  const int16_t *src, ptrdiff_t sstep,
283  int16_t *dst, ptrdiff_t dstep)
284 {
285  if (c->order == 2) {
286  FILTER_O2(int16_t, S16)
287  } else if (c->order == 4) {
288  FILTER_BW_O4(int16_t, S16)
289  } else {
290  FILTER_DIRECT_FORM_II(int16_t, S16)
291  }
292 }
293 
294 /**
295  * Perform IIR filtering on floating-point input samples.
296  *
297  * @param coeffs pointer to filter coefficients
298  * @param state pointer to filter state
299  * @param size input length
300  * @param src source samples
301  * @param sstep source stride
302  * @param dst filtered samples (destination may be the same as input)
303  * @param dstep destination stride
304  */
305 static void iir_filter_flt(const struct FFIIRFilterCoeffs *c,
306  struct FFIIRFilterState *s, int size,
307  const float *src, ptrdiff_t sstep,
308  float *dst, ptrdiff_t dstep)
309 {
310  if (c->order == 2) {
311  FILTER_O2(float, FLT)
312  } else if (c->order == 4) {
313  FILTER_BW_O4(float, FLT)
314  } else {
315  FILTER_DIRECT_FORM_II(float, FLT)
316  }
317 }
318 
320 {
321  av_freep(state);
322 }
323 
325 {
326  struct FFIIRFilterCoeffs *coeffs = *coeffsp;
327  if (coeffs) {
328  av_freep(&coeffs->cx);
329  av_freep(&coeffs->cy);
330  }
331  av_freep(coeffsp);
332 }
333 
335  f->filter_flt = iir_filter_flt;
336 
337 #if HAVE_MIPSFPU
339 #endif
340 }
FFIIRFilterState
IIR filter state.
Definition: iirfilter.c:51
iirfilter.h
FFIIRFilterCoeffs::gain
float gain
Definition: iirfilter.c:43
MAXORDER
#define MAXORDER
maximum supported filter order
Definition: iirfilter.c:56
FILTER_BW_O4
#define FILTER_BW_O4(type, fmt)
Definition: iirfilter.c:229
FFIIRFilterContext
Definition: iirfilter.h:51
ff_iir_filter
void ff_iir_filter(const struct FFIIRFilterCoeffs *c, struct FFIIRFilterState *s, int size, const int16_t *src, ptrdiff_t sstep, int16_t *dst, ptrdiff_t dstep)
Perform IIR filtering on signed 16-bit input samples.
Definition: iirfilter.c:280
FFIIRFilterCoeffs::cx
int * cx
Definition: iirfilter.c:44
av_malloc
#define av_malloc(s)
Definition: tableprint_vlc.h:30
FFIIRFilterCoeffs
IIR filter global parameters.
Definition: iirfilter.c:41
AV_LOG_ERROR
#define AV_LOG_ERROR
Something went wrong and cannot losslessly be recovered.
Definition: log.h:180
av_cold
#define av_cold
Definition: attributes.h:90
s
#define s(width, name)
Definition: cbs_vp9.c:198
IIRFilterMode
IIRFilterMode
Definition: iirfilter.h:44
iir_filter_flt
static void iir_filter_flt(const struct FFIIRFilterCoeffs *c, struct FFIIRFilterState *s, int size, const float *src, ptrdiff_t sstep, float *dst, ptrdiff_t dstep)
Perform IIR filtering on floating-point input samples.
Definition: iirfilter.c:305
ff_iir_filter_init
void ff_iir_filter_init(FFIIRFilterContext *f)
Initialize FFIIRFilterContext.
Definition: iirfilter.c:334
NULL
#define NULL
Definition: coverity.c:32
ff_iir_filter_init_coeffs
av_cold struct FFIIRFilterCoeffs * ff_iir_filter_init_coeffs(void *avc, enum IIRFilterType filt_type, enum IIRFilterMode filt_mode, int order, float cutoff_ratio, float stopband, float ripple)
Initialize filter coefficients.
Definition: iirfilter.c:166
FF_FILTER_MODE_LOWPASS
@ FF_FILTER_MODE_LOWPASS
Definition: iirfilter.h:45
FF_FILTER_MODE_HIGHPASS
@ FF_FILTER_MODE_HIGHPASS
Definition: iirfilter.h:46
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
FILTER_O2
#define FILTER_O2(type, fmt)
Definition: iirfilter.c:264
f
f
Definition: af_crystalizer.c:121
size
int size
Definition: twinvq_data.h:10344
biquad_init_coeffs
static av_cold int biquad_init_coeffs(void *avc, struct FFIIRFilterCoeffs *c, enum IIRFilterMode filt_mode, int order, float cutoff_ratio, float stopband)
Definition: iirfilter.c:123
attributes.h
a0
#define a0
Definition: regdef.h:46
M_PI
#define M_PI
Definition: mathematics.h:67
th
#define th
Definition: regdef.h:75
IIRFilterType
IIRFilterType
Definition: iirfilter.h:36
log.h
lrintf
#define lrintf(x)
Definition: libm_mips.h:72
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:255
FF_FILTER_TYPE_BIQUAD
@ FF_FILTER_TYPE_BIQUAD
Definition: iirfilter.h:38
FFIIRFilterCoeffs::cy
float * cy
Definition: iirfilter.c:45
ff_iir_filter_free_coeffsp
av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp)
Free filter coefficients.
Definition: iirfilter.c:324
common.h
av_mallocz
void * av_mallocz(size_t size)
Allocate a memory block with alignment suitable for all memory accesses (including vectors if availab...
Definition: mem.c:254
ff_iir_filter_init_mips
void ff_iir_filter_init_mips(FFIIRFilterContext *f)
Definition: iirfilter_mips.c:203
FFIIRFilterState::x
float x[1]
Definition: iirfilter.c:52
FILTER_DIRECT_FORM_II
#define FILTER_DIRECT_FORM_II(type, fmt)
Definition: iirfilter.c:242
ret
ret
Definition: filter_design.txt:187
ff_iir_filter_init_state
av_cold struct FFIIRFilterState * ff_iir_filter_init_state(int order)
Create new filter state.
Definition: iirfilter.c:205
butterworth_init_coeffs
static av_cold int butterworth_init_coeffs(void *avc, struct FFIIRFilterCoeffs *c, enum IIRFilterMode filt_mode, int order, float cutoff_ratio, float stopband)
Definition: iirfilter.c:58
state
static struct @362 state
FFIIRFilterCoeffs::order
int order
Definition: iirfilter.c:42
mem.h
ff_iir_filter_free_statep
av_cold void ff_iir_filter_free_statep(struct FFIIRFilterState **state)
Free and zero filter state.
Definition: iirfilter.c:319
FF_FILTER_TYPE_BUTTERWORTH
@ FF_FILTER_TYPE_BUTTERWORTH
Definition: iirfilter.h:39
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:34
src
INIT_CLIP pixel * src
Definition: h264pred_template.c:418
av_log
#define av_log(a,...)
Definition: tableprint_vlc.h:27