changeset 141:54c9e0811ae7

More resampler fixes (particularly to latency calculation) and tests

--- a/dsp/rateconversion/Resampler.cpp	Mon Oct 14 08:19:45 2013 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Mon Oct 14 16:15:32 2013 +0100
@@ -11,6 +11,8 @@
 
 using std::vector;
 
+//#define DEBUG_RESAMPLER 1
+
 Resampler::Resampler(int sourceRate, int targetRate) :
     m_sourceRate(sourceRate),
     m_targetRate(targetRate)
@@ -52,9 +54,12 @@
     int inputSpacing = m_targetRate / m_gcd;
     int outputSpacing = m_sourceRate / m_gcd;
 
-    m_latency = int(ceil((m_filterLength / 2.0) / outputSpacing));
-
-    int bufferLength = 0;
+#ifdef DEBUG_RESAMPLER
+    std::cerr << "resample " << m_sourceRate << " -> " << m_targetRate
+	      << ": inputSpacing " << inputSpacing << ", outputSpacing "
+	      << outputSpacing << ": filter length " << m_filterLength
+	      << std::endl;
+#endif
 
     m_phaseData = new Phase[inputSpacing];
 
@@ -66,16 +71,11 @@
 	while (p.nextPhase < 0) p.nextPhase += inputSpacing;
 	p.nextPhase %= inputSpacing;
 	
-	p.drop = int(ceil(std::max(0, outputSpacing - phase) / inputSpacing));
-	p.take = int((outputSpacing +
-		      ((m_filterLength - 1 - phase) % inputSpacing))
-		     / outputSpacing);
+	p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase))
+			  / inputSpacing));
 
-	int filtZipLength = int(ceil((m_filterLength - phase) / inputSpacing));
-	if (filtZipLength > bufferLength) {
-	    bufferLength = filtZipLength;
-	}
-	
+	int filtZipLength = int(ceil(double(m_filterLength - phase)
+				     / inputSpacing));
 	for (int i = 0; i < filtZipLength; ++i) {
 	    p.filter.push_back(filter[i * inputSpacing + phase]);
 	}
@@ -83,6 +83,24 @@
 	m_phaseData[phase] = p;
     }
 
+#ifdef DEBUG_RESAMPLER
+    for (int phase = 0; phase < inputSpacing; ++phase) {
+	std::cerr << "filter for phase " << phase << " of " << inputSpacing << " (with length " << m_phaseData[phase].filter.size() << "):";
+	for (int i = 0; i < m_phaseData[phase].filter.size(); ++i) {
+	    if (i % 4 == 0) {
+		std::cerr << std::endl << i << ": ";
+	    }
+	    float v = m_phaseData[phase].filter[i];
+	    if (v == 1) {
+		std::cerr << " *** " << v << " ***  ";
+	    } else {
+		std::cerr << v << " ";
+	    }
+	}
+	std::cerr << std::endl;
+    }
+#endif
+
     delete[] filter;
 
     // The May implementation of this uses a pull model -- we ask the
@@ -102,55 +120,58 @@
     // else have a lengthy declared latency on the output. We do the
     // latter. (What do other implementations do?)
 
-    m_phase = m_filterLength % inputSpacing;
-    m_buffer = vector<double>(bufferLength, 0);
+    m_phase = (m_filterLength/2) % inputSpacing;
+    
+    m_buffer = vector<double>(m_phaseData[0].filter.size(), 0);
+
+    m_latency =
+	((m_buffer.size() * inputSpacing) - (m_filterLength/2)) / outputSpacing
+	+ m_phase;
+
+#ifdef DEBUG_RESAMPLER
+    std::cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
+	      << ", latency " << m_latency << std::endl;
+#endif
 }
 
 double
-Resampler::reconstructOne(const double *src)
+Resampler::reconstructOne()
 {
     Phase &pd = m_phaseData[m_phase];
     double *filt = pd.filter.data();
+    double v = 0.0;
     int n = pd.filter.size();
-    double v = 0.0;
     for (int i = 0; i < n; ++i) {
 	v += m_buffer[i] * filt[i];
     }
     m_buffer = vector<double>(m_buffer.begin() + pd.drop, m_buffer.end());
-    for (int i = 0; i < pd.take; ++i) {
-	m_buffer.push_back(src[i]);
-    }
+    m_phase = pd.nextPhase;
     return v;
 }
 
 int
-Resampler::process(const double *src, double *dst, int remaining)
+Resampler::process(const double *src, double *dst, int n)
 {
-    int m = 0;
-    int offset = 0;
-
-    while (remaining >= m_phaseData[m_phase].take) {
-//	std::cerr << "remaining = " << remaining << ", m = " << m << ", take = " << m_phaseData[m_phase].take << std::endl;
-	int advance = m_phaseData[m_phase].take;
-	dst[m] = reconstructOne(src + offset);
-	offset += advance;
-	remaining -= advance;
-	m_phase = m_phaseData[m_phase].nextPhase;
-//	std::cerr << "remaining -> " << remaining << ", new phase has advance " << m_phaseData[m_phase].take << std::endl;
-	++m;
+    for (int i = 0; i < n; ++i) {
+	m_buffer.push_back(src[i]);
     }
 
-//    if (remaining > 0) {
-//	std::cerr << "have " << remaining << " spare, pushing to buffer" << std::endl;
-//    }	
+    int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate));
+    int outidx = 0;
 
-    for (int i = 0; i < remaining; ++i) {
-	m_buffer.push_back(src[offset + i]);
+#ifdef DEBUG_RESAMPLER
+    std::cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << std::endl;
+#endif
+
+    while (outidx < maxout &&
+	   m_buffer.size() >= m_phaseData[m_phase].filter.size()) {
+	dst[outidx] = reconstructOne();
+	outidx++;
     }
-
-    return m;
+    
+    return outidx;
 }
-
+    
 std::vector<double>
 Resampler::resample(int sourceRate, int targetRate, const double *data, int n)
 {
@@ -158,9 +179,9 @@
 
     int latency = r.getLatency();
 
-    int m = int(ceil((n * targetRate) / sourceRate));
+    int m = int(ceil(double(n * targetRate) / sourceRate));
     int m1 = m + latency;
-    int n1 = int((m1 * sourceRate) / targetRate);
+    int n1 = int(double(m1 * sourceRate) / targetRate);
 
     vector<double> pad(n1 - n, 0.0);
     vector<double> out(m1, 0.0);
@@ -168,6 +189,15 @@
     int got = r.process(data, out.data(), n);
     got += r.process(pad.data(), out.data() + got, pad.size());
 
+#ifdef DEBUG_RESAMPLER
+    std::cerr << "resample: " << n << " in, " << got << " out" << std::endl;
+    for (int i = 0; i < got; ++i) {
+	if (i % 5 == 0) std::cout << std::endl << i << "... ";
+	std::cout << (float) out[i] << " ";
+    }
+    std::cout << std::endl;
+#endif
+
     return vector<double>(out.begin() + latency, out.begin() + got);
 }
 

--- a/dsp/rateconversion/Resampler.h	Mon Oct 14 08:19:45 2013 +0100
+++ b/dsp/rateconversion/Resampler.h	Mon Oct 14 16:15:32 2013 +0100
@@ -50,7 +50,6 @@
         int nextPhase;
         std::vector<double> filter;
         int drop;
-        int take;
     };
 
     Phase *m_phaseData;
@@ -58,7 +57,7 @@
     std::vector<double> m_buffer;
 
     void initialise();
-    double reconstructOne(const double *);
+    double reconstructOne();
 };
 
 #endif

--- a/dsp/rateconversion/TestResampler.cpp	Mon Oct 14 08:19:45 2013 +0100
+++ b/dsp/rateconversion/TestResampler.cpp	Mon Oct 14 16:15:32 2013 +0100
@@ -23,13 +23,16 @@
 		     int n,
 		     double *in,
 		     int m,
-		     double *expected)
+		     double *expected,
+		     int skip)
 {
     vector<double> resampled = Resampler::resample(sourceRate, targetRate,
 						   in, n);
-    BOOST_CHECK_EQUAL(resampled.size(), m);
+    if (skip == 0) {
+	BOOST_CHECK_EQUAL(resampled.size(), m);
+    }
     for (int i = 0; i < m; ++i) {
-	BOOST_CHECK_SMALL(resampled[i] - expected[i], 1e-8);
+	BOOST_CHECK_SMALL(resampled[i + skip] - expected[i], 1e-8);
     }
 }
 
@@ -45,7 +48,6 @@
 
     Resampler r(sourceRate, targetRate);
     int latency = r.getLatency();
-    std::cerr << "latency = " << latency << std::endl;
 
     int m1 = m + latency;
     int n1 = int((m1 * sourceRate) / targetRate);
@@ -67,32 +69,24 @@
     int i = 0;
 
     while (true) {
-	std::cerr << "i = " << i << ", n1 = " << n1 << ", chunkSize = " << chunkSize << std::endl;
 	got += r.process(inPadded + i, outPadded + got, chunkSize);
 	i = i + chunkSize;
 	chunkSize = chunkSize + 1;
-	if (i + 1 >= n1) {
+	if (i >= n1) {
 	    break;
 	} else if (i + chunkSize >= n1) {
 	    chunkSize = n1 - i;
+	} else if (chunkSize > 15) {
+	    chunkSize = 1;
 	}
     }
 
-//    int got = r.process(inPadded, outPadded, n1);
-    std::cerr << i << " in, " << got << " out" << std::endl;
+    BOOST_CHECK_EQUAL(got, m1);
 
-    BOOST_CHECK_EQUAL(got, m1);
-/*
-    std::cerr << "results including latency padding:" << std::endl;
-    for (int i = 0; i < m1; ++i) {
-	std::cerr << outPadded[i] << " ";
-	if (i % 6 == 5) std::cerr << "\n";
-    }
-    std::cerr << "\n";
-*/
     for (int i = latency; i < m1; ++i) {
 	BOOST_CHECK_SMALL(outPadded[i] - expected[i-latency], 1e-8);
     }
+
     delete[] outPadded;
     delete[] inPadded;
 }
@@ -100,7 +94,7 @@
 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);
+    testResamplerOneShot(4, 4, 10, d, 10, d, 0);
 }
 
 BOOST_AUTO_TEST_CASE(sameRate)
@@ -109,5 +103,32 @@
     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)
+{
+    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, 400);
+}
+
 BOOST_AUTO_TEST_SUITE_END()