Mercurial > hg > qm-dsp

changeset 398:333e27d1efa1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixes to resampler frequency tests
author Chris Cannam <c.cannam@qmul.ac.uk>
date Mon, 12 May 2014 17:56:08 +0100
parents fd207df9432e
children e48b3f641038
files dsp/rateconversion/Resampler.cpp tests/TestResampler.cpp
diffstat 2 files changed, 91 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- a/dsp/rateconversion/Resampler.cpp	Sat May 10 13:41:06 2014 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Mon May 12 17:56:08 2014 +0100
@@ -26,9 +26,11 @@
 
 using std::vector;
 using std::map;
+using std::cerr;
+using std::endl;
 
 //#define DEBUG_RESAMPLER 1
-//#define DEBUG_RESAMPLER_VERBOSE
+//#define DEBUG_RESAMPLER_VERBOSE 1
 
 Resampler::Resampler(int sourceRate, int targetRate) :
     m_sourceRate(sourceRate),
@@ -106,10 +108,10 @@
     int outputSpacing = m_sourceRate / m_gcd;
 
 #ifdef DEBUG_RESAMPLER
-    std::cerr << "resample " << m_sourceRate << " -> " << m_targetRate
+    cerr << "resample " << m_sourceRate << " -> " << m_targetRate
 	      << ": inputSpacing " << inputSpacing << ", outputSpacing "
 	      << outputSpacing << ": filter length " << m_filterLength
-	      << std::endl;
+	      << endl;
 #endif
 
     // Now we have a filter of (odd) length flen in which the lower
@@ -205,6 +207,19 @@
 	m_phaseData[phase] = p;
     }
 
+#ifdef DEBUG_RESAMPLER
+    int cp = 0;
+    int totDrop = 0;
+    for (int i = 0; i < inputSpacing; ++i) {
+        cerr << "phase = " << cp << ", drop = " << m_phaseData[cp].drop
+             << ", filter length = " << m_phaseData[cp].filter.size()
+             << ", next phase = " << m_phaseData[cp].nextPhase << endl;
+        totDrop += m_phaseData[cp].drop;
+        cp = m_phaseData[cp].nextPhase;
+    }
+    cerr << "total drop = " << totDrop << endl;
+#endif
+
     // The May implementation of this uses a pull model -- we ask the
     // resampler for a certain number of output samples, and it asks
     // its source stream for as many as it needs to calculate
@@ -266,8 +281,8 @@
     m_bufferOrigin = 0;
 
 #ifdef DEBUG_RESAMPLER
-    std::cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
-	      << ", latency " << m_latency << std::endl;
+    cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
+	      << ", latency " << m_latency << endl;
 #endif
 }
 
@@ -283,12 +298,12 @@
     const double *const __restrict__ buf = m_buffer.data() + m_bufferOrigin;
     const double *const __restrict__ filt = pd.filter.data();
 
-//    std::cerr << "phase = " << m_phase << ", drop = " << pd.drop << ", buffer for reconstruction starts...";
+//    cerr << "phase = " << m_phase << ", drop = " << pd.drop << ", buffer for reconstruction starts...";
 //    for (int i = 0; i < 20; ++i) {
-//        if (i % 5 == 0) std::cerr << "\n" << i << " ";
-//        std::cerr << buf[i] << " ";
+//        if (i % 5 == 0) cerr << "\n" << i << " ";
+//        cerr << buf[i] << " ";
 //    }
-//    std::cerr << std::endl;
+//    cerr << endl;
 
     for (int i = 0; i < n; ++i) {
 	// NB gcc can only vectorize this with -ffast-math
@@ -311,7 +326,7 @@
     int outidx = 0;
 
 #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;
+    cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << endl;
 #endif
 
     double scaleFactor = (double(m_targetRate) / m_gcd) / m_peakToPole;
@@ -328,18 +343,18 @@
     return outidx;
 }
     
-std::vector<double>
+vector<double>
 Resampler::process(const double *src, int n)
 {
     int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate));
-    std::vector<double> out(maxout, 0.0);
+    vector<double> out(maxout, 0.0);
     int got = process(src, out.data(), n);
     assert(got <= maxout);
     if (got < maxout) out.resize(got);
     return out;
 }
 
-std::vector<double>
+vector<double>
 Resampler::resample(int sourceRate, int targetRate, const double *data, int n)
 {
     Resampler r(sourceRate, targetRate);
@@ -364,24 +379,28 @@
 
     int m = int(ceil((double(n) * targetRate) / sourceRate));
     
-//    std::cerr << "n = " << n << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", m = " << m << ", latency = " << latency << ", inputPad = " << inputPad << ", m1 = " << m1 << ", n1 = " << n1 << ", n1 - n = " << n1 - n << std::endl;
+#ifdef DEBUG_RESAMPLER
+    cerr << "n = " << n << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", m = " << m << ", latency = " << latency << ", inputPad = " << inputPad << ", m1 = " << m1 << ", n1 = " << n1 << ", n1 - n = " << n1 - n << endl;
+#endif
 
     vector<double> pad(n1 - n, 0.0);
     vector<double> out(m1 + 1, 0.0);
 
-    int got = r.process(data, out.data(), n);
-    got += r.process(pad.data(), out.data() + got, pad.size());
-
+    int gotData = r.process(data, out.data(), n);
+    int gotPad = r.process(pad.data(), out.data() + gotData, pad.size());
+    int got = gotData + gotPad;
+    
 #ifdef DEBUG_RESAMPLER
-    std::cerr << "resample: " << n << " in, " << got << " out" << std::endl;
+    cerr << "resample: " << n << " in, " << pad.size() << " padding, " << got << " out (" << gotData << " data, " << gotPad << " padding, latency = " << latency << ")" << endl;
 #endif
 #ifdef DEBUG_RESAMPLER_VERBOSE
-    std::cerr << "first 10 in:" << std::endl;
-    for (int i = 0; i < 10; ++i) {
-        std::cerr << data[i] << " ";
-        if (i == 5) std::cerr << std::endl;
+    int printN = 50;
+    cerr << "first " << printN << " in:" << endl;
+    for (int i = 0; i < printN && i < n; ++i) {
+	if (i % 5 == 0) cerr << endl << i << "... ";
+        cerr << data[i] << " ";
     }
-    std::cerr << std::endl;
+    cerr << endl;
 #endif
 
     int toReturn = got - latency;
@@ -391,12 +410,12 @@
 			  out.begin() + latency + toReturn);
 
 #ifdef DEBUG_RESAMPLER_VERBOSE
-    std::cerr << "all out (after latency compensation), length " << sliced.size() << ":";
-    for (int i = 0; i < sliced.size(); ++i) {
-	if (i % 5 == 0) std::cerr << std::endl << i << "... ";
-	std::cerr << sliced[i] << " ";
+    cerr << "first " << printN << " out (after latency compensation), length " << sliced.size() << ":";
+    for (int i = 0; i < printN && i < sliced.size(); ++i) {
+	if (i % 5 == 0) cerr << endl << i << "... ";
+	cerr << sliced[i] << " ";
     }
-    std::cerr << std::endl;
+    cerr << endl;
 #endif
 
     return sliced;
--- a/tests/TestResampler.cpp	Sat May 10 13:41:06 2014 +0100
+++ b/tests/TestResampler.cpp	Mon May 12 17:56:08 2014 +0100
@@ -154,22 +154,26 @@
 measureSinFreq(const vector<double> &v, int rate, int countCycles)
 {
     int n = v.size();
-    int firstCrossing = -1;
-    int lastCrossing = -1;
-    int nCrossings = 0;
-    // count -ve -> +ve transitions
-    for (int i = 0; i + 1 < n; ++i) {
-        if (v[i] <= 0.0 && v[i+1] > 0.0) {
-            if (firstCrossing < 0) firstCrossing = i;
-            lastCrossing = i;
-            ++nCrossings;
-            if (nCrossings == countCycles) break;
+    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 = nCrossings - 1;
+    int nCycles = nPeaks - 1;
     if (nCycles <= 0) return 0.0;
-    cout << "lastCrossing = " << lastCrossing << ", firstCrossing = " << firstCrossing << ", dist = " << lastCrossing - firstCrossing << ", nCycles = " << nCycles << endl;
-    double cycle = double(lastCrossing - firstCrossing) / nCycles;
+    double cycle = double(lastPeak - firstPeak) / nCycles;
+//    cout << "lastPeak = " << lastPeak << ", firstPeak = " << firstPeak << ", dist = " << lastPeak - firstPeak << ", nCycles = " << nCycles << ", cycle = " << cycle << endl;
     return rate / cycle;
 }
 
@@ -182,10 +186,10 @@
     // sinusoid of the same frequency as we started with. Let's start
     // with a few thousand cycles of it
 
-    int nCycles = 5000;
+    int nCycles = 500;
 
     int duration = int(nCycles * float(sourceRate) / float(freq));
-    cout << "freq = " << freq << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", duration = " << duration << endl;
+//    cout << "freq = " << freq << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", duration = " << duration << endl;
 
     vector<double> in(duration, 0);
     for (int i = 0; i < duration; ++i) {
@@ -200,36 +204,46 @@
 
     BOOST_CHECK_EQUAL(in.size(), back.size());
 
-    double inFreq = measureSinFreq(in, sourceRate, nCycles - 2);
-    double backFreq = measureSinFreq(back, sourceRate, nCycles - 2);
+    double inFreq = measureSinFreq(in, sourceRate, nCycles / 2);
+    double backFreq = measureSinFreq(back, sourceRate, nCycles / 2);
     
-    cout << "inFreq = " << inFreq << ", backFreq = " << backFreq << endl;
+    BOOST_CHECK_SMALL(inFreq - backFreq, 1e-8);
+}
 
-    BOOST_CHECK_SMALL(inFreq - backFreq, 1e-8);
-
-    //    for (int i = 0; i < int(in.size()); ++i) {
-//	BOOST_CHECK_SMALL(in[i] - back[i], 1e-6);
-//    }
-}
+// 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(440, 44100, 22050);
+    testSinFrequency(441, 44100, 22050);
+}
+
+BOOST_AUTO_TEST_CASE(downUp5)
+{
+    testSinFrequency(300, 15000, 3000);
 }
 
 BOOST_AUTO_TEST_CASE(downUp16)
 {
-    testSinFrequency(440, 48000, 3000);
+    testSinFrequency(300, 48000, 3000);
 }
 
 BOOST_AUTO_TEST_CASE(upDown2)
 {
-    testSinFrequency(440, 44100, 88200);
+    testSinFrequency(441, 44100, 88200);
+}
+
+BOOST_AUTO_TEST_CASE(upDown5)
+{
+    testSinFrequency(300, 3000, 15000);
 }
 
 BOOST_AUTO_TEST_CASE(upDown16)
 {
-    testSinFrequency(440, 3000, 48000);
+    testSinFrequency(300, 3000, 48000);
 }
 
 vector<double>
@@ -287,7 +301,7 @@
 	BOOST_CHECK_SMALL(inSpectrum[i] - outSpectrum[i], 1e-7);
     }
 }
-
+/*
 BOOST_AUTO_TEST_CASE(spectrum)
 {
     int rates[] = { 8000, 22050, 44100, 48000 };
@@ -297,6 +311,6 @@
 	}
     }
 }
-
+*/
 BOOST_AUTO_TEST_SUITE_END()