← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/iaxmodem/intrepid

« back to all changes in this revision

Viewing changes to lib/spandsp/src/filter_tools.c

  • Committer: Bazaar Package Importer
  • Author(s): Julien BLACHE
  • Date: 2008-02-12 15:29:42 UTC
  • mfrom: (1.1.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20080212152942-28cxxstfy8iujm0p
Tags: 1.1.0~dfsg-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/spandsp/src/filter_tools.c
 
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 General Public License version 2, as
 
17
 * 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 General Public License for more details.
 
23
 *
 
24
 * You should have received a copy of the GNU General Public License
 
25
 * 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.2 2008/01/10 14:06:21 steveu Exp $
 
29
 */
 
30
 
 
31
#ifdef HAVE_CONFIG_H
 
32
#include "config.h"
 
33
#endif
 
34
 
 
35
#include <inttypes.h>
 
36
#include <stdlib.h>
 
37
#include <unistd.h>
 
38
#if defined(HAVE_TGMATH_H)
 
39
#include <tgmath.h>
 
40
#endif
 
41
#if defined(HAVE_MATH_H)
 
42
#include <math.h>
 
43
#endif
 
44
#include <string.h>
 
45
#include <stdio.h>
 
46
#include <time.h>
 
47
#include <fcntl.h>
 
48
 
 
49
#include "spandsp.h"
 
50
#include "filter_tools.h"
 
51
 
 
52
#define FALSE 0
 
53
#define TRUE (!FALSE)
 
54
 
 
55
#define MAXPZ       8192
 
56
#define SEQLEN      8192
 
57
 
 
58
static complex_t circle[SEQLEN/2];
 
59
 
 
60
static __inline__ complex_t expj(double theta)
 
61
{
 
62
    return complex_set(cos(theta), sin(theta));
 
63
}
 
64
/*- End of function --------------------------------------------------------*/
 
65
 
 
66
static __inline__ double fix(double x)
 
67
{
 
68
    /* Nearest integer */
 
69
    return (x >= 0.0)  ?  floor(0.5 + x)  :  -floor(0.5 - x);
 
70
}
 
71
/*- End of function --------------------------------------------------------*/
 
72
 
 
73
static void fftx(complex_t data[], complex_t temp[], int n)
 
74
{
 
75
    int i;
 
76
    int h;
 
77
    int p;
 
78
    int t;
 
79
    int i2;
 
80
    complex_t wkt;
 
81
 
 
82
    if (n > 1)
 
83
    {
 
84
        h = n/2;
 
85
        for (i = 0;  i < h;  i++)
 
86
        {
 
87
            i2 = i*2;
 
88
            temp[i] = data[i2];         /* Even */
 
89
            temp[h + i] = data[i2 + 1]; /* Odd */
 
90
        }
 
91
        fftx(&temp[0], &data[0], h);
 
92
        fftx(&temp[h], &data[h], h);
 
93
        p = 0;
 
94
        t = SEQLEN/n;
 
95
        for (i = 0;  i < h;  i++)
 
96
        {
 
97
            wkt = complex_mul(&circle[p], &temp[h + i]);
 
98
            data[i] = complex_add(&temp[i], &wkt);
 
99
            data[h + i] = complex_sub(&temp[i], &wkt);
 
100
            p += t;
 
101
        }
 
102
    }
 
103
}
 
104
/*- End of function --------------------------------------------------------*/
 
105
 
 
106
void ifft(complex_t data[], int len)
 
107
{
 
108
    int i;
 
109
    double x;
 
110
    complex_t temp[SEQLEN];
 
111
 
 
112
    /* A very slow and clunky FFT, that's just fine for filter design. */
 
113
    for (i = 0;  i < SEQLEN/2;  i++)
 
114
    {
 
115
        x = (2.0*3.1415926535*i)/(double) SEQLEN;
 
116
        circle[i] = expj(x);
 
117
    }
 
118
    fftx(data, temp, len);
 
119
}
 
120
/*- End of function --------------------------------------------------------*/
 
121
 
 
122
void compute_raised_cosine_filter(double coeffs[],
 
123
                                  int len,
 
124
                                  int root,
 
125
                                  int sinc_compensate,
 
126
                                  double alpha,
 
127
                                  double beta)
 
128
{
 
129
    double f;
 
130
    double x;
 
131
    double f1;
 
132
    double f2;
 
133
    double tau;
 
134
    complex_t vec[SEQLEN];
 
135
    int i;
 
136
    int j;
 
137
    int h;
 
138
    
 
139
    f1 = (1.0 - beta)*alpha;
 
140
    f2 = (1.0 + beta)*alpha;
 
141
    tau = 0.5/alpha;
 
142
    /* (Root) raised cosine */
 
143
    for (i = 0;  i <= SEQLEN/2;  i++)
 
144
    {
 
145
        f = (double) i/(double) SEQLEN;
 
146
        if (f <= f1)
 
147
            vec[i] = complex_set(1.0, 0.0);
 
148
        else if (f <= f2)
 
149
            vec[i] = complex_set(0.5*(1.0 + cos((3.1415926535*tau/beta) * (f - f1))), 0.0);
 
150
        else
 
151
            vec[i] = complex_set(0.0, 0.0);
 
152
    }
 
153
    if (root)
 
154
    {
 
155
        for (i = 0;  i <= SEQLEN/2;  i++)
 
156
            vec[i].re = sqrt(vec[i].re);
 
157
    }
 
158
    if (sinc_compensate)
 
159
    {
 
160
        for (i = 1;  i <= SEQLEN/2;  i++)
 
161
            {
 
162
            x = 3.1415926535*(double) i/(double) SEQLEN;
 
163
                vec[i].re *= (x/sin(x));
 
164
            }
 
165
    }
 
166
    for (i = 0;  i <= SEQLEN/2;  i++)
 
167
        vec[i].re *= tau;
 
168
    for (i = 1;  i < SEQLEN/2;  i++)
 
169
        vec[SEQLEN - i] = vec[i];
 
170
    ifft(vec, SEQLEN);
 
171
    h = (len - 1)/2;
 
172
    for (i = 0;  i < len;  i++)
 
173
    {
 
174
        j = (SEQLEN - h + i)%SEQLEN;
 
175
        coeffs[i] = vec[j].re/(double) SEQLEN;
 
176
    }
 
177
}
 
178
/*- End of function --------------------------------------------------------*/
 
179
 
 
180
void compute_hilbert_transform(double coeffs[], int len)
 
181
{
 
182
    double x;
 
183
    int i;
 
184
    int h;
 
185
 
 
186
    h = (len - 1)/2;
 
187
    coeffs[h] = 0.0;
 
188
    for (i = 1;  i <= h;  i++)
 
189
    {
 
190
        if ((i & 1))
 
191
        {
 
192
            x = 1.0/(double) i;
 
193
            coeffs[h + i] = -x;
 
194
            coeffs[h - i] = x;
 
195
        }
 
196
    }
 
197
}
 
198
/*- End of function --------------------------------------------------------*/
 
199
 
 
200
void apply_hamming_window(double coeffs[], int len)
 
201
{
 
202
    double w;
 
203
    int i;
 
204
    int h;
 
205
 
 
206
    h = (len - 1)/2;
 
207
    for (i = 1;  i <= h;  i++)
 
208
    {
 
209
        w = 0.53836 - 0.46164*cos(2.0*3.1415926535*(double) (h + i)/(double) (len - 1.0));
 
210
        coeffs[h + i] *= w;
 
211
        coeffs[h - i] *= w;
 
212
    }
 
213
}
 
214
/*- End of function --------------------------------------------------------*/
 
215
 
 
216
void truncate_coeffs(double coeffs[], int len, int bits, int hilbert)
 
217
{
 
218
    double x;
 
219
    double fac;
 
220
    double max;
 
221
    double scale;
 
222
    int h;
 
223
    int i;
 
224
 
 
225
    fac = pow(2.0, (double) (bits - 1.0));
 
226
    h = (len - 1)/2;
 
227
    max = (hilbert)  ?  coeffs[h - 1]  :  coeffs[h];    /* Max coeff */
 
228
    scale = (fac - 1.0)/(fac*max);
 
229
    for (i = 0;  i < len;  i++)
 
230
    {
 
231
        x = coeffs[i]*scale;           /* Scale coeffs so max is (fac - 1.0)/fac */
 
232
        coeffs[i] = fix(x*fac)/fac;    /* Truncate */
 
233
    }
 
234
}
 
235
/*- End of function --------------------------------------------------------*/
 
236
/*- End of file ------------------------------------------------------------*/