source: src/fft.c @ ac5f22c

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

src/fft.c: fix norm computation when HAVE_COMPLEX_H is undefined

  • Property mode set to 100644
File size: 5.8 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  /* create plans */
58#ifdef HAVE_COMPLEX_H
59  s->fft_size = size/2+1;
60  s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size);
61  s->pfw = fftw_plan_dft_r2c_1d(size, s->in,  s->specdata, FFTW_ESTIMATE);
62  s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE);
63#else
64  s->fft_size = size;
65  s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size);
66  s->pfw = fftw_plan_r2r_1d(size, s->in,  s->specdata, FFTW_R2HC, FFTW_ESTIMATE);
67  s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE);
68#endif
69  return s;
70}
71
72void del_aubio_fft(aubio_fft_t * s) {
73  /* destroy data */
74  fftw_destroy_plan(s->pfw);
75  fftw_destroy_plan(s->pbw);
76  fftw_free(s->specdata);
77  AUBIO_FREE(s->out);
78  AUBIO_FREE(s->in );
79  AUBIO_FREE(s);
80}
81
82void aubio_fft_do(const aubio_fft_t * s, 
83    const smpl_t * data, fft_data_t * spectrum, const uint_t size) {
84  uint_t i;
85  for (i=0;i<size;i++) s->in[i] = data[i];
86  fftw_execute(s->pfw);
87  for (i=0; i < s->fft_size; i++) spectrum[i] = s->specdata[i];
88}
89
90void aubio_fft_rdo(const aubio_fft_t * s, 
91    const fft_data_t * spectrum, smpl_t * data, const uint_t size) {
92  uint_t i;
93  const smpl_t renorm = 1./(smpl_t)size;
94  for (i=0; i < s->fft_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}
105
106void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
107  uint_t i;
108  for (i=0;i<size/2+1;i++) phas[i] = ARGC(spectrum[i]);
109}
110
111void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
112  uint_t j;
113  for (j=0; j<size/2+1; j++) {
114    spectrum[j]  = CEXPC(I*phas[j]);
115    spectrum[j] *= norm[j];
116  }
117}
118
119#else
120
121void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
122  uint_t i;
123  norm[0] = spectrum[0];
124  for (i=1;i<size/2;i++) norm[i] = SQRT((SQR(spectrum[i]) + SQR(spectrum[size-i])));
125  norm[size/2] = spectrum[size/2];
126}
127
128void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
129  uint_t i;
130  phas[0] = 0;
131  for (i=1;i<size/2+1;i++) phas[i] = atan2f(spectrum[size-i] , spectrum[i]);
132  phas[size/2] = 0;
133}
134
135void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
136  uint_t j;
137  for (j=0; j<size/2+1; j++) {
138    spectrum[j]       = norm[j]*COS(phas[j]);
139  }
140  for (j=1; j<size/2+1; j++) {
141    spectrum[size-j]  = norm[j]*SIN(phas[j]);
142  }
143}
144
145#endif
146
147/* new interface aubio_mfft */
148struct _aubio_mfft_t {
149        aubio_fft_t * fft;      /* fftw interface */
150        fft_data_t ** spec;     /* complex spectral data */
151        uint_t winsize;
152        uint_t channels;
153};
154
155aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels){
156  uint_t i;
157  aubio_mfft_t * fft = AUBIO_NEW(aubio_mfft_t);
158  fft->winsize       = winsize;
159  fft->channels      = channels;
160  fft->fft           = new_aubio_fft(winsize);
161  fft->spec          = AUBIO_ARRAY(fft_data_t*,channels);
162  for (i=0; i < channels; i++)
163    fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize);
164  return fft;
165}
166
167/* execute stft */
168void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){
169  uint_t i=0;
170  /* execute stft */
171  for (i=0; i < fft->channels; i++) {
172    aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize);
173    /* put norm and phase into fftgrain */
174    aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize);
175    aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize);
176  }
177}
178
179/* execute inverse fourier transform */
180void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){
181  uint_t i=0;
182  for (i=0; i < fft->channels; i++) {
183    aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize);
184    aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize);
185  }
186}
187
188void del_aubio_mfft(aubio_mfft_t * fft) {
189  uint_t i;
190  for (i=0; i < fft->channels; i++)
191    AUBIO_FREE(fft->spec[i]);
192  AUBIO_FREE(fft->spec);
193  del_aubio_fft(fft->fft);
194  AUBIO_FREE(fft);       
195}
Note: See TracBrowser for help on using the repository browser.