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