FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavcodec/iirfilter.c
Date: 2023-09-24 13:02:57
Exec Total Coverage
Lines: 70 116 60.3%
Functions: 7 9 77.8%
Branches: 26 73 35.6%
Line Branch Exec Source
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;
46 } FFIIRFilterCoeffs;
47
48 /**
49 * IIR filter state
50 */
51 typedef struct FFIIRFilterState {
52 float x[1];
53 } FFIIRFilterState;
54
55 /// maximum supported filter order
56 #define MAXORDER 30
57
58 1 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 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 1 wa = 2 * tan(M_PI * 0.5 * cutoff_ratio);
80
81 1 c->cx[0] = 1;
82
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (i = 1; i < (order >> 1) + 1; i++)
83 2 c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i;
84
85 1 p[0][0] = 1.0;
86 1 p[0][1] = 0.0;
87
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (i = 1; i <= order; i++)
88 4 p[i][0] = p[i][1] = 0.0;
89
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (i = 0; i < order; i++) {
90 double zp[2];
91 4 double th = (i + (order >> 1) + 0.5) * M_PI / order;
92 double a_re, a_im, c_re, c_im;
93 4 zp[0] = cos(th) * wa;
94 4 zp[1] = sin(th) * wa;
95 4 a_re = zp[0] + 2.0;
96 4 c_re = zp[0] - 2.0;
97 4 a_im =
98 4 c_im = zp[1];
99 4 zp[0] = (a_re * c_re + a_im * c_im) / (c_re * c_re + c_im * c_im);
100 4 zp[1] = (a_im * c_re - a_re * c_im) / (c_re * c_re + c_im * c_im);
101
102
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 20 for (j = order; j >= 1; j--) {
103 16 a_re = p[j][0];
104 16 a_im = p[j][1];
105 16 p[j][0] = a_re * zp[0] - a_im * zp[1] + p[j - 1][0];
106 16 p[j][1] = a_re * zp[1] + a_im * zp[0] + p[j - 1][1];
107 }
108 4 a_re = p[0][0] * zp[0] - p[0][1] * zp[1];
109 4 p[0][1] = p[0][0] * zp[1] + p[0][1] * zp[0];
110 4 p[0][0] = a_re;
111 }
112 1 c->gain = p[order][0];
113
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (i = 0; i < order; i++) {
114 4 c->gain += p[i][0];
115 4 c->cy[i] = (-p[i][0] * p[order][0] + -p[i][1] * p[order][1]) /
116 4 (p[order][0] * p[order][0] + p[order][1] * p[order][1]);
117 }
118 1 c->gain /= 1 << order;
119
120 1 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
166 1 av_cold struct FFIIRFilterCoeffs *ff_iir_filter_init_coeffs(void *avc,
167 enum IIRFilterType filt_type,
168 enum IIRFilterMode filt_mode,
169 int order, float cutoff_ratio,
170 float stopband, float ripple)
171 {
172 FFIIRFilterCoeffs *c;
173 1 int ret = 0;
174
175
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 1 if (order <= 0 || order > MAXORDER || cutoff_ratio >= 1.0)
176 ✗ return NULL;
177
178
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (!(c = av_mallocz(sizeof(*c))) ||
179
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 !(c->cx = av_malloc (sizeof(c->cx[0]) * ((order >> 1) + 1))) ||
180
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 !(c->cy = av_malloc (sizeof(c->cy[0]) * order)))
181 ✗ goto free;
182 1 c->order = order;
183
184
1/3
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1 switch (filt_type) {
185 1 case FF_FILTER_TYPE_BUTTERWORTH:
186 1 ret = butterworth_init_coeffs(avc, c, filt_mode, order, cutoff_ratio,
187 stopband);
188 1 break;
189 ✗ case FF_FILTER_TYPE_BIQUAD:
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
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (!ret)
199 1 return c;
200 ✗ free:
201 ✗ ff_iir_filter_free_coeffsp(&c);
202 ✗ return NULL;
203 }
204
205 1 av_cold struct FFIIRFilterState *ff_iir_filter_init_state(int order)
206 {
207 1 FFIIRFilterState *s = av_mallocz(sizeof(FFIIRFilterState) + sizeof(s->x[0]) * (order - 1));
208 1 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 1 void ff_iir_filter(const struct FFIIRFilterCoeffs *c,
281 struct FFIIRFilterState *s, int size,
282 const int16_t *src, ptrdiff_t sstep,
283 int16_t *dst, ptrdiff_t dstep)
284 {
285
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (c->order == 2) {
286 ✗ FILTER_O2(int16_t, S16)
287
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (c->order == 4) {
288
2/2
✓ Branch 0 taken 256 times.
✓ Branch 1 taken 1 times.
 257 FILTER_BW_O4(int16_t, S16)
289 } else {
290 ✗ FILTER_DIRECT_FORM_II(int16_t, S16)
291 }
292 1 }
293
294 /**
295 * Perform IIR filtering on floating-point input samples.
296 *
297 * @param coeffs pointer to filter coefficients
298 * @param state pointer to filter state
299 * @param size input length
300 * @param src source samples
301 * @param sstep source stride
302 * @param dst filtered samples (destination may be the same as input)
303 * @param dstep destination stride
304 */
305 ✗ static void iir_filter_flt(const struct FFIIRFilterCoeffs *c,
306 struct FFIIRFilterState *s, int size,
307 const float *src, ptrdiff_t sstep,
308 float *dst, ptrdiff_t dstep)
309 {
310 ✗ if (c->order == 2) {
311 ✗ FILTER_O2(float, FLT)
312 ✗ } else if (c->order == 4) {
313 ✗ FILTER_BW_O4(float, FLT)
314 } else {
315 ✗ FILTER_DIRECT_FORM_II(float, FLT)
316 }
317 ✗ }
318
319 1 av_cold void ff_iir_filter_free_statep(struct FFIIRFilterState **state)
320 {
321 1 av_freep(state);
322 1 }
323
324 12 av_cold void ff_iir_filter_free_coeffsp(struct FFIIRFilterCoeffs **coeffsp)
325 {
326 12 struct FFIIRFilterCoeffs *coeffs = *coeffsp;
327
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 11 times.
 12 if (coeffs) {
328 1 av_freep(&coeffs->cx);
329 1 av_freep(&coeffs->cy);
330 }
331 12 av_freep(coeffsp);
332 12 }
333
334 11 void ff_iir_filter_init(FFIIRFilterContext *f) {
335 11 f->filter_flt = iir_filter_flt;
336
337 #if HAVE_MIPSFPU
338 ff_iir_filter_init_mips(f);
339 #endif
340 11 }
341