c@362: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@362: 
c@362: #include "Resampler.h"
c@362: 
c@368: #include "qm-dsp/base/Window.h"
c@368: #include "qm-dsp/dsp/transforms/FFT.h"
c@368: 
c@362: #include <iostream>
c@362: 
c@362: #include <cmath>
c@362: 
c@362: #define BOOST_TEST_DYN_LINK
c@362: #define BOOST_TEST_MAIN
c@362: 
c@362: #include <boost/test/unit_test.hpp>
c@362: 
c@362: BOOST_AUTO_TEST_SUITE(TestResampler)
c@362: 
c@362: using std::cout;
c@362: using std::endl;
c@363: using std::vector;
c@363: 
c@363: void
c@363: testResamplerOneShot(int sourceRate,
c@363: 		     int targetRate,
c@363: 		     int n,
c@363: 		     double *in,
c@363: 		     int m,
c@366: 		     double *expected,
c@366: 		     int skip)
c@363: {
c@363:     vector<double> resampled = Resampler::resample(sourceRate, targetRate,
c@363: 						   in, n);
c@366:     if (skip == 0) {
c@366: 	BOOST_CHECK_EQUAL(resampled.size(), m);
c@366:     }
c@363:     for (int i = 0; i < m; ++i) {
c@373: 	BOOST_CHECK_SMALL(resampled[i + skip] - expected[i], 1e-6);
c@363:     }
c@363: }
c@362: 
c@362: void
c@362: testResampler(int sourceRate,
c@362: 	      int targetRate,
c@362: 	      int n,
c@362: 	      double *in,
c@362: 	      int m,
c@362: 	      double *expected)
c@362: {
c@364:     // Here we provide the input in chunks (of varying size)
c@363: 
c@362:     Resampler r(sourceRate, targetRate);
c@362:     int latency = r.getLatency();
c@362: 
c@362:     int m1 = m + latency;
c@362:     int n1 = int((m1 * sourceRate) / targetRate);
c@362: 
c@362:     double *inPadded = new double[n1];
c@362:     double *outPadded = new double[m1];
c@362: 
c@362:     for (int i = 0; i < n1; ++i) {
c@362: 	if (i < n) inPadded[i] = in[i];
c@362: 	else inPadded[i] = 0.0;
c@362:     }
c@362:     
c@362:     for (int i = 0; i < m1; ++i) {
c@362: 	outPadded[i] = -999.0;
c@362:     }
c@362: 
c@364:     int chunkSize = 1;
c@364:     int got = 0;
c@364:     int i = 0;
c@362: 
c@364:     while (true) {
c@364: 	got += r.process(inPadded + i, outPadded + got, chunkSize);
c@364: 	i = i + chunkSize;
c@364: 	chunkSize = chunkSize + 1;
c@366: 	if (i >= n1) {
c@364: 	    break;
c@364: 	} else if (i + chunkSize >= n1) {
c@364: 	    chunkSize = n1 - i;
c@366: 	} else if (chunkSize > 15) {
c@366: 	    chunkSize = 1;
c@364: 	}
c@364:     }
c@364: 
c@366:     BOOST_CHECK_EQUAL(got, m1);
c@362: 
c@362:     for (int i = latency; i < m1; ++i) {
c@363: 	BOOST_CHECK_SMALL(outPadded[i] - expected[i-latency], 1e-8);
c@362:     }
c@366: 
c@362:     delete[] outPadded;
c@362:     delete[] inPadded;
c@362: }
c@373: 
c@365: BOOST_AUTO_TEST_CASE(sameRateOneShot)
c@365: {
c@365:     double d[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
c@366:     testResamplerOneShot(4, 4, 10, d, 10, d, 0);
c@365: }
c@365: 
c@362: BOOST_AUTO_TEST_CASE(sameRate)
c@362: {
c@362:     double d[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
c@362:     testResampler(4, 4, 10, d, 10, d);
c@362: }
c@362: 
c@366: BOOST_AUTO_TEST_CASE(interpolatedMisc)
c@366: {
c@366:     // Interpolating any signal by N should give a signal in which
c@366:     // every Nth sample is the original signal
c@366:     double in[] = { 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 };
c@366:     int n = sizeof(in)/sizeof(in[0]);
c@366:     for (int factor = 2; factor < 10; ++factor) {
c@366: 	vector<double> out = Resampler::resample(6, 6 * factor, in, n);
c@366: 	for (int i = 0; i < n; ++i) {
c@366: 	    BOOST_CHECK_SMALL(out[i * factor] - in[i], 1e-5);
c@366: 	}
c@366:     }
c@366: }
c@366: 
c@366: BOOST_AUTO_TEST_CASE(interpolatedSine)
c@366: {
c@367:     // Interpolating a sinusoid should give us a sinusoid, once we've
c@367:     // dropped the first few samples
c@366:     double in[1000];
c@366:     double out[2000];
c@366:     for (int i = 0; i < 1000; ++i) {
c@366: 	in[i] = sin(i * M_PI / 2.0);
c@366:     }
c@366:     for (int i = 0; i < 2000; ++i) {
c@366: 	out[i] = sin(i * M_PI / 4.0);
c@366:     }
c@367:     testResamplerOneShot(8, 16, 1000, in, 200, out, 512);
c@367: }
c@367: 
c@367: BOOST_AUTO_TEST_CASE(decimatedSine)
c@367: {
c@367:     // Decimating a sinusoid should give us a sinusoid, once we've
c@367:     // dropped the first few samples
c@367:     double in[2000];
c@367:     double out[1000];
c@367:     for (int i = 0; i < 2000; ++i) {
c@367: 	in[i] = sin(i * M_PI / 8.0);
c@367:     }
c@367:     for (int i = 0; i < 1000; ++i) {
c@367: 	out[i] = sin(i * M_PI / 4.0);
c@367:     }
c@367:     testResamplerOneShot(16, 8, 2000, in, 200, out, 256);
c@366: }
c@373: 
c@368: vector<double>
c@369: squareWave(int rate, double freq, int n)
c@368: {
c@368:     //!!! todo: hoist, test
c@368:     vector<double> v(n, 0.0);
c@368:     for (int h = 0; h < (rate/4)/freq; ++h) {
c@368: 	double m = h * 2 + 1;
c@369: 	double scale = 1.0 / m;
c@368: 	for (int i = 0; i < n; ++i) {
c@369: 	    double s = scale * sin((i * 2.0 * M_PI * m * freq) / rate);
c@369: 	    v[i] += s;
c@368: 	}
c@368:     }
c@368:     return v;
c@368: }
c@368: 
c@368: void
c@368: testSpectrum(int inrate, int outrate)
c@368: {
c@368:     // One second of a square wave
c@368:     int freq = 500;
c@368: 
c@368:     vector<double> square =
c@368: 	squareWave(inrate, freq, inrate);
c@368: 
c@368:     vector<double> maybeSquare = 
c@368: 	Resampler::resample(inrate, outrate, square.data(), square.size());
c@368: 
c@368:     BOOST_CHECK_EQUAL(maybeSquare.size(), outrate);
c@368: 
c@368:     Window<double>(HanningWindow, inrate).cut(square.data());
c@368:     Window<double>(HanningWindow, outrate).cut(maybeSquare.data());
c@368: 
c@368:     // forward magnitude with size inrate, outrate
c@368: 
c@368:     vector<double> inSpectrum(inrate, 0.0);
c@368:     FFTReal(inrate).forwardMagnitude(square.data(), inSpectrum.data());
c@373:     for (int i = 0; i < (int)inSpectrum.size(); ++i) {
c@369: 	inSpectrum[i] /= inrate;
c@369:     }
c@368: 
c@368:     vector<double> outSpectrum(outrate, 0.0);
c@368:     FFTReal(outrate).forwardMagnitude(maybeSquare.data(), outSpectrum.data());
c@373:     for (int i = 0; i < (int)outSpectrum.size(); ++i) {
c@369: 	outSpectrum[i] /= outrate;
c@369:     }
c@368: 
c@368:     // Don't compare bins any higher than 99% of Nyquist freq of lower sr
c@368:     int lengthOfInterest = (inrate < outrate ? inrate : outrate) / 2;
c@368:     lengthOfInterest = lengthOfInterest - (lengthOfInterest / 100);
c@368: 
c@368:     for (int i = 0; i < lengthOfInterest; ++i) {
c@368: 	BOOST_CHECK_SMALL(inSpectrum[i] - outSpectrum[i], 1e-7);
c@368:     }
c@368: }
c@368: 
c@368: BOOST_AUTO_TEST_CASE(spectrum)
c@368: {
c@368:     int rates[] = { 8000, 22050, 44100, 48000 };
c@373:     for (int i = 0; i < (int)(sizeof(rates)/sizeof(rates[0])); ++i) {
c@373: 	    for (int j = 0; j < (int)(sizeof(rates)/sizeof(rates[0])); ++j) {
c@368: 	    testSpectrum(rates[i], rates[j]);
c@368: 	}
c@368:     }
c@368: }
c@368: 
c@362: BOOST_AUTO_TEST_SUITE_END()
c@362: