changeset 143:a4aa37f7af28

Some fixes, and start on spectrum test

--- a/dsp/rateconversion/Resampler.cpp	Mon Oct 14 16:20:00 2013 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Tue Oct 15 18:27:19 2013 +0100
@@ -43,8 +43,10 @@
     
     m_filterLength = params.length;
     
+    std::cerr << "making filter... ";
     KaiserWindow kw(params);
     SincWindow sw(m_filterLength, peakToPole * 2);
+    std::cerr << "done" << std::endl;
 
     double *filter = new double[m_filterLength];
     for (int i = 0; i < m_filterLength; ++i) filter[i] = 1.0;
@@ -138,11 +140,12 @@
 Resampler::reconstructOne()
 {
     Phase &pd = m_phaseData[m_phase];
-    double *filt = pd.filter.data();
     double v = 0.0;
     int n = pd.filter.size();
+    const double *buf = m_buffer.data();
+    const double *filt = pd.filter.data();
     for (int i = 0; i < n; ++i) {
-	v += m_buffer[i] * filt[i];
+	v += buf[i] * filt[i]; //!!! gcc can't vectorize: why?
     }
     m_buffer = vector<double>(m_buffer.begin() + pd.drop, m_buffer.end());
     m_phase = pd.nextPhase;
@@ -168,6 +171,8 @@
 	scaleFactor = double(m_targetRate) / double(m_sourceRate);
     }
 
+    std::cerr << "maxout = " << maxout << std::endl;
+
     while (outidx < maxout &&
 	   m_buffer.size() >= m_phaseData[m_phase].filter.size()) {
 	dst[outidx] = scaleFactor * reconstructOne();
@@ -184,12 +189,30 @@
 
     int latency = r.getLatency();
 
+    // latency is the output latency. We need to provide enough
+    // padding input samples at the end of input to guarantee at
+    // *least* the latency's worth of output samples. that is,
+
+    int inputPad = int(ceil(double(latency * sourceRate) / targetRate));
+
+    std::cerr << "latency = " << latency << ", inputPad = " << inputPad << std::endl;
+
+    // that means we are providing this much input in total:
+
+    int n1 = n + inputPad;
+
+    // and obtaining this much output in total:
+
+    int m1 = int(ceil(double(n1 * targetRate) / sourceRate));
+
+    // in order to return this much output to the user:
+
     int m = int(ceil(double(n * targetRate) / sourceRate));
-    int m1 = m + latency;
-    int n1 = int(double(m1 * sourceRate) / targetRate);
+    
+    std::cerr << "n = " << n << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", m = " << m << ", latency = " << latency << ", m1 = " << m1 << ", n1 = " << n1 << ", n1 - n = " << n1 - n << std::endl;
 
     vector<double> pad(n1 - n, 0.0);
-    vector<double> out(m1, 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());
@@ -203,6 +226,10 @@
     std::cout << std::endl;
 #endif
 
-    return vector<double>(out.begin() + latency, out.begin() + got);
+    int toReturn = got - latency;
+    if (toReturn > m) toReturn = m;
+
+    return vector<double>(out.begin() + latency, 
+			  out.begin() + latency + toReturn);
 }
 

--- a/dsp/rateconversion/TestResampler.cpp	Mon Oct 14 16:20:00 2013 +0100
+++ b/dsp/rateconversion/TestResampler.cpp	Tue Oct 15 18:27:19 2013 +0100
@@ -2,6 +2,9 @@
 
 #include "Resampler.h"
 
+#include "qm-dsp/base/Window.h"
+#include "qm-dsp/dsp/transforms/FFT.h"
+
 #include <iostream>
 
 #include <cmath>
@@ -147,5 +150,68 @@
     testResamplerOneShot(16, 8, 2000, in, 200, out, 256);
 }
 
+vector<double>
+squareWave(int rate, int 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 / m;
+	for (int i = 0; i < n; ++i) {
+	    v[i] += scale * sin(i * 2 * M_PI * freq / rate);
+	}
+    }
+    return v;
+}
+
+void
+testSpectrum(int inrate, int outrate)
+{
+    // One second of a square wave
+    int freq = 500;
+
+    std::cerr << "inrate = " << inrate << ", outrate = " << outrate << ", freq * outrate / inrate = " << (double(freq) * double(outrate)) / double(inrate) << std::endl;
+
+    std::cerr << "making square wave... ";
+    vector<double> square =
+	squareWave(inrate, freq, inrate);
+    std::cerr << "done" << std::endl;
+
+    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());
+
+    vector<double> outSpectrum(outrate, 0.0);
+    FFTReal(outrate).forwardMagnitude(maybeSquare.data(), outSpectrum.data());
+
+    // 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 < sizeof(rates)/sizeof(rates[0]); ++i) {
+	for (int j = 0; j < sizeof(rates)/sizeof(rates[0]); ++j) {
+	    testSpectrum(rates[i], rates[j]);
+	}
+    }
+}
+
 BOOST_AUTO_TEST_SUITE_END()