Commit | Line | Data |
---|---|---|
2ba45a60 DM |
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 "iirfilter.h" | |
28 | #include <math.h> | |
29 | #include "libavutil/attributes.h" | |
30 | #include "libavutil/common.h" | |
31 | ||
32 | /** | |
33 | * IIR filter global parameters | |
34 | */ | |
35 | typedef struct FFIIRFilterCoeffs{ | |
36 | int order; | |
37 | float gain; | |
38 | int *cx; | |
39 | float *cy; | |
40 | }FFIIRFilterCoeffs; | |
41 | ||
42 | /** | |
43 | * IIR filter state | |
44 | */ | |
45 | typedef struct FFIIRFilterState{ | |
46 | float x[1]; | |
47 | }FFIIRFilterState; | |
48 | ||
49 | /// maximum supported filter order | |
50 | #define MAXORDER 30 | |
51 | ||
52 | static av_cold int butterworth_init_coeffs(void *avc, | |
53 | struct FFIIRFilterCoeffs *c, | |
54 | enum IIRFilterMode filt_mode, | |
55 | int order, float cutoff_ratio, | |
56 | float stopband) | |
57 | { | |
58 | int i, j; | |
59 | double wa; | |
60 | double p[MAXORDER + 1][2]; | |
61 | ||
62 | if (filt_mode != FF_FILTER_MODE_LOWPASS) { | |
63 | av_log(avc, AV_LOG_ERROR, "Butterworth filter currently only supports " | |
64 | "low-pass filter mode\n"); | |
65 | return -1; | |
66 | } | |
67 | if (order & 1) { | |
68 | av_log(avc, AV_LOG_ERROR, "Butterworth filter currently only supports " | |
69 | "even filter orders\n"); | |
70 | return -1; | |
71 | } | |
72 | ||
73 | wa = 2 * tan(M_PI * 0.5 * cutoff_ratio); | |
74 | ||
75 | c->cx[0] = 1; | |
76 | for(i = 1; i < (order >> 1) + 1; i++) | |
77 | c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i; | |
78 | ||
79 | p[0][0] = 1.0; | |
80 | p[0][1] = 0.0; | |
81 | for(i = 1; i <= order; i++) | |
82 | p[i][0] = p[i][1] = 0.0; | |
83 | for(i = 0; i < order; i++){ | |
84 | double zp[2]; | |
85 | double th = (i + (order >> 1) + 0.5) * M_PI / order; | |
86 | double a_re, a_im, c_re, c_im; | |
87 | zp[0] = cos(th) * wa; | |
88 | zp[1] = sin(th) * wa; | |
89 | a_re = zp[0] + 2.0; | |
90 | c_re = zp[0] - 2.0; | |
91 | a_im = | |
92 | c_im = zp[1]; | |
93 | zp[0] = (a_re * c_re + a_im * c_im) / (c_re * c_re + c_im * c_im); | |
94 | zp[1] = (a_im * c_re - a_re * c_im) / (c_re * c_re + c_im * c_im); | |
95 | ||
96 | for(j = order; j >= 1; j--) | |
97 | { | |
98 | a_re = p[j][0]; | |
99 | a_im = p[j][1]; | |
100 | p[j][0] = a_re*zp[0] - a_im*zp[1] + p[j-1][0]; | |
101 | p[j][1] = a_re*zp[1] + a_im*zp[0] + p[j-1][1]; | |
102 | } | |
103 | a_re = p[0][0]*zp[0] - p[0][1]*zp[1]; | |
104 | p[0][1] = p[0][0]*zp[1] + p[0][1]*zp[0]; | |
105 | p[0][0] = a_re; | |
106 | } | |
107 | c->gain = p[order][0]; | |
108 | for(i = 0; i < order; i++){ | |
109 | c->gain += p[i][0]; | |
110 | c->cy[i] = (-p[i][0] * p[order][0] + -p[i][1] * p[order][1]) / | |
111 | (p[order][0] * p[order][0] + p[order][1] * p[order][1]); | |
112 | } | |
113 | c->gain /= 1 << order; | |
114 | ||
115 | return 0; | |
116 | } | |
117 | ||
118 | static av_cold int biquad_init_coeffs(void *avc, struct FFIIRFilterCoeffs *c, | |
119 | enum IIRFilterMode filt_mode, int order, | |
120 | float cutoff_ratio, float stopband) | |
121 | { | |
122 | double cos_w0, sin_w0; | |
123 | double a0, x0, x1; | |
124 | ||
125 | if (filt_mode != FF_FILTER_MODE_HIGHPASS && | |
126 | filt_mode != FF_FILTER_MODE_LOWPASS) { | |
127 | av_log(avc, AV_LOG_ERROR, "Biquad filter currently only supports " | |
128 | "high-pass and low-pass filter modes\n"); | |
129 | return -1; | |
130 | } | |
131 | if (order != 2) { | |
132 | av_log(avc, AV_LOG_ERROR, "Biquad filter must have order of 2\n"); | |
133 | return -1; | |
134 | } | |
135 | ||
136 | cos_w0 = cos(M_PI * cutoff_ratio); | |
137 | sin_w0 = sin(M_PI * cutoff_ratio); | |
138 | ||
139 | a0 = 1.0 + (sin_w0 / 2.0); | |
140 | ||
141 | if (filt_mode == FF_FILTER_MODE_HIGHPASS) { | |
142 | c->gain = ((1.0 + cos_w0) / 2.0) / a0; | |
143 | x0 = ((1.0 + cos_w0) / 2.0) / a0; | |
144 | x1 = (-(1.0 + cos_w0)) / a0; | |
145 | } else { // FF_FILTER_MODE_LOWPASS | |
146 | c->gain = ((1.0 - cos_w0) / 2.0) / a0; | |
147 | x0 = ((1.0 - cos_w0) / 2.0) / a0; | |
148 | x1 = (1.0 - cos_w0) / a0; | |
149 | } | |
150 | c->cy[0] = (-1.0 + (sin_w0 / 2.0)) / a0; | |
151 | c->cy[1] = (2.0 * cos_w0) / a0; | |
152 | ||
153 | // divide by gain to make the x coeffs integers. | |
154 | // during filtering, the delay state will include the gain multiplication | |
155 | c->cx[0] = lrintf(x0 / c->gain); | |
156 | c->cx[1] = lrintf(x1 / c->gain); | |
157 | ||
158 | return 0; | |
159 | } | |
160 | ||
161 | av_cold struct FFIIRFilterCoeffs* ff_iir_filter_init_coeffs(void *avc, | |
162 | enum IIRFilterType filt_type, | |
163 | enum IIRFilterMode filt_mode, | |
164 | int order, float cutoff_ratio, | |
165 | float stopband, float ripple) | |
166 | { | |
167 | FFIIRFilterCoeffs *c; | |
168 | int ret = 0; | |
169 | ||
170 | if (order <= 0 || order > MAXORDER || cutoff_ratio >= 1.0) | |
171 | return NULL; | |
172 | ||
173 | FF_ALLOCZ_OR_GOTO(avc, c, sizeof(FFIIRFilterCoeffs), | |
174 | init_fail); | |
175 | FF_ALLOC_OR_GOTO (avc, c->cx, sizeof(c->cx[0]) * ((order >> 1) + 1), | |
176 | init_fail); | |
177 | FF_ALLOC_OR_GOTO (avc, c->cy, sizeof(c->cy[0]) * order, | |
178 | init_fail); | |
179 | c->order = order; | |
180 | ||
181 | switch (filt_type) { | |
182 | case FF_FILTER_TYPE_BUTTERWORTH: | |
183 | ret = butterworth_init_coeffs(avc, c, filt_mode, order, cutoff_ratio, | |
184 | stopband); | |
185 | break; | |
186 | case FF_FILTER_TYPE_BIQUAD: | |
187 | ret = biquad_init_coeffs(avc, c, filt_mode, order, cutoff_ratio, | |
188 | stopband); | |
189 | break; | |
190 | default: | |
191 | av_log(avc, AV_LOG_ERROR, "filter type is not currently implemented\n"); | |
192 | goto init_fail; | |
193 | } | |
194 | ||
195 | if (!ret) | |
196 | return c; | |
197 | ||
198 | init_fail: | |
f6fa7814 | 199 | ff_iir_filter_free_coeffsp(&c); |
2ba45a60 DM |
200 | return NULL; |
201 | } | |
202 | ||
203 | av_cold struct FFIIRFilterState* ff_iir_filter_init_state(int order) | |
204 | { | |
205 | FFIIRFilterState* s = av_mallocz(sizeof(FFIIRFilterState) + sizeof(s->x[0]) * (order - 1)); | |
206 | return s; | |
207 | } | |
208 | ||
209 | #define CONV_S16(dest, source) dest = av_clip_int16(lrintf(source)); | |
210 | ||
211 | #define CONV_FLT(dest, source) dest = source; | |
212 | ||
213 | #define FILTER_BW_O4_1(i0, i1, i2, i3, fmt) \ | |
214 | in = *src0 * c->gain \ | |
215 | + c->cy[0]*s->x[i0] + c->cy[1]*s->x[i1] \ | |
216 | + c->cy[2]*s->x[i2] + 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, int sstep, int16_t *dst, int dstep) | |
279 | { | |
280 | if (c->order == 2) { | |
281 | FILTER_O2(int16_t, S16) | |
282 | } else if (c->order == 4) { | |
283 | FILTER_BW_O4(int16_t, S16) | |
284 | } else { | |
285 | FILTER_DIRECT_FORM_II(int16_t, S16) | |
286 | } | |
287 | } | |
288 | ||
289 | void ff_iir_filter_flt(const struct FFIIRFilterCoeffs *c, | |
290 | struct FFIIRFilterState *s, int size, | |
291 | const float *src, int sstep, float *dst, int dstep) | |
292 | { | |
293 | if (c->order == 2) { | |
294 | FILTER_O2(float, FLT) | |
295 | } else if (c->order == 4) { | |
296 | FILTER_BW_O4(float, FLT) | |
297 | } else { | |
298 | FILTER_DIRECT_FORM_II(float, FLT) | |
299 | } | |
300 | } | |
301 | ||
f6fa7814 | 302 | av_cold void ff_iir_filter_free_statep(struct FFIIRFilterState **state) |
2ba45a60 | 303 | { |
f6fa7814 | 304 | av_freep(state); |
2ba45a60 DM |
305 | } |
306 | ||
f6fa7814 | 307 | av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp) |
2ba45a60 | 308 | { |
f6fa7814 | 309 | struct FFIIRFilterCoeffs *coeffs = *coeffsp; |
2ba45a60 | 310 | if(coeffs){ |
f6fa7814 DM |
311 | av_freep(&coeffs->cx); |
312 | av_freep(&coeffs->cy); | |
2ba45a60 | 313 | } |
f6fa7814 | 314 | av_freep(coeffsp); |
2ba45a60 DM |
315 | } |
316 | ||
317 | void ff_iir_filter_init(FFIIRFilterContext *f) { | |
318 | f->filter_flt = ff_iir_filter_flt; | |
319 | ||
320 | if (HAVE_MIPSFPU) | |
321 | ff_iir_filter_init_mips(f); | |
322 | } | |
323 | ||
324 | #ifdef TEST | |
325 | #include <stdio.h> | |
326 | ||
327 | #define FILT_ORDER 4 | |
328 | #define SIZE 1024 | |
329 | int main(void) | |
330 | { | |
331 | struct FFIIRFilterCoeffs *fcoeffs = NULL; | |
332 | struct FFIIRFilterState *fstate = NULL; | |
333 | float cutoff_coeff = 0.4; | |
334 | int16_t x[SIZE], y[SIZE]; | |
335 | int i; | |
336 | ||
337 | fcoeffs = ff_iir_filter_init_coeffs(NULL, FF_FILTER_TYPE_BUTTERWORTH, | |
338 | FF_FILTER_MODE_LOWPASS, FILT_ORDER, | |
339 | cutoff_coeff, 0.0, 0.0); | |
340 | fstate = ff_iir_filter_init_state(FILT_ORDER); | |
341 | ||
342 | for (i = 0; i < SIZE; i++) { | |
343 | x[i] = lrint(0.75 * INT16_MAX * sin(0.5*M_PI*i*i/SIZE)); | |
344 | } | |
345 | ||
346 | ff_iir_filter(fcoeffs, fstate, SIZE, x, 1, y, 1); | |
347 | ||
348 | for (i = 0; i < SIZE; i++) | |
349 | printf("%6d %6d\n", x[i], y[i]); | |
350 | ||
f6fa7814 DM |
351 | ff_iir_filter_free_coeffsp(&fcoeffs); |
352 | ff_iir_filter_free_statep(&fstate); | |
2ba45a60 DM |
353 | return 0; |
354 | } | |
355 | #endif /* TEST */ |