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 /**
281  * Perform IIR filtering on floating-point input samples.
282  *
283  * @param coeffs pointer to filter coefficients
284  * @param state pointer to filter state
285  * @param size input length
286  * @param src source samples
287  * @param sstep source stride
288  * @param dst filtered samples (destination may be the same as input)
289  * @param dstep destination stride
290  */
291 static void iir_filter_flt(const struct FFIIRFilterCoeffs *c,
292  struct FFIIRFilterState *s, int size,
293  const float *src, ptrdiff_t sstep,
294  float *dst, ptrdiff_t dstep)
295 {
296  if (c->order == 2) {
297  FILTER_O2(float, FLT)
298  } else if (c->order == 4) {
299  FILTER_BW_O4(float, FLT)
300  } else {
301  FILTER_DIRECT_FORM_II(float, FLT)
302  }
303 }
304 
306 {
307  av_freep(state);
308 }
309 
311 {
312  struct FFIIRFilterCoeffs *coeffs = *coeffsp;
313  if (coeffs) {
314  av_freep(&coeffs->cx);
315  av_freep(&coeffs->cy);
316  }
317  av_freep(coeffsp);
318 }
319 
321  f->filter_flt = iir_filter_flt;
322 
323 #if HAVE_MIPSFPU
325 #endif
326 }
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:50
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:209
av_cold
#define av_cold
Definition: attributes.h:90
s
#define s(width, name)
Definition: cbs_vp9.c:198
IIRFilterMode
IIRFilterMode
Definition: iirfilter.h:43
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:291
ff_iir_filter_init
void ff_iir_filter_init(FFIIRFilterContext *f)
Initialize FFIIRFilterContext.
Definition: iirfilter.c:320
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:44
FF_FILTER_MODE_HIGHPASS
@ FF_FILTER_MODE_HIGHPASS
Definition: iirfilter.h:45
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:122
dst
uint8_t ptrdiff_t const uint8_t ptrdiff_t int intptr_t intptr_t int int16_t * dst
Definition: dsp.h:83
size
int size
Definition: twinvq_data.h:10344
a0
static double a0(void *priv, double x, double y)
Definition: vf_xfade.c:2028
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
M_PI
#define M_PI
Definition: mathematics.h:67
IIRFilterType
IIRFilterType
Definition: iirfilter.h:35
log.h
lrintf
#define lrintf(x)
Definition: libm_mips.h:72
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
FF_FILTER_TYPE_BIQUAD
@ FF_FILTER_TYPE_BIQUAD
Definition: iirfilter.h:37
FFIIRFilterCoeffs::cy
float * cy
Definition: iirfilter.c:45
state
static struct @457 state
ff_iir_filter_free_coeffsp
av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp)
Free filter coefficients.
Definition: iirfilter.c:310
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:256
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
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:305
FF_FILTER_TYPE_BUTTERWORTH
@ FF_FILTER_TYPE_BUTTERWORTH
Definition: iirfilter.h:38
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:34
av_log
#define av_log(a,...)
Definition: tableprint_vlc.h:27
src
#define src
Definition: vp8dsp.c:248