Mercurial > libavcodec.hg

comparison iirfilter.c @ 9945:12cf93d4b282 libavcodec

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Eliminate use of complex.h from iirfilter.c
author alexc
date Fri, 10 Jul 2009 20:45:13 +0000
parents 2cd0d1447bd3
children 0de8fdd06303
comparison
equal deleted inserted replaced
9944:c5ca5e520fe1 9945:12cf93d4b282
23 * @file libavcodec/iirfilter.c 23 * @file libavcodec/iirfilter.c
24 * different IIR filters implementation 24 * different IIR filters implementation
25 */ 25 */
26 26
27 #include "iirfilter.h" 27 #include "iirfilter.h"
28 #include <complex.h>
29 #include <math.h> 28 #include <math.h>
30 29
31 /** 30 /**
32 * IIR filter global parameters 31 * IIR filter global parameters
33 */ 32 */
54 float stopband, float ripple) 53 float stopband, float ripple)
55 { 54 {
56 int i, j, size; 55 int i, j, size;
57 FFIIRFilterCoeffs *c; 56 FFIIRFilterCoeffs *c;
58 double wa; 57 double wa;
59 double complex p[MAXORDER + 1]; 58 double p[MAXORDER + 1][2];
60 59
61 if(filt_type != FF_FILTER_TYPE_BUTTERWORTH || filt_mode != FF_FILTER_MODE_LOWPASS) 60 if(filt_type != FF_FILTER_TYPE_BUTTERWORTH || filt_mode != FF_FILTER_MODE_LOWPASS)
62 return NULL; 61 return NULL;
63 if(order <= 1 || (order & 1) || order > MAXORDER || cutoff_ratio >= 1.0) 62 if(order <= 1 || (order & 1) || order > MAXORDER || cutoff_ratio >= 1.0)
64 return NULL; 63 return NULL;
72 71
73 c->cx[0] = 1; 72 c->cx[0] = 1;
74 for(i = 1; i < (order >> 1) + 1; i++) 73 for(i = 1; i < (order >> 1) + 1; i++)
75 c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i; 74 c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i;
76 75
77 p[0] = 1.0; 76 p[0][0] = 1.0;
77 p[0][1] = 0.0;
78 for(i = 1; i <= order; i++) 78 for(i = 1; i <= order; i++)
79 p[i] = 0.0; 79 p[i][0] = p[i][1] = 0.0;
80 for(i = 0; i < order; i++){ 80 for(i = 0; i < order; i++){
81 double complex zp; 81 double zp[2];
82 double th = (i + (order >> 1) + 0.5) * M_PI / order; 82 double th = (i + (order >> 1) + 0.5) * M_PI / order;
83 zp = cexp(I*th) * wa; 83 double a_re, a_im, c_re, c_im;
84 zp = (zp + 2.0) / (zp - 2.0); 84 zp[0] = cos(th) * wa;
85 zp[1] = sin(th) * wa;
86 a_re = zp[0] + 2.0;
87 c_re = zp[0] - 2.0;
88 a_im =
89 c_im = zp[1];
90 zp[0] = (a_re * c_re + a_im * c_im) / (c_re * c_re + c_im * c_im);
91 zp[1] = (a_im * c_re - a_re * c_im) / (c_re * c_re + c_im * c_im);
85 92
86 for(j = order; j >= 1; j--) 93 for(j = order; j >= 1; j--)
87 p[j] = zp*p[j] + p[j - 1]; 94 {
88 p[0] *= zp; 95 a_re = p[j][0];
96 a_im = p[j][1];
97 p[j][0] = a_re*zp[0] - a_im*zp[1] + p[j-1][0];
98 p[j][1] = a_re*zp[1] + a_im*zp[0] + p[j-1][1];
99 }
100 a_re = p[0][0]*zp[0] - p[0][1]*zp[1];
101 p[0][1] = p[0][0]*zp[1] + p[0][1]*zp[0];
102 p[0][0] = a_re;
89 } 103 }
90 c->gain = creal(p[order]); 104 c->gain = p[order][0];
91 for(i = 0; i < order; i++){ 105 for(i = 0; i < order; i++){
92 c->gain += creal(p[i]); 106 c->gain += p[i][0];
93 c->cy[i] = creal(-p[i] / p[order]); 107 c->cy[i] = (-p[i][0] * p[order][0] + -p[i][1] * p[order][1]) /
108 (p[order][0] * p[order][0] + p[order][1] * p[order][1]);
94 } 109 }
95 c->gain /= 1 << order; 110 c->gain /= 1 << order;
96 111
97 return c; 112 return c;
98 } 113 }