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 
291  struct FFIIRFilterState *s, int size,
292  const float *src, ptrdiff_t sstep,
293  float *dst, ptrdiff_t dstep)
294 {
295  if (c->order == 2) {
296  FILTER_O2(float, FLT)
297  } else if (c->order == 4) {
298  FILTER_BW_O4(float, FLT)
299  } else {
300  FILTER_DIRECT_FORM_II(float, FLT)
301  }
302 }
303 
305 {
306  av_freep(state);
307 }
308 
310 {
311  struct FFIIRFilterCoeffs *coeffs = *coeffsp;
312  if (coeffs) {
313  av_freep(&coeffs->cx);
314  av_freep(&coeffs->cy);
315  }
316  av_freep(coeffsp);
317 }
318 
321 
322  if (HAVE_MIPSFPU)
324 }
#define NULL
Definition: coverity.c:32
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
#define a0
Definition: regdef.h:46
IIRFilterMode
Definition: iirfilter.h:44
void * av_mallocz(size_t size)
Allocate a memory block with alignment suitable for all memory accesses (including vectors if availab...
Definition: mem.c:237
av_cold struct FFIIRFilterState * ff_iir_filter_init_state(int order)
Create new filter state.
Definition: iirfilter.c:201
Macro definitions for various function/variable attributes.
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
#define av_cold
Definition: attributes.h:88
#define av_malloc(s)
void(* filter_flt)(const struct FFIIRFilterCoeffs *coeffs, struct FFIIRFilterState *state, int size, const float *src, ptrdiff_t sstep, float *dst, ptrdiff_t dstep)
Perform IIR filtering on floating-point input samples.
Definition: iirfilter.h:63
#define f(width, name)
Definition: cbs_vp9.c:255
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
#define lrintf(x)
Definition: libm_mips.h:70
ptrdiff_t size
Definition: opengl_enc.c:100
#define av_log(a,...)
#define FILTER_DIRECT_FORM_II(type, fmt)
Definition: iirfilter.c:238
#define src
Definition: vp8dsp.c:254
#define AV_LOG_ERROR
Something went wrong and cannot losslessly be recovered.
Definition: log.h:194
av_cold void ff_iir_filter_free_statep(struct FFIIRFilterState **state)
Free and zero filter state.
Definition: iirfilter.c:304
IIR filter state.
Definition: iirfilter.c:47
#define th
Definition: regdef.h:75
#define FILTER_O2(type, fmt)
Definition: iirfilter.c:260
#define MAXORDER
maximum supported filter order
Definition: iirfilter.c:52
#define FILTER_BW_O4(type, fmt)
Definition: iirfilter.c:225
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
static struct @318 state
#define s(width, name)
Definition: cbs_vp9.c:257
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
float x[1]
Definition: iirfilter.c:48
void ff_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:290
IIRFilterType
Definition: iirfilter.h:36
IIR filter global parameters.
Definition: iirfilter.c:37
void ff_iir_filter_init_mips(FFIIRFilterContext *f)
void ff_iir_filter_init(FFIIRFilterContext *f)
Initialize FFIIRFilterContext.
Definition: iirfilter.c:319
common internal and external API header
IIR filter interface.
#define av_freep(p)
#define M_PI
Definition: mathematics.h:52
av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp)
Free filter coefficients.
Definition: iirfilter.c:309
int i
Definition: input.c:407