From c04af18902322cea1bd05b6f01b7860ca7bfa432 Mon Sep 17 00:00:00 2001 From: Piotr Pawlow Date: Thu, 13 Mar 2014 20:51:31 +0100 Subject: - import version 0.9.1 from upstream --- src/lingot-fft.c | 167 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 src/lingot-fft.c (limited to 'src/lingot-fft.c') diff --git a/src/lingot-fft.c b/src/lingot-fft.c new file mode 100644 index 0000000..4d21b0a --- /dev/null +++ b/src/lingot-fft.c @@ -0,0 +1,167 @@ +/* + * lingot, a musical instrument tuner. + * + * Copyright (C) 2004-2011 Ibán Cereijo Graña, Jairo Chapela Martínez. + * + * This file is part of lingot. + * + * lingot is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * lingot is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with lingot; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include +#include + +#include "lingot-fft.h" +#include "lingot-config.h" + +/* + DFT functions. + */ + +#ifndef LIB_FFTW + +#include "lingot-complex.h" + +// phase factor table, for FFT optimization. +LingotComplex* wn = NULL; + +// creates the phase factor table. +void lingot_fft_create_phase_factors(LingotConfig* conf) { + FLT alpha; + unsigned int i; + + wn = (LingotComplex*) malloc((conf->fft_size >> 1) * sizeof(LingotComplex)); + + for (i = 0; i < (conf->fft_size >> 1); i++) { + alpha = -2.0 * i * M_PI / conf->fft_size; + wn[i].r = cos(alpha); + wn[i].i = sin(alpha); + } +} + +void lingot_fft_destroy_phase_factors() { + if (wn != NULL) { + free(wn); + } +} + +//------------------------------------------------------------------------ + +// Fast Fourier Transform. +void _lingot_fft_fft(FLT* in, LingotComplex* out, unsigned long int N, + unsigned long int offset, unsigned long int d1, unsigned long int step) { + LingotComplex X1, X2; + unsigned long int Np2 = (N >> 1); // N/2 + register unsigned long int a, b, c, q; + + if (N == 2) { // butterfly for N = 2; + + X1.r = in[offset]; + X1.i = 0.0; + X2.r = in[offset + step]; + X2.i = 0.0; + + lingot_complex_add(&X1, &X2, &out[d1]); + lingot_complex_sub(&X1, &X2, &out[d1 + Np2]); + + return; + } + + _lingot_fft_fft(in, out, Np2, offset, d1, step << 1); + _lingot_fft_fft(in, out, Np2, offset + step, d1 + Np2, step << 1); + + for (q = 0, c = 0; q < (N >> 1); q++, c += step) { + + a = q + d1; + b = a + Np2; + + X1 = out[a]; + lingot_complex_mul(&out[b], &wn[c], &X2); + lingot_complex_add(&X1, &X2, &out[a]); + lingot_complex_sub(&X1, &X2, &out[b]); + } +} + +void lingot_fft_fft(FLT* in, LingotComplex* out, unsigned long int N) { + _lingot_fft_fft(in, out, N, 0, 0, 1); +} + +#endif + +//------------------------------------------------------------------------ + + +/* Spectral Power Distribution esteem, selectively in frequency, by DFT. + transforms signal in of N1 samples from frequency wi, with sample + separation of dw rads, storing the result on buffer out with N2 samples. */ +void lingot_fft_spd(FLT* in, int N1, FLT wi, FLT dw, FLT* out, int N2) { + FLT Xr, Xi; + FLT wn; + FLT N1_2 = N1 * N1; + register int i, n; + + for (i = 0; i < N2; i++) { + + Xr = 0.0; + Xi = 0.0; + + for (n = 0; n < N1; n++) { // O(n1*n2) :( + + wn = (wi + dw * i) * n; + Xr = Xr + cos(wn) * in[n]; + Xi = Xi - sin(wn) * in[n]; + } + + out[i] = (Xr * Xr + Xi * Xi) / N1_2; // normalized squared module. + } +} + +//------------------------------------------------------------------------ + +/* + Evaluates 1st and 2nd derivatives from SPD at frequency w. + */ +void lingot_fft_spd_diffs(FLT* in, int N, FLT w, FLT* out_d1, FLT* out_d2) { + FLT x_cos_wn; + FLT x_sin_wn; + register int n; + + FLT SUM_x_sin_wn = 0.0; + FLT SUM_x_cos_wn = 0.0; + FLT SUM_x_n_sin_wn = 0.0; + FLT SUM_x_n_cos_wn = 0.0; + FLT SUM_x_n2_sin_wn = 0.0; + FLT SUM_x_n2_cos_wn = 0.0; + + for (n = 0; n < N; n++) { + + x_cos_wn = in[n] * cos(w * n); + x_sin_wn = in[n] * sin(w * n); + + SUM_x_sin_wn += x_sin_wn; + SUM_x_cos_wn += x_cos_wn; + SUM_x_n_sin_wn += x_sin_wn * n; + SUM_x_n_cos_wn += x_cos_wn * n; + SUM_x_n2_sin_wn += x_sin_wn * n * n; + SUM_x_n2_cos_wn += x_cos_wn * n * n; + } + + FLT N_2 = N * N; + *out_d1 = 2.0 * (SUM_x_sin_wn * SUM_x_n_cos_wn - SUM_x_cos_wn + * SUM_x_n_sin_wn) / N_2; + *out_d2 = 2.0 * (SUM_x_n_cos_wn * SUM_x_n_cos_wn - SUM_x_sin_wn + * SUM_x_n2_sin_wn + SUM_x_n_sin_wn * SUM_x_n_sin_wn - SUM_x_cos_wn + * SUM_x_n2_cos_wn) / N_2; +} -- cgit v1.2.3