/* $Id$ */ /* * Copyright (c) 2010 Dimitri Sokolyuk * * Permission to use, copy, modify, and distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ #include #include #include #include #include #include /* must be prior fftw3 */ #include #include "fft.h" static fftw_plan plan; static fftw_complex *out; static double *in; static size_t sz; static double *window; static double *sq; #define ALPHA 0.53836 #define BETA (1.0 - ALPHA) static double * hamming(size_t n) { double *p; int i; p = calloc(n, sizeof(double)); assert(p); for (i = 0; i < n; i++) { p[i] = ALPHA - BETA * cos(2 * M_PI * i / (n - 1)); p[i] /= INT16_MAX * M_SQRT1_2; } return p; } static double * squares(size_t n) { double *p; int i; p = calloc(n, sizeof(double)); assert(p); for (i = 0; i < n; i++) p[i] = sqrt(i + 1); return p; } void resize_fft(size_t n) { fftw_plan oldplan = plan; double *oldwindow = window; double *oldsq = sq; plan = fftw_plan_dft_r2c_1d(n, in, out, FFTW_MEASURE); window = hamming(n); sq = squares(n / 2); sz = n; fftw_destroy_plan(oldplan); free(oldwindow); free(oldsq); } void init_fft(size_t n) { in = fftw_malloc(n * sizeof(double)); out = fftw_malloc(n * sizeof(fftw_complex) / 2); assert(in && out); resize_fft(n); } void exec_fft(int *io, size_t n) { int i; if (n != sz) resize_fft(n); for (i = 0; i < sz; i++) in[i] = window[i] * io[i]; fftw_execute(plan); for (i = 0; i < sz / 2; i++) io[i] = sq[i] * cabs(out[i]); } void free_fft(void) { fftw_destroy_plan(plan); fftw_free(in); fftw_free(out); free(window); free(sq); }