source: src/fft.c @ 237f632

feature/autosinkfeature/constantqfeature/pitchshiftfeature/pydocstringsfeature/timestretchpitchshiftsamplertimestretchyinfft+
Last change on this file since 237f632 was 237f632, checked in by Paul Brossier <piem@altern.org>, 14 years ago

use replacement if complex.h not found
use replacement if complex.h not found

  • Property mode set to 100644
File size: 5.9 KB
Line 
1/*
2   Copyright (C) 2003 Paul Brossier
3
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.
8
9   This program is distributed in the hope that it will be useful,
10   but WITHOUT ANY WARRANTY; without even the implied warranty of
11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12   GNU General Public License for more details.
13
14   You should have received a copy of the GNU General Public License
15   along with this program; if not, write to the Free Software
16   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17
18*/
19
20#include "aubio_priv.h"
21#include "sample.h"
22#include "mathutils.h"
23#include "fft.h"
24
25#if FFTW3F_SUPPORT
26#define fftw_malloc             fftwf_malloc
27#define fftw_free               fftwf_free
28#define fftw_execute            fftwf_execute
29#define fftw_plan_dft_r2c_1d    fftwf_plan_dft_r2c_1d
30#define fftw_plan_dft_c2r_1d    fftwf_plan_dft_c2r_1d
31#define fftw_plan_r2r_1d      fftwf_plan_r2r_1d
32#define fftw_plan               fftwf_plan
33#define fftw_destroy_plan       fftwf_destroy_plan
34#endif
35
36#if FFTW3F_SUPPORT
37#define real_t smpl_t
38#else
39#define real_t lsmp_t
40#endif
41
42struct _aubio_fft_t {
43        uint_t fft_size;
44        uint_t channels;
45        real_t          *in, *out;
46        fft_data_t      *specdata;
47        fftw_plan       pfw, pbw;
48};
49
50static void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size);
51
52aubio_fft_t * new_aubio_fft(uint_t size) {
53        aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);
54        /* allocate memory */
55        s->in       = AUBIO_ARRAY(real_t,size);
56        s->out      = AUBIO_ARRAY(real_t,size);
57        s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*size);
58        /* create plans */
59#ifdef HAVE_COMPLEX_H
60        s->pfw = fftw_plan_dft_r2c_1d(size, s->in,  s->specdata, FFTW_ESTIMATE);
61        s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE);
62#else
63        s->pfw = fftw_plan_r2r_1d(size, s->in,  s->specdata, FFTW_R2HC, FFTW_ESTIMATE);
64        s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE);
65#endif
66        return s;
67}
68
69void del_aubio_fft(aubio_fft_t * s) {
70        /* destroy data */
71        fftw_destroy_plan(s->pfw);
72        fftw_destroy_plan(s->pbw);
73        fftw_free(s->specdata);
74        AUBIO_FREE(s->out);
75        AUBIO_FREE(s->in );
76        AUBIO_FREE(s);
77}
78
79void aubio_fft_do(const aubio_fft_t * s, 
80                const smpl_t * data, fft_data_t * spectrum, 
81                const uint_t size) {
82        uint_t i;
83        for (i=0;i<size;i++) s->in[i] = data[i];
84        fftw_execute(s->pfw);
85        for (i=0;i<size;i++) spectrum[i] = s->specdata[i];
86}
87
88void aubio_fft_rdo(const aubio_fft_t * s, 
89                const fft_data_t * spectrum, 
90                smpl_t * data, 
91                const uint_t size) {
92        uint_t i;
93        const smpl_t renorm = 1./(smpl_t)size;
94        for (i=0;i<size;i++) s->specdata[i] = spectrum[i];
95        fftw_execute(s->pbw);
96        for (i=0;i<size;i++) data[i] = s->out[i]*renorm;
97}
98
99#ifdef HAVE_COMPLEX_H
100
101void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
102        uint_t i;
103        for (i=0;i<size/2+1;i++) norm[i] = ABSC(spectrum[i]);
104        //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]);
105}
106
107void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
108        uint_t i;
109        for (i=0;i<size/2+1;i++) phas[i] = ARGC(spectrum[i]);
110        //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]);
111}
112
113void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
114  uint_t j;
115  for (j=0; j<size/2+1; j++) {
116    spectrum[j]  = CEXPC(I*phas[j]);
117    spectrum[j] *= norm[j];
118  }
119}
120
121#else
122
123void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
124        uint_t i;
125  norm[0] = -spectrum[0];
126        for (i=1;i<size/2+1;i++) norm[i] = SQRT(SQR(spectrum[i]) + SQR(spectrum[size-i]));
127        //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]);
128}
129
130void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
131        uint_t i;
132  phas[0] = PI;
133        for (i=1;i<size/2+1;i++) phas[i] = atan2f(spectrum[size-i] , spectrum[i]);
134        //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]);
135}
136
137void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
138  uint_t j;
139  for (j=0; j<size/2+1; j++) {
140    spectrum[j]       = norm[j]*COS(phas[j]);
141  }
142  for (j=1; j<size/2+1; j++) {
143    spectrum[size-j]  = norm[j]*SIN(phas[j]);
144  }
145}
146
147#endif
148
149/* new interface aubio_mfft */
150struct _aubio_mfft_t {
151        aubio_fft_t * fft;      /* fftw interface */
152        fft_data_t ** spec;     /* complex spectral data */
153        uint_t winsize;
154        uint_t channels;
155};
156
157aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels){
158        uint_t i;
159        aubio_mfft_t * fft = AUBIO_NEW(aubio_mfft_t);
160        fft->winsize       = winsize;
161        fft->channels      = channels;
162        fft->fft           = new_aubio_fft(winsize);
163        fft->spec          = AUBIO_ARRAY(fft_data_t*,channels);
164        for (i=0; i < channels; i++)
165                fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize);
166        return fft;
167}
168
169/* execute stft */
170void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){
171        uint_t i=0;
172        /* execute stft */
173        for (i=0; i < fft->channels; i++) {
174                aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize);
175                /* put norm and phase into fftgrain */
176                aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize);
177                aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize);
178        }
179}
180
181/* execute inverse fourier transform */
182void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){
183        uint_t i=0;
184        for (i=0; i < fft->channels; i++) {
185                aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize);
186                aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize);
187        }
188}
189
190void del_aubio_mfft(aubio_mfft_t * fft) {
191        uint_t i;
192        for (i=0; i < fft->channels; i++)
193                AUBIO_FREE(fft->spec[i]);
194        AUBIO_FREE(fft->spec);
195        del_aubio_fft(fft->fft);
196        AUBIO_FREE(fft);       
197}
Note: See TracBrowser for help on using the repository browser.