Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
6147 | serge | 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 |
||
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: |
||
199 | ff_iir_filter_free_coeffsp(&c); |
||
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 | |||
302 | av_cold void ff_iir_filter_free_statep(struct FFIIRFilterState **state) |
||
303 | { |
||
304 | av_freep(state); |
||
305 | } |
||
306 | |||
307 | av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp) |
||
308 | { |
||
309 | struct FFIIRFilterCoeffs *coeffs = *coeffsp; |
||
310 | if(coeffs){ |
||
311 | av_freep(&coeffs->cx); |
||
312 | av_freep(&coeffs->cy); |
||
313 | } |
||
314 | av_freep(coeffsp); |
||
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 |
||
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 | |||
351 | ff_iir_filter_free_coeffsp(&fcoeffs); |
||
352 | ff_iir_filter_free_statep(&fstate); |
||
353 | return 0; |
||
354 | } |
||
355 | #endif /* TEST */>>>>>>>>=>><>>>=>> |