← Back to branch summary

~ubuntu-branches/ubuntu/saucy/solfege/saucy

« back to all changes in this revision

Viewing changes to soundcard/fft.c

  • Committer: Bazaar Package Importer
  • Author(s): Tom Cato Amundsen
  • Date: 2010-03-28 06:34:28 UTC
  • mfrom: (1.1.10 upstream) (2.1.7 sid)
  • Revision ID: james.westby@ubuntu.com-20100328063428-wg2bqvoce2aq4xfb
Tags: 3.15.9-1
* New upstream release.
* Redo packaging. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
soundcard/fft.c
1
 
/* 
2
 
 * fft.c: fft and frequency calculation.
3
 
 * Copyright (C) 1999 Richard Boulton <richard@tartarus.org>
4
 
 * Convolution stuff by Ralph Loader <suckfish@ihug.co.nz>
5
 
 *
6
 
 * August 2000: almost completely rewritten by Fabio Checconi <fchecconi@libero.it>
7
 
 *      ( see src/fft.old.c for the original one, used by xmms-0.9.5.1 ).
8
 
 *
9
 
 *  This program is free software; you can redistribute it and/or modify
10
 
 *  it under the terms of the GNU General Public License as published by
11
 
 *  the Free Software Foundation; either version 2 of the License, or
12
 
 *  (at your option) any later version.
13
 
 *
14
 
 *  This program is distributed in the hope that it will be useful,
15
 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
 
 *  GNU General Public License for more details.
18
 
 *
19
 
 *  You should have received a copy of the GNU General Public License
20
 
 *  along with this program; if not, write to the Free Software
21
 
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22
 
 */
23
 
 
24
 
#ifdef HAVE_CONFIG_H
25
 
#include <config.h>
26
 
#endif
27
 
#ifdef HAVE_FFTW
28
 
#include <rfftw.h>
29
 
#endif
30
 
 
31
 
#include <stdlib.h>
32
 
#include <math.h>
33
 
#include "xmalloc.h"
34
 
#include "fft.h"
35
 
 
36
 
#ifdef M_PI
37
 
#define PI      M_PI
38
 
#else
39
 
#define PI      3.14159265358979323846
40
 
#endif
41
 
 
42
 
unsigned long nn = 0;
43
 
double *window = NULL;
44
 
 
45
 
int     _fft_init               ( unsigned long points );
46
 
 
47
 
#ifndef HAVE_FFTW
48
 
 
49
 
int     _fft_reverse_bits       ( unsigned int initial );
50
 
void    _fft_calculate          ( void );
51
 
 
52
 
double *rout = NULL, *iout = NULL, *costable = NULL, *sintable = NULL;
53
 
unsigned long lognn;
54
 
 
55
 
int *bit_reverse = NULL;
56
 
 
57
 
#else
58
 
 
59
 
int have_plan = 0;
60
 
fftw_real *rin = NULL, *rout = NULL;
61
 
rfftw_plan p;
62
 
 
63
 
#endif
64
 
 
65
 
#ifndef HAVE_FFTW
66
 
 
67
 
int
68
 
_fft_reverse_bits( unsigned int initial )
69
 
{
70
 
        unsigned int reversed = 0, loop;
71
 
 
72
 
        for( loop = 0; loop < lognn; loop++) {
73
 
                reversed <<= 1;
74
 
                reversed += (initial & 1);
75
 
                initial >>= 1;
76
 
        }
77
 
 
78
 
        return reversed;
79
 
}
80
 
 
81
 
void 
82
 
_fft_calculate( void )
83
 
{
84
 
        unsigned int i, j, k;
85
 
        unsigned int exchanges;
86
 
        float fact_real, fact_imag;
87
 
        float tmp_real, tmp_imag;
88
 
        unsigned int factfact;
89
 
    
90
 
        exchanges = 1;
91
 
        factfact = nn / 2;
92
 
 
93
 
        for( i = lognn; i != 0; i-- ) {
94
 
                for( j = 0; j != exchanges; j++ ) {
95
 
                        fact_real = costable[j * factfact];
96
 
                        fact_imag = sintable[j * factfact];
97
 
            
98
 
                        for( k = j; k < nn; k += exchanges << 1) {
99
 
                                int k1 = k + exchanges;
100
 
                                tmp_real = fact_real * rout[k1] - fact_imag * iout[k1];
101
 
                                tmp_imag = fact_real * iout[k1] + fact_imag * rout[k1];
102
 
                                rout[k1] = rout[k] - tmp_real;
103
 
                                iout[k1] = iout[k] - tmp_imag;
104
 
                                rout[k]  += tmp_real;
105
 
                                iout[k]  += tmp_imag;
106
 
                        }
107
 
                }
108
 
                exchanges <<= 1;
109
 
                factfact >>= 1;
110
 
        }
111
 
}
112
 
 
113
 
#define __fft_calculate _fft_calculate();
114
 
 
115
 
#else
116
 
 
117
 
#define __fft_calculate rfftw_one( p, rin, rout );
118
 
 
119
 
#endif
120
 
 
121
 
int
122
 
_fft_init( unsigned long points )
123
 
{
124
 
        long c;
125
 
#ifndef HAVE_FFTW
126
 
 
127
 
        unsigned int i = points, k = 1;
128
 
        float j;
129
 
 
130
 
        while( i > 2 ) {
131
 
                if( i & 0x0001 )
132
 
                        return -1;
133
 
                i >>= 1;
134
 
                k++;
135
 
        }
136
 
 
137
 
        nn = points;
138
 
        lognn = k;
139
 
 
140
 
        if( iout != NULL )
141
 
                free( iout );
142
 
        if( rout != NULL )
143
 
                free( rout );
144
 
        if( sintable != NULL )
145
 
                free( sintable );
146
 
        if( costable != NULL )
147
 
                free( costable );
148
 
        if( bit_reverse != NULL )
149
 
                free( bit_reverse );
150
 
 
151
 
        bit_reverse = (int*) xmalloc( sizeof( int ) * nn );
152
 
        costable = (double*) xmalloc( sizeof( double ) * nn / 2 );
153
 
        sintable = (double*) xmalloc( sizeof( double ) * nn / 2 );
154
 
 
155
 
        rout = (double*) xmalloc ( sizeof( double ) * nn );
156
 
        iout = (double*) xmalloc ( sizeof( double ) * nn );
157
 
 
158
 
        for( i = 0; i < nn; i++ ) {
159
 
                bit_reverse[i] = _fft_reverse_bits( i );
160
 
        }
161
 
        
162
 
        for( i = 0; i < nn / 2; i++ ) {
163
 
                j = 2 * PI * i / nn;
164
 
                costable[i] = cos( j );
165
 
                sintable[i] = sin( j );
166
 
        }
167
 
 
168
 
#else
169
 
 
170
 
        nn = points;
171
 
 
172
 
        if( rout )
173
 
                free( rout );
174
 
        if( rin )
175
 
                free( rin );
176
 
        if( have_plan )
177
 
                rfftw_destroy_plan( p );
178
 
 
179
 
        rin = (fftw_real*) xmalloc( sizeof( fftw_real ) * nn );
180
 
        rout = (fftw_real*) xmalloc( sizeof( fftw_real ) * nn );
181
 
 
182
 
        p = rfftw_create_plan( nn, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE );
183
 
        have_plan = 1;
184
 
 
185
 
#endif
186
 
 
187
 
        if( window )
188
 
                free( window );
189
 
 
190
 
        window = (double*) xmalloc( sizeof( double ) * nn );
191
 
        for ( c = 0; c < nn / 2; c++ )
192
 
                window[c] = window[nn - c - 1] = 0.54 - 0.46 * cos( 2 * PI * c / nn );
193
 
 
194
 
        return 0;
195
 
}
196
 
 
197
 
int
198
 
fft( void *snd_sample, unsigned long size, double *outbuf, int typ )
199
 
{
200
 
        unsigned long n;
201
 
 
202
 
        if( nn != size ) {
203
 
                if( _fft_init( size ) != 0 )
204
 
                        return -1;
205
 
        }
206
 
 
207
 
        if( typ == BYTE ) {
208
 
                for( n = 0; n < nn; n++ ) {
209
 
#ifndef HAVE_FFTW
210
 
                        rout[bit_reverse[n]] = (((unsigned char*)snd_sample)[n] - 128) * window[n];
211
 
                        iout[bit_reverse[n]] = 0;
212
 
#else
213
 
                        rin[n] = (((unsigned char*)snd_sample)[n] - 128) * window[n];
214
 
#endif
215
 
                }
216
 
        } else {
217
 
                for( n = 0; n < nn; n++ ) {
218
 
#ifndef HAVE_FFTW
219
 
                        rout[bit_reverse[n]] = ((unsigned int*)snd_sample)[n];
220
 
                        iout[bit_reverse[n]] = 0;
221
 
#else
222
 
                        rin[n] = ((int*)snd_sample)[n];
223
 
#endif
224
 
                }
225
 
        }
226
 
 
227
 
        __fft_calculate;
228
 
 
229
 
        for( n = 0; n < size / 2; n ++ )
230
 
#ifndef HAVE_FFTW
231
 
                outbuf[n] = rout[n] * rout[n] + iout[n] * iout[n];
232
 
#else
233
 
                outbuf[n] = rout[n] * rout[n] + rout[size - n] * rout[size - n];
234
 
#endif
235
 
 
236
 
        return 0;
237
 
}