Mercurial > hg > audiostuff

comparison spandsp-0.0.6pre17/src/filter_tools.c @ 4:26cd8f1ef0b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import spandsp-0.0.6pre17
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 15:50:58 +0200
parents
children
comparison
equal deleted inserted replaced
3:c6c5a16ce2f2 4:26cd8f1ef0b1
1 /*
2 * SpanDSP - a series of DSP components for telephony
3 *
4 * filter_tools.c - A collection of routines used for filter design.
5 *
6 * Written by Steve Underwood <steveu@coppice.org>
7 *
8 * Copyright (C) 2008 Steve Underwood
9 *
10 * This includes some elements based on the mkfilter package by
11 * A.J. Fisher, University of York <fisher@minster.york.ac.uk>, November 1996
12 *
13 * All rights reserved.
14 *
15 * This program is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU Lesser General Public License version 2.1,
17 * as published by the Free Software Foundation.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Lesser General Public License for more details.
23 *
24 * You should have received a copy of the GNU Lesser General Public
25 * License along with this program; if not, write to the Free Software
26 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
27 *
28 * $Id: filter_tools.c,v 1.11 2009/10/05 16:33:25 steveu Exp $
29 */
30
31 #if defined(HAVE_CONFIG_H)
32 #include "config.h"
33 #endif
34
35 #include <inttypes.h>
36 #include <stdlib.h>
37 #if defined(HAVE_TGMATH_H)
38 #include <tgmath.h>
39 #endif
40 #if defined(HAVE_MATH_H)
41 #include <math.h>
42 #endif
43 #include "floating_fudge.h"
44 #include <string.h>
45 #include <stdio.h>
46 #include <time.h>
47 #include <fcntl.h>
48
49 #include "spandsp/telephony.h"
50 #include "spandsp/complex.h"
51 #include "filter_tools.h"
52
53 #if !defined(FALSE)
54 #define FALSE 0
55 #endif
56 #if !defined(TRUE)
57 #define TRUE (!FALSE)
58 #endif
59
60 #define MAXPZ 8192
61 #define SEQ_LEN 8192
62 #define MAX_FFT_LEN SEQ_LEN
63
64 static complex_t circle[MAX_FFT_LEN/2];
65
66 static __inline__ complex_t expj(double theta)
67 {
68 return complex_set(cos(theta), sin(theta));
69 }
70 /*- End of function --------------------------------------------------------*/
71
72 static __inline__ double fix(double x)
73 {
74 /* Nearest integer */
75 return (x >= 0.0) ? floor(0.5 + x) : -floor(0.5 - x);
76 }
77 /*- End of function --------------------------------------------------------*/
78
79 static void fftx(complex_t data[], complex_t temp[], int n)
80 {
81 int i;
82 int h;
83 int p;
84 int t;
85 int i2;
86 complex_t wkt;
87
88 if (n > 1)
89 {
90 h = n/2;
91 for (i = 0; i < h; i++)
92 {
93 i2 = i*2;
94 temp[i] = data[i2]; /* Even */
95 temp[h + i] = data[i2 + 1]; /* Odd */
96 }
97 fftx(&temp[0], &data[0], h);
98 fftx(&temp[h], &data[h], h);
99 p = 0;
100 t = MAX_FFT_LEN/n;
101 for (i = 0; i < h; i++)
102 {
103 wkt = complex_mul(&circle[p], &temp[h + i]);
104 data[i] = complex_add(&temp[i], &wkt);
105 data[h + i] = complex_sub(&temp[i], &wkt);
106 p += t;
107 }
108 }
109 }
110 /*- End of function --------------------------------------------------------*/
111
112 void ifft(complex_t data[], int len)
113 {
114 int i;
115 double x;
116 complex_t temp[MAX_FFT_LEN];
117
118 /* A very slow and clunky FFT, that's just fine for filter design. */
119 for (i = 0; i < MAX_FFT_LEN/2; i++)
120 {
121 x = (2.0*3.1415926535*i)/(double) MAX_FFT_LEN;
122 circle[i] = expj(x);
123 }
124 fftx(data, temp, len);
125 }
126 /*- End of function --------------------------------------------------------*/
127
128 void compute_raised_cosine_filter(double coeffs[],
129 int len,
130 int root,
131 int sinc_compensate,
132 double alpha,
133 double beta)
134 {
135 double f;
136 double x;
137 double f1;
138 double f2;
139 double tau;
140 complex_t vec[SEQ_LEN];
141 int i;
142 int j;
143 int h;
144
145 f1 = (1.0 - beta)*alpha;
146 f2 = (1.0 + beta)*alpha;
147 tau = 0.5/alpha;
148 /* (Root) raised cosine */
149 for (i = 0; i <= SEQ_LEN/2; i++)
150 {
151 f = (double) i/(double) SEQ_LEN;
152 if (f <= f1)
153 vec[i] = complex_set(1.0, 0.0);
154 else if (f <= f2)
155 vec[i] = complex_set(0.5*(1.0 + cos((3.1415926535*tau/beta)*(f - f1))), 0.0);
156 else
157 vec[i] = complex_set(0.0, 0.0);
158 }
159 if (root)
160 {
161 for (i = 0; i <= SEQ_LEN/2; i++)
162 vec[i].re = sqrt(vec[i].re);
163 }
164 if (sinc_compensate)
165 {
166 for (i = 1; i <= SEQ_LEN/2; i++)
167 {
168 x = 3.1415926535*(double) i/(double) SEQ_LEN;
169 vec[i].re *= (x/sin(x));
170 }
171 }
172 for (i = 0; i <= SEQ_LEN/2; i++)
173 vec[i].re *= tau;
174 for (i = 1; i < SEQ_LEN/2; i++)
175 vec[SEQ_LEN - i] = vec[i];
176 ifft(vec, SEQ_LEN);
177 h = (len - 1)/2;
178 for (i = 0; i < len; i++)
179 {
180 j = (SEQ_LEN - h + i)%SEQ_LEN;
181 coeffs[i] = vec[j].re/(double) SEQ_LEN;
182 }
183 }
184 /*- End of function --------------------------------------------------------*/
185
186 void compute_hilbert_transform(double coeffs[], int len)
187 {
188 double x;
189 int i;
190 int h;
191
192 h = (len - 1)/2;
193 coeffs[h] = 0.0;
194 for (i = 1; i <= h; i++)
195 {
196 if ((i & 1))
197 {
198 x = 1.0/(double) i;
199 coeffs[h + i] = -x;
200 coeffs[h - i] = x;
201 }
202 }
203 }
204 /*- End of function --------------------------------------------------------*/
205
206 void apply_hamming_window(double coeffs[], int len)
207 {
208 double w;
209 int i;
210 int h;
211
212 h = (len - 1)/2;
213 for (i = 1; i <= h; i++)
214 {
215 w = 0.53836 - 0.46164*cos(2.0*3.1415926535*(double) (h + i)/(double) (len - 1.0));
216 coeffs[h + i] *= w;
217 coeffs[h - i] *= w;
218 }
219 }
220 /*- End of function --------------------------------------------------------*/
221
222 void truncate_coeffs(double coeffs[], int len, int bits, int hilbert)
223 {
224 double x;
225 double fac;
226 double max;
227 double scale;
228 int h;
229 int i;
230
231 fac = pow(2.0, (double) (bits - 1.0));
232 h = (len - 1)/2;
233 max = (hilbert) ? coeffs[h - 1] : coeffs[h]; /* Max coeff */
234 scale = (fac - 1.0)/(fac*max);
235 for (i = 0; i < len; i++)
236 {
237 x = coeffs[i]*scale; /* Scale coeffs so max is (fac - 1.0)/fac */
238 coeffs[i] = fix(x*fac)/fac; /* Truncate */
239 }
240 }
241 /*- End of function --------------------------------------------------------*/
242 /*- End of file ------------------------------------------------------------*/

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.