diff options
Diffstat (limited to 'src/lingot-fft.c')
| -rw-r--r-- | src/lingot-fft.c | 177 |
1 files changed, 177 insertions, 0 deletions
diff --git a/src/lingot-fft.c b/src/lingot-fft.c new file mode 100644 index 0000000..18de134 --- /dev/null +++ b/src/lingot-fft.c @@ -0,0 +1,177 @@ +//-*- C++ -*-
+/*
+ * lingot, a musical instrument tuner.
+ *
+ * Copyright (C) 2004-2007 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 <math.h>
+#include <stdlib.h>
+
+#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;
+
+// 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()
+ {
+ 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;
+ }
|