# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1383307628 0
# Node ID 88971211795c00cdcf393e319f7b96f48453a0a7
# Parent  1494e18439aedb6a5d0a581ff37500578869c13e
Lower filter cutoff to below target Nyquist when downsampling

diff -r 1494e18439ae -r 88971211795c dsp/rateconversion/Resampler.cpp
--- a/dsp/rateconversion/Resampler.cpp	Tue Oct 22 08:58:28 2013 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Fri Nov 01 12:07:08 2013 +0000
@@ -50,7 +50,7 @@
 }
 
 // peakToPole -> length -> beta -> window
-static map<int, map<int, map<double, vector<double> > > >
+static map<double, map<int, map<double, vector<double> > > >
 knownFilters;
 
 static Mutex
@@ -63,11 +63,15 @@
     int lower = std::min(m_sourceRate, m_targetRate);
 
     m_gcd = MathUtilities::gcd(lower, higher);
+    m_peakToPole = higher / m_gcd;
 
-    int peakToPole = higher / m_gcd;
+    if (m_targetRate < m_sourceRate) {
+        // antialiasing filter, should be slightly below nyquist
+        m_peakToPole = m_peakToPole / (1.0 - bandwidth/2.0);
+    }
 
     KaiserWindow::Parameters params =
-	KaiserWindow::parametersForBandwidth(snr, bandwidth, peakToPole);
+	KaiserWindow::parametersForBandwidth(snr, bandwidth, higher / m_gcd);
 
     params.length =
 	(params.length % 2 == 0 ? params.length + 1 : params.length);
@@ -80,21 +84,21 @@
     vector<double> filter;
     knownFilterMutex.lock();
 
-    if (knownFilters[peakToPole][m_filterLength].find(params.beta) ==
-	knownFilters[peakToPole][m_filterLength].end()) {
+    if (knownFilters[m_peakToPole][m_filterLength].find(params.beta) ==
+	knownFilters[m_peakToPole][m_filterLength].end()) {
 
 	KaiserWindow kw(params);
-	SincWindow sw(m_filterLength, peakToPole * 2);
+	SincWindow sw(m_filterLength, m_peakToPole * 2);
 
 	filter = vector<double>(m_filterLength, 0.0);
 	for (int i = 0; i < m_filterLength; ++i) filter[i] = 1.0;
 	sw.cut(filter.data());
 	kw.cut(filter.data());
 
-	knownFilters[peakToPole][m_filterLength][params.beta] = filter;
+	knownFilters[m_peakToPole][m_filterLength][params.beta] = filter;
     }
 
-    filter = knownFilters[peakToPole][m_filterLength][params.beta];
+    filter = knownFilters[m_peakToPole][m_filterLength][params.beta];
     knownFilterMutex.unlock();
 
     int inputSpacing = m_targetRate / m_gcd;
@@ -309,10 +313,7 @@
     std::cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << std::endl;
 #endif
 
-    double scaleFactor = 1.0;
-    if (m_targetRate < m_sourceRate) {
-	scaleFactor = double(m_targetRate) / double(m_sourceRate);
-    }
+    double scaleFactor = (double(m_targetRate) / m_gcd) / m_peakToPole;
 
     while (outidx < maxout &&
 	   m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) {
diff -r 1494e18439ae -r 88971211795c dsp/rateconversion/Resampler.h
--- a/dsp/rateconversion/Resampler.h	Tue Oct 22 08:58:28 2013 +0100
+++ b/dsp/rateconversion/Resampler.h	Fri Nov 01 12:07:08 2013 +0000
@@ -75,6 +75,7 @@
     int m_filterLength;
     int m_bufferLength;
     int m_latency;
+    double m_peakToPole;
     
     struct Phase {
         int nextPhase;
diff -r 1494e18439ae -r 88971211795c tests/TestResampler.cpp
--- a/tests/TestResampler.cpp	Tue Oct 22 08:58:28 2013 +0100
+++ b/tests/TestResampler.cpp	Fri Nov 01 12:07:08 2013 +0000
@@ -197,9 +197,9 @@
 	outSpectrum[i] /= outrate;
     }
 
-    // Don't compare bins any higher than 99% of Nyquist freq of lower sr
+    // Don't compare bins any higher than 96% of Nyquist freq of lower sr
     int lengthOfInterest = (inrate < outrate ? inrate : outrate) / 2;
-    lengthOfInterest = lengthOfInterest - (lengthOfInterest / 100);
+    lengthOfInterest = lengthOfInterest - (lengthOfInterest / 25);
 
     for (int i = 0; i < lengthOfInterest; ++i) {
 	BOOST_CHECK_SMALL(inSpectrum[i] - outSpectrum[i], 1e-7);