changeset 103:790832f18419

Merge

--- a/Makefile.inc	Thu May 08 09:35:10 2014 +0100
+++ b/Makefile.inc	Tue May 13 13:31:24 2014 +0100
@@ -16,12 +16,14 @@
 LDFLAGS := $(LDFLAGS) 
 PLUGIN_LDFLAGS := $(LDFLAGS) $(PLUGIN_LDFLAGS)
 TEST_LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework
+PF_LDFLAGS := $(LDFLAGS) -lsndfile
 
 PLUGIN	:= cqvamp$(PLUGIN_EXT)
 TEST	:= $(LIB_DIR)/test
+PF	:= $(LIB_DIR)/processfile
 
-LIB_HEADERS := $(LIB_DIR)/CQKernel.h $(LIB_DIR)/ConstantQ.h $(LIB_DIR)/CQInterpolated.h
-LIB_SOURCES := $(LIB_DIR)/CQKernel.cpp $(LIB_DIR)/ConstantQ.cpp $(LIB_DIR)/CQInterpolated.cpp
+LIB_HEADERS := $(LIB_DIR)/CQBase.h $(LIB_DIR)/CQKernel.h $(LIB_DIR)/ConstantQ.h $(LIB_DIR)/CQSpectrogram.h $(LIB_DIR)/CQInverse.h
+LIB_SOURCES := $(LIB_DIR)/CQKernel.cpp $(LIB_DIR)/ConstantQ.cpp $(LIB_DIR)/CQSpectrogram.cpp $(LIB_DIR)/CQInverse.cpp
 
 VAMP_HEADERS := $(VAMP_DIR)/CQVamp.h
 VAMP_SOURCES := $(VAMP_DIR)/CQVamp.cpp $(VAMP_DIR)/libmain.cpp
@@ -30,12 +32,15 @@
 SOURCES	     := $(LIB_SOURCES) $(VAMP_SOURCES)
 OBJECTS	     := $(SOURCES:.cpp=.o)
 
-TEST_SOURCES := $(LIB_DIR)/test.cpp $(SOURCES)
-TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o)
+TEST_SOURCES := $(LIB_DIR)/test.cpp 
+TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o) $(OBJECTS)
+
+PF_SOURCES := $(LIB_DIR)/processfile.cpp
+PF_OBJECTS := $(PF_SOURCES:.cpp=.o) $(OBJECTS)
 
 LIBS	:= $(QMDSP_DIR)/libqm-dsp.a $(VAMPSDK_DIR)/libvamp-sdk.a -lpthread
 
-all: $(PLUGIN) $(TEST)
+all: $(PLUGIN) $(TEST) $(PF)
 
 $(PLUGIN):	$(OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(PLUGIN_LDFLAGS)
@@ -43,24 +48,40 @@
 $(TEST):	$(TEST_OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(TEST_LDFLAGS)
 
+$(PF):	$(PF_OBJECTS)
+	$(CXX) -o $@ $^ $(LIBS) $(PF_LDFLAGS)
+
 clean:		
-	rm -f $(OBJECTS)
+	rm -f $(OBJECTS) $(TEST_OBJECTS) $(PF_OBJECTS)
 
 distclean:	clean
 	rm -f $(PLUGIN)
 
 depend:
-	makedepend -Y -fMakefile.inc $(SOURCES) $(HEADERS)
+	makedepend -Y -fMakefile.inc $(SOURCES) $(TEST_SOURCES) $(PF_SOURCES) $(HEADERS)
 
 # DO NOT DELETE
 
 cpp-qm-dsp/CQKernel.o: cpp-qm-dsp/CQKernel.h
-cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQKernel.h
-vamp/CQVamp.o: vamp/CQVamp.h cpp-qm-dsp/CQInterpolated.h
-vamp/CQVamp.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQKernel.h
-vamp/libmain.o: vamp/CQVamp.h cpp-qm-dsp/CQInterpolated.h
-vamp/libmain.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
 cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQKernel.h
-cpp-qm-dsp/CQInterpolated.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQKernel.h
-vamp/CQVamp.o: cpp-qm-dsp/CQInterpolated.h cpp-qm-dsp/ConstantQ.h
+cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
+cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQBase.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQKernel.h
+vamp/CQVamp.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h
+vamp/CQVamp.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
 vamp/CQVamp.o: cpp-qm-dsp/CQKernel.h
+vamp/libmain.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h
+vamp/libmain.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
+vamp/libmain.o: cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/test.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
+cpp-qm-dsp/test.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/processfile.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
+cpp-qm-dsp/processfile.o: cpp-qm-dsp/CQKernel.h cpp-qm-dsp/CQInverse.h
+cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
+cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
+vamp/CQVamp.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
+vamp/CQVamp.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h

--- a/Makefile.linux	Thu May 08 09:35:10 2014 +0100
+++ b/Makefile.linux	Tue May 13 13:31:24 2014 +0100
@@ -1,6 +1,6 @@
 
-#CFLAGS := -Wall -O3 -fPIC -I../vamp-plugin-sdk/
-CFLAGS := -g -fPIC -I../vamp-plugin-sdk
+CFLAGS := -Wall -O3 -fPIC -I../vamp-plugin-sdk/
+#CFLAGS := -g -fPIC -I../vamp-plugin-sdk
 
 CXXFLAGS := $(CFLAGS)
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQBase.h	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,60 @@
+/* -*- 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 CQBASE_H
+#define CQBASE_H
+
+#include <vector>
+#include <complex>
+
+class CQBase // interface class
+{
+public:
+    typedef std::complex<double> Complex;
+    typedef std::vector<double> RealSequence;
+    typedef std::vector<double> RealColumn;
+    typedef std::vector<Complex> ComplexSequence;
+    typedef std::vector<Complex> ComplexColumn;
+    typedef std::vector<RealColumn> RealBlock;
+    typedef std::vector<ComplexColumn> ComplexBlock;
+
+    virtual double getSampleRate() const = 0;
+    virtual int getBinsPerOctave() const = 0;
+    virtual int getOctaves() const = 0; 
+    virtual int getTotalBins() const = 0;
+    virtual int getColumnHop() const = 0;
+    virtual int getLatency() const = 0;
+    virtual double getMaxFrequency() const = 0;
+    virtual double getMinFrequency() const = 0; // actual min, not that passed to ctor
+    virtual double getBinFrequency(int bin) const = 0;
+};
+
+#endif

--- a/cpp-qm-dsp/CQInterpolated.cpp	Thu May 08 09:35:10 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,247 +0,0 @@
-/* -*- 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 "CQInterpolated.h"
-
-#include <iostream>
-#include <stdexcept>
-
-using std::vector;
-
-using std::cerr;
-using std::endl;
-
-CQInterpolated::CQInterpolated(double sampleRate,
-			       double minFreq, double maxFreq,
-			       int binsPerOctave,
-			       Interpolation interpolation) :
-    m_cq(sampleRate, minFreq, maxFreq, binsPerOctave),
-    m_interpolation(interpolation)
-{
-}
-
-CQInterpolated::~CQInterpolated()
-{
-}
-
-vector<vector<double> >
-CQInterpolated::process(const vector<double> &td)
-{
-    return postProcess(m_cq.process(td), false);
-}
-
-vector<vector<double> >
-CQInterpolated::getRemainingBlocks()
-{
-    return postProcess(m_cq.getRemainingBlocks(), true);
-}
-
-vector<vector<double> >
-CQInterpolated::postProcess(const vector<vector<double> > &cq, bool insist)
-{
-    if (m_interpolation == None) {
-	return cq;
-    }
-
-    int width = cq.size();
-
-    for (int i = 0; i < width; ++i) {
-	m_buffer.push_back(cq[i]);
-    }
-    
-    if (m_interpolation == Hold) {
-	return fetchHold(insist);
-    } else {
-	return fetchLinear(insist);
-    }
-}
-
-vector<vector<double> >
-CQInterpolated::fetchHold(bool)
-{
-    Grid out;
-    
-    int width = m_buffer.size();
-    int height = getTotalBins();
-
-    for (int i = 0; i < width; ++i) {
-	
-	vector<double> col = m_buffer[i];
-
-	int thisHeight = col.size();
-	int prevHeight = m_prevColumn.size();
-
-	for (int j = thisHeight; j < height; ++j) {
-	    if (j < prevHeight) {
-		col.push_back(m_prevColumn[j]);
-	    } else {
-		col.push_back(0.0);
-	    }
-	}
-
-	m_prevColumn = col;
-	out.push_back(col);
-    }
-
-    m_buffer.clear();
-
-    return out;
-}
-
-vector<vector<double> >
-CQInterpolated::fetchLinear(bool insist)
-{
-    Grid out;
-
-    //!!! This is surprisingly messy. I must be missing something.
-
-    // We can only return any data when we have at least one column
-    // that has the full height in the buffer, that is not the first
-    // column.
-    //
-    // If the first col has full height, and there is another one
-    // later that also does, then we can interpolate between those, up
-    // to but not including the second full height column.  Then we
-    // drop and return the columns we interpolated, leaving the second
-    // full-height col as the first col in the buffer. And repeat as
-    // long as enough columns are available.
-    //
-    // If the first col does not have full height, then (so long as
-    // we're following the logic above) we must simply have not yet
-    // reached the first full-height column in the CQ output, and we
-    // can interpolate nothing.
-    
-    int width = m_buffer.size();
-    int height = getTotalBins();
-
-    if (width == 0) return out;
-    
-    int firstFullHeight = -1;
-    int secondFullHeight = -1;
-
-    for (int i = 0; i < width; ++i) {
-	if ((int)m_buffer[i].size() == height) {
-	    if (firstFullHeight == -1) {
-		firstFullHeight = i;
-	    } else if (secondFullHeight == -1) {
-		secondFullHeight = i;
-		break;
-	    }
-	}
-    }
-
-    if (firstFullHeight < 0) {
-	if (insist) {
-	    out = m_buffer;
-	    m_buffer.clear();
-	    return out;
-	} else {
-	    return out;
-	}
-    } else if (firstFullHeight > 0) {
-	// can interpolate nothing, stash up to first full height & recurse
-	out = Grid(m_buffer.begin(), m_buffer.begin() + firstFullHeight);
-	m_buffer = Grid(m_buffer.begin() + firstFullHeight, m_buffer.end());
-	Grid more = fetchLinear(insist);
-	out.insert(out.end(), more.begin(), more.end());
-	return out;
-    } else if (secondFullHeight < 0) {
-	// firstFullHeight == 0, but there is no second full height --
-	// wait for it unless insist flag is set
-	if (insist) {
-	    out = m_buffer;
-	    m_buffer.clear();
-	    return out;
-	} else {
-	    return out;
-	}
-    } else {
-	// firstFullHeight == 0 and secondFullHeight also valid. Can interpolate
-	out = linearInterpolated(m_buffer, 0, secondFullHeight);
-	m_buffer = Grid(m_buffer.begin() + secondFullHeight, m_buffer.end());
-	Grid more = fetchLinear(insist);
-	out.insert(out.end(), more.begin(), more.end());
-	return out;
-    }
-}
-
-vector<vector<double> >
-CQInterpolated::linearInterpolated(const Grid &g, int x0, int x1)
-{
-    // g must be a grid with full-height columns at x0 and x1
-
-    if (x0 >= x1) {
-	throw std::logic_error("x0 >= x1");
-    }
-    if (x1 >= (int)g.size()) {
-	throw std::logic_error("x1 >= g.size()");
-    }
-    if (g[x0].size() != g[x1].size()) {
-	throw std::logic_error("x0 and x1 are not the same height");
-    }
-
-    int height = g[x0].size();
-    int width = x1 - x0;
-
-    Grid out(g.begin() + x0, g.begin() + x1);
-
-    for (int y = 0; y < height; ++y) {
-
-	int spacing = width;
-	for (int i = 1; i < width; ++i) {
-	    int thisHeight = g[x0 + i].size();
-	    if (thisHeight > height) {
-		throw std::logic_error("First column not full-height");
-	    }
-	    if (thisHeight > y) {
-		spacing = i;
-		break;
-	    }
-	}
-
-	if (spacing < 2) continue;
-
-	for (int i = 0; i + spacing <= width; i += spacing) {
-	    for (int j = 1; j < spacing; ++j) {
-		double proportion = double(j)/double(spacing);
-		double interpolated = 
-		    g[x0 + i][y] * (1.0 - proportion) +
-		    g[x0 + i + spacing][y] * proportion;
-		out[i + j].push_back(interpolated);
-	    }
-	}
-    }
-
-    return out;
-}
-
-
-	

--- a/cpp-qm-dsp/CQInterpolated.h	Thu May 08 09:35:10 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-/* -*- 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 CQ_INTERPOLATED_H
-#define CQ_INTERPOLATED_H
-
-#include <vector>
-#include "ConstantQ.h"
-
-class CQInterpolated
-{
-public:
-    enum Interpolation {
-	None, // leave empty cells empty
-	Hold, // repeat prior cell
-	Linear, // linear interpolation between consecutive time cells
-    };
-
-    CQInterpolated(double sampleRate,
-		   double minFreq, double maxFreq,
-		   int binsPerOctave,
-		   Interpolation interpolation);
-    ~CQInterpolated();
-
-    double getSampleRate() const { return m_cq.getSampleRate(); }
-    int getBinsPerOctave() const { return m_cq.getBinsPerOctave(); }
-    int getOctaves() const { return m_cq.getOctaves(); }
-    int getTotalBins() const { return m_cq.getTotalBins(); }
-    int getColumnHop() const { return m_cq.getColumnHop(); }
-    int getLatency() const { return m_cq.getLatency(); } 
-    double getMaxFrequency() const { return m_cq.getMaxFrequency(); }
-    double getMinFrequency() const { return m_cq.getMinFrequency(); }
-    double getBinFrequency(int bin) const { return m_cq.getBinFrequency(bin); }
-
-    std::vector<std::vector<double> > process(const std::vector<double> &);
-    std::vector<std::vector<double> > getRemainingBlocks();
-
-private:
-    ConstantQ m_cq;
-    Interpolation m_interpolation;
-
-    typedef std::vector<std::vector<double> > Grid;
-    Grid m_buffer;
-    Grid postProcess(const Grid &, bool insist);
-    Grid fetchHold(bool insist);
-    Grid fetchLinear(bool insist);
-    Grid linearInterpolated(const Grid &, int, int);
-    std::vector<double> m_prevColumn;
-};
-
-#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQInverse.cpp	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,409 @@
+/* -*- 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 "CQInverse.h"
+
+#include "dsp/rateconversion/Resampler.h"
+#include "maths/MathUtilities.h"
+#include "dsp/transforms/FFT.h"
+
+#include <algorithm>
+#include <iostream>
+#include <stdexcept>
+
+using std::vector;
+using std::cerr;
+using std::endl;
+
+//#define DEBUG_CQ 1
+
+CQInverse::CQInverse(double sampleRate,
+                     double minFreq,
+                     double maxFreq,
+                     int binsPerOctave) :
+    m_sampleRate(sampleRate),
+    m_maxFrequency(maxFreq),
+    m_minFrequency(minFreq),
+    m_binsPerOctave(binsPerOctave),
+    m_fft(0)
+{
+    if (minFreq <= 0.0 || maxFreq <= 0.0) {
+        throw std::invalid_argument("Frequency extents must be positive");
+    }
+
+    initialise();
+}
+
+CQInverse::~CQInverse()
+{
+    delete m_fft;
+    for (int i = 0; i < (int)m_upsamplers.size(); ++i) {
+        delete m_upsamplers[i];
+    }
+    delete m_kernel;
+}
+
+double
+CQInverse::getMinFrequency() const
+{
+    return m_p.minFrequency / pow(2.0, m_octaves - 1);
+}
+
+double
+CQInverse::getBinFrequency(int bin) const
+{
+    return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave()));
+}
+
+void
+CQInverse::initialise()
+{
+    m_octaves = int(ceil(log2(m_maxFrequency / m_minFrequency)));
+    m_kernel = new CQKernel(m_sampleRate, m_maxFrequency, m_binsPerOctave);
+    m_p = m_kernel->getProperties();
+    
+    // Use exact powers of two for resampling rates. They don't have
+    // to be related to our actual samplerate: the resampler only
+    // cares about the ratio, but it only accepts integer source and
+    // target rates, and if we start from the actual samplerate we
+    // risk getting non-integer rates for lower octaves
+
+    int sourceRate = pow(2, m_octaves);
+    vector<int> latencies;
+
+    // top octave, no resampling
+    latencies.push_back(0);
+    m_upsamplers.push_back(0);
+
+    for (int i = 1; i < m_octaves; ++i) {
+
+        int factor = pow(2, i);
+
+        Resampler *r = new Resampler
+            (sourceRate / factor, sourceRate, 60, 0.02);
+
+#ifdef DEBUG_CQ
+        cerr << "inverse: octave " << i << ": resample from " << sourceRate/factor << " to " << sourceRate << endl;
+#endif
+
+	// See ConstantQ.cpp for discussion on latency -- output
+	// latency here is at target rate which, this way around, is
+	// what we want
+
+        latencies.push_back(r->getLatency());
+        m_upsamplers.push_back(r);
+    }
+
+    // additionally we will have fftHop latency at individual octave
+    // rate (before upsampling) for the overlap-add in each octave
+    for (int i = 0; i < m_octaves; ++i) {
+        latencies[i] += m_p.fftHop * pow(2, i);
+    }
+
+    // Now reverse the drop adjustment made in ConstantQ to align the
+    // atom centres across different octaves (but this time at output
+    // sample rate)
+
+    int emptyHops = m_p.firstCentre / m_p.atomSpacing;
+
+    vector<int> pushes;
+    for (int i = 0; i < m_octaves; ++i) {
+	int factor = pow(2, i);
+	int pushHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops;
+	int push = ((pushHops * m_p.fftHop) * factor) / m_p.atomsPerFrame;
+	pushes.push_back(push);
+    }
+
+    int maxLatLessPush = 0;
+    for (int i = 0; i < m_octaves; ++i) {
+	int latLessPush = latencies[i] - pushes[i];
+	if (latLessPush > maxLatLessPush) maxLatLessPush = latLessPush;
+    }
+
+    int totalLatency = maxLatLessPush + 10;
+    if (totalLatency < 0) totalLatency = 0;
+
+    m_outputLatency = totalLatency + m_p.firstCentre * pow(2, m_octaves-1);
+
+#ifdef DEBUG_CQ
+    cerr << "totalLatency = " << totalLatency << ", m_outputLatency = " << m_outputLatency << endl;
+#endif
+
+    for (int i = 0; i < m_octaves; ++i) {
+
+	// Calculate the difference between the total latency applied
+	// across all octaves, and the existing latency due to the
+	// upsampler for this octave.
+
+        int latencyPadding = totalLatency - latencies[i] + pushes[i];
+
+#ifdef DEBUG_CQ
+        cerr << "octave " << i << ": push " << pushes[i] << ", resampler latency inc overlap space " << latencies[i] << ", latencyPadding = " << latencyPadding << " (/factor = " << latencyPadding / pow(2, i) << ")" << endl;
+#endif
+
+        m_buffers.push_back(RealSequence(latencyPadding, 0.0));
+    }
+
+    for (int i = 0; i < m_octaves; ++i) {
+        // Fixed-size buffer for IFFT overlap-add
+        m_olaBufs.push_back(RealSequence(m_p.fftSize, 0.0));
+    }
+
+    m_fft = new FFTReal(m_p.fftSize);
+}
+
+CQInverse::RealSequence
+CQInverse::process(const ComplexBlock &block)
+{
+    // The input data is of the form produced by ConstantQ::process --
+    // an unknown number N of columns of varying height. We assert
+    // that N is a multiple of atomsPerFrame * 2^(octaves-1), as must
+    // be the case for data that came directly from our ConstantQ
+    // implementation.
+
+    int widthProvided = block.size();
+
+    if (widthProvided == 0) {
+        return drawFromBuffers();
+    }
+
+    int blockWidth = m_p.atomsPerFrame * int(pow(2, m_octaves - 1));
+
+    if (widthProvided % blockWidth != 0) {
+        cerr << "ERROR: CQInverse::process: Input block size ("
+             << widthProvided
+             << ") must be a multiple of processing block width "
+             << "(atoms-per-frame * 2^(octaves-1) = "
+             << m_p.atomsPerFrame << " * 2^(" << m_octaves << "-1) = "
+             << blockWidth << ")" << endl;
+        throw std::invalid_argument
+            ("Input block size must be a multiple of processing block width");
+    }
+
+    // Procedure:
+    // 
+    // 1. Slice the list of columns into a set of lists of columns,
+    // one per octave, each of width N / (2^octave-1) and height
+    // binsPerOctave, containing the values present in that octave
+    //
+    // 2. Group each octave list by atomsPerFrame columns at a time,
+    // and stack these so as to achieve a list, for each octave, of
+    // taller columns of height binsPerOctave * atomsPerFrame
+    //
+    // 3. For each taller column, take the product with the inverse CQ
+    // kernel (which is the conjugate of the forward kernel) and
+    // perform an inverse FFT
+    //
+    // 4. Overlap-add each octave's resynthesised blocks (unwindowed)
+    //
+    // 5. Resample each octave's overlap-add stream to the original
+    // rate
+    //
+    // 6. Sum the resampled streams and return
+    
+    for (int i = 0; i < m_octaves; ++i) {
+        
+        // Step 1
+
+        ComplexBlock oct;
+
+        for (int j = 0; j < widthProvided; ++j) {
+            int h = block[j].size();
+            if (h < m_binsPerOctave * (i+1)) {
+                continue;
+            }
+            ComplexColumn col(block[j].begin() + m_binsPerOctave * i,
+                              block[j].begin() + m_binsPerOctave * (i+1));
+            oct.push_back(col);
+        }
+
+        // Steps 2, 3, 4, 5
+        processOctave(i, oct);
+    }
+    
+    // Step 6
+    return drawFromBuffers();
+}
+
+CQInverse::RealSequence
+CQInverse::drawFromBuffers()
+{
+    // 6. Sum the resampled streams and return
+
+    int available = 0;
+
+    for (int i = 0; i < m_octaves; ++i) {
+        if (i == 0 || int(m_buffers[i].size()) < available) {
+            available = m_buffers[i].size();
+        }
+    }
+
+    RealSequence result(available, 0);
+
+    if (available == 0) {
+        return result;
+    }
+
+    for (int i = 0; i < m_octaves; ++i) {
+        for (int j = 0; j < available; ++j) {
+            result[j] += m_buffers[i][j];
+        }
+        m_buffers[i] = RealSequence(m_buffers[i].begin() + available,
+                                    m_buffers[i].end());
+    }
+
+    return result;
+}
+
+CQInverse::RealSequence
+CQInverse::getRemainingOutput()
+{
+    for (int i = 0; i < (m_p.fftSize - m_p.fftHop) / m_p.fftHop; ++i) {
+        for (int j = 0; j < m_octaves; ++j) {
+            overlapAddAndResample(j, RealSequence(m_olaBufs[j].size(), 0));
+        }
+    }
+
+    return drawFromBuffers();
+}
+
+void
+CQInverse::processOctave(int octave, const ComplexBlock &columns)
+{
+    // 2. Group each octave list by atomsPerFrame columns at a time,
+    // and stack these so as to achieve a list, for each octave, of
+    // taller columns of height binsPerOctave * atomsPerFrame
+
+    int ncols = columns.size();
+
+    if (ncols % m_p.atomsPerFrame != 0) {
+        cerr << "ERROR: CQInverse::process: Number of columns ("
+             << ncols
+             << ") in octave " << octave
+             << " must be a multiple of atoms-per-frame ("
+             << m_p.atomsPerFrame << ")" << endl;
+        throw std::invalid_argument
+            ("Columns in octave must be a multiple of atoms per frame");
+    }
+
+    for (int i = 0; i < ncols; i += m_p.atomsPerFrame) {
+
+        ComplexColumn tallcol;
+        for (int b = 0; b < m_binsPerOctave; ++b) {
+            for (int a = 0; a < m_p.atomsPerFrame; ++a) {
+                tallcol.push_back(columns[i + a][m_binsPerOctave - b - 1]);
+            }
+        }
+        
+        processOctaveColumn(octave, tallcol);
+    }
+}
+
+void
+CQInverse::processOctaveColumn(int octave, const ComplexColumn &column)
+{
+    // 3. For each taller column, take the product with the inverse CQ
+    // kernel (which is the conjugate of the forward kernel) and
+    // perform an inverse FFT
+
+    if ((int)column.size() != m_p.atomsPerFrame * m_binsPerOctave) {
+        cerr << "ERROR: CQInverse::processOctaveColumn: Height of column ("
+             << column.size() << ") in octave " << octave
+             << " must be atoms-per-frame * bins-per-octave ("
+             << m_p.atomsPerFrame << " * " << m_binsPerOctave << " = "
+             << m_p.atomsPerFrame * m_binsPerOctave << ")" << endl;
+        throw std::invalid_argument
+            ("Column height must match atoms-per-frame * bins-per-octave");
+    }
+
+    ComplexSequence transformed = m_kernel->processInverse(column);
+
+    int halfLen = m_p.fftSize/2 + 1;
+
+    RealSequence ri(halfLen, 0);
+    RealSequence ii(halfLen, 0);
+
+    for (int i = 0; i < halfLen; ++i) {
+        ri[i] = transformed[i].real();
+        ii[i] = transformed[i].imag();
+    }
+
+    RealSequence timeDomain(m_p.fftSize, 0);
+
+    m_fft->inverse(ri.data(), ii.data(), timeDomain.data());
+
+    overlapAddAndResample(octave, timeDomain);
+}
+
+void
+CQInverse::overlapAddAndResample(int octave, const RealSequence &seq)
+{
+    // 4. Overlap-add each octave's resynthesised blocks (unwindowed)
+    //
+    // and
+    //
+    // 5. Resample each octave's overlap-add stream to the original
+    // rate
+
+    if (seq.size() != m_olaBufs[octave].size()) {
+        cerr << "ERROR: CQInverse::overlapAdd: input sequence length ("
+             << seq.size() << ") is expected to match OLA buffer size ("
+             << m_olaBufs[octave].size() << ")" << endl;
+        throw std::invalid_argument
+            ("Input sequence length should match OLA buffer size");
+    }
+
+    RealSequence toResample(m_olaBufs[octave].begin(),
+                            m_olaBufs[octave].begin() + m_p.fftHop);
+
+    RealSequence resampled = 
+        octave > 0 ?
+        m_upsamplers[octave]->process(toResample.data(), toResample.size()) :
+        toResample;
+
+    m_buffers[octave].insert(m_buffers[octave].end(),
+                             resampled.begin(),
+                             resampled.end());
+    
+    m_olaBufs[octave] = RealSequence(m_olaBufs[octave].begin() + m_p.fftHop,
+                                     m_olaBufs[octave].end());
+    
+    RealSequence pad(m_p.fftHop, 0);
+
+    m_olaBufs[octave].insert(m_olaBufs[octave].end(),
+                             pad.begin(),
+                             pad.end());
+
+    for (int i = 0; i < m_p.fftSize; ++i) {
+        m_olaBufs[octave][i] += seq[i];
+    }
+}
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQInverse.h	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,91 @@
+/* -*- 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 CQINVERSE_H
+#define CQINVERSE_H
+
+#include "CQBase.h"
+#include "CQKernel.h"
+
+class Resampler;
+class FFTReal;
+
+class CQInverse : public CQBase
+{
+public:
+    CQInverse(double sampleRate,
+	      double minFreq, double maxFreq,
+	      int binsPerOctave);
+    virtual ~CQInverse();
+
+    virtual double getSampleRate() const { return m_sampleRate; }
+    virtual int getBinsPerOctave() const { return m_binsPerOctave; }
+    virtual int getOctaves() const { return m_octaves; }
+    virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; }
+    virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; }
+    virtual int getLatency() const { return m_outputLatency; } 
+    virtual double getMaxFrequency() const { return m_p.maxFrequency; }
+    virtual double getMinFrequency() const; // actual min, not that passed to ctor
+    virtual double getBinFrequency(int bin) const;
+
+    // Input is the format produced by ConstantQ class,
+    // i.e. uninterpolated complex, not the real-valued stuff produced
+    // by CQSpectrogram
+
+    RealSequence process(const ComplexBlock &);
+    RealSequence getRemainingOutput();
+
+private:
+    double m_sampleRate;
+    double m_maxFrequency;
+    double m_minFrequency;
+    int m_binsPerOctave;
+    int m_octaves;
+
+    CQKernel *m_kernel;
+    CQKernel::Properties m_p;
+
+    std::vector<Resampler *> m_upsamplers;
+    std::vector<RealSequence> m_buffers;
+    std::vector<RealSequence> m_olaBufs; // fixed-length, for overlap-add
+    
+    int m_outputLatency;
+
+    FFTReal *m_fft;
+    
+    void initialise();
+    void processOctave(int octave, const ComplexBlock &block);
+    void processOctaveColumn(int octave, const ComplexColumn &column);
+    void overlapAddAndResample(int octave, const RealSequence &);
+    RealSequence drawFromBuffers();
+};
+
+#endif

--- a/cpp-qm-dsp/CQKernel.cpp	Thu May 08 09:35:10 2014 +0100
+++ b/cpp-qm-dsp/CQKernel.cpp	Tue May 13 13:31:24 2014 +0100
@@ -241,9 +241,10 @@
     
     cerr << "weight = " << weight << endl;
 
-    // apply normalisation weight, make sparse, and store conjugates
-    // (our multiplication order means we will effectively be using
-    // the adjoint or conjugate transpose of the kernel matrix)
+    // apply normalisation weight, make sparse, and store conjugate
+    // (we use the adjoint or conjugate transpose of the kernel matrix
+    // for the forward transform, the plain kernel for the inverse
+    // which we expect to be less common)
 
     KernelMatrix sk;
 
@@ -274,14 +275,14 @@
 }
 
 vector<C>
-CQKernel::process(const vector<C> &cv)
+CQKernel::processForward(const vector<C> &cv)
 {
-    // matrix multiply m_kernel.data by in, converting in to complex
-    // as we go
+    // straightforward matrix multiply (taking into account m_kernel's
+    // slightly-sparse representation)
 
     int nrows = m_p.binsPerOctave * m_p.atomsPerFrame;
 
-    vector<C> rv(nrows, C(0, 0));
+    vector<C> rv(nrows, C());
 
     for (int i = 0; i < nrows; ++i) {
         for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) {
@@ -292,4 +293,28 @@
     return rv;
 }
 
+vector<C>
+CQKernel::processInverse(const vector<C> &cv)
+{
+    // matrix multiply by conjugate transpose of m_kernel. This is
+    // actually the original kernel as calculated, we just stored the
+    // conjugate-transpose of the kernel because we expect to be doing
+    // more forward transforms than inverse ones.
 
+    int ncols = m_p.binsPerOctave * m_p.atomsPerFrame;
+    int nrows = m_p.fftSize;
+
+    vector<C> rv(nrows, C());
+
+    for (int j = 0; j < ncols; ++j) {
+        int i0 = m_kernel.origin[j];
+        int i1 = i0 + m_kernel.data[j].size();
+        for (int i = i0; i < i1; ++i) {
+            rv[i] += cv[j] * conj(m_kernel.data[j][i - i0]);
+        }
+    }
+
+    return rv;
+}
+
+

--- a/cpp-qm-dsp/CQKernel.h	Thu May 08 09:35:10 2014 +0100
+++ b/cpp-qm-dsp/CQKernel.h	Tue May 13 13:31:24 2014 +0100
@@ -59,7 +59,10 @@
 
     Properties getProperties() const { return m_p; }
 
-    std::vector<std::complex<double> > process
+    std::vector<std::complex<double> > processForward
+        (const std::vector<std::complex<double> > &);
+
+    std::vector<std::complex<double> > processInverse
         (const std::vector<std::complex<double> > &);
 
 private:

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQSpectrogram.cpp	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,263 @@
+/* -*- 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 "CQSpectrogram.h"
+
+#include <iostream>
+#include <stdexcept>
+
+using std::cerr;
+using std::endl;
+
+CQSpectrogram::CQSpectrogram(double sampleRate,
+			       double minFreq, double maxFreq,
+			       int binsPerOctave,
+			       Interpolation interpolation) :
+    m_cq(sampleRate, minFreq, maxFreq, binsPerOctave),
+    m_interpolation(interpolation)
+{
+}
+
+CQSpectrogram::~CQSpectrogram()
+{
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::process(const RealSequence &td)
+{
+    return postProcess(m_cq.process(td), false);
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::getRemainingOutput()
+{
+    return postProcess(m_cq.getRemainingOutput(), true);
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::postProcess(const ComplexBlock &cq, bool insist)
+{
+    int width = cq.size();
+
+    // convert to magnitudes
+    RealBlock spec;
+    for (int i = 0; i < width; ++i) {
+        int height = cq[i].size();
+        RealColumn col(height, 0);
+        for (int j = 0; j < height; ++j) {
+            col[j] = abs(cq[i][j]);
+        }
+        spec.push_back(col);
+    }
+
+    if (m_interpolation == InterpolateZeros) {
+        for (int i = 0; i < width; ++i) {
+            int sh = spec[i].size();
+            int fh = getTotalBins();
+            for (int j = sh; j < fh; ++j) {
+                spec[i].push_back(0);
+            }
+        }
+	return spec;
+    }
+
+    for (int i = 0; i < width; ++i) {
+	m_buffer.push_back(spec[i]);
+    }
+    
+    if (m_interpolation == InterpolateHold) {
+	return fetchHold(insist);
+    } else {
+	return fetchLinear(insist);
+    }
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::fetchHold(bool)
+{
+    RealBlock out;
+    
+    int width = m_buffer.size();
+    int height = getTotalBins();
+
+    for (int i = 0; i < width; ++i) {
+	
+	RealColumn col = m_buffer[i];
+
+	int thisHeight = col.size();
+	int prevHeight = m_prevColumn.size();
+
+	for (int j = thisHeight; j < height; ++j) {
+	    if (j < prevHeight) {
+		col.push_back(m_prevColumn[j]);
+	    } else {
+		col.push_back(0.0);
+	    }
+	}
+
+	m_prevColumn = col;
+	out.push_back(col);
+    }
+
+    m_buffer.clear();
+
+    return out;
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::fetchLinear(bool insist)
+{
+    RealBlock out;
+
+    //!!! This is surprisingly messy. I must be missing something.
+
+    // We can only return any data when we have at least one column
+    // that has the full height in the buffer, that is not the first
+    // column.
+    //
+    // If the first col has full height, and there is another one
+    // later that also does, then we can interpolate between those, up
+    // to but not including the second full height column.  Then we
+    // drop and return the columns we interpolated, leaving the second
+    // full-height col as the first col in the buffer. And repeat as
+    // long as enough columns are available.
+    //
+    // If the first col does not have full height, then (so long as
+    // we're following the logic above) we must simply have not yet
+    // reached the first full-height column in the CQ output, and we
+    // can interpolate nothing.
+    
+    int width = m_buffer.size();
+    int height = getTotalBins();
+
+    if (width == 0) return out;
+    
+    int firstFullHeight = -1;
+    int secondFullHeight = -1;
+
+    for (int i = 0; i < width; ++i) {
+	if ((int)m_buffer[i].size() == height) {
+	    if (firstFullHeight == -1) {
+		firstFullHeight = i;
+	    } else if (secondFullHeight == -1) {
+		secondFullHeight = i;
+		break;
+	    }
+	}
+    }
+
+    if (firstFullHeight < 0) {
+	if (insist) {
+	    out = m_buffer;
+	    m_buffer.clear();
+	    return out;
+	} else {
+	    return out;
+	}
+    } else if (firstFullHeight > 0) {
+	// can interpolate nothing, stash up to first full height & recurse
+	out = RealBlock(m_buffer.begin(), m_buffer.begin() + firstFullHeight);
+	m_buffer = RealBlock(m_buffer.begin() + firstFullHeight, m_buffer.end());
+	RealBlock more = fetchLinear(insist);
+	out.insert(out.end(), more.begin(), more.end());
+	return out;
+    } else if (secondFullHeight < 0) {
+	// firstFullHeight == 0, but there is no second full height --
+	// wait for it unless insist flag is set
+	if (insist) {
+	    out = m_buffer;
+	    m_buffer.clear();
+	    return out;
+	} else {
+	    return out;
+	}
+    } else {
+	// firstFullHeight == 0 and secondFullHeight also valid. Can interpolate
+	out = linearInterpolated(m_buffer, 0, secondFullHeight);
+	m_buffer = RealBlock(m_buffer.begin() + secondFullHeight, m_buffer.end());
+	RealBlock more = fetchLinear(insist);
+	out.insert(out.end(), more.begin(), more.end());
+	return out;
+    }
+}
+
+CQSpectrogram::RealBlock
+CQSpectrogram::linearInterpolated(const RealBlock &g, int x0, int x1)
+{
+    // g must be a grid with full-height columns at x0 and x1
+
+    if (x0 >= x1) {
+	throw std::logic_error("x0 >= x1");
+    }
+    if (x1 >= (int)g.size()) {
+	throw std::logic_error("x1 >= g.size()");
+    }
+    if (g[x0].size() != g[x1].size()) {
+	throw std::logic_error("x0 and x1 are not the same height");
+    }
+
+    int height = g[x0].size();
+    int width = x1 - x0;
+
+    RealBlock out(g.begin() + x0, g.begin() + x1);
+
+    for (int y = 0; y < height; ++y) {
+
+	int spacing = width;
+	for (int i = 1; i < width; ++i) {
+	    int thisHeight = g[x0 + i].size();
+	    if (thisHeight > height) {
+		throw std::logic_error("First column not full-height");
+	    }
+	    if (thisHeight > y) {
+		spacing = i;
+		break;
+	    }
+	}
+
+	if (spacing < 2) continue;
+
+	for (int i = 0; i + spacing <= width; i += spacing) {
+	    for (int j = 1; j < spacing; ++j) {
+		double proportion = double(j)/double(spacing);
+		double interpolated = 
+		    g[x0 + i][y] * (1.0 - proportion) +
+		    g[x0 + i + spacing][y] * proportion;
+		out[i + j].push_back(interpolated);
+	    }
+	}
+    }
+
+    return out;
+}
+
+
+	

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQSpectrogram.h	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,77 @@
+/* -*- 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 CQSPECTROGRAM_H
+#define CQSPECTROGRAM_H
+
+#include "ConstantQ.h"
+
+class CQSpectrogram : public CQBase
+{
+public:
+    enum Interpolation {
+	InterpolateZeros, // leave empty cells as zero
+	InterpolateHold, // repeat prior cell
+	InterpolateLinear, // linear interpolation between consecutive time cells
+    };
+
+    CQSpectrogram(double sampleRate,
+                  double minFreq, double maxFreq,
+                  int binsPerOctave,
+                  Interpolation interpolation);
+    virtual ~CQSpectrogram();
+
+    virtual double getSampleRate() const { return m_cq.getSampleRate(); }
+    virtual int getBinsPerOctave() const { return m_cq.getBinsPerOctave(); }
+    virtual int getOctaves() const { return m_cq.getOctaves(); }
+    virtual int getTotalBins() const { return m_cq.getTotalBins(); }
+    virtual int getColumnHop() const { return m_cq.getColumnHop(); }
+    virtual int getLatency() const { return m_cq.getLatency(); } 
+    virtual double getMaxFrequency() const { return m_cq.getMaxFrequency(); }
+    virtual double getMinFrequency() const { return m_cq.getMinFrequency(); }
+    virtual double getBinFrequency(int bin) const { return m_cq.getBinFrequency(bin); }
+
+    RealBlock process(const RealSequence &);
+    RealBlock getRemainingOutput();
+
+private:
+    ConstantQ m_cq;
+    Interpolation m_interpolation;
+
+    RealBlock m_buffer;
+    RealBlock postProcess(const ComplexBlock &, bool insist);
+    RealBlock fetchHold(bool insist);
+    RealBlock fetchLinear(bool insist);
+    RealBlock linearInterpolated(const RealBlock &, int, int);
+    RealColumn m_prevColumn;
+};
+
+#endif

--- a/cpp-qm-dsp/ConstantQ.cpp	Thu May 08 09:35:10 2014 +0100
+++ b/cpp-qm-dsp/ConstantQ.cpp	Tue May 13 13:31:24 2014 +0100
@@ -38,16 +38,14 @@
 #include "dsp/transforms/FFT.h"
 
 #include <algorithm>
-#include <complex>
 #include <iostream>
 #include <stdexcept>
 
 using std::vector;
-using std::complex;
 using std::cerr;
 using std::endl;
 
-typedef std::complex<double> C;
+//#define DEBUG_CQ 1
 
 ConstantQ::ConstantQ(double sampleRate,
                      double minFreq,
@@ -114,6 +112,10 @@
         Resampler *r = new Resampler
             (sourceRate, sourceRate / factor, 60, 0.02);
 
+#ifdef DEBUG_CQ
+        cerr << "forward: octave " << i << ": resample from " << sourceRate << " to " << sourceRate / factor << endl;
+#endif
+
         // We need to adapt the latencies so as to get the first input
         // sample to be aligned, in time, at the decimator output
         // across all octaves.
@@ -172,26 +174,37 @@
 	if (latPlusDrop > maxLatPlusDrop) maxLatPlusDrop = latPlusDrop;
     }
 
-    // we want to design totalLatency such that totalLatency -
-    // latencies[0] - drops[0] is a multiple of m_p.fftHop, so that we
-    // can get identical results in octave 0 to our reference
-    // implementation, making for easier testing (though other octaves
-    // will differ because of different resampler implementations)
+    int totalLatency = maxLatPlusDrop;
 
-    int totalLatency = maxLatPlusDrop;
     int lat0 = totalLatency - latencies[0] - drops[0];
     totalLatency = ceil(double(lat0 / m_p.fftHop) * m_p.fftHop)
 	+ latencies[0] + drops[0];
 
-//    cerr << "total latency = " << totalLatency << endl;
+    // We want (totalLatency - latencies[i]) to be a multiple of 2^i
+    // for each octave i, so that we do not end up with fractional
+    // octave latencies below. In theory this is hard, in practice if
+    // we ensure it for the last octave we should be OK.
+    double finalOctLat = latencies[m_octaves-1];
+    double finalOctFact = pow(2, m_octaves-1);
+    totalLatency =
+        int(round(finalOctLat +
+                  finalOctFact *
+                  ceil((totalLatency - finalOctLat) / finalOctFact)));
+
+#ifdef DEBUG_CQ
+    cerr << "total latency = " << totalLatency << endl;
+#endif
 
     // Padding as in the reference (will be introduced with the
     // latency compensation in the loop below)
     m_outputLatency = totalLatency + m_bigBlockSize
 	- m_p.firstCentre * pow(2, m_octaves-1);
 
-//    cerr << "m_bigBlockSize = " << m_bigBlockSize << ", firstCentre = "
-//	 << m_p.firstCentre << ", m_octaves = " << m_octaves << ", so m_outputLatency = " << m_outputLatency << endl;
+#ifdef DEBUG_CQ
+    cerr << "m_bigBlockSize = " << m_bigBlockSize << ", firstCentre = "
+	 << m_p.firstCentre << ", m_octaves = " << m_octaves
+         << ", so m_outputLatency = " << m_outputLatency << endl;
+#endif
 
     for (int i = 0; i < m_octaves; ++i) {
 
@@ -208,24 +221,39 @@
 	    double(totalLatency - latencies[i] - drops[i]
 		   + m_bigBlockSize) / factor;
 
+#ifdef DEBUG_CQ
+        cerr << "octave " << i << ": resampler latency = " << latencies[i]
+             << ", drop " << drops[i] << " (/factor = " << drops[i]/factor
+             << "), octaveLatency = " << octaveLatency << " -> "
+             << int(round(octaveLatency)) << " (diff * factor = "
+             << (octaveLatency - round(octaveLatency)) << " * "
+             << factor << " = "
+             << (octaveLatency - round(octaveLatency)) * factor << ")" << endl;
+
+        cerr << "double(" << totalLatency << " - " 
+             << latencies[i] << " - " << drops[i] << " + " 
+             << m_bigBlockSize << ") / " << factor << " = " 
+             << octaveLatency << endl;
+#endif
+
         m_buffers.push_back
-            (vector<double>(int(round(octaveLatency)), 0.0));
+            (RealSequence(int(round(octaveLatency)), 0.0));
     }
 
     m_fft = new FFTReal(m_p.fftSize);
 }
 
-vector<vector<double> > 
-ConstantQ::process(const vector<double> &td)
+ConstantQ::ComplexBlock
+ConstantQ::process(const RealSequence &td)
 {
     m_buffers[0].insert(m_buffers[0].end(), td.begin(), td.end());
 
     for (int i = 1; i < m_octaves; ++i) {
-        vector<double> dec = m_decimators[i]->process(td.data(), td.size());
+        RealSequence dec = m_decimators[i]->process(td.data(), td.size());
         m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end());
     }
 
-    vector<vector<double> > out;
+    ComplexBlock out;
 
     while (true) {
 
@@ -236,8 +264,6 @@
 	bool enough = true;
 	for (int i = 0; i < m_octaves; ++i) {
 	    int required = m_p.fftSize * pow(2, m_octaves - i - 1);
-//	    cerr << "for octave " << i << ", buf len =  "<< m_buffers[i].size() << " (need " << required << ")" << endl;
-
 	    if ((int)m_buffers[i].size() < required) {
 		enough = false;
 	    }
@@ -247,7 +273,7 @@
         int base = out.size();
         int totalColumns = pow(2, m_octaves - 1) * m_p.atomsPerFrame;
         for (int i = 0; i < totalColumns; ++i) {
-            out.push_back(vector<double>());
+            out.push_back(ComplexColumn());
         }
 
         for (int octave = 0; octave < m_octaves; ++octave) {
@@ -255,7 +281,7 @@
             int blocksThisOctave = pow(2, (m_octaves - octave - 1));
 
             for (int b = 0; b < blocksThisOctave; ++b) {
-                vector<vector<double> > block = processOctaveBlock(octave);
+                ComplexBlock block = processOctaveBlock(octave);
                 
                 for (int j = 0; j < m_p.atomsPerFrame; ++j) {
 
@@ -266,7 +292,7 @@
 
                     while (int(out[target].size()) < 
                            m_p.binsPerOctave * (octave + 1)) {
-                        out[target].push_back(0.0);
+                        out[target].push_back(Complex());
                     }
                     
                     for (int i = 0; i < m_p.binsPerOctave; ++i) {
@@ -281,42 +307,42 @@
     return out;
 }
 
-vector<vector<double> >
-ConstantQ::getRemainingBlocks()
+ConstantQ::ComplexBlock
+ConstantQ::getRemainingOutput()
 {
     // Same as padding added at start, though rounded up
     int pad = ceil(double(m_outputLatency) / m_bigBlockSize) * m_bigBlockSize;
-    vector<double> zeros(pad, 0.0);
+    RealSequence zeros(pad, 0.0);
     return process(zeros);
 }
 
-vector<vector<double> >
+ConstantQ::ComplexBlock
 ConstantQ::processOctaveBlock(int octave)
 {
-    vector<double> ro(m_p.fftSize, 0.0);
-    vector<double> io(m_p.fftSize, 0.0);
+    RealSequence ro(m_p.fftSize, 0.0);
+    RealSequence io(m_p.fftSize, 0.0);
 
     m_fft->forward(m_buffers[octave].data(), ro.data(), io.data());
 
-    vector<double> shifted;
+    RealSequence shifted;
     shifted.insert(shifted.end(), 
                    m_buffers[octave].begin() + m_p.fftHop,
                    m_buffers[octave].end());
     m_buffers[octave] = shifted;
 
-    vector<C> cv;
+    ComplexSequence cv;
     for (int i = 0; i < m_p.fftSize; ++i) {
-        cv.push_back(C(ro[i], io[i]));
+        cv.push_back(Complex(ro[i], io[i]));
     }
 
-    vector<C> cqrowvec = m_kernel->process(cv);
+    ComplexSequence cqrowvec = m_kernel->processForward(cv);
 
-    // Reform into a column matrix and use only the magnitude
-    vector<vector<double> > cqblock;
+    // Reform into a column matrix
+    ComplexBlock cqblock;
     for (int j = 0; j < m_p.atomsPerFrame; ++j) {
-        cqblock.push_back(vector<double>());
+        cqblock.push_back(ComplexColumn());
         for (int i = 0; i < m_p.binsPerOctave; ++i) {
-            cqblock[j].push_back(abs(cqrowvec[i * m_p.atomsPerFrame + j]));
+            cqblock[j].push_back(cqrowvec[i * m_p.atomsPerFrame + j]);
         }
     }
 

--- a/cpp-qm-dsp/ConstantQ.h	Thu May 08 09:35:10 2014 +0100
+++ b/cpp-qm-dsp/ConstantQ.h	Tue May 13 13:31:24 2014 +0100
@@ -32,33 +32,64 @@
 #ifndef CONSTANTQ_H
 #define CONSTANTQ_H
 
+#include "CQBase.h"
 #include "CQKernel.h"
 
-#include <vector>
-
 class Resampler;
 class FFTReal;
 
-class ConstantQ
+/**
+ * Calculate a complex sparse constant-Q representation from
+ * time-domain input.
+ *
+ * For a real (magnitude-only) or interpolated representation, see
+ * CQSpectrogram.
+ */
+class ConstantQ : public CQBase
 {
 public:
     ConstantQ(double sampleRate, 
 	      double minFreq, double maxFreq, 
 	      int binsPerOctave);
-    ~ConstantQ();
+    virtual ~ConstantQ();
 
-    double getSampleRate() const { return m_sampleRate; }
-    int getBinsPerOctave() const { return m_binsPerOctave; }
-    int getOctaves() const { return m_octaves; }
-    int getTotalBins() const { return m_octaves * m_binsPerOctave; }
-    int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; }
-    int getLatency() const { return m_outputLatency; } 
-    double getMaxFrequency() const { return m_p.maxFrequency; }
-    double getMinFrequency() const; // actual min, not that passed to ctor
-    double getBinFrequency(int bin) const;
+    virtual double getSampleRate() const { return m_sampleRate; }
+    virtual int getBinsPerOctave() const { return m_binsPerOctave; }
+    virtual int getOctaves() const { return m_octaves; }
+    virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; }
+    virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; }
+    virtual int getLatency() const { return m_outputLatency; } 
+    virtual double getMaxFrequency() const { return m_p.maxFrequency; }
+    virtual double getMinFrequency() const;
+    virtual double getBinFrequency(int bin) const;
 
-    std::vector<std::vector<double> > process(const std::vector<double> &);
-    std::vector<std::vector<double> > getRemainingBlocks();
+    /**
+     * Given a series of time-domain samples, return a series of
+     * constant-Q columns. Any samples left over (that did not fit
+     * into a constant-Q processing block) are saved for the next call
+     * to process or getRemainingBlocks.
+     *
+     * Each column contains a series of constant-Q bin values ordered
+     * from highest to lowest frequency.
+     *
+     * Columns are of variable height: each will contain at least
+     * getBinsPerOctave() values, because the highest-frequency octave
+     * is always present, but a second octave (if requested) will
+     * appear only in alternate columns, a third octave only in every
+     * fourth column, and so on.
+     *
+     * If you need a format in which all columns are of equal height
+     * and every bin contains a value, use CQInterpolated instead of
+     * ConstantQ.
+     */
+    ComplexBlock process(const RealSequence &);
+
+    /**
+     * Return the remaining constant-Q columns following the end of
+     * processing. Any buffered input is padded so as to ensure that
+     * all input provided to process() will have been returned.
+     */
+    ComplexBlock getRemainingOutput();
 
 private:
     double m_sampleRate;
@@ -72,14 +103,14 @@
     int m_bigBlockSize;
 
     std::vector<Resampler *> m_decimators;
-    std::vector<std::vector<double> > m_buffers;
+    std::vector<RealSequence> m_buffers;
 
     int m_outputLatency;
 
     FFTReal *m_fft;
 
     void initialise();
-    std::vector<std::vector<double> > processOctaveBlock(int octave);
+    ComplexBlock processOctaveBlock(int octave);
 };
 
 #endif

--- a/cpp-qm-dsp/Makefile	Thu May 08 09:35:10 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,38 +0,0 @@
-
-DEFINES := -DUSE_PTHREADS
-
-CFLAGS := -I../.. $(CFLAGS) $(DEFINES)
-
-#CXXFLAGS := -I../.. -fPIC -Wall -g $(CXXFLAGS) $(DEFINES)
-CXXFLAGS := -I../.. -fPIC -Wall -O3 -ffast-math -ftree-vectorize $(CXXFLAGS) $(DEFINES)
-
-LDFLAGS := $(LDFLAGS)
-
-#VG	:= valgrind
-
-LIBS 		:= ../../qm-dsp/libqm-dsp.a -lpthread
-PROGRAM_LIBS	:= -lsndfile
-
-SOURCES	:= CQKernel.cpp ConstantQ.cpp test.cpp
-
-OBJECTS := $(SOURCES:.cpp=.o)
-OBJECTS := $(OBJECTS:.c=.o)
-
-PROGRAM	:= test
-
-all: $(PROGRAM)
-
-test:	$(OBJECTS)
-	$(CXX) -o $@ $^ $(LDFLAGS) $(LIBS) $(PROGRAM_LIBS)
-
-clean: 
-	rm -f *.o 
-
-depend:
-	makedepend -Y $(SOURCES)
-
-# DO NOT DELETE
-
-CQKernel.o: CQKernel.h
-ConstantQ.o: ConstantQ.h CQKernel.h
-test.o: ConstantQ.h CQKernel.h

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/processfile.cpp	Tue May 13 13:31:24 2014 +0100
@@ -0,0 +1,218 @@
+
+#include "ConstantQ.h"
+#include "CQInverse.h"
+
+#include <sndfile.h>
+
+#include <iostream>
+
+using std::vector;
+using std::cerr;
+using std::endl;
+
+#include <cstring>
+
+#include <unistd.h>
+#include <sys/time.h>
+
+int main(int argc, char **argv)
+{
+    if (argc < 3 || argc > 4) {
+	cerr << "Usage: " << argv[0] << " infile.wav outfile.wav [differencefile.wav]" << endl;
+	return 2;
+    }
+
+    char *fileName = strdup(argv[1]);
+    char *fileNameOut = strdup(argv[2]);
+    char *diffFileName = (argc == 4 ? strdup(argv[3]) : 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];
+
+    ConstantQ cq(sfinfo.samplerate,
+		 100, sfinfo.samplerate / 3,
+		 60);
+
+    CQInverse cqi(sfinfo.samplerate,
+		 100, sfinfo.samplerate / 3,
+		 60);
+
+    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;
+
+    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) {
+		    cqout[i] -= buffer[outframe + i - latency];
+		}
+	    }
+	    
+	    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) {
+		r[i] -= buffer[outframe + i - latency];
+	    }
+	}
+	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;
+
+    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/cpp-qm-dsp/test.cpp	Thu May 08 09:35:10 2014 +0100
+++ b/cpp-qm-dsp/test.cpp	Tue May 13 13:31:24 2014 +0100
@@ -28,7 +28,7 @@
     authorization.
 */
 
-#include "ConstantQ.h"
+#include "CQSpectrogram.h"
 
 #include <iostream>
 #include <vector>
@@ -50,10 +50,10 @@
 	in.push_back(sin(i * M_PI / 2.0));
     }
 
-    ConstantQ k(8, 1, 4, 4);
+    CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros);
 
     vector<vector<double> > out = k.process(in);
-    vector<vector<double> > rest = k.getRemainingBlocks();
+    vector<vector<double> > rest = k.getRemainingOutput();
 
     out.insert(out.end(), rest.begin(), rest.end());
 

--- a/vamp/CQVamp.cpp	Thu May 08 09:35:10 2014 +0100
+++ b/vamp/CQVamp.cpp	Tue May 13 13:31:24 2014 +0100
@@ -47,7 +47,7 @@
     m_maxMIDIPitch(84),
     m_tuningFrequency(440),
     m_bpo(24),
-    m_interpolation(CQInterpolated::Linear),
+    m_interpolation(CQSpectrogram::InterpolateLinear),
     m_cq(0),
     m_maxFrequency(inputSampleRate/2),
     m_minFrequency(46),
@@ -155,7 +155,7 @@
     desc.defaultValue = 2;
     desc.isQuantized = true;
     desc.quantizeStep = 1;
-    desc.valueNames.push_back("None, leave empty");
+    desc.valueNames.push_back("None, leave as zero");
     desc.valueNames.push_back("None, repeat prior value");
     desc.valueNames.push_back("Linear interpolation");
     list.push_back(desc);
@@ -198,7 +198,7 @@
     } else if (param == "bpo") {
         m_bpo = lrintf(value);
     } else if (param == "interpolation") {
-        m_interpolation = (CQInterpolated::Interpolation)lrintf(value);
+        m_interpolation = (CQSpectrogram::Interpolation)lrintf(value);
     } else {
         std::cerr << "WARNING: CQVamp::setParameter: unknown parameter \""
                   << param << "\"" << std::endl;
@@ -224,7 +224,7 @@
     m_maxFrequency = Pitch::getFrequencyForPitch
         (m_maxMIDIPitch, 0, m_tuningFrequency);
 
-    m_cq = new CQInterpolated
+    m_cq = new CQSpectrogram
 	(m_inputSampleRate, m_minFrequency, m_maxFrequency, m_bpo,
          m_interpolation);
 
@@ -236,7 +236,7 @@
 {
     if (m_cq) {
 	delete m_cq;
-	m_cq = new CQInterpolated
+	m_cq = new CQSpectrogram
 	    (m_inputSampleRate, m_minFrequency, m_maxFrequency, m_bpo,
              m_interpolation);
     }
@@ -272,7 +272,7 @@
 
     if (m_cq) {
         char name[20];
-        for (int i = 0; i < d.binCount; ++i) {
+        for (int i = 0; i < (int)d.binCount; ++i) {
             float freq = m_cq->getBinFrequency(i);
             sprintf(name, "%.1f Hz", freq);
             d.binNames.push_back(name);
@@ -314,7 +314,7 @@
 CQVamp::FeatureSet
 CQVamp::getRemainingFeatures()
 {
-    vector<vector<double> > cqout = m_cq->getRemainingBlocks();
+    vector<vector<double> > cqout = m_cq->getRemainingOutput();
     return convertToFeatures(cqout);
 }
 

--- a/vamp/CQVamp.h	Thu May 08 09:35:10 2014 +0100
+++ b/vamp/CQVamp.h	Tue May 13 13:31:24 2014 +0100
@@ -34,7 +34,7 @@
 
 #include <vamp-sdk/Plugin.h>
 
-#include "cpp-qm-dsp/CQInterpolated.h"
+#include "cpp-qm-dsp/CQSpectrogram.h"
 
 class ConstantQ;
 
@@ -75,9 +75,9 @@
     int m_maxMIDIPitch;
     float m_tuningFrequency;
     int m_bpo;
-    CQInterpolated::Interpolation m_interpolation;
+    CQSpectrogram::Interpolation m_interpolation;
 
-    CQInterpolated *m_cq;
+    CQSpectrogram *m_cq;
     float m_maxFrequency;
     float m_minFrequency;
     int m_stepSize;