# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1400491311 -3600
# Node ID 6b13f9c694a83a5c101b2b8cf0f403dc547e2cb6
# Parent  1e33f719dde181e2d0cfc651ecd8b63b0da52d94
Start introducing the tests

diff -r 1e33f719dde1 -r 6b13f9c694a8 .hgignore
--- a/.hgignore	Fri May 16 13:48:57 2014 +0100
+++ b/.hgignore	Mon May 19 10:21:51 2014 +0100
@@ -2,5 +2,10 @@
 *~
 *.o
 *.so
-cpp-qm-dsp/test
+*.a
+*.dylib
+*.dll
 *.class
+
+syntax: re
+^test/Test[^.]*$
diff -r 1e33f719dde1 -r 6b13f9c694a8 Makefile.inc
--- a/Makefile.inc	Fri May 16 13:48:57 2014 +0100
+++ b/Makefile.inc	Mon May 19 10:21:51 2014 +0100
@@ -11,8 +11,10 @@
 
 CXX	?= g++
 CC	?= gcc
+VALGRIND	?= valgrind -q
 
 GENERAL_FLAGS	:= -I. -I$(VAMPSDK_DIR) -I$(INC_DIR) -I$(LIB_DIR) -I$(KFFT_DIR) -I$(KFFT_DIR)/tools -Dkiss_fft_scalar=double
+
 CFLAGS := $(GENERAL_FLAGS) $(CFLAGS) 
 CXXFLAGS := $(GENERAL_FLAGS) $(CXXFLAGS)
 
@@ -23,7 +25,6 @@
 
 LIB	:= libcq.a
 PLUGIN	:= cqvamp$(PLUGIN_EXT)
-TEST	:= $(TEST_DIR)/test
 PF	:= $(TEST_DIR)/processfile
 
 LIB_HEADERS	:= \
@@ -66,29 +67,35 @@
 	$(VAMP_DIR)/libmain.cpp \
 	$(VAMP_DIR)/Pitch.cpp
 
+TEST_SOURCES	:= \
+	$(TEST_DIR)/TestCQFrequency.cpp \
+	$(TEST_DIR)/TestFFT.cpp \
+	$(TEST_DIR)/TestMathUtilities.cpp \
+	$(TEST_DIR)/TestResampler.cpp \
+	$(TEST_DIR)/TestWindow.cpp
+
 HEADERS	     := $(LIB_HEADERS) $(VAMP_HEADERS)
 SOURCES	     := $(LIB_SOURCES) $(VAMP_SOURCES)
+
 LIB_OBJECTS  := $(LIB_SOURCES:.cpp=.o)
 LIB_OBJECTS  := $(LIB_OBJECTS:.c=.o)
+
 OBJECTS	     := $(SOURCES:.cpp=.o)
 OBJECTS	     := $(OBJECTS:.c=.o)
 
-TEST_SOURCES := $(TEST_DIR)/test.cpp 
-TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o) $(OBJECTS)
+TEST_TARGETS	:= $(TEST_SOURCES:.cpp=)
 
 PF_SOURCES := $(TEST_DIR)/processfile.cpp
 PF_OBJECTS := $(PF_SOURCES:.cpp=.o) $(OBJECTS)
 
 LIBS	:= $(VAMPSDK_DIR)/libvamp-sdk.a -lpthread
 
-all: $(LIB) $(PLUGIN) $(TEST) $(PF)
+all: $(LIB) $(PLUGIN) $(TEST_TARGETS) $(PF)
+	for t in $(TEST_TARGETS); do echo; echo "Running $$t"; $(VALGRIND) ./"$$t" || exit 1; done
 
 $(PLUGIN):	$(OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(PLUGIN_LDFLAGS)
 
-$(TEST):	$(TEST_OBJECTS)
-	$(CXX) -o $@ $^ $(LIBS) $(TEST_LDFLAGS)
-
 $(PF):	$(PF_OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(PF_LDFLAGS)
 
@@ -96,6 +103,9 @@
 	ar cr $@ $^
 	ranlib $@
 
+Test%:	Test%.o
+	$(CXX) -o $@ $^ $(LIB) $(LIBS) $(TEST_LDFLAGS)
+
 clean:		
 	rm -f $(OBJECTS) $(TEST_OBJECTS) $(PF_OBJECTS)
 
diff -r 1e33f719dde1 -r 6b13f9c694a8 test/TestCQFrequency.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestCQFrequency.cpp	Mon May 19 10:21:51 2014 +0100
@@ -0,0 +1,74 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "cq/CQSpectrogram.h"
+
+#include "dsp/Window.h"
+
+#include <cmath>
+#include <vector>
+
+using std::vector;
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestCQFrequency)
+
+// The principle here is to feed a single windowed sinusoid into a
+// small CQ transform and check that the output has its peak bin at
+// the correct frequency. We can repeat for different frequencies both
+// inside and outside the frequency range supported by the CQ. We
+// should also repeat for CQSpectrogram outputs as well as the raw CQ.
+
+// Set up fs/2 = 50, frequency range 10 -> 40 i.e. 2 octaves, fixed
+// duration of 2 seconds
+static const double sampleRate = 100;
+static const double cqmin = 10;
+static const double cqmax = 40;
+static const double bpo = 4;
+static const int duration = sampleRate * 2;
+
+// Fairly arbitrary max value for CQ bins other than the "correct" one
+static const double threshold = 0.08;
+
+void
+checkCQFreqOutput(const CQSpectrogram::RealBlock &output, double freq)
+{
+    
+}
+
+void
+testCQFrequency(double freq)
+{
+    CQParameters params(sampleRate, cqmin, cqmax, bpo);
+    CQSpectrogram cq(params, CQSpectrogram::InterpolateLinear);
+
+    vector<double> input;
+    for (int i = 0; i < duration; ++i) {
+        input.push_back(sin((i * 2 * M_PI * freq) / sampleRate));
+    }
+    Window<double>(HanningWindow, duration).cut(input.data());
+
+    CQSpectrogram::RealBlock output = cq.process(input);
+    CQSpectrogram::RealBlock rest = cq.getRemainingOutput();
+    output.insert(output.end(), rest.begin(), rest.end());
+
+    checkCQFreqOutput(output, freq);
+}
+
+BOOST_AUTO_TEST_CASE(freq_0) { testCQFrequency(0); }
+BOOST_AUTO_TEST_CASE(freq_5) { testCQFrequency(5); }
+BOOST_AUTO_TEST_CASE(freq_10) { testCQFrequency(10); }
+BOOST_AUTO_TEST_CASE(freq_15) { testCQFrequency(15); }
+BOOST_AUTO_TEST_CASE(freq_20) { testCQFrequency(20); }
+BOOST_AUTO_TEST_CASE(freq_25) { testCQFrequency(25); }
+BOOST_AUTO_TEST_CASE(freq_30) { testCQFrequency(30); }
+BOOST_AUTO_TEST_CASE(freq_35) { testCQFrequency(35); }
+BOOST_AUTO_TEST_CASE(freq_40) { testCQFrequency(40); }
+BOOST_AUTO_TEST_CASE(freq_45) { testCQFrequency(45); }
+BOOST_AUTO_TEST_CASE(freq_50) { testCQFrequency(50); }
+
+BOOST_AUTO_TEST_SUITE_END()
+
diff -r 1e33f719dde1 -r 6b13f9c694a8 test/TestFFT.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestFFT.cpp	Mon May 19 10:21:51 2014 +0100
@@ -0,0 +1,466 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/FFT.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+#include <stdexcept>
+
+BOOST_AUTO_TEST_SUITE(TestFFT)
+
+#define COMPARE_CONST(a, n) \
+    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
+        BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14);				\
+    }
+
+#define COMPARE_ARRAY(a, b)						\
+    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
+        BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14);			\
+    }
+
+//!!! need at least one test with complex time-domain signal
+
+BOOST_AUTO_TEST_CASE(forwardArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double re[] = { 999, 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re+1, im+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(re[0], 999.0);
+    BOOST_CHECK_EQUAL(im[0], 999.0);
+    BOOST_CHECK_EQUAL(re[5], 999.0);
+    BOOST_CHECK_EQUAL(im[5], 999.0);
+    delete[] in;
+}
+
+BOOST_AUTO_TEST_CASE(r_forwardArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double re[] = { 999, 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re+1, im+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(re[0], 999.0);
+    BOOST_CHECK_EQUAL(im[0], 999.0);
+    BOOST_CHECK_EQUAL(re[5], 999.0);
+    BOOST_CHECK_EQUAL(im[5], 999.0);
+    delete[] in;
+}
+
+BOOST_AUTO_TEST_CASE(inverseArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *re = new double[4];
+    double *im = new double[4];
+    re[0] = 0;
+    re[1] = 1;
+    re[2] = 0;
+    re[3] = 1;
+    im[0] = 0;
+    im[1] = -2;
+    im[2] = 0;
+    im[3] = 2;
+    double outre[] = { 999, 999, 999, 999, 999, 999 };
+    double outim[] = { 999, 999, 999, 999, 999, 999 };
+    FFT(4).process(true, re, im, outre+1, outim+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(outre[0], 999.0);
+    BOOST_CHECK_EQUAL(outim[0], 999.0);
+    BOOST_CHECK_EQUAL(outre[5], 999.0);
+    BOOST_CHECK_EQUAL(outim[5], 999.0);
+    delete[] re;
+    delete[] im;
+}
+
+BOOST_AUTO_TEST_CASE(r_inverseArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *re = new double[3];
+    double *im = new double[3];
+    re[0] = 0;
+    re[1] = 1;
+    re[2] = 0;
+    im[0] = 0;
+    im[1] = -2;
+    im[2] = 0;
+    double outre[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(4).inverse(re, im, outre+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(outre[0], 999.0);
+    BOOST_CHECK_EQUAL(outre[5], 999.0);
+    delete[] re;
+    delete[] im;
+}
+
+BOOST_AUTO_TEST_CASE(dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double in[] = { 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double in[] = { 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(c_dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double rin[] = { 1, 1, 1, 1 };
+    double iin[] = { 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, rin, iin, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    BOOST_CHECK_EQUAL(im[0], 4.0);
+    BOOST_CHECK_EQUAL(im[1], 0.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, rin);
+    COMPARE_ARRAY(backim, iin);
+}
+
+BOOST_AUTO_TEST_CASE(sine)
+{
+    // Sine. Output is purely imaginary
+    double in[] = { 0, 1, 0, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    COMPARE_CONST(re, 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], -2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 2.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_sine)
+{
+    // Sine. Output is purely imaginary
+    double in[] = { 0, 1, 0, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    COMPARE_CONST(re, 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], -2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 2.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(cosine)
+{
+    // Cosine. Output is purely real
+    double in[] = { 1, 0, -1, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 2.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_cosine)
+{
+    // Cosine. Output is purely real
+    double in[] = { 1, 0, -1, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 2.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(c_cosine)
+{
+    // Cosine. Output is purely real
+    double rin[] = { 1, 0, -1, 0 };
+    double iin[] = { 1, 0, -1, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, rin, iin, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 2.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], 2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 2.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, rin);
+    COMPARE_ARRAY(backim, iin);
+}
+	
+BOOST_AUTO_TEST_CASE(sineCosine)
+{
+    // Sine and cosine mixed
+    double in[] = { 0.5, 1, -0.5, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+	
+BOOST_AUTO_TEST_CASE(r_sineCosine)
+{
+    // Sine and cosine mixed
+    double in[] = { 0.5, 1, -0.5, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+	
+BOOST_AUTO_TEST_CASE(c_sineCosine)
+{
+    double rin[] = { 1, 0, -1, 0 };
+    double iin[] = { 0, 1, 0, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, rin, iin, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 4.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, rin);
+    COMPARE_ARRAY(backim, iin);
+}
+
+BOOST_AUTO_TEST_CASE(nyquist)
+{
+    double in[] = { 1, -1, 1, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 4.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_nyquist)
+{
+    double in[] = { 1, -1, 1, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 4.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 1.0);
+    BOOST_CHECK_EQUAL(re[1], 1.0);
+    BOOST_CHECK_EQUAL(re[2], 1.0);
+    BOOST_CHECK_EQUAL(re[3], 1.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 1.0);
+    BOOST_CHECK_EQUAL(re[1], 1.0);
+    BOOST_CHECK_EQUAL(re[2], 1.0);
+    BOOST_CHECK_EQUAL(re[3], 1.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(sizes) 
+{
+    // Complex supports any size. A single test with an odd size
+    // will do here, without getting too much into our expectations
+    // about supported butterflies etc
+
+    double in[] = { 1, 1, 1 };
+    double re[] = { 999, 999, 999 };
+    double im[] = { 999, 999, 999 };
+    FFT(3).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 3.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[3];
+    double backim[3];
+    FFT(3).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_sizes)
+{
+    // Real supports any even size, but not odd ones
+
+    BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument);
+
+    double in[] = { 1, 1, 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(6).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 6.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    BOOST_CHECK_EQUAL(re[4], 0.0);
+    BOOST_CHECK_EQUAL(re[5], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[6];
+    FFTReal(6).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
diff -r 1e33f719dde1 -r 6b13f9c694a8 test/TestMathUtilities.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestMathUtilities.cpp	Mon May 19 10:21:51 2014 +0100
@@ -0,0 +1,153 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/MathUtilities.h"
+
+#include <cmath>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestMathUtilities)
+
+BOOST_AUTO_TEST_CASE(round)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::round(0.5), 1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(0.49), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(0.99), 1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(0.01), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(0.0), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(100.0), 100.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-0.2), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-0.5), -1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-0.99), -1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-1.0), -1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-1.1), -1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::round(-1.5), -2.0);
+}
+
+BOOST_AUTO_TEST_CASE(mean)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::mean(0, 0), 0);
+    double d0[] = { 0, 4, 3, -1 };
+    BOOST_CHECK_EQUAL(MathUtilities::mean(d0, 4), 1.5);
+    double d1[] = { -2.6 };
+    BOOST_CHECK_EQUAL(MathUtilities::mean(d1, 1), -2.6);
+    std::vector<double> v;
+    v.push_back(0);
+    v.push_back(4);
+    v.push_back(3);
+    v.push_back(-1);
+    BOOST_CHECK_EQUAL(MathUtilities::mean(v, 0, 4), 1.5);
+    BOOST_CHECK_EQUAL(MathUtilities::mean(v, 1, 2), 3.5);
+    BOOST_CHECK_EQUAL(MathUtilities::mean(v, 3, 1), -1);
+    BOOST_CHECK_EQUAL(MathUtilities::mean(v, 3, 0), 0);
+}
+
+BOOST_AUTO_TEST_CASE(sum)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::sum(0, 0), 0);
+    double d0[] = { 0, 4, 3, -1 };
+    BOOST_CHECK_EQUAL(MathUtilities::sum(d0, 4), 6);
+    double d1[] = { -2.6 };
+    BOOST_CHECK_EQUAL(MathUtilities::sum(d1, 1), -2.6);
+}
+
+BOOST_AUTO_TEST_CASE(median)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::median(0, 0), 0);
+    double d0[] = { 0, 4, 3, -1 };
+    BOOST_CHECK_EQUAL(MathUtilities::median(d0, 4), 1.5);
+    double d1[] = { 0, 4, 3, -1, -1 };
+    BOOST_CHECK_EQUAL(MathUtilities::median(d1, 5), 0);
+    double d2[] = { 1.0, -2.0 };
+    BOOST_CHECK_EQUAL(MathUtilities::median(d2, 2), -0.5);
+    double d3[] = { -2.6 };
+    BOOST_CHECK_EQUAL(MathUtilities::median(d3, 1), -2.6);
+}
+
+BOOST_AUTO_TEST_CASE(princarg)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(M_PI), M_PI);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(-M_PI), M_PI);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(2 * M_PI), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(5 * M_PI), M_PI);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(1.0), 1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(-1.0), -1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::princarg(-10.0), -10.0 + 4 * M_PI);
+}
+
+BOOST_AUTO_TEST_CASE(isPowerOfTwo)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(0), false);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(1), true);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(-2), false);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(2), true);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(3), false);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(12), false);
+    BOOST_CHECK_EQUAL(MathUtilities::isPowerOfTwo(16), true);
+}
+
+BOOST_AUTO_TEST_CASE(nextPowerOfTwo)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(0), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(1), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(-2), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(2), 2);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(3), 4);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(12), 16);
+    BOOST_CHECK_EQUAL(MathUtilities::nextPowerOfTwo(16), 16);
+}
+
+BOOST_AUTO_TEST_CASE(previousPowerOfTwo)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(0), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(1), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(-2), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(2), 2);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(3), 2);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(12), 8);
+    BOOST_CHECK_EQUAL(MathUtilities::previousPowerOfTwo(16), 16);
+}
+
+BOOST_AUTO_TEST_CASE(nearestPowerOfTwo)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(0), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(1), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(-2), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(2), 2);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(3), 4);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(11), 8);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(12), 16);
+    BOOST_CHECK_EQUAL(MathUtilities::nearestPowerOfTwo(16), 16);
+}
+
+BOOST_AUTO_TEST_CASE(factorial)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(-10), 0.0);
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(0), 1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(1), 1.0);
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(2), 2.0);
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(3), 6.0);
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(4), 24.0);
+
+    // Too big for an int, hence double return value from factorial
+    BOOST_CHECK_EQUAL(MathUtilities::factorial(20), 2432902008176640000.0);
+}
+
+BOOST_AUTO_TEST_CASE(gcd)
+{
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(1, 1), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(2, 1), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(2, 3), 1);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(4, 2), 2);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(18, 24), 6);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(27, 18), 9);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(18, 36), 18);
+    BOOST_CHECK_EQUAL(MathUtilities::gcd(37, 18), 1);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
diff -r 1e33f719dde1 -r 6b13f9c694a8 test/TestResampler.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestResampler.cpp	Mon May 19 10:21:51 2014 +0100
@@ -0,0 +1,315 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/Resampler.h"
+
+#include "dsp/Window.h"
+#include "dsp/FFT.h"
+
+#include <iostream>
+
+#include <cmath>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestResampler)
+
+using std::cout;
+using std::endl;
+using std::vector;
+
+void
+testResamplerOneShot(int sourceRate,
+		     int targetRate,
+		     int n,
+		     double *in,
+		     int m,
+		     double *expected,
+		     int skip)
+{
+    vector<double> resampled = Resampler::resample(sourceRate, targetRate,
+						   in, n);
+    if (skip == 0) {
+	BOOST_CHECK_EQUAL(resampled.size(), m);
+    }
+    for (int i = 0; i < m; ++i) {
+	BOOST_CHECK_SMALL(resampled[i + skip] - expected[i], 1e-6);
+    }
+}
+
+void
+testResampler(int sourceRate,
+	      int targetRate,
+	      int n,
+	      double *in,
+	      int m,
+	      double *expected)
+{
+    // Here we provide the input in chunks (of varying size)
+
+    Resampler r(sourceRate, targetRate);
+    int latency = r.getLatency();
+
+    int m1 = m + latency;
+    int n1 = int((m1 * sourceRate) / targetRate);
+
+    double *inPadded = new double[n1];
+    double *outPadded = new double[m1];
+
+    for (int i = 0; i < n1; ++i) {
+	if (i < n) inPadded[i] = in[i];
+	else inPadded[i] = 0.0;
+    }
+    
+    for (int i = 0; i < m1; ++i) {
+	outPadded[i] = -999.0;
+    }
+
+    int chunkSize = 1;
+    int got = 0;
+    int i = 0;
+
+    while (true) {
+	got += r.process(inPadded + i, outPadded + got, chunkSize);
+	i = i + chunkSize;
+	chunkSize = chunkSize + 1;
+	if (i >= n1) {
+	    break;
+	} else if (i + chunkSize >= n1) {
+	    chunkSize = n1 - i;
+	} else if (chunkSize > 15) {
+	    chunkSize = 1;
+	}
+    }
+
+    BOOST_CHECK_EQUAL(got, m1);
+
+    for (int i = latency; i < m1; ++i) {
+	BOOST_CHECK_SMALL(outPadded[i] - expected[i-latency], 1e-8);
+    }
+
+    delete[] outPadded;
+    delete[] inPadded;
+}
+
+BOOST_AUTO_TEST_CASE(sameRateOneShot)
+{
+    double d[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
+    testResamplerOneShot(4, 4, 10, d, 10, d, 0);
+}
+
+BOOST_AUTO_TEST_CASE(sameRate)
+{
+    double d[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
+    testResampler(4, 4, 10, d, 10, d);
+}
+
+BOOST_AUTO_TEST_CASE(interpolatedMisc)
+{
+    // Interpolating any signal by N should give a signal in which
+    // every Nth sample is the original signal
+    double in[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
+    int n = sizeof(in)/sizeof(in[0]);
+    for (int factor = 2; factor < 10; ++factor) {
+	vector<double> out = Resampler::resample(6, 6 * factor, in, n);
+	for (int i = 0; i < n; ++i) {
+	    BOOST_CHECK_SMALL(out[i * factor] - in[i], 1e-5);
+	}
+    }
+}
+
+BOOST_AUTO_TEST_CASE(interpolatedSine)
+{
+    // Interpolating a sinusoid should give us a sinusoid, once we've
+    // dropped the first few samples
+    double in[1000];
+    double out[2000];
+    for (int i = 0; i < 1000; ++i) {
+	in[i] = sin(i * M_PI / 2.0);
+    }
+    for (int i = 0; i < 2000; ++i) {
+	out[i] = sin(i * M_PI / 4.0);
+    }
+    testResamplerOneShot(8, 16, 1000, in, 200, out, 512);
+}
+
+BOOST_AUTO_TEST_CASE(decimatedSine)
+{
+    // Decimating a sinusoid should give us a sinusoid, once we've
+    // dropped the first few samples
+    double in[2000];
+    double out[1000];
+    for (int i = 0; i < 2000; ++i) {
+	in[i] = sin(i * M_PI / 8.0);
+    }
+    for (int i = 0; i < 1000; ++i) {
+	out[i] = sin(i * M_PI / 4.0);
+    }
+    testResamplerOneShot(16, 8, 2000, in, 200, out, 256);
+}
+
+double
+measureSinFreq(const vector<double> &v, int rate, int countCycles)
+{
+    int n = v.size();
+    int firstPeak = -1;
+    int lastPeak = -1;
+    int nPeaks = 0;
+    // count +ve peaks
+    for (int i = v.size()/4; i + 1 < n; ++i) {
+        // allow some fuzz
+        int x0 = int(10000 * v[i-1]);
+        int x1 = int(10000 * v[i]);
+        int x2 = int(10000 * v[i+1]);
+        if (x1 > 0 && x1 > x0 && x1 >= x2) {
+            if (firstPeak < 0) firstPeak = i;
+            lastPeak = i;
+            ++nPeaks;
+            if (nPeaks == countCycles) break;
+        }
+    }
+    int nCycles = nPeaks - 1;
+    if (nCycles <= 0) return 0.0;
+    double cycle = double(lastPeak - firstPeak) / nCycles;
+//    cout << "lastPeak = " << lastPeak << ", firstPeak = " << firstPeak << ", dist = " << lastPeak - firstPeak << ", nCycles = " << nCycles << ", cycle = " << cycle << endl;
+    return rate / cycle;
+}
+
+void
+testSinFrequency(int freq,
+                 int sourceRate,
+                 int targetRate)
+{
+    // Resampling a sinusoid and then resampling back should give us a
+    // sinusoid of the same frequency as we started with
+
+    int nCycles = 500;
+
+    int duration = int(nCycles * float(sourceRate) / float(freq));
+//    cout << "freq = " << freq << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", duration = " << duration << endl;
+
+    vector<double> in(duration, 0);
+    for (int i = 0; i < duration; ++i) {
+        in[i] = sin(i * M_PI * 2.0 * freq / sourceRate);
+    }
+
+    vector<double> out = Resampler::resample(sourceRate, targetRate,
+                                             in.data(), in.size());
+
+    vector<double> back = Resampler::resample(targetRate, sourceRate,
+                                              out.data(), out.size());
+
+    BOOST_CHECK_EQUAL(in.size(), back.size());
+
+    double inFreq = measureSinFreq(in, sourceRate, nCycles / 2);
+    double backFreq = measureSinFreq(back, sourceRate, nCycles / 2);
+    
+    BOOST_CHECK_SMALL(inFreq - backFreq, 1e-8);
+}
+
+// In each of the following we use a frequency that has an exact cycle
+// length in samples at the lowest sample rate, so that we can easily
+// rule out errors in measuring the cycle length after resampling. If
+// the resampler gets its input or output rate wrong, that will show
+// up no matter what the test signal's initial frequency is.
+
+BOOST_AUTO_TEST_CASE(downUp2)
+{
+    testSinFrequency(441, 44100, 22050);
+}
+
+BOOST_AUTO_TEST_CASE(downUp5)
+{
+    testSinFrequency(300, 15000, 3000);
+}
+
+BOOST_AUTO_TEST_CASE(downUp16)
+{
+    testSinFrequency(300, 48000, 3000);
+}
+
+BOOST_AUTO_TEST_CASE(upDown2)
+{
+    testSinFrequency(441, 44100, 88200);
+}
+
+BOOST_AUTO_TEST_CASE(upDown5)
+{
+    testSinFrequency(300, 3000, 15000);
+}
+
+BOOST_AUTO_TEST_CASE(upDown16)
+{
+    testSinFrequency(300, 3000, 48000);
+}
+
+vector<double>
+squareWave(int rate, double freq, int n)
+{
+    //!!! todo: hoist, test
+    vector<double> v(n, 0.0);
+    for (int h = 0; h < (rate/4)/freq; ++h) {
+	double m = h * 2 + 1;
+	double scale = 1.0 / m;
+	for (int i = 0; i < n; ++i) {
+	    double s = scale * sin((i * 2.0 * M_PI * m * freq) / rate);
+	    v[i] += s;
+	}
+    }
+    return v;
+}
+
+void
+testSpectrum(int inrate, int outrate)
+{
+    // One second of a square wave
+    int freq = 500;
+
+    vector<double> square =
+	squareWave(inrate, freq, inrate);
+
+    vector<double> maybeSquare = 
+	Resampler::resample(inrate, outrate, square.data(), square.size());
+
+    BOOST_CHECK_EQUAL(maybeSquare.size(), outrate);
+
+    Window<double>(HanningWindow, inrate).cut(square.data());
+    Window<double>(HanningWindow, outrate).cut(maybeSquare.data());
+
+    // forward magnitude with size inrate, outrate
+
+    vector<double> inSpectrum(inrate, 0.0);
+    FFTReal(inrate).forwardMagnitude(square.data(), inSpectrum.data());
+    for (int i = 0; i < (int)inSpectrum.size(); ++i) {
+	inSpectrum[i] /= inrate;
+    }
+
+    vector<double> outSpectrum(outrate, 0.0);
+    FFTReal(outrate).forwardMagnitude(maybeSquare.data(), outSpectrum.data());
+    for (int i = 0; i < (int)outSpectrum.size(); ++i) {
+	outSpectrum[i] /= outrate;
+    }
+
+    // Don't compare bins any higher than 96% of Nyquist freq of lower sr
+    int lengthOfInterest = (inrate < outrate ? inrate : outrate) / 2;
+    lengthOfInterest = lengthOfInterest - (lengthOfInterest / 25);
+
+    for (int i = 0; i < lengthOfInterest; ++i) {
+	BOOST_CHECK_SMALL(inSpectrum[i] - outSpectrum[i], 1e-7);
+    }
+}
+/*
+BOOST_AUTO_TEST_CASE(spectrum)
+{
+    int rates[] = { 8000, 22050, 44100, 48000 };
+    for (int i = 0; i < (int)(sizeof(rates)/sizeof(rates[0])); ++i) {
+	    for (int j = 0; j < (int)(sizeof(rates)/sizeof(rates[0])); ++j) {
+	    testSpectrum(rates[i], rates[j]);
+	}
+    }
+}
+*/
+BOOST_AUTO_TEST_SUITE_END()
+
diff -r 1e33f719dde1 -r 6b13f9c694a8 test/TestWindow.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestWindow.cpp	Mon May 19 10:21:51 2014 +0100
@@ -0,0 +1,202 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/Window.h"
+#include "dsp/KaiserWindow.h"
+#include "dsp/SincWindow.h"
+
+#include <iostream>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestWindow)
+
+using std::cout;
+using std::endl;
+
+#define COMPARE_ARRAY(a, b)						\
+    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
+        BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-4);			\
+    }
+
+void
+testSymmetric(double *d, int n)
+{
+    for (int i = 0; i <= n/2; ++i) {
+	BOOST_CHECK_CLOSE(d[i], d[n-i-1], 1e-10);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(periodic)
+{
+    // We can't actually test whether a function is periodic, given
+    // only one cycle of it! But we can make sure that all but the
+    // first sample is symmetric, which is what a symmetric window
+    // becomes when generated in periodic mode
+    double d[9];
+    for (int n = 8; n <= 9; ++n) {
+	for (int wt = (int)FirstWindow; wt <= (int)LastWindow; ++wt) {
+	    for (int i = 0; i < n; ++i) d[i] = 1.0;
+	    Window<double> w((WindowType)wt, n);
+	    w.cut(d);
+	    testSymmetric(d + 1, n - 1);
+	}
+    }
+}
+
+template <int N>
+void testWindow(WindowType type, const double expected[N])
+{
+    double d[N];
+    for (int i = 0; i < N; ++i) d[i] = 1.0;
+    Window<double> w(type, N);
+    w.cut(d);
+    COMPARE_ARRAY(d, expected);
+
+    double d0[N], d1[N];
+    for (int i = 0; i < N; ++i) d0[i] = 0.5 + (1.0 / (N * 2)) * (i + 1);
+    w.cut(d0, d1);
+    for (int i = 0; i < N; ++i) {
+	BOOST_CHECK_SMALL(d1[i] - d0[i] * expected[i], 1e-4);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(bartlett)
+{
+    double e1[] = { 1 };
+    testWindow<1>(BartlettWindow, e1);
+
+    double e2[] = { 0, 0 };
+    testWindow<2>(BartlettWindow, e2);
+
+    double e3[] = { 0, 2./3., 2./3. };
+    testWindow<3>(BartlettWindow, e3);
+
+    double e4[] = { 0, 1./2., 1., 1./2. };
+    testWindow<4>(BartlettWindow, e4);
+
+    double e5[] = { 0, 1./2., 1., 1., 1./2. };
+    testWindow<5>(BartlettWindow, e5);
+
+    double e6[] = { 0, 1./3., 2./3., 1., 2./3., 1./3. };
+    testWindow<6>(BartlettWindow, e6);
+}
+    
+BOOST_AUTO_TEST_CASE(hamming)
+{
+    double e1[] = { 1 };
+    testWindow<1>(HammingWindow, e1);
+
+    double e10[] = {
+	0.0800, 0.1679, 0.3979, 0.6821, 0.9121,
+        1.0000, 0.9121, 0.6821, 0.3979, 0.1679
+    };
+    testWindow<10>(HammingWindow, e10);
+}
+    
+BOOST_AUTO_TEST_CASE(hann)
+{
+    double e1[] = { 1 };
+    testWindow<1>(HanningWindow, e1);
+
+    double e10[] = {
+        0, 0.0955, 0.3455, 0.6545, 0.9045,
+        1.0000, 0.9045, 0.6545, 0.3455, 0.0955,
+    };
+    testWindow<10>(HanningWindow, e10);
+}
+    
+BOOST_AUTO_TEST_CASE(blackman)
+{
+    double e1[] = { 1 };
+    testWindow<1>(BlackmanWindow, e1);
+
+    double e10[] = {
+        0, 0.0402, 0.2008, 0.5098, 0.8492,
+        1.0000, 0.8492, 0.5098, 0.2008, 0.0402,
+    };
+    testWindow<10>(BlackmanWindow, e10);
+}
+    
+BOOST_AUTO_TEST_CASE(blackmanHarris)
+{
+    double e1[] = { 1 };
+    testWindow<1>(BlackmanHarrisWindow, e1);
+
+    double e10[] = {
+        0.0001, 0.0110, 0.1030, 0.3859, 0.7938,
+        1.0000, 0.7938, 0.3859, 0.1030, 0.0110,
+    };
+    testWindow<10>(BlackmanHarrisWindow, e10);
+}
+
+BOOST_AUTO_TEST_CASE(kaiser_4_10)
+{
+    double d[10];
+    for (int i = 0; i < 10; ++i) d[i] = 1.0;
+    double e[] = { 
+        0.0885, 0.2943, 0.5644, 0.8216, 0.9789,
+        0.9789, 0.8216, 0.5644, 0.2943, 0.0885
+    };
+    KaiserWindow::Parameters p;
+    p.length = 10;
+    p.beta = 4;
+    KaiserWindow k(p);
+    k.cut(d);
+    COMPARE_ARRAY(d, e);
+}
+    
+BOOST_AUTO_TEST_CASE(kaiser_2p5_11)
+{
+    double d[11];
+    for (int i = 0; i < 11; ++i) d[i] = 1.0;
+    double e[] = { 
+        0.3040, 0.5005, 0.6929, 0.8546, 0.9622,
+        1.0000, 0.9622, 0.8546, 0.6929, 0.5005, 0.3040
+    };
+    KaiserWindow::Parameters p;
+    p.length = 11;
+    p.beta = 2.5;
+    KaiserWindow k(p);
+    k.cut(d);
+    COMPARE_ARRAY(d, e);
+}
+
+//!!! todo: tests for kaiser with attenuation and bandwidth parameters
+
+template <int N>
+void testSinc(double p, const double expected[N])
+{
+    double d[N];
+    for (int i = 0; i < N; ++i) d[i] = 1.0;
+    SincWindow w(N, p);
+    w.cut(d);
+    COMPARE_ARRAY(d, expected);
+
+    double d0[N], d1[N];
+    for (int i = 0; i < N; ++i) d0[i] = 0.5 + (1.0 / (N * 2)) * (i + 1);
+    w.cut(d0, d1);
+    for (int i = 0; i < N; ++i) {
+	BOOST_CHECK_SMALL(d1[i] - d0[i] * expected[i], 1e-4);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(sinc)
+{
+    double e1[] = { 0, 0, 1, 0, 0 };
+    testSinc<5>(1, e1);
+
+    double e2[] = { 0, 0, 1, 0, 0 };
+    testSinc<5>(2, e2);
+
+    double e3[] = { -0.2122, 0.0, 0.6366, 1.0, 0.6366, 0, -0.2122 };
+    testSinc<7>(4, e3);
+
+    double e4[] = { -0.2122, 0, 0.6366, 1, 0.6366, 0 };
+    testSinc<6>(4, e4);
+}
+        
+BOOST_AUTO_TEST_SUITE_END()
+