# HG changeset patch
# User Chris Cannam
# Date 1171020754 0
# Node ID a867be73b6388cea2461da53c0077c29b6e387d0
# Parent  185454896a76eb75e08f756c12fbded0ca3c917f
* Add non-fftw3 fft alternative

diff -r 185454896a76 -r a867be73b638 data/data.pro
--- a/data/data.pro	Wed Feb 07 14:21:14 2007 +0000
+++ b/data/data.pro	Fri Feb 09 11:32:34 2007 +0000
@@ -14,7 +14,8 @@
 MOC_DIR = tmp_moc
 
 # Input
-HEADERS += fft/FFTCache.h \
+HEADERS += fft/FFTapi.h \
+           fft/FFTCache.h \
            fft/FFTDataServer.h \
            fft/FFTFileCache.h \
            fft/FFTMemoryCache.h \
@@ -51,7 +52,8 @@
            model/TextModel.h \
            model/WaveFileModel.h \
            model/WritableWaveFileModel.h
-SOURCES += fft/FFTDataServer.cpp \
+SOURCES += fft/FFTapi.cpp \
+           fft/FFTDataServer.cpp \
            fft/FFTFileCache.cpp \
            fft/FFTMemoryCache.cpp \
            fileio/AudioFileReader.cpp \
diff -r 185454896a76 -r a867be73b638 data/fft/FFTDataServer.cpp
--- a/data/fft/FFTDataServer.cpp	Wed Feb 07 14:21:14 2007 +0000
+++ b/data/fft/FFTDataServer.cpp	Fri Feb 09 11:32:34 2007 +0000
@@ -583,21 +583,21 @@
     }
 
     m_fftInput = (fftsample *)
-        fftwf_malloc(fftSize * sizeof(fftsample));
+        fftf_malloc(fftSize * sizeof(fftsample));
 
-    m_fftOutput = (fftwf_complex *)
-        fftwf_malloc((fftSize/2 + 1) * sizeof(fftwf_complex));
+    m_fftOutput = (fftf_complex *)
+        fftf_malloc((fftSize/2 + 1) * sizeof(fftf_complex));
 
     m_workbuffer = (float *)
-        fftwf_malloc((fftSize+2) * sizeof(float));
+        fftf_malloc((fftSize+2) * sizeof(float));
 
-    m_fftPlan = fftwf_plan_dft_r2c_1d(m_fftSize,
+    m_fftPlan = fftf_plan_dft_r2c_1d(m_fftSize,
                                       m_fftInput,
                                       m_fftOutput,
                                       FFTW_ESTIMATE);
 
     if (!m_fftPlan) {
-        std::cerr << "ERROR: fftwf_plan_dft_r2c_1d(" << m_windowSize << ") failed!" << std::endl;
+        std::cerr << "ERROR: fftf_plan_dft_r2c_1d(" << m_windowSize << ") failed!" << std::endl;
         throw(0);
     }
 
@@ -642,10 +642,10 @@
     std::cerr << "FFTDataServer(" << this << " [" << (void *)QThread::currentThreadId() << "]): deleteProcessingData" << std::endl;
 #endif
     if (m_fftInput) {
-        fftwf_destroy_plan(m_fftPlan);
-        fftwf_free(m_fftInput);
-        fftwf_free(m_fftOutput);
-        fftwf_free(m_workbuffer);
+        fftf_destroy_plan(m_fftPlan);
+        fftf_free(m_fftInput);
+        fftf_free(m_fftOutput);
+        fftf_free(m_workbuffer);
     }
     m_fftInput = 0;
 }
@@ -1041,7 +1041,7 @@
 	m_fftInput[i + m_fftSize/2] = temp;
     }
 
-    fftwf_execute(m_fftPlan);
+    fftf_execute(m_fftPlan);
 
     fftsample factor = 0.0;
 
diff -r 185454896a76 -r a867be73b638 data/fft/FFTDataServer.h
--- a/data/fft/FFTDataServer.h	Wed Feb 07 14:21:14 2007 +0000
+++ b/data/fft/FFTDataServer.h	Fri Feb 09 11:32:34 2007 +0000
@@ -19,7 +19,7 @@
 #include "base/Window.h"
 #include "base/Thread.h"
 
-#include <fftw3.h>
+#include "FFTapi.h"
 
 #include <QMutex>
 #include <QWaitCondition>
@@ -159,9 +159,9 @@
     QWaitCondition m_condition;
 
     fftsample *m_fftInput;
-    fftwf_complex *m_fftOutput;
+    fftf_complex *m_fftOutput;
     float *m_workbuffer;
-    fftwf_plan m_fftPlan;
+    fftf_plan m_fftPlan;
 
     class FillThread : public Thread
     {
diff -r 185454896a76 -r a867be73b638 data/fft/FFTapi.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/data/fft/FFTapi.cpp	Fri Feb 09 11:32:34 2007 +0000
@@ -0,0 +1,231 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Sonic Visualiser
+    An audio file viewer and annotation editor.
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Chris Cannam and QMUL.
+    FFT code from Don Cross's public domain FFT implementation.
+    
+    This program 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.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#include "FFTapi.h"
+
+#ifndef HAVE_FFTW3F
+
+#include <cmath>
+
+void
+fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
+{
+    if (!ri || !ro || !io) return;
+
+    unsigned int bits;
+    unsigned int i, j, k, m;
+    unsigned int blockSize, blockEnd;
+
+    double tr, ti;
+
+    if (n < 2) return;
+    if (n & (n-1)) return;
+
+    double angle = 2.0 * M_PI;
+    if (inverse) angle = -angle;
+
+    for (i = 0; ; ++i) {
+	if (n & (1 << i)) {
+	    bits = i;
+	    break;
+	}
+    }
+
+    static unsigned int tableSize = 0;
+    static int *table = 0;
+
+    if (tableSize != n) {
+
+	delete[] table;
+
+	table = new int[n];
+
+	for (i = 0; i < n; ++i) {
+	
+	    m = i;
+
+	    for (j = k = 0; j < bits; ++j) {
+		k = (k << 1) | (m & 1);
+		m >>= 1;
+	    }
+
+	    table[i] = k;
+	}
+
+	tableSize = n;
+    }
+
+    if (ii) {
+	for (i = 0; i < n; ++i) {
+	    ro[table[i]] = ri[i];
+	    io[table[i]] = ii[i];
+	}
+    } else {
+	for (i = 0; i < n; ++i) {
+	    ro[table[i]] = ri[i];
+	    io[table[i]] = 0.0;
+	}
+    }
+
+    blockEnd = 1;
+
+    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
+
+	double delta = angle / (double)blockSize;
+	double sm2 = -sin(-2 * delta);
+	double sm1 = -sin(-delta);
+	double cm2 = cos(-2 * delta);
+	double cm1 = cos(-delta);
+	double w = 2 * cm1;
+	double ar[3], ai[3];
+
+	for (i = 0; i < n; i += blockSize) {
+
+	    ar[2] = cm2;
+	    ar[1] = cm1;
+
+	    ai[2] = sm2;
+	    ai[1] = sm1;
+
+	    for (j = i, m = 0; m < blockEnd; j++, m++) {
+
+		ar[0] = w * ar[1] - ar[2];
+		ar[2] = ar[1];
+		ar[1] = ar[0];
+
+		ai[0] = w * ai[1] - ai[2];
+		ai[2] = ai[1];
+		ai[1] = ai[0];
+
+		k = j + blockEnd;
+		tr = ar[0] * ro[k] - ai[0] * io[k];
+		ti = ar[0] * io[k] + ai[0] * ro[k];
+
+		ro[k] = ro[j] - tr;
+		io[k] = io[j] - ti;
+
+		ro[j] += tr;
+		io[j] += ti;
+	    }
+	}
+
+	blockEnd = blockSize;
+    }
+
+    if (inverse) {
+
+	double denom = (double)n;
+
+	for (i = 0; i < n; i++) {
+	    ro[i] /= denom;
+	    io[i] /= denom;
+	}
+    }
+}
+
+struct fftf_plan_ {
+    int size;
+    int inverse;
+    float *real;
+    fftf_complex *cplx;
+};
+
+fftf_plan
+fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned)
+{
+    if (n < 2) return 0;
+    if (n & (n-1)) return 0;
+    
+    fftf_plan_ *plan = new fftf_plan_;
+    plan->size = n;
+    plan->inverse = 0;
+    plan->real = in;
+    plan->cplx = out;
+    return plan;
+}
+
+fftf_plan
+fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned)
+{
+    if (n < 2) return 0;
+    if (n & (n-1)) return 0;
+    
+    fftf_plan_ *plan = new fftf_plan_;
+    plan->size = n;
+    plan->inverse = 1;
+    plan->real = out;
+    plan->cplx = in;
+    return plan;
+}
+
+void
+fftf_destroy_plan(fftf_plan p)
+{
+    delete p;
+}
+
+void
+fftf_execute(const fftf_plan p)
+{
+    //!!! doesn't work yet for inverse transforms
+
+    float *real = p->real;
+    fftf_complex *cplx = p->cplx;
+    int n = p->size;
+    int inv = p->inverse;
+
+    double *ri = new double[n];
+    double *ro = new double[n];
+    double *io = new double[n];
+
+    double *ii = 0;
+    if (inv) ii = new double[n];
+
+    if (!inv) {
+        for (int i = 0; i < n; ++i) {
+            ri[i] = real[i];
+        }
+    } else {
+        for (int i = 0; i < n/2+1; ++i) {
+            ri[i] = cplx[i][0];
+            ii[i] = cplx[i][1];
+            if (i > 0) {
+                ri[n-i] = ri[i];
+                ii[n-i] = -ii[i];
+            }
+        }
+    }
+
+    fft(n, inv, ri, ii, ro, io);
+
+    if (!inv) {
+        for (int i = 0; i < n/2+1; ++i) {
+            cplx[i][0] = ro[i];
+            cplx[i][1] = io[i];
+        }
+    } else {
+        for (int i = 0; i < n; ++i) {
+            real[i] = ro[i];
+        }
+    }
+
+    delete[] ri;
+    delete[] ro;
+    delete[] io;
+    delete[] ii;
+}
+
+#endif
diff -r 185454896a76 -r a867be73b638 data/fft/FFTapi.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/data/fft/FFTapi.h	Fri Feb 09 11:32:34 2007 +0000
@@ -0,0 +1,53 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Sonic Visualiser
+    An audio file viewer and annotation editor.
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Chris Cannam and QMUL.
+    
+    This program 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.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef _FFT_API_H_
+#define _FFT_API_H_
+
+#ifdef HAVE_FFTW3F
+
+#include <fftw3.h>
+
+#define fftf_complex fftwf_complex
+#define fftf_malloc fftwf_malloc
+#define fftf_free fftwf_free
+#define fftf_plan fftwf_plan
+#define fftf_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
+#define fftf_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
+#define fftf_execute fftwf_execute
+#define fftf_destroy_plan fftwf_destroy_plan
+
+#else
+
+// Provide a fallback FFT implementation if FFTW3f is not available.
+
+typedef float fftf_complex[2];
+#define fftf_malloc malloc
+#define fftf_free free
+
+struct fftf_plan_;
+typedef fftf_plan_ *fftf_plan;
+
+fftf_plan fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned);
+fftf_plan fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned);
+void fftf_execute(const fftf_plan p);
+void fftf_destroy_plan(fftf_plan p);
+
+#define FFTW_ESTIMATE 0
+
+#endif
+
+#endif
+