Mercurial > hg > constant-q-cpp

changeset 119:a38d6940f8fb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Bring in dsp dependencies
author Chris Cannam <c.cannam@qmul.ac.uk>
date Thu, 15 May 2014 12:24:11 +0100
parents 36bfbc606642
children 3fa7e938e7d4
files .hgsub .hgsubstate COPYING src/dsp/KaiserWindow.cpp src/dsp/KaiserWindow.h src/dsp/MathUtilities.cpp src/dsp/MathUtilities.h src/dsp/Resampler.cpp src/dsp/Resampler.h src/dsp/SincWindow.cpp src/dsp/SincWindow.h src/processfile.cpp src/test.cpp test/processfile.cpp test/test.cpp
diffstat 15 files changed, 1879 insertions(+), 377 deletions(-) [+]
line wrap: on
line diff
--- a/.hgsub	Thu May 15 12:06:31 2014 +0100
+++ b/.hgsub	Thu May 15 12:24:11 2014 +0100
@@ -1,1 +1,2 @@
+src/ext/kissfft = http://hg.code.sf.net/p/kissfft/code
 misc/matlab = https://code.soundsoftware.ac.uk/hg/constant-q-toolbox
--- a/.hgsubstate	Thu May 15 12:06:31 2014 +0100
+++ b/.hgsubstate	Thu May 15 12:24:11 2014 +0100
@@ -1,2 +1,3 @@
 2b03ca77abcc91140dbb37827f06034769c53547 matlab
 2b03ca77abcc91140dbb37827f06034769c53547 misc/matlab
+fbe1bb0bc7b94ec252842b8b7e3f3347ec75d92f src/ext/kissfft
--- a/COPYING	Thu May 15 12:06:31 2014 +0100
+++ b/COPYING	Thu May 15 12:24:11 2014 +0100
@@ -1,12 +1,6 @@
 
-The code in this library is provided under an MIT-style licence as
-follows. However, note that other code which it depends on may be
-licensed differently (for example, qm-dsp is under the GPL). This
-licence applies only to the source files actually contained in this
-repository, and not to any of its subrepositories. See individual
-files for more details.
+This library is provided to you under the following licence:
 
-    Constant-Q library
     Copyright (c) 2013-2014 Queen Mary, University of London
 
     Permission is hereby granted, free of charge, to any person
@@ -34,4 +28,38 @@
     use or other dealings in this Software without prior written
     authorization.
 
+The code in src/ext/kissfft is under the following licence:
 
+    Copyright (c) 2003-2010 Mark Borgerding
+
+    All rights reserved.
+
+    Redistribution and use in source and binary forms, with or without
+    modification, are permitted provided that the following conditions
+    are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the author nor the names of any contributors may be used
+      to endorse or promote products derived from this software
+      without specific prior written permission.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
+    CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+    INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+    MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+    DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+    TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+    ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
+    TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
+    THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+    SUCH DAMAGE.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/KaiserWindow.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,83 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#include "KaiserWindow.h"
+
+#include "maths/MathUtilities.h"
+
+KaiserWindow::Parameters
+KaiserWindow::parametersForTransitionWidth(double attenuation,
+					   double transition)
+{
+    Parameters p;
+    p.length = 1 + (attenuation > 21.0 ?
+		    ceil((attenuation - 7.95) / (2.285 * transition)) :
+		    ceil(5.79 / transition));
+    p.beta = (attenuation > 50.0 ? 
+	      0.1102 * (attenuation - 8.7) :
+	      attenuation > 21.0 ? 
+	      0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) :
+	      0);
+    return p;
+}
+
+static double besselTerm(double x, int i)
+{
+    if (i == 0) {
+	return 1;
+    } else {
+	double f = MathUtilities::factorial(i);
+	return pow(x/2, i*2) / (f*f);
+    }
+}
+
+static double bessel0(double x)
+{
+    double b = 0.0;
+    for (int i = 0; i < 20; ++i) {
+	b += besselTerm(x, i);
+    }
+    return b;
+}
+
+void
+KaiserWindow::init()
+{
+    double denominator = bessel0(m_beta);
+    bool even = (m_length % 2 == 0);
+    for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) {
+	double k = double(2*i) / double(m_length-1) - 1.0;
+	m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator);
+    }
+    for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) {
+        m_window.push_back(m_window[int(m_length/2) - i - 1]);
+    }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/KaiserWindow.h	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,123 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#ifndef KAISER_WINDOW_H
+#define KAISER_WINDOW_H
+
+#include <vector>
+#include <cmath>
+
+/**
+ * Kaiser window: A windower whose bandwidth and sidelobe height
+ * (signal-noise ratio) can be specified. These parameters are traded
+ * off against the window length.
+ */
+class KaiserWindow
+{
+public:
+    struct Parameters {
+	int length;
+	double beta;
+    };
+
+    /**
+     * Construct a Kaiser windower with the given length and beta
+     * parameter.
+     */
+    KaiserWindow(Parameters p) : m_length(p.length), m_beta(p.beta) { init(); }
+
+    /**
+     * Construct a Kaiser windower with the given attenuation in dB
+     * and transition width in samples.
+     */
+    static KaiserWindow byTransitionWidth(double attenuation,
+					  double transition) {
+	return KaiserWindow
+	    (parametersForTransitionWidth(attenuation, transition));
+    }
+
+    /**
+     * Construct a Kaiser windower with the given attenuation in dB
+     * and transition bandwidth in Hz for the given samplerate.
+     */
+    static KaiserWindow byBandwidth(double attenuation,
+				    double bandwidth,
+				    double samplerate) {
+	return KaiserWindow
+	    (parametersForBandwidth(attenuation, bandwidth, samplerate));
+    }
+
+    /**
+     * Obtain the parameters necessary for a Kaiser window of the
+     * given attenuation in dB and transition width in samples.
+     */
+    static Parameters parametersForTransitionWidth(double attenuation,
+						   double transition);
+
+    /**
+     * Obtain the parameters necessary for a Kaiser window of the
+     * given attenuation in dB and transition bandwidth in Hz for the
+     * given samplerate.
+     */
+    static Parameters parametersForBandwidth(double attenuation,
+					     double bandwidth,
+					     double samplerate) {
+	return parametersForTransitionWidth
+	    (attenuation, (bandwidth * 2 * M_PI) / samplerate);
+    } 
+
+    int getLength() const {
+	return m_length;
+    }
+
+    const double *getWindow() const { 
+	return m_window.data();
+    }
+
+    void cut(double *src) const { 
+	cut(src, src); 
+    }
+
+    void cut(const double *src, double *dst) const {
+	for (int i = 0; i < m_length; ++i) {
+	    dst[i] = src[i] * m_window[i];
+	}
+    }
+
+private:
+    int m_length;
+    double m_beta;
+    std::vector<double> m_window;
+
+    void init();
+};
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/MathUtilities.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,419 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#include "MathUtilities.h"
+
+#include <iostream>
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+
+double MathUtilities::mod(double x, double y)
+{
+    double a = floor( x / y );
+
+    double b = x - ( y * a );
+    return b;
+}
+
+double MathUtilities::princarg(double ang)
+{
+    double ValOut;
+
+    ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
+
+    return ValOut;
+}
+
+void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
+{
+    unsigned int i;
+    double temp = 0.0;
+    double a=0.0;
+	
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+		
+	a  += ::pow( fabs(temp), double(alpha) );
+    }
+    a /= ( double )len;
+    a = ::pow( a, ( 1.0 / (double) alpha ) );
+
+    *ANorm = a;
+}
+
+double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
+{
+    unsigned int i;
+    unsigned int len = data.size();
+    double temp = 0.0;
+    double a=0.0;
+	
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+		
+	a  += ::pow( fabs(temp), double(alpha) );
+    }
+    a /= ( double )len;
+    a = ::pow( a, ( 1.0 / (double) alpha ) );
+
+    return a;
+}
+
+double MathUtilities::round(double x)
+{
+    if (x < 0) {
+        return -floor(-x + 0.5);
+    } else {
+        return floor(x + 0.5);
+    }
+}
+
+double MathUtilities::median(const double *src, unsigned int len)
+{
+    if (len == 0) return 0;
+    
+    std::vector<double> scratch;
+    for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
+    std::sort(scratch.begin(), scratch.end());
+
+    int middle = len/2;
+    if (len % 2 == 0) {
+        return (scratch[middle] + scratch[middle - 1]) / 2;
+    } else {
+        return scratch[middle];
+    }
+}
+
+double MathUtilities::sum(const double *src, unsigned int len)
+{
+    unsigned int i ;
+    double retVal =0.0;
+
+    for(  i = 0; i < len; i++)
+    {
+	retVal += src[ i ];
+    }
+
+    return retVal;
+}
+
+double MathUtilities::mean(const double *src, unsigned int len)
+{
+    double retVal =0.0;
+
+    if (len == 0) return 0;
+
+    double s = sum( src, len );
+	
+    retVal =  s  / (double)len;
+
+    return retVal;
+}
+
+double MathUtilities::mean(const std::vector<double> &src,
+                           unsigned int start,
+                           unsigned int count)
+{
+    double sum = 0.;
+	
+    if (count == 0) return 0;
+    
+    for (int i = 0; i < (int)count; ++i)
+    {
+        sum += src[start + i];
+    }
+
+    return sum / count;
+}
+
+void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
+{
+    unsigned int i;
+    double temp = 0.0;
+
+    if (len == 0) {
+        *min = *max = 0;
+        return;
+    }
+	
+    *min = data[0];
+    *max = data[0];
+
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+
+	if( temp < *min )
+	{
+	    *min =  temp ;
+	}
+	if( temp > *max )
+	{
+	    *max =  temp ;
+	}
+		
+    }
+}
+
+int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
+{
+	unsigned int index = 0;
+	unsigned int i;
+	double temp = 0.0;
+	
+	double max = pData[0];
+
+	for( i = 0; i < Length; i++)
+	{
+		temp = pData[ i ];
+
+		if( temp > max )
+		{
+			max =  temp ;
+			index = i;
+		}
+		
+   	}
+
+	if (pMax) *pMax = max;
+
+
+	return index;
+}
+
+int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
+{
+	unsigned int index = 0;
+	unsigned int i;
+	double temp = 0.0;
+	
+	double max = data[0];
+
+	for( i = 0; i < data.size(); i++)
+	{
+		temp = data[ i ];
+
+		if( temp > max )
+		{
+			max =  temp ;
+			index = i;
+		}
+		
+   	}
+
+	if (pMax) *pMax = max;
+
+
+	return index;
+}
+
+void MathUtilities::circShift( double* pData, int length, int shift)
+{
+	shift = shift % length;
+	double temp;
+	int i,n;
+
+	for( i = 0; i < shift; i++)
+	{
+		temp=*(pData + length - 1);
+
+		for( n = length-2; n >= 0; n--)
+		{
+			*(pData+n+1)=*(pData+n);
+		}
+
+        *pData = temp;
+    }
+}
+
+int MathUtilities::compareInt (const void * a, const void * b)
+{
+  return ( *(int*)a - *(int*)b );
+}
+
+void MathUtilities::normalise(double *data, int length, NormaliseType type)
+{
+    switch (type) {
+
+    case NormaliseNone: return;
+
+    case NormaliseUnitSum:
+    {
+        double sum = 0.0;
+        for (int i = 0; i < length; ++i) {
+            sum += data[i];
+        }
+        if (sum != 0.0) {
+            for (int i = 0; i < length; ++i) {
+                data[i] /= sum;
+            }
+        }
+    }
+    break;
+
+    case NormaliseUnitMax:
+    {
+        double max = 0.0;
+        for (int i = 0; i < length; ++i) {
+            if (fabs(data[i]) > max) {
+                max = fabs(data[i]);
+            }
+        }
+        if (max != 0.0) {
+            for (int i = 0; i < length; ++i) {
+                data[i] /= max;
+            }
+        }
+    }
+    break;
+
+    }
+}
+
+void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
+{
+    switch (type) {
+
+    case NormaliseNone: return;
+
+    case NormaliseUnitSum:
+    {
+        double sum = 0.0;
+        for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
+        if (sum != 0.0) {
+            for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
+        }
+    }
+    break;
+
+    case NormaliseUnitMax:
+    {
+        double max = 0.0;
+        for (int i = 0; i < (int)data.size(); ++i) {
+            if (fabs(data[i]) > max) max = fabs(data[i]);
+        }
+        if (max != 0.0) {
+            for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
+        }
+    }
+    break;
+
+    }
+}
+
+void MathUtilities::adaptiveThreshold(std::vector<double> &data)
+{
+    int sz = int(data.size());
+    if (sz == 0) return;
+
+    std::vector<double> smoothed(sz);
+	
+    int p_pre = 8;
+    int p_post = 7;
+
+    for (int i = 0; i < sz; ++i) {
+
+        int first = std::max(0,      i - p_pre);
+        int last  = std::min(sz - 1, i + p_post);
+
+        smoothed[i] = mean(data, first, last - first + 1);
+    }
+
+    for (int i = 0; i < sz; i++) {
+        data[i] -= smoothed[i];
+        if (data[i] < 0.0) data[i] = 0.0;
+    }
+}
+
+bool
+MathUtilities::isPowerOfTwo(int x)
+{
+    if (x < 1) return false;
+    if (x & (x-1)) return false;
+    return true;
+}
+
+int
+MathUtilities::nextPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    if (x < 1) return 1;
+    int n = 1;
+    while (x) { x >>= 1; n <<= 1; }
+    return n;
+}
+
+int
+MathUtilities::previousPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    if (x < 1) return 1;
+    int n = 1;
+    x >>= 1;
+    while (x) { x >>= 1; n <<= 1; }
+    return n;
+}
+
+int
+MathUtilities::nearestPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
+    if (x - n0 < n1 - x) return n0;
+    else return n1;
+}
+
+double
+MathUtilities::factorial(int x)
+{
+    if (x < 0) return 0;
+    double f = 1;
+    for (int i = 1; i <= x; ++i) {
+	f = f * i;
+    }
+    return f;
+}
+
+int
+MathUtilities::gcd(int a, int b)
+{
+    int c = a % b;
+    if (c == 0) {
+        return b;
+    } else {
+        return gcd(b, c);
+    }
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/MathUtilities.h	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,153 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#ifndef MATHUTILITIES_H
+#define MATHUTILITIES_H
+
+#include <vector>
+
+#include "nan-inf.h"
+
+/**
+ * Static helper functions for simple mathematical calculations.
+ */
+class MathUtilities  
+{
+public:	
+    /**
+     * Round x to the nearest integer.
+     */
+    static double round( double x );
+
+    /**
+     * Return through min and max pointers the highest and lowest
+     * values in the given array of the given length.
+     */
+    static void	  getFrameMinMax( const double* data, unsigned int len,  double* min, double* max );
+
+    /**
+     * Return the mean of the given array of the given length.
+     */
+    static double mean( const double* src, unsigned int len );
+
+    /**
+     * Return the mean of the subset of the given vector identified by
+     * start and count.
+     */
+    static double mean( const std::vector<double> &data,
+                        unsigned int start, unsigned int count );
+    
+    /**
+     * Return the sum of the values in the given array of the given
+     * length.
+     */
+    static double sum( const double* src, unsigned int len );
+
+    /**
+     * Return the median of the values in the given array of the given
+     * length. If the array is even in length, the returned value will
+     * be half-way between the two values adjacent to median.
+     */
+    static double median( const double* src, unsigned int len );
+
+    /**
+     * The principle argument function. Map the phase angle ang into
+     * the range [-pi,pi).
+     */
+    static double princarg( double ang );
+
+    /**
+     * Floating-point division modulus: return x % y.
+     */
+    static double mod( double x, double y);
+
+    static void	  getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm);
+    static double getAlphaNorm(const std::vector <double> &data, unsigned int alpha );
+
+    static void   circShift( double* data, int length, int shift);
+
+    static int	  getMax( double* data, unsigned int length, double* max = 0 );
+    static int	  getMax( const std::vector<double> &data, double* max = 0 );
+    static int    compareInt(const void * a, const void * b);
+
+    enum NormaliseType {
+        NormaliseNone,
+        NormaliseUnitSum,
+        NormaliseUnitMax
+    };
+
+    static void normalise(double *data, int length,
+                          NormaliseType n = NormaliseUnitMax);
+
+    static void normalise(std::vector<double> &data,
+                          NormaliseType n = NormaliseUnitMax);
+
+    /**
+     * Threshold the input/output vector data against a moving-mean
+     * average filter.
+     */
+    static void adaptiveThreshold(std::vector<double> &data);
+
+    /** 
+     * Return true if x is 2^n for some integer n >= 0.
+     */
+    static bool isPowerOfTwo(int x);
+
+    /**
+     * Return the next higher integer power of two from x, e.g. 1300
+     * -> 2048, 2048 -> 2048.
+     */
+    static int nextPowerOfTwo(int x);
+
+    /**
+     * Return the next lower integer power of two from x, e.g. 1300 ->
+     * 1024, 2048 -> 2048.
+     */
+    static int previousPowerOfTwo(int x);
+
+    /**
+     * Return the nearest integer power of two to x, e.g. 1300 -> 1024,
+     * 12 -> 16 (not 8; if two are equidistant, the higher is returned).
+     */
+    static int nearestPowerOfTwo(int x);
+
+    /**
+     * Return x!
+     */
+    static double factorial(int x); // returns double in case it is large
+
+    /**
+     * Return the greatest common divisor of natural numbers a and b.
+     */
+    static int gcd(int a, int b);
+};
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/Resampler.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,433 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#include "Resampler.h"
+
+#include "maths/MathUtilities.h"
+#include "base/KaiserWindow.h"
+#include "base/SincWindow.h"
+#include "thread/Thread.h"
+
+#include <iostream>
+#include <vector>
+#include <map>
+#include <cassert>
+
+using std::vector;
+using std::map;
+using std::cerr;
+using std::endl;
+
+//#define DEBUG_RESAMPLER 1
+//#define DEBUG_RESAMPLER_VERBOSE 1
+
+Resampler::Resampler(int sourceRate, int targetRate) :
+    m_sourceRate(sourceRate),
+    m_targetRate(targetRate)
+{
+    initialise(100, 0.02);
+}
+
+Resampler::Resampler(int sourceRate, int targetRate, 
+                     double snr, double bandwidth) :
+    m_sourceRate(sourceRate),
+    m_targetRate(targetRate)
+{
+    initialise(snr, bandwidth);
+}
+
+Resampler::~Resampler()
+{
+    delete[] m_phaseData;
+}
+
+// peakToPole -> length -> beta -> window
+static map<double, map<int, map<double, vector<double> > > >
+knownFilters;
+
+static Mutex
+knownFilterMutex;
+
+void
+Resampler::initialise(double snr, double bandwidth)
+{
+    int higher = std::max(m_sourceRate, m_targetRate);
+    int lower = std::min(m_sourceRate, m_targetRate);
+
+    m_gcd = MathUtilities::gcd(lower, higher);
+    m_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, higher / m_gcd);
+
+    params.length =
+	(params.length % 2 == 0 ? params.length + 1 : params.length);
+    
+    params.length =
+        (params.length > 200001 ? 200001 : params.length);
+
+    m_filterLength = params.length;
+
+    vector<double> filter;
+    knownFilterMutex.lock();
+
+    if (knownFilters[m_peakToPole][m_filterLength].find(params.beta) ==
+	knownFilters[m_peakToPole][m_filterLength].end()) {
+
+	KaiserWindow kw(params);
+	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[m_peakToPole][m_filterLength][params.beta] = filter;
+    }
+
+    filter = knownFilters[m_peakToPole][m_filterLength][params.beta];
+    knownFilterMutex.unlock();
+
+    int inputSpacing = m_targetRate / m_gcd;
+    int outputSpacing = m_sourceRate / m_gcd;
+
+#ifdef DEBUG_RESAMPLER
+    cerr << "resample " << m_sourceRate << " -> " << m_targetRate
+         << ": inputSpacing " << inputSpacing << ", outputSpacing "
+         << outputSpacing << ": filter length " << m_filterLength
+         << endl;
+#endif
+
+    // Now we have a filter of (odd) length flen in which the lower
+    // sample rate corresponds to every n'th point and the higher rate
+    // to every m'th where n and m are higher and lower rates divided
+    // by their gcd respectively. So if x coordinates are on the same
+    // scale as our filter resolution, then source sample i is at i *
+    // (targetRate / gcd) and target sample j is at j * (sourceRate /
+    // gcd).
+
+    // To reconstruct a single target sample, we want a buffer (real
+    // or virtual) of flen values formed of source samples spaced at
+    // intervals of (targetRate / gcd), in our example case 3.  This
+    // is initially formed with the first sample at the filter peak.
+    //
+    // 0  0  0  0  a  0  0  b  0
+    //
+    // and of course we have our filter
+    //
+    // f1 f2 f3 f4 f5 f6 f7 f8 f9
+    //
+    // We take the sum of products of non-zero values from this buffer
+    // with corresponding values in the filter
+    //
+    // a * f5 + b * f8
+    //
+    // Then we drop (sourceRate / gcd) values, in our example case 4,
+    // from the start of the buffer and fill until it has flen values
+    // again
+    //
+    // a  0  0  b  0  0  c  0  0
+    //
+    // repeat to reconstruct the next target sample
+    //
+    // a * f1 + b * f4 + c * f7
+    //
+    // and so on.
+    //
+    // Above I said the buffer could be "real or virtual" -- ours is
+    // virtual. We don't actually store all the zero spacing values,
+    // except for padding at the start; normally we store only the
+    // values that actually came from the source stream, along with a
+    // phase value that tells us how many virtual zeroes there are at
+    // the start of the virtual buffer.  So the two examples above are
+    //
+    // 0 a b  [ with phase 1 ]
+    // a b c  [ with phase 0 ]
+    //
+    // Having thus broken down the buffer so that only the elements we
+    // need to multiply are present, we can also unzip the filter into
+    // every-nth-element subsets at each phase, allowing us to do the
+    // filter multiplication as a simply vector multiply. That is, rather
+    // than store
+    //
+    // f1 f2 f3 f4 f5 f6 f7 f8 f9
+    // 
+    // we store separately
+    //
+    // f1 f4 f7
+    // f2 f5 f8
+    // f3 f6 f9
+    //
+    // Each time we complete a multiply-and-sum, we need to work out
+    // how many (real) samples to drop from the start of our buffer,
+    // and how many to add at the end of it for the next multiply.  We
+    // know we want to drop enough real samples to move along by one
+    // computed output sample, which is our outputSpacing number of
+    // virtual buffer samples. Depending on the relationship between
+    // input and output spacings, this may mean dropping several real
+    // samples, one real sample, or none at all (and simply moving to
+    // a different "phase").
+
+    m_phaseData = new Phase[inputSpacing];
+
+    for (int phase = 0; phase < inputSpacing; ++phase) {
+
+	Phase p;
+
+	p.nextPhase = phase - outputSpacing;
+	while (p.nextPhase < 0) p.nextPhase += inputSpacing;
+	p.nextPhase %= inputSpacing;
+	
+	p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase))
+			  / inputSpacing));
+
+	int filtZipLength = int(ceil(double(m_filterLength - phase)
+				     / inputSpacing));
+
+	for (int i = 0; i < filtZipLength; ++i) {
+	    p.filter.push_back(filter[i * inputSpacing + phase]);
+	}
+
+	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
+    // those. This means (among other things) that the source stream
+    // can be asked for enough samples up-front to fill the buffer
+    // before the first output sample is generated.
+    // 
+    // In this implementation we're using a push model in which a
+    // certain number of source samples is provided and we're asked
+    // for as many output samples as that makes available. But we
+    // can't return any samples from the beginning until half the
+    // filter length has been provided as input. This means we must
+    // either return a very variable number of samples (none at all
+    // until the filter fills, then half the filter length at once) or
+    // else have a lengthy declared latency on the output. We do the
+    // latter. (What do other implementations do?)
+    //
+    // We want to make sure the first "real" sample will eventually be
+    // aligned with the centre sample in the filter (it's tidier, and
+    // easier to do diagnostic calculations that way). So we need to
+    // pick the initial phase and buffer fill accordingly.
+    // 
+    // Example: if the inputSpacing is 2, outputSpacing is 3, and
+    // filter length is 7,
+    // 
+    //    x     x     x     x     a     b     c ... input samples
+    // 0  1  2  3  4  5  6  7  8  9 10 11 12 13 ... 
+    //          i        j        k        l    ... output samples
+    // [--------|--------] <- filter with centre mark
+    //
+    // Let h be the index of the centre mark, here 3 (generally
+    // int(filterLength/2) for odd-length filters).
+    //
+    // The smallest n such that h + n * outputSpacing > filterLength
+    // is 2 (that is, ceil((filterLength - h) / outputSpacing)), and
+    // (h + 2 * outputSpacing) % inputSpacing == 1, so the initial
+    // phase is 1.
+    //
+    // To achieve our n, we need to pre-fill the "virtual" buffer with
+    // 4 zero samples: the x's above. This is int((h + n *
+    // outputSpacing) / inputSpacing). It's the phase that makes this
+    // buffer get dealt with in such a way as to give us an effective
+    // index for sample a of 9 rather than 8 or 10 or whatever.
+    //
+    // This gives us output latency of 2 (== n), i.e. output samples i
+    // and j will appear before the one in which input sample a is at
+    // the centre of the filter.
+
+    int h = int(m_filterLength / 2);
+    int n = ceil(double(m_filterLength - h) / outputSpacing);
+    
+    m_phase = (h + n * outputSpacing) % inputSpacing;
+
+    int fill = (h + n * outputSpacing) / inputSpacing;
+    
+    m_latency = n;
+
+    m_buffer = vector<double>(fill, 0);
+    m_bufferOrigin = 0;
+
+#ifdef DEBUG_RESAMPLER
+    cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
+	      << ", latency " << m_latency << endl;
+#endif
+}
+
+double
+Resampler::reconstructOne()
+{
+    Phase &pd = m_phaseData[m_phase];
+    double v = 0.0;
+    int n = pd.filter.size();
+
+    assert(n + m_bufferOrigin <= (int)m_buffer.size());
+
+    const double *const __restrict__ buf = m_buffer.data() + m_bufferOrigin;
+    const double *const __restrict__ filt = pd.filter.data();
+
+    for (int i = 0; i < n; ++i) {
+	// NB gcc can only vectorize this with -ffast-math
+	v += buf[i] * filt[i];
+    }
+
+    m_bufferOrigin += pd.drop;
+    m_phase = pd.nextPhase;
+    return v;
+}
+
+int
+Resampler::process(const double *src, double *dst, int n)
+{
+    for (int i = 0; i < n; ++i) {
+	m_buffer.push_back(src[i]);
+    }
+
+    int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate));
+    int outidx = 0;
+
+#ifdef DEBUG_RESAMPLER
+    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;
+
+    while (outidx < maxout &&
+	   m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) {
+	dst[outidx] = scaleFactor * reconstructOne();
+	outidx++;
+    }
+
+    m_buffer = vector<double>(m_buffer.begin() + m_bufferOrigin, m_buffer.end());
+    m_bufferOrigin = 0;
+    
+    return outidx;
+}
+    
+vector<double>
+Resampler::process(const double *src, int n)
+{
+    int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate));
+    vector<double> out(maxout, 0.0);
+    int got = process(src, out.data(), n);
+    assert(got <= maxout);
+    if (got < maxout) out.resize(got);
+    return out;
+}
+
+vector<double>
+Resampler::resample(int sourceRate, int targetRate, const double *data, int n)
+{
+    Resampler r(sourceRate, targetRate);
+
+    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));
+
+    // 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));
+    
+#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 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
+    cerr << "resample: " << n << " in, " << pad.size() << " padding, " << got << " out (" << gotData << " data, " << gotPad << " padding, latency = " << latency << ")" << endl;
+#endif
+#ifdef DEBUG_RESAMPLER_VERBOSE
+    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] << " ";
+    }
+    cerr << endl;
+#endif
+
+    int toReturn = got - latency;
+    if (toReturn > m) toReturn = m;
+
+    vector<double> sliced(out.begin() + latency, 
+			  out.begin() + latency + toReturn);
+
+#ifdef DEBUG_RESAMPLER_VERBOSE
+    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] << " ";
+    }
+    cerr << endl;
+#endif
+
+    return sliced;
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/Resampler.h	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,119 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#ifndef RESAMPLER_H
+#define RESAMPLER_H
+
+#include <vector>
+
+/**
+ * Resampler resamples a stream from one integer sample rate to
+ * another (arbitrary) rate, using a kaiser-windowed sinc filter.  The
+ * results and performance are pretty similar to libraries such as
+ * libsamplerate, though this implementation does not support
+ * time-varying ratios (the ratio is fixed on construction).
+ *
+ * See also Decimator, which is faster and rougher but supports only
+ * power-of-two downsampling factors.
+ */
+class Resampler
+{
+public:
+    /**
+     * Construct a Resampler to resample from sourceRate to
+     * targetRate.
+     */
+    Resampler(int sourceRate, int targetRate);
+
+    /**
+     * Construct a Resampler to resample from sourceRate to
+     * targetRate, using the given filter parameters.
+     */
+    Resampler(int sourceRate, int targetRate,
+              double snr, double bandwidth);
+
+    virtual ~Resampler();
+
+    /**
+     * Read n input samples from src and write resampled data to
+     * dst. The return value is the number of samples written, which
+     * will be no more than ceil((n * targetRate) / sourceRate). The
+     * caller must ensure the dst buffer has enough space for the
+     * samples returned.
+     */
+    int process(const double *src, double *dst, int n);
+
+    /**
+     * Read n input samples from src and return resampled data by
+     * value.
+     */
+    std::vector<double> process(const double *src, int n);
+
+    /**
+     * Return the number of samples of latency at the output due by
+     * the filter. (That is, the output will be delayed by this number
+     * of samples relative to the input.)
+     */
+    int getLatency() const { return m_latency; }
+
+    /**
+     * Carry out a one-off resample of a single block of n
+     * samples. The output is latency-compensated.
+     */
+    static std::vector<double> resample
+    (int sourceRate, int targetRate, const double *data, int n);
+
+private:
+    int m_sourceRate;
+    int m_targetRate;
+    int m_gcd;
+    int m_filterLength;
+    int m_bufferLength;
+    int m_latency;
+    double m_peakToPole;
+    
+    struct Phase {
+        int nextPhase;
+        std::vector<double> filter;
+        int drop;
+    };
+
+    Phase *m_phaseData;
+    int m_phase;
+    std::vector<double> m_buffer;
+    int m_bufferOrigin;
+
+    void initialise(double, double);
+    double reconstructOne();
+};
+
+#endif
+    
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/SincWindow.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,63 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#include "SincWindow.h"
+
+#include <cmath>
+
+void
+SincWindow::init()
+{
+    if (m_length < 1) {
+	return;
+    } else if (m_length < 2) {
+	m_window.push_back(1);
+	return;
+    } else {
+
+	int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2);
+	int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2);
+	double m = 2 * M_PI / m_p;
+
+	for (int i = 0; i < n0; ++i) {
+	    double x = ((m_length / 2) - i) * m;
+	    m_window.push_back(sin(x) / x);
+	}
+
+	m_window.push_back(1.0);
+
+	for (int i = 1; i < n1; ++i) {
+	    double x = i * m;
+	    m_window.push_back(sin(x) / x);
+	}
+    }
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dsp/SincWindow.h	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,79 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#ifndef SINC_WINDOW_H
+#define SINC_WINDOW_H
+
+#include <vector>
+
+/**
+ * A window containing values of the sinc function, i.e. sin(x)/x with
+ * sinc(0) == 1, with x == 0 at the centre.
+ */
+class SincWindow
+{
+public:
+    /**
+     * Construct a windower of the given length, containing the values
+     * of sinc(x) with x=0 in the middle, i.e. at sample (length-1)/2
+     * for odd or (length/2)+1 for even length, such that the distance
+     * from -pi to pi (the nearest zero crossings either side of the
+     * peak) is p samples.
+     */
+    SincWindow(int length, double p) : m_length(length), m_p(p) { init(); }
+
+    int getLength() const {
+	return m_length;
+    }
+
+    const double *getWindow() const { 
+	return m_window.data();
+    }
+
+    void cut(double *src) const { 
+	cut(src, src); 
+    }
+
+    void cut(const double *src, double *dst) const {
+	for (int i = 0; i < m_length; ++i) {
+	    dst[i] = src[i] * m_window[i];
+	}
+    }
+
+private:
+    int m_length;
+    double m_p;
+    std::vector<double> m_window;
+
+    void init();
+};
+
+#endif
--- a/src/processfile.cpp	Thu May 15 12:06:31 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,294 +0,0 @@
-
-#include "ConstantQ.h"
-#include "CQInverse.h"
-
-#include <sndfile.h>
-
-#include <iostream>
-
-using std::vector;
-using std::cerr;
-using std::endl;
-
-#include <cstring>
-
-#include <getopt.h>
-#include <unistd.h>
-#include <sys/time.h>
-#include <cstdlib>
-
-int main(int argc, char **argv)
-{
-    double maxFreq = 0;
-    double minFreq = 0;
-    int bpo = 0;
-    bool help = false;
-    
-    int c;
-
-    while (1) {
-	int optionIndex = 0;
-	
-	static struct option longOpts[] = {
-	    { "help", 0, 0, 'h', },
-	    { "maxfreq", 1, 0, 'x', },
-	    { "minfreq", 1, 0, 'n', },
-	    { "bpo", 1, 0, 'b' },
-	    { 0, 0, 0, 0 },
-	};
-
-	c = getopt_long(argc, argv,
-			"hx:n:b:",
-			longOpts, &optionIndex);
-	if (c == -1) break;
-
-	switch (c) {
-	case 'h': help = true; break;
-	case 'x': maxFreq = atof(optarg); break;
-	case 'n': minFreq = atof(optarg); break;
-	case 'b': bpo = atoi(optarg); break;
-	default: help = true; break;
-	}
-    }
-
-    if (help || (optind + 2 != argc && optind + 3 != argc)) {
-	cerr << endl;
-	cerr << "Usage: " << argv[0] << " [options] infile.wav outfile.wav [differencefile.wav]" << endl;
-	cerr << endl;
-	cerr << "Options:" << endl;
-	cerr << "  -x<X>, --maxfreq <X>  Maximum frequency (default = sample rate / 3)" << endl;
-	cerr << "  -n<X>, --minfreq <X>  Minimum frequency (default = 100, actual min may vary)" << endl;
-	cerr << "  -b<X>, --bpo <X>      Bins per octave   (default = 60)" << endl;
-	cerr << "  -h, --help            Print this help" << endl;
-	cerr << endl;
-	cerr << "This rather useless program simply performs a forward Constant-Q transform with" << endl;
-	cerr << "the requested parameters, followed by its inverse, and writes the result to the" << endl;
-	cerr << "output file. If a diff file name is provided, the difference between input and" << endl;
-	cerr << "output signals is also written to that. All this accomplishes is to produce a" << endl;
-	cerr << "signal that approximates the input: it's intended for test purposes only." << endl;
-	cerr << endl;
-	cerr << "(Want to calculate and obtain a Constant-Q spectrogram? Use the CQVamp plugin" << endl;
-	cerr << "in a Vamp plugin host.)" << endl;
-	cerr << endl;
-	return 2;
-    }
-
-    char *fileName = strdup(argv[optind++]);
-    char *fileNameOut = strdup(argv[optind++]);
-    char *diffFileName = (optind < argc ? strdup(argv[optind++]) : 0);
-    bool doDiff = (diffFileName != 0);
-
-    SNDFILE *sndfile;
-    SNDFILE *sndfileOut;
-    SNDFILE *sndDiffFile = 0;
-    SF_INFO sfinfo;
-    SF_INFO sfinfoOut;
-    SF_INFO sfinfoDiff;
-    memset(&sfinfo, 0, sizeof(SF_INFO));
-
-    sndfile = sf_open(fileName, SFM_READ, &sfinfo);
-    if (!sndfile) {
-	cerr << "ERROR: Failed to open input file \"" << fileName << "\": "
-	     << sf_strerror(sndfile) << endl;
-	return 1;
-    }
-
-    sfinfoOut.channels = 1;
-    sfinfoOut.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
-    sfinfoOut.frames = sfinfo.frames;
-    sfinfoOut.samplerate = sfinfo.samplerate;
-    sfinfoOut.sections = sfinfo.sections;
-    sfinfoOut.seekable = sfinfo.seekable;
-
-    sndfileOut = sf_open(fileNameOut, SFM_WRITE, &sfinfoOut) ;
-    if (!sndfileOut) {
-	cerr << "ERROR: Failed to open output file \"" << fileNameOut << "\" for writing: "
-	     << sf_strerror(sndfileOut) << endl;
-	return 1;
-    }
-
-    if (doDiff) {
-	sfinfoDiff = sfinfoOut;
-	sndDiffFile = sf_open(diffFileName, SFM_WRITE, &sfinfoDiff);
-	if (!sndDiffFile) {
-	    cerr << "ERROR: Failed to open diff output file \"" << diffFileName << "\" for writing: "
-		 << sf_strerror(sndDiffFile) << endl;
-	    return 1;
-	}
-    }
-
-    int ibs = 1024;
-    int channels = sfinfo.channels;
-    float *fbuf = new float[channels * ibs];
-
-    if (maxFreq == 0.0) maxFreq = sfinfo.samplerate / 3;
-    if (minFreq == 0.0) minFreq = 100;
-    if (bpo == 0) bpo = 60;
-
-    ConstantQ cq(sfinfo.samplerate, minFreq, maxFreq, bpo);
-    CQInverse cqi(sfinfo.samplerate, minFreq, maxFreq, bpo);
-
-    cerr << "max freq = " << cq.getMaxFrequency() << ", min freq = "
-	 << cq.getMinFrequency() << ", octaves = " << cq.getOctaves() << endl;
-
-    cerr << "octave boundaries: ";
-    for (int i = 0; i < cq.getOctaves(); ++i) {
-	cerr << cq.getMaxFrequency() / pow(2, i) << " ";
-    }
-    cerr << endl;
-
-    int inframe = 0;
-    int outframe = 0;
-    int latency = cq.getLatency() + cqi.getLatency();
-
-    vector<double> buffer;
-
-    double maxdiff = 0.0;
-    int maxdiffidx = 0;
-
-    cerr << "forward latency = " << cq.getLatency() << ", inverse latency = " 
-	 << cqi.getLatency() << ", total = " << latency << endl;
-
-    timeval tv;
-    (void)gettimeofday(&tv, 0);
-
-    while (inframe < sfinfo.frames) {
-
-        int count = -1;
-	
-	if ((count = sf_readf_float(sndfile, fbuf, ibs)) < 0) {
-	    break;
-	}
-
-	vector<double> cqin;
-	for (int i = 0; i < count; ++i) {
-	    double v = fbuf[i * channels];
-	    if (channels > 1) {
-		for (int c = 1; c < channels; ++c) {
-		    v += fbuf[i * channels + c];
-		}
-		v /= channels;
-	    }
-	    cqin.push_back(v);
-	}
-
-	if (doDiff) {
-	    buffer.insert(buffer.end(), cqin.begin(), cqin.end());
-	}
-	
-	vector<double> cqout = cqi.process(cq.process(cqin));
-
-	for (int i = 0; i < int(cqout.size()); ++i) {
-	    if (cqout[i] > 1.0) cqout[i] = 1.0;
-	    if (cqout[i] < -1.0) cqout[i] = -1.0;
-	}
-
-	if (outframe >= latency) {
-
-	    sf_writef_double(sndfileOut, 
-			     cqout.data(), 
-			     cqout.size());
-
-	} else if (outframe + (int)cqout.size() >= latency) {
-
-	    int offset = latency - outframe;
-	    sf_writef_double(sndfileOut, 
-			     cqout.data() + offset,
-			     cqout.size() - offset);
-	}
-
-	if (doDiff) {
-	    for (int i = 0; i < (int)cqout.size(); ++i) {
-		if (outframe + i >= latency) {
-		    int dframe = outframe + i - latency;
-		    if (dframe >= (int)buffer.size()) cqout[i] = 0;
-		    else cqout[i] -= buffer[dframe];
-		    if (fabs(cqout[i]) > maxdiff &&
-			dframe > sfinfo.samplerate && // ignore first/last sec
-			dframe + sfinfo.samplerate < sfinfo.frames) {
-			maxdiff = fabs(cqout[i]);
-			maxdiffidx = dframe;
-		    }
-		}
-	    }
-	    
-	    if (outframe >= latency) {
-
-		sf_writef_double(sndDiffFile, 
-				 cqout.data(), 
-				 cqout.size());
-
-	    } else if (outframe + (int)cqout.size() >= latency) {
-
-		int offset = latency - outframe;
-		sf_writef_double(sndDiffFile, 
-				 cqout.data() + offset,
-				 cqout.size() - offset);
-	    }
-	}
-
-	inframe += count;
-	outframe += cqout.size();
-    }
-
-    vector<double> r = cqi.process(cq.getRemainingOutput());
-    vector<double> r2 = cqi.getRemainingOutput();
-
-    r.insert(r.end(), r2.begin(), r2.end());
-
-    for (int i = 0; i < int(r.size()); ++i) {
-	if (r[i] > 1.0) r[i] = 1.0;
-	if (r[i] < -1.0) r[i] = -1.0;
-    }
-
-    sf_writef_double(sndfileOut, r.data(), r.size());
-    if (doDiff) {
-	for (int i = 0; i < (int)r.size(); ++i) {
-	    if (outframe + i >= latency) {
-		int dframe = outframe + i - latency;
-		if (dframe >= (int)buffer.size()) r[i] = 0;
-		else r[i] -= buffer[dframe];
-		if (fabs(r[i]) > maxdiff &&
-		    dframe > sfinfo.samplerate && // ignore first/last sec
-		    dframe + sfinfo.samplerate < sfinfo.frames) {
-		    maxdiff = fabs(r[i]);
-		    maxdiffidx = dframe;
-		}
-	    }
-	}
-	sf_writef_double(sndDiffFile, r.data(), r.size());
-    }
-    outframe += r.size();
-
-    sf_close(sndfile);
-    sf_close(sndfileOut);
-
-    if (doDiff) {
-	sf_close(sndDiffFile);
-    }
-
-    cerr << "in: " << inframe << ", out: " << outframe - latency << endl;
-
-    if (doDiff) {
-	double db = 10 * log10(maxdiff);
-	cerr << "max diff [excluding first and last second of audio] is "
-	     << maxdiff << " (" << db << " dBFS)"
-	     << " at sample index " << maxdiffidx << endl;
-    }
-    
-    timeval etv;
-    (void)gettimeofday(&etv, 0);
-        
-    etv.tv_sec -= tv.tv_sec;
-    if (etv.tv_usec < tv.tv_usec) {
-	etv.tv_usec += 1000000;
-	etv.tv_sec -= 1;
-    }
-    etv.tv_usec -= tv.tv_usec;
-        
-    double sec = double(etv.tv_sec) + (double(etv.tv_usec) / 1000000.0);
-    cerr << "elapsed time (not counting init): " << sec << " sec, frames/sec at input: " << inframe/sec << endl;
-
-    return 0;
-}
-
--- a/src/test.cpp	Thu May 15 12:06:31 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,76 +0,0 @@
-/*
-    Constant-Q library
-    Copyright (c) 2013-2014 Queen Mary, University of London
-
-    Permission is hereby granted, free of charge, to any person
-    obtaining a copy of this software and associated documentation
-    files (the "Software"), to deal in the Software without
-    restriction, including without limitation the rights to use, copy,
-    modify, merge, publish, distribute, sublicense, and/or sell copies
-    of the Software, and to permit persons to whom the Software is
-    furnished to do so, subject to the following conditions:
-
-    The above copyright notice and this permission notice shall be
-    included in all copies or substantial portions of the Software.
-
-    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
-    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
-    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
-    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-
-    Except as contained in this notice, the names of the Centre for
-    Digital Music; Queen Mary, University of London; and Chris Cannam
-    shall not be used in advertising or otherwise to promote the sale,
-    use or other dealings in this Software without prior written
-    authorization.
-*/
-
-#include "CQSpectrogram.h"
-
-#include <iostream>
-#include <vector>
-
-using std::vector;
-using std::cerr;
-using std::cout;
-using std::endl;
-
-#include <cstdio>
-
-int main(int argc, char **argv)
-{
-    vector<double> in;
-
-    for (int i = 0; i < 64; ++i) {
-//	if (i == 0) in.push_back(1);
-//	else in.push_back(0);
-	in.push_back(sin(i * M_PI / 2.0));
-    }
-
-    CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros);
-
-    vector<vector<double> > out = k.process(in);
-    vector<vector<double> > rest = k.getRemainingOutput();
-
-    out.insert(out.end(), rest.begin(), rest.end());
-
-    cerr << "got " << out.size() << " back (" << out[0].size() << " in each?)" << endl;
-
-    for (int b = 0; b < (int)out.size() / 8; ++b) {
-	printf("\nColumns %d to %d:\n\n", b * 8, b * 8 + 7);
-	for (int j = int(out[0].size()) - 1; j >= 0; --j) {
-	    for (int i = 0; i < 8; ++i) {
-		if (i + b * 8 < (int)out.size()) {
-		    double v = out[i + b * 8][j];
-		    if (v < 0.0001) printf("  0      ");
-		    else printf("  %.4f ", out[i + b * 8][j]);
-		}
-	    }
-	    printf("\n");
-	}
-    }
-}
-
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/processfile.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,294 @@
+
+#include "ConstantQ.h"
+#include "CQInverse.h"
+
+#include <sndfile.h>
+
+#include <iostream>
+
+using std::vector;
+using std::cerr;
+using std::endl;
+
+#include <cstring>
+
+#include <getopt.h>
+#include <unistd.h>
+#include <sys/time.h>
+#include <cstdlib>
+
+int main(int argc, char **argv)
+{
+    double maxFreq = 0;
+    double minFreq = 0;
+    int bpo = 0;
+    bool help = false;
+    
+    int c;
+
+    while (1) {
+	int optionIndex = 0;
+	
+	static struct option longOpts[] = {
+	    { "help", 0, 0, 'h', },
+	    { "maxfreq", 1, 0, 'x', },
+	    { "minfreq", 1, 0, 'n', },
+	    { "bpo", 1, 0, 'b' },
+	    { 0, 0, 0, 0 },
+	};
+
+	c = getopt_long(argc, argv,
+			"hx:n:b:",
+			longOpts, &optionIndex);
+	if (c == -1) break;
+
+	switch (c) {
+	case 'h': help = true; break;
+	case 'x': maxFreq = atof(optarg); break;
+	case 'n': minFreq = atof(optarg); break;
+	case 'b': bpo = atoi(optarg); break;
+	default: help = true; break;
+	}
+    }
+
+    if (help || (optind + 2 != argc && optind + 3 != argc)) {
+	cerr << endl;
+	cerr << "Usage: " << argv[0] << " [options] infile.wav outfile.wav [differencefile.wav]" << endl;
+	cerr << endl;
+	cerr << "Options:" << endl;
+	cerr << "  -x<X>, --maxfreq <X>  Maximum frequency (default = sample rate / 3)" << endl;
+	cerr << "  -n<X>, --minfreq <X>  Minimum frequency (default = 100, actual min may vary)" << endl;
+	cerr << "  -b<X>, --bpo <X>      Bins per octave   (default = 60)" << endl;
+	cerr << "  -h, --help            Print this help" << endl;
+	cerr << endl;
+	cerr << "This rather useless program simply performs a forward Constant-Q transform with" << endl;
+	cerr << "the requested parameters, followed by its inverse, and writes the result to the" << endl;
+	cerr << "output file. If a diff file name is provided, the difference between input and" << endl;
+	cerr << "output signals is also written to that. All this accomplishes is to produce a" << endl;
+	cerr << "signal that approximates the input: it's intended for test purposes only." << endl;
+	cerr << endl;
+	cerr << "(Want to calculate and obtain a Constant-Q spectrogram? Use the CQVamp plugin" << endl;
+	cerr << "in a Vamp plugin host.)" << endl;
+	cerr << endl;
+	return 2;
+    }
+
+    char *fileName = strdup(argv[optind++]);
+    char *fileNameOut = strdup(argv[optind++]);
+    char *diffFileName = (optind < argc ? strdup(argv[optind++]) : 0);
+    bool doDiff = (diffFileName != 0);
+
+    SNDFILE *sndfile;
+    SNDFILE *sndfileOut;
+    SNDFILE *sndDiffFile = 0;
+    SF_INFO sfinfo;
+    SF_INFO sfinfoOut;
+    SF_INFO sfinfoDiff;
+    memset(&sfinfo, 0, sizeof(SF_INFO));
+
+    sndfile = sf_open(fileName, SFM_READ, &sfinfo);
+    if (!sndfile) {
+	cerr << "ERROR: Failed to open input file \"" << fileName << "\": "
+	     << sf_strerror(sndfile) << endl;
+	return 1;
+    }
+
+    sfinfoOut.channels = 1;
+    sfinfoOut.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
+    sfinfoOut.frames = sfinfo.frames;
+    sfinfoOut.samplerate = sfinfo.samplerate;
+    sfinfoOut.sections = sfinfo.sections;
+    sfinfoOut.seekable = sfinfo.seekable;
+
+    sndfileOut = sf_open(fileNameOut, SFM_WRITE, &sfinfoOut) ;
+    if (!sndfileOut) {
+	cerr << "ERROR: Failed to open output file \"" << fileNameOut << "\" for writing: "
+	     << sf_strerror(sndfileOut) << endl;
+	return 1;
+    }
+
+    if (doDiff) {
+	sfinfoDiff = sfinfoOut;
+	sndDiffFile = sf_open(diffFileName, SFM_WRITE, &sfinfoDiff);
+	if (!sndDiffFile) {
+	    cerr << "ERROR: Failed to open diff output file \"" << diffFileName << "\" for writing: "
+		 << sf_strerror(sndDiffFile) << endl;
+	    return 1;
+	}
+    }
+
+    int ibs = 1024;
+    int channels = sfinfo.channels;
+    float *fbuf = new float[channels * ibs];
+
+    if (maxFreq == 0.0) maxFreq = sfinfo.samplerate / 3;
+    if (minFreq == 0.0) minFreq = 100;
+    if (bpo == 0) bpo = 60;
+
+    ConstantQ cq(sfinfo.samplerate, minFreq, maxFreq, bpo);
+    CQInverse cqi(sfinfo.samplerate, minFreq, maxFreq, bpo);
+
+    cerr << "max freq = " << cq.getMaxFrequency() << ", min freq = "
+	 << cq.getMinFrequency() << ", octaves = " << cq.getOctaves() << endl;
+
+    cerr << "octave boundaries: ";
+    for (int i = 0; i < cq.getOctaves(); ++i) {
+	cerr << cq.getMaxFrequency() / pow(2, i) << " ";
+    }
+    cerr << endl;
+
+    int inframe = 0;
+    int outframe = 0;
+    int latency = cq.getLatency() + cqi.getLatency();
+
+    vector<double> buffer;
+
+    double maxdiff = 0.0;
+    int maxdiffidx = 0;
+
+    cerr << "forward latency = " << cq.getLatency() << ", inverse latency = " 
+	 << cqi.getLatency() << ", total = " << latency << endl;
+
+    timeval tv;
+    (void)gettimeofday(&tv, 0);
+
+    while (inframe < sfinfo.frames) {
+
+        int count = -1;
+	
+	if ((count = sf_readf_float(sndfile, fbuf, ibs)) < 0) {
+	    break;
+	}
+
+	vector<double> cqin;
+	for (int i = 0; i < count; ++i) {
+	    double v = fbuf[i * channels];
+	    if (channels > 1) {
+		for (int c = 1; c < channels; ++c) {
+		    v += fbuf[i * channels + c];
+		}
+		v /= channels;
+	    }
+	    cqin.push_back(v);
+	}
+
+	if (doDiff) {
+	    buffer.insert(buffer.end(), cqin.begin(), cqin.end());
+	}
+	
+	vector<double> cqout = cqi.process(cq.process(cqin));
+
+	for (int i = 0; i < int(cqout.size()); ++i) {
+	    if (cqout[i] > 1.0) cqout[i] = 1.0;
+	    if (cqout[i] < -1.0) cqout[i] = -1.0;
+	}
+
+	if (outframe >= latency) {
+
+	    sf_writef_double(sndfileOut, 
+			     cqout.data(), 
+			     cqout.size());
+
+	} else if (outframe + (int)cqout.size() >= latency) {
+
+	    int offset = latency - outframe;
+	    sf_writef_double(sndfileOut, 
+			     cqout.data() + offset,
+			     cqout.size() - offset);
+	}
+
+	if (doDiff) {
+	    for (int i = 0; i < (int)cqout.size(); ++i) {
+		if (outframe + i >= latency) {
+		    int dframe = outframe + i - latency;
+		    if (dframe >= (int)buffer.size()) cqout[i] = 0;
+		    else cqout[i] -= buffer[dframe];
+		    if (fabs(cqout[i]) > maxdiff &&
+			dframe > sfinfo.samplerate && // ignore first/last sec
+			dframe + sfinfo.samplerate < sfinfo.frames) {
+			maxdiff = fabs(cqout[i]);
+			maxdiffidx = dframe;
+		    }
+		}
+	    }
+	    
+	    if (outframe >= latency) {
+
+		sf_writef_double(sndDiffFile, 
+				 cqout.data(), 
+				 cqout.size());
+
+	    } else if (outframe + (int)cqout.size() >= latency) {
+
+		int offset = latency - outframe;
+		sf_writef_double(sndDiffFile, 
+				 cqout.data() + offset,
+				 cqout.size() - offset);
+	    }
+	}
+
+	inframe += count;
+	outframe += cqout.size();
+    }
+
+    vector<double> r = cqi.process(cq.getRemainingOutput());
+    vector<double> r2 = cqi.getRemainingOutput();
+
+    r.insert(r.end(), r2.begin(), r2.end());
+
+    for (int i = 0; i < int(r.size()); ++i) {
+	if (r[i] > 1.0) r[i] = 1.0;
+	if (r[i] < -1.0) r[i] = -1.0;
+    }
+
+    sf_writef_double(sndfileOut, r.data(), r.size());
+    if (doDiff) {
+	for (int i = 0; i < (int)r.size(); ++i) {
+	    if (outframe + i >= latency) {
+		int dframe = outframe + i - latency;
+		if (dframe >= (int)buffer.size()) r[i] = 0;
+		else r[i] -= buffer[dframe];
+		if (fabs(r[i]) > maxdiff &&
+		    dframe > sfinfo.samplerate && // ignore first/last sec
+		    dframe + sfinfo.samplerate < sfinfo.frames) {
+		    maxdiff = fabs(r[i]);
+		    maxdiffidx = dframe;
+		}
+	    }
+	}
+	sf_writef_double(sndDiffFile, r.data(), r.size());
+    }
+    outframe += r.size();
+
+    sf_close(sndfile);
+    sf_close(sndfileOut);
+
+    if (doDiff) {
+	sf_close(sndDiffFile);
+    }
+
+    cerr << "in: " << inframe << ", out: " << outframe - latency << endl;
+
+    if (doDiff) {
+	double db = 10 * log10(maxdiff);
+	cerr << "max diff [excluding first and last second of audio] is "
+	     << maxdiff << " (" << db << " dBFS)"
+	     << " at sample index " << maxdiffidx << endl;
+    }
+    
+    timeval etv;
+    (void)gettimeofday(&etv, 0);
+        
+    etv.tv_sec -= tv.tv_sec;
+    if (etv.tv_usec < tv.tv_usec) {
+	etv.tv_usec += 1000000;
+	etv.tv_sec -= 1;
+    }
+    etv.tv_usec -= tv.tv_usec;
+        
+    double sec = double(etv.tv_sec) + (double(etv.tv_usec) / 1000000.0);
+    cerr << "elapsed time (not counting init): " << sec << " sec, frames/sec at input: " << inframe/sec << endl;
+
+    return 0;
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test.cpp	Thu May 15 12:24:11 2014 +0100
@@ -0,0 +1,76 @@
+/*
+    Constant-Q library
+    Copyright (c) 2013-2014 Queen Mary, University of London
+
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+    CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    Except as contained in this notice, the names of the Centre for
+    Digital Music; Queen Mary, University of London; and Chris Cannam
+    shall not be used in advertising or otherwise to promote the sale,
+    use or other dealings in this Software without prior written
+    authorization.
+*/
+
+#include "CQSpectrogram.h"
+
+#include <iostream>
+#include <vector>
+
+using std::vector;
+using std::cerr;
+using std::cout;
+using std::endl;
+
+#include <cstdio>
+
+int main(int argc, char **argv)
+{
+    vector<double> in;
+
+    for (int i = 0; i < 64; ++i) {
+//	if (i == 0) in.push_back(1);
+//	else in.push_back(0);
+	in.push_back(sin(i * M_PI / 2.0));
+    }
+
+    CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros);
+
+    vector<vector<double> > out = k.process(in);
+    vector<vector<double> > rest = k.getRemainingOutput();
+
+    out.insert(out.end(), rest.begin(), rest.end());
+
+    cerr << "got " << out.size() << " back (" << out[0].size() << " in each?)" << endl;
+
+    for (int b = 0; b < (int)out.size() / 8; ++b) {
+	printf("\nColumns %d to %d:\n\n", b * 8, b * 8 + 7);
+	for (int j = int(out[0].size()) - 1; j >= 0; --j) {
+	    for (int i = 0; i < 8; ++i) {
+		if (i + b * 8 < (int)out.size()) {
+		    double v = out[i + b * 8][j];
+		    if (v < 0.0001) printf("  0      ");
+		    else printf("  %.4f ", out[i + b * 8][j]);
+		}
+	    }
+	    printf("\n");
+	}
+    }
+}
+