2 * Copyright © 2009 Keith Packard <keithp@keithp.com>
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
9 * This program is distributed in the hope that it will be useful, but
10 * WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * General Public License for more details.
14 * You should have received a copy of the GNU General Public License along
15 * with this program; if not, write to the Free Software Foundation, Inc.,
16 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
24 static inline double sqr (double x) { return x * x; }
27 * Kaiser Window digital filter
31 /* not used in this program */
32 static double highpass (double n, double m, double wc)
39 return (M_PI/2 - wc) / M_PI;
40 return -sin(dist * (M_PI/2-wc)) / (M_PI * dist);
44 static double lowpass (double n, double m, double wc)
51 return sin (wc * dist) / (M_PI * dist);
54 static double kaiser (double n, double m, double beta)
57 return i0 (beta * sqrt (1 - sqr((n - alpha) / alpha))) / i0(beta);
60 static double beta (double A)
63 return 0.1102 * (A - 8.7);
65 return 0.5842 * pow((A - 21), 0.4) + 0.07886 * (A - 21);
70 static int M (double A, double delta_omega)
73 return ceil ((A - 7.95) / (2.285 * delta_omega));
75 return ceil(5.79 / delta_omega);
86 static struct filter_param
87 filter (double omega_pass, double omega_stop, double error)
89 struct filter_param p;
90 p.omega_pass = omega_pass;
91 p.delta_omega = omega_stop - omega_pass;
92 p.A = -20 * log10 (error);
94 p.M = M (p.A, p.delta_omega);
101 make_low_pass_filter(double omega_pass, double omega_stop, double error, int *length_p)
103 struct filter_param p = filter(omega_pass, omega_stop, error);
109 lpf = calloc (length, sizeof(double));
110 for (n = 0; n < length; n++)
111 lpf[n] = lowpass(n, p.M, omega_pass) * kaiser(n, p.M, p.beta);
117 convolve(double *d, int d_len, double *e, int e_len)
119 int w = (e_len - 1) / 2;
121 double *con = calloc (d_len, sizeof (double));
123 for (n = 0; n < d_len; n++) {
126 for (o = -w; o <= w; o++) {
128 double sample = p < 0 ? d[0] : p >= d_len ? d[d_len-1] : d[p];
129 v += sample * e[o + w];
137 cc_low_pass(double *data, int data_len, double omega_pass, double omega_stop, double error)
140 double *fir = make_low_pass_filter(omega_pass, omega_stop, error, &fir_len);
143 result = convolve(data, data_len, fir, fir_len);