Mercurial > hg > qm-dsp

changeset 375:ad21307eaf99

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Integrate resampler and tests into build system etc
author Chris Cannam <c.cannam@qmul.ac.uk>
date Mon, 21 Oct 2013 09:40:22 +0100
parents 3e5f13ac984f
children 8e32cd2af755
files dsp/rateconversion/Decimator.h dsp/rateconversion/Resampler.cpp dsp/rateconversion/Resampler.h dsp/rateconversion/TestResampler.cpp qm-dsp.pro tests/Makefile tests/TestResampler.cpp
diffstat 7 files changed, 271 insertions(+), 225 deletions(-) [+]
line wrap: on
line diff
--- a/dsp/rateconversion/Decimator.h	Fri Oct 18 14:57:48 2013 +0100
+++ b/dsp/rateconversion/Decimator.h	Mon Oct 21 09:40:22 2013 +0100
@@ -15,6 +15,12 @@
 #ifndef DECIMATOR_H
 #define DECIMATOR_H
 
+/**
+ * Decimator carries out a fast downsample by a power-of-two
+ * factor. Only a limited number of factors are supported, from two to
+ * whatever getHighestSupportedFactor() returns. This is much faster
+ * than Resampler but has a worse signal-noise ratio.
+ */
 class Decimator  
 {
 public:
--- a/dsp/rateconversion/Resampler.cpp	Fri Oct 18 14:57:48 2013 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Mon Oct 21 09:40:22 2013 +0100
@@ -1,11 +1,23 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file by Chris Cannam.
+
+    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 "Resampler.h"
 
-#include "qm-dsp/maths/MathUtilities.h"
-#include "qm-dsp/base/KaiserWindow.h"
-#include "qm-dsp/base/SincWindow.h"
-#include "qm-dsp/thread/Thread.h"
+#include "maths/MathUtilities.h"
+#include "base/KaiserWindow.h"
+#include "base/SincWindow.h"
+#include "thread/Thread.h"
 
 #include <iostream>
 #include <vector>
--- a/dsp/rateconversion/Resampler.h	Fri Oct 18 14:57:48 2013 +0100
+++ b/dsp/rateconversion/Resampler.h	Mon Oct 21 09:40:22 2013 +0100
@@ -1,10 +1,32 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file by Chris Cannam.
+
+    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 RESAMPLER_H
 #define RESAMPLER_H
 
 #include <vector>
 
+/**
+ * Resampler resamples a stream from one integer sample rate to
+ * another (arbitrary) rate, using a kaiser-windowed sinc filter.  The
+ * results and performance are pretty similar to libraries such as
+ * libsamplerate, though this implementation does not support
+ * time-varying ratios (the ratio is fixed on construction).
+ *
+ * See also Decimator, which is faster and rougher but supports only
+ * power-of-two downsampling factors.
+ */
 class Resampler
 {
 public:
--- a/dsp/rateconversion/TestResampler.cpp	Fri Oct 18 14:57:48 2013 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,220 +0,0 @@
-/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
-#include "Resampler.h"
-
-#include "qm-dsp/base/Window.h"
-#include "qm-dsp/dsp/transforms/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);
-}
-
-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 99% of Nyquist freq of lower sr
-    int lengthOfInterest = (inrate < outrate ? inrate : outrate) / 2;
-    lengthOfInterest = lengthOfInterest - (lengthOfInterest / 100);
-
-    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()
-
--- a/qm-dsp.pro	Fri Oct 18 14:57:48 2013 +0100
+++ b/qm-dsp.pro	Mon Oct 21 09:40:22 2013 +0100
@@ -60,6 +60,7 @@
            dsp/onsets/PeakPicking.h \
            dsp/phasevocoder/PhaseVocoder.h \
            dsp/rateconversion/Decimator.h \
+           dsp/rateconversion/Resampler.h \
            dsp/rhythm/BeatSpectrum.h \
            dsp/segmentation/cluster_melt.h \
            dsp/segmentation/ClusterMeltSegmenter.h \
@@ -102,6 +103,7 @@
            dsp/onsets/PeakPicking.cpp \
            dsp/phasevocoder/PhaseVocoder.cpp \
            dsp/rateconversion/Decimator.cpp \
+           dsp/rateconversion/Resampler.cpp \
            dsp/rhythm/BeatSpectrum.cpp \
            dsp/segmentation/cluster_melt.c \
            dsp/segmentation/ClusterMeltSegmenter.cpp \
--- a/tests/Makefile	Fri Oct 18 14:57:48 2013 +0100
+++ b/tests/Makefile	Mon Oct 21 09:40:22 2013 +0100
@@ -5,7 +5,7 @@
 LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework
 LIBS := ../libqm-dsp.a
 
-TESTS	:= test-mathutilities test-window test-fft test-pvoc
+TESTS	:= test-mathutilities test-window test-fft test-pvoc test-resampler
 
 #VG	:= valgrind
 
@@ -24,10 +24,14 @@
 test-pvoc: TestPhaseVocoder.o $(LIBS)
 	$(CXX) -o $@ $^ $(LDFLAGS)
 
+test-resampler: TestResampler.o $(LIBS)
+	$(CXX) -o $@ $^ $(LDFLAGS)
+
 TestMathUtilities.o: $(LIBS)
 TestWindow.o: $(LIBS)
 TestFFT.o: $(LIBS)
 TestPhaseVocoder.o: $(LIBS)
+TestResampler.o: $(LIBS)
 
 clean: 
 	rm -f *.o $(TESTS)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/TestResampler.cpp	Mon Oct 21 09:40:22 2013 +0100
@@ -0,0 +1,220 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/rateconversion/Resampler.h"
+
+#include "base/Window.h"
+#include "dsp/transforms/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);
+}
+
+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 99% of Nyquist freq of lower sr
+    int lengthOfInterest = (inrate < outrate ? inrate : outrate) / 2;
+    lengthOfInterest = lengthOfInterest - (lengthOfInterest / 100);
+
+    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()
+