diff options
Diffstat (limited to 'src/lingot-fft.c')
| -rw-r--r-- | src/lingot-fft.c | 167 |
1 files changed, 167 insertions, 0 deletions
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 <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 = 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;
+}
|