changeset 116:6deec2a51d13

Moving to standalone library layout

--- a/cpp-qm-dsp/CQBase.h	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,60 +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 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/CQInverse.cpp	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,411 +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 "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, 50, 0.05);
-
-#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 j = 0; j < m_octaves; ++j) {
-        int factor = pow(2, j);
-        int latency = (j > 0 ? m_upsamplers[j]->getLatency() : 0) / factor;
-        for (int i = 0; i < (latency + m_p.fftSize) / m_p.fftHop; ++i) {
-            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];
-    }
-}
-

--- a/cpp-qm-dsp/CQInverse.h	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,91 +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 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 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,337 +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 "CQKernel.h"
-
-#include "maths/MathUtilities.h"
-#include "dsp/transforms/FFT.h"
-#include "base/Window.h"
-
-#include <cmath>
-#include <cassert>
-#include <vector>
-#include <iostream>
-#include <algorithm>
-
-using std::vector;
-using std::complex;
-using std::cerr;
-using std::endl;
-
-typedef std::complex<double> C;
-
-CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave) :
-    m_fft(0)
-{
-    m_p.sampleRate = sampleRate;
-    m_p.maxFrequency = maxFreq;
-    m_p.binsPerOctave = binsPerOctave;
-    generateKernel();
-}
-
-CQKernel::~CQKernel()
-{
-    delete m_fft;
-}
-
-void
-CQKernel::generateKernel()
-{
-    double q = 1;
-    double atomHopFactor = 0.25;
-    double thresh = 0.0005;
-
-    double bpo = m_p.binsPerOctave;
-
-    m_p.minFrequency = (m_p.maxFrequency / 2) * pow(2, 1.0/bpo);
-    m_p.Q = q / (pow(2, 1.0/bpo) - 1.0);
-
-    double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency);
-    double minNK = round
-        (m_p.Q * m_p.sampleRate /
-         (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo)));
-
-    if (minNK == 0 || maxNK == 0) {
-        // most likely pathological parameters of some sort
-        cerr << "WARNING: CQKernel::generateKernel: minNK or maxNK is zero (minNK == " << minNK << ", maxNK == " << maxNK << "), not generating a kernel" << endl;
-        m_p.atomSpacing = 0;
-        m_p.firstCentre = 0;
-        m_p.fftSize = 0;
-        m_p.atomsPerFrame = 0;
-        m_p.lastCentre = 0;
-        m_p.fftHop = 0;
-        return;
-    }
-
-    m_p.atomSpacing = round(minNK * atomHopFactor);
-    m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing);
-    m_p.fftSize = MathUtilities::nextPowerOfTwo
-        (m_p.firstCentre + ceil(maxNK / 2.0));
-
-    m_p.atomsPerFrame = floor
-        (1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing);
-
-    cerr << "atomsPerFrame = " << m_p.atomsPerFrame << " (atomHopFactor = " << atomHopFactor << ", atomSpacing = " << m_p.atomSpacing << ", fftSize = " << m_p.fftSize << ", maxNK = " << maxNK << ", firstCentre = " << m_p.firstCentre << ")" << endl;
-
-    m_p.lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing;
-
-    m_p.fftHop = (m_p.lastCentre + m_p.atomSpacing) - m_p.firstCentre;
-
-    cerr << "fftHop = " << m_p.fftHop << endl;
-
-    m_fft = new FFT(m_p.fftSize);
-
-    for (int k = 1; k <= m_p.binsPerOctave; ++k) {
-        
-        int nk = round(m_p.Q * m_p.sampleRate /
-                       (m_p.minFrequency * pow(2, ((k-1.0) / bpo))));
-
-        // The MATLAB version uses a symmetric window, but our windows
-        // are periodic. A symmetric window of size N is a periodic
-        // one of size N-1 with the first element stuck on the end
-        Window<double> w(BlackmanHarrisWindow, nk-1);
-        vector<double> win = w.getWindowData();
-        win.push_back(win[0]);
-
-        for (int i = 0; i < (int)win.size(); ++i) {
-            win[i] = sqrt(win[i]) / nk;
-        }
-
-        double fk = m_p.minFrequency * pow(2, ((k-1.0) / bpo));
-
-        vector<double> reals, imags;
-        
-        for (int i = 0; i < nk; ++i) {
-            double arg = (2.0 * M_PI * fk * i) / m_p.sampleRate;
-            reals.push_back(win[i] * cos(arg));
-            imags.push_back(win[i] * sin(arg));
-        }
-
-        int atomOffset = m_p.firstCentre - int(ceil(nk/2.0));
-
-        for (int i = 0; i < m_p.atomsPerFrame; ++i) {
-
-            int shift = atomOffset + (i * m_p.atomSpacing);
-
-            vector<double> rin(m_p.fftSize, 0.0);
-            vector<double> iin(m_p.fftSize, 0.0);
-
-            for (int j = 0; j < nk; ++j) {
-                rin[j + shift] = reals[j];
-                iin[j + shift] = imags[j];
-            }
-
-            vector<double> rout(m_p.fftSize, 0.0);
-            vector<double> iout(m_p.fftSize, 0.0);
-
-            m_fft->process(false,
-                           rin.data(), iin.data(),
-                           rout.data(), iout.data());
-
-            // Keep this dense for the moment (until after
-            // normalisation calculations)
-
-            vector<C> row;
-
-            for (int j = 0; j < m_p.fftSize; ++j) {
-                if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) {
-                    row.push_back(C(0, 0));
-                } else {
-                    row.push_back(C(rout[j] / m_p.fftSize,
-                                    iout[j] / m_p.fftSize));
-                }
-            }
-
-            m_kernel.origin.push_back(0);
-            m_kernel.data.push_back(row);
-        }
-    }
-
-    assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
-
-    // print density as diagnostic
-
-    int nnz = 0;
-    for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
-        for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) {
-            if (m_kernel.data[i][j] != C(0, 0)) {
-                ++nnz;
-            }
-        }
-    }
-
-    cerr << "size = " << m_kernel.data.size() << "*" << m_kernel.data[0].size() << " (fft size = " << m_p.fftSize << ")" << endl;
-
-    assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
-    assert((int)m_kernel.data[0].size() == m_p.fftSize);
-
-    cerr << "density = " << double(nnz) / double(m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize) << " (" << nnz << " of " << m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize << ")" << endl;
-
-    finaliseKernel();
-}
-
-static bool ccomparator(C &c1, C &c2)
-{
-    return abs(c1) < abs(c2);
-}
-
-static int maxidx(vector<C> &v)
-{
-    return std::max_element(v.begin(), v.end(), ccomparator) - v.begin();
-}
-
-void
-CQKernel::finaliseKernel()
-{
-    // calculate weight for normalisation
-
-    int wx1 = maxidx(m_kernel.data[0]);
-    int wx2 = maxidx(m_kernel.data[m_kernel.data.size()-1]);
-
-    vector<vector<C> > subset(m_kernel.data.size());
-    for (int j = wx1; j <= wx2; ++j) {
-        for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
-            subset[i].push_back(m_kernel.data[i][j]);
-        }
-    }
-
-    int nrows = subset.size();
-    int ncols = subset[0].size();
-    vector<vector<C> > square(ncols); // conjugate transpose of subset * subset
-
-    for (int i = 0; i < nrows; ++i) {
-        assert((int)subset[i].size() == ncols);
-    }
-
-    for (int j = 0; j < ncols; ++j) {
-        for (int i = 0; i < ncols; ++i) {
-            C v(0, 0);
-            for (int k = 0; k < nrows; ++k) {
-                v += subset[k][i] * conj(subset[k][j]);
-            }
-            square[i].push_back(v);
-        }
-    }
-
-    vector<double> wK;
-    double q = 1; //!!! duplicated from constructor
-    for (int i = round(1.0/q); i < ncols - round(1.0/q) - 2; ++i) {
-        wK.push_back(abs(square[i][i]));
-    }
-
-    double weight = double(m_p.fftHop) / m_p.fftSize;
-    weight /= MathUtilities::mean(wK.data(), wK.size());
-    weight = sqrt(weight);
-    
-    cerr << "weight = " << weight << endl;
-
-    // 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;
-
-    for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
-
-        sk.origin.push_back(0);
-        sk.data.push_back(vector<C>());
-
-        int lastNZ = 0;
-        for (int j = (int)m_kernel.data[i].size()-1; j >= 0; --j) {
-            if (abs(m_kernel.data[i][j]) != 0.0) {
-                lastNZ = j;
-                break;
-            }
-        }
-
-        bool haveNZ = false;
-        for (int j = 0; j <= lastNZ; ++j) {
-            if (haveNZ || abs(m_kernel.data[i][j]) != 0.0) {
-                if (!haveNZ) sk.origin[i] = j;
-                haveNZ = true;
-                sk.data[i].push_back(conj(m_kernel.data[i][j]) * weight);
-            }
-        }
-    }
-
-    m_kernel = sk;
-}
-
-vector<C>
-CQKernel::processForward(const vector<C> &cv)
-{
-    // straightforward matrix multiply (taking into account m_kernel's
-    // slightly-sparse representation)
-
-    if (m_kernel.data.empty()) return vector<C>();
-
-    int nrows = m_p.binsPerOctave * m_p.atomsPerFrame;
-
-    vector<C> rv(nrows, C());
-
-    for (int i = 0; i < nrows; ++i) {
-        int len = m_kernel.data[i].size();
-        for (int j = 0; j < len; ++j) {
-            rv[i] += cv[j + m_kernel.origin[i]] * m_kernel.data[i][j];
-        }
-    }
-
-    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.
-
-    if (m_kernel.data.empty()) return vector<C>();
-
-    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 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,82 +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_KERNEL_H
-#define CQ_KERNEL_H
-
-#include <vector>
-#include <complex>
-
-class FFT;
-
-class CQKernel
-{
-public:
-    CQKernel(double sampleRate, double maxFreq, int binsPerOctave);
-    ~CQKernel();
-    
-    struct Properties {
-        double sampleRate;
-        double maxFrequency;
-        double minFrequency;
-        int binsPerOctave;
-        int fftSize;
-        int fftHop;
-        int atomsPerFrame;
-        int atomSpacing;
-        int firstCentre;
-        int lastCentre;
-        double Q;
-    };
-
-    Properties getProperties() const { return m_p; }
-
-    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:
-    Properties m_p;
-    FFT *m_fft;
-
-    struct KernelMatrix {
-        std::vector<int> origin;
-        std::vector<std::vector<std::complex<double> > > data;
-    };
-    KernelMatrix m_kernel;
-
-    void generateKernel();
-    void finaliseKernel();
-};
-
-#endif

--- a/cpp-qm-dsp/CQSpectrogram.cpp	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,261 +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 "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;
-	    }
-	}
-    }
-
-//    cerr << "fetchLinear: firstFullHeight = " << firstFullHeight << ", secondFullHeight = " << secondFullHeight << endl;
-
-    if (firstFullHeight < 0) {
-	if (insist) {
-            return fetchHold(true);
-	} 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) {
-            return fetchHold(true);
-	} 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;
-}
-
-
-	

--- a/cpp-qm-dsp/CQSpectrogram.h	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +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 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 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,352 +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 "ConstantQ.h"
-
-#include "CQKernel.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
-
-ConstantQ::ConstantQ(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();
-}
-
-ConstantQ::~ConstantQ()
-{
-    delete m_fft;
-    for (int i = 0; i < (int)m_decimators.size(); ++i) {
-        delete m_decimators[i];
-    }
-    delete m_kernel;
-}
-
-double
-ConstantQ::getMinFrequency() const
-{
-    return m_p.minFrequency / pow(2.0, m_octaves - 1);
-}
-
-double
-ConstantQ::getBinFrequency(int bin) const
-{
-    return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave()));
-}
-
-void
-ConstantQ::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_decimators.push_back(0);
-
-    for (int i = 1; i < m_octaves; ++i) {
-
-        int factor = pow(2, i);
-
-        Resampler *r = new Resampler
-            (sourceRate, sourceRate / factor, 50, 0.05);
-
-#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.
-        // 
-        // Our decimator uses a linear phase filter, but being causal
-        // it is not zero phase: it has a latency that depends on the
-        // decimation factor. Those latencies have been calculated
-        // per-octave and are available to us in the latencies
-        // array. Left to its own devices, the first input sample will
-        // appear at output sample 0 in the highest octave (where no
-        // decimation is needed), sample number latencies[1] in the
-        // next octave down, latencies[2] in the next one, etc. We get
-        // to apply some artificial per-octave latency after the
-        // decimator in the processing chain, in order to compensate
-        // for the differing latencies associated with different
-        // decimation factors. How much should we insert?
-        //
-        // The outputs of the decimators are at different rates (in
-        // terms of the relation between clock time and samples) and
-        // we want them aligned in terms of time. So, for example, a
-        // latency of 10 samples with a decimation factor of 2 is
-        // equivalent to a latency of 20 with no decimation -- they
-        // both result in the first output sample happening at the
-        // same equivalent time in milliseconds.
-	// 
-	// So here we record the latency added by the decimator, in
-	// terms of the sample rate of the undecimated signal. Then we
-	// use that to compensate in a moment, when we've discovered
-	// what the longest latency across all octaves is.
-
-        latencies.push_back(r->getLatency() * factor);
-        m_decimators.push_back(r);
-    }
-
-    m_bigBlockSize = m_p.fftSize * pow(2, m_octaves - 1);
-
-    // Now add in the extra padding and compensate for hops that must
-    // be dropped in order to align the atom centres across
-    // octaves. Again this is a bit trickier because we are doing it
-    // at input rather than output and so must work in per-octave
-    // sample rates rather than output blocks
-
-    int emptyHops = m_p.firstCentre / m_p.atomSpacing;
-
-    vector<int> drops;
-    for (int i = 0; i < m_octaves; ++i) {
-	int factor = pow(2, i);
-	int dropHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops;
-	int drop = ((dropHops * m_p.fftHop) * factor) / m_p.atomsPerFrame;
-	drops.push_back(drop);
-    }
-
-    int maxLatPlusDrop = 0;
-    for (int i = 0; i < m_octaves; ++i) {
-	int latPlusDrop = latencies[i] + drops[i];
-	if (latPlusDrop > maxLatPlusDrop) maxLatPlusDrop = latPlusDrop;
-    }
-
-    int totalLatency = maxLatPlusDrop;
-
-    int lat0 = totalLatency - latencies[0] - drops[0];
-    totalLatency = ceil(double(lat0 / m_p.fftHop) * m_p.fftHop)
-	+ latencies[0] + drops[0];
-
-    // 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);
-
-#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) {
-
-	double factor = pow(2, i);
-
-	// Calculate the difference between the total latency applied
-	// across all octaves, and the existing latency due to the
-	// decimator for this octave, and then convert it back into
-	// the sample rate appropriate for the output latency of this
-	// decimator -- including one additional big block of padding
-	// (as in the reference).
-
-	double octaveLatency =
-	    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
-            (RealSequence(int(round(octaveLatency)), 0.0));
-    }
-
-    m_fft = new FFTReal(m_p.fftSize);
-}
-
-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) {
-        RealSequence dec = m_decimators[i]->process(td.data(), td.size());
-        m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end());
-    }
-
-    ComplexBlock out;
-
-    while (true) {
-
-	// We could have quite different remaining sample counts in
-	// different octaves, because (apart from the predictable
-	// added counts for decimator output on each block) we also
-	// have variable additional latency per octave
-	bool enough = true;
-	for (int i = 0; i < m_octaves; ++i) {
-	    int required = m_p.fftSize * pow(2, m_octaves - i - 1);
-	    if ((int)m_buffers[i].size() < required) {
-		enough = false;
-	    }
-	}
-	if (!enough) break;
-
-        int base = out.size();
-        int totalColumns = pow(2, m_octaves - 1) * m_p.atomsPerFrame;
-        for (int i = 0; i < totalColumns; ++i) {
-            out.push_back(ComplexColumn());
-        }
-
-        for (int octave = 0; octave < m_octaves; ++octave) {
-
-            int blocksThisOctave = pow(2, (m_octaves - octave - 1));
-
-            for (int b = 0; b < blocksThisOctave; ++b) {
-                ComplexBlock block = processOctaveBlock(octave);
-                
-                for (int j = 0; j < m_p.atomsPerFrame; ++j) {
-
-                    int target = base +
-			    (b * (totalColumns / blocksThisOctave) + 
-			     (j * ((totalColumns / blocksThisOctave) /
-				   m_p.atomsPerFrame)));
-
-                    while (int(out[target].size()) < 
-                           m_p.binsPerOctave * (octave + 1)) {
-                        out[target].push_back(Complex());
-                    }
-                    
-                    for (int i = 0; i < m_p.binsPerOctave; ++i) {
-                        out[target][m_p.binsPerOctave * octave + i] = 
-                            block[j][m_p.binsPerOctave - i - 1];
-                    }
-                }
-            }
-        }
-    }
-
-    return out;
-}
-
-ConstantQ::ComplexBlock
-ConstantQ::getRemainingOutput()
-{
-    // Same as padding added at start, though rounded up
-    int pad = ceil(double(m_outputLatency) / m_bigBlockSize) * m_bigBlockSize;
-    RealSequence zeros(pad, 0.0);
-    return process(zeros);
-}
-
-ConstantQ::ComplexBlock
-ConstantQ::processOctaveBlock(int octave)
-{
-    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());
-
-    RealSequence shifted;
-    shifted.insert(shifted.end(), 
-                   m_buffers[octave].begin() + m_p.fftHop,
-                   m_buffers[octave].end());
-    m_buffers[octave] = shifted;
-
-    ComplexSequence cv;
-    for (int i = 0; i < m_p.fftSize; ++i) {
-        cv.push_back(Complex(ro[i], io[i]));
-    }
-
-    ComplexSequence cqrowvec = m_kernel->processForward(cv);
-
-    // Reform into a column matrix
-    ComplexBlock cqblock;
-    for (int j = 0; j < m_p.atomsPerFrame; ++j) {
-        cqblock.push_back(ComplexColumn());
-        for (int i = 0; i < m_p.binsPerOctave; ++i) {
-            cqblock[j].push_back(cqrowvec[i * m_p.atomsPerFrame + j]);
-        }
-    }
-
-    return cqblock;
-}
-
-

--- a/cpp-qm-dsp/ConstantQ.h	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,117 +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 CONSTANTQ_H
-#define CONSTANTQ_H
-
-#include "CQBase.h"
-#include "CQKernel.h"
-
-class Resampler;
-class FFTReal;
-
-/**
- * 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);
-    virtual ~ConstantQ();
-
-    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;
-
-    /**
-     * 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;
-    double m_maxFrequency;
-    double m_minFrequency;
-    int m_binsPerOctave;
-    int m_octaves;
-
-    CQKernel *m_kernel;
-    CQKernel::Properties m_p;
-    int m_bigBlockSize;
-
-    std::vector<Resampler *> m_decimators;
-    std::vector<RealSequence> m_buffers;
-
-    int m_outputLatency;
-
-    FFTReal *m_fft;
-
-    void initialise();
-    ComplexBlock processOctaveBlock(int octave);
-};
-
-#endif
-

--- a/cpp-qm-dsp/processfile.cpp	Thu May 15 11:59:11 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/cpp-qm-dsp/test.cpp	Thu May 15 11:59:11 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/cq/CQBase.h	Thu May 15 12:04:00 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cq/CQInverse.h	Thu May 15 12:04:00 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cq/CQKernel.h	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,82 @@
+/* -*- 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_KERNEL_H
+#define CQ_KERNEL_H
+
+#include <vector>
+#include <complex>
+
+class FFT;
+
+class CQKernel
+{
+public:
+    CQKernel(double sampleRate, double maxFreq, int binsPerOctave);
+    ~CQKernel();
+    
+    struct Properties {
+        double sampleRate;
+        double maxFrequency;
+        double minFrequency;
+        int binsPerOctave;
+        int fftSize;
+        int fftHop;
+        int atomsPerFrame;
+        int atomSpacing;
+        int firstCentre;
+        int lastCentre;
+        double Q;
+    };
+
+    Properties getProperties() const { return m_p; }
+
+    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:
+    Properties m_p;
+    FFT *m_fft;
+
+    struct KernelMatrix {
+        std::vector<int> origin;
+        std::vector<std::vector<std::complex<double> > > data;
+    };
+    KernelMatrix m_kernel;
+
+    void generateKernel();
+    void finaliseKernel();
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cq/CQSpectrogram.h	Thu May 15 12:04:00 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cq/ConstantQ.h	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,117 @@
+/* -*- 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 CONSTANTQ_H
+#define CONSTANTQ_H
+
+#include "CQBase.h"
+#include "CQKernel.h"
+
+class Resampler;
+class FFTReal;
+
+/**
+ * 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);
+    virtual ~ConstantQ();
+
+    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;
+
+    /**
+     * 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;
+    double m_maxFrequency;
+    double m_minFrequency;
+    int m_binsPerOctave;
+    int m_octaves;
+
+    CQKernel *m_kernel;
+    CQKernel::Properties m_p;
+    int m_bigBlockSize;
+
+    std::vector<Resampler *> m_decimators;
+    std::vector<RealSequence> m_buffers;
+
+    int m_outputLatency;
+
+    FFTReal *m_fft;
+
+    void initialise();
+    ComplexBlock processOctaveBlock(int octave);
+};
+
+#endif
+

Binary file misc/paper/smc2010.pdf has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/Makefile	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,6 @@
+
+MAY_DIR	:= ../../may
+MAY_YC	:= $(MAY_DIR)/bin/yc
+
+test:
+	$(MAY_YC) ./test.yeti

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/build.xml	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,74 @@
+<project name="cqt" default="test" basedir=".">
+
+  <property name="maydir" value="${basedir}/../../may"/>
+  <property name="yetidir" value="${basedir}/../../yeti"/>
+
+  <property name="jardir" value="${maydir}/ext/jar"/>
+
+  <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/jzy3d-swt-0.9.1.jar:${jardir}/jzy3d-api-0.9.1.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/>
+
+  <condition property="archtag" value="linux32">
+    <os family="unix" arch="i386"/>
+  </condition>
+  <condition property="archtag" value="linux64">
+    <os family="unix" arch="amd64"/>
+  </condition>
+  <condition property="archtag" value="win32">
+    <os family="windows" arch="x86"/>
+  </condition>
+  <condition property="archtag" value="win64">
+    <os family="windows" arch="amd64"/>
+  </condition>
+  <condition property="archtag" value="osx">
+    <os family="mac"/>
+  </condition>
+
+  <target name="taskdef">
+    <taskdef name="yetic" classname="yeti.lang.compiler.YetiTask"
+	     classpath="${yetidir}/yeti.jar:${maydir}/may.jar:${extjars}" />
+  </target>
+
+  <target name="prepare">
+    <mkdir dir="${basedir}/classes"/>
+  </target>
+
+  <target name="yeticlasses" depends="taskdef">
+    <yetic srcdir="${basedir}" 
+	   destdir="${basedir}/classes"
+	   includes="**/*.yeti"
+	   preload="yeti/lang/std:yeti/lang/io"/>
+  </target>
+
+  <target name="classes" depends="prepare,yeticlasses"/>
+
+  <target name="jar" depends="classes,taskdef">
+    <jar jarfile="${basedir}/cqt.jar">
+      <fileset dir="${basedir}/classes" 
+	       includes="**/*.class" 
+	       excludes="**/test*.class"/>
+    </jar>
+  </target>
+
+  <target name="testjar" depends="classes,taskdef">
+    <jar jarfile="${basedir}/test.jar">
+      <fileset dir="${basedir}/classes" 
+	       includes="**/test*.class"/>
+    </jar>
+  </target>
+
+  <target name="test" depends="jar,testjar,taskdef">
+    <java classpath="${basedir}/test.jar:${basedir}/cqt.jar:${maydir}/may.jar:${yetidir}/yeti.jar:${extjars}"
+	  classname="test"
+	  fork="true" failonerror="true">
+      <sysproperty key="java.library.path" path="${maydir}/ext/native/${archtag}"/>
+    </java>
+  </target>
+
+  <target name="clean">
+    <delete dir="${basedir}/classes"/>
+  </target>
+
+  <target name="rebuild" depends="clean,jar"/>
+
+</project>
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/cqt.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,260 @@
+/*
+    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.
+*/
+
+module cqt;
+
+cqtkernel = load cqtkernel;
+resample = load may.stream.resample;
+manipulate = load may.stream.manipulate;
+mat = load may.matrix;
+cm = load may.matrix.complex;
+framer = load may.stream.framer;
+cplx = load may.complex;
+fft = load may.transform.fft;
+vec = load may.vector;
+ch = load may.stream.channels;
+
+{ pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc;
+
+cqt { maxFreq, minFreq, binsPerOctave } str =
+   (sampleRate = str.sampleRate;
+    octaves = ceil (log2 (maxFreq / minFreq));
+//    actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave));
+
+    kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave };
+
+//    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave), fftSize = \(kdata.fftSize), hop = \(kdata.fftHop)";
+
+//    eprintln "atomsPerFrame = \(kdata.atomsPerFrame)";
+
+    padding = (kdata.fftSize * (pow 2 (octaves-1)));
+
+//    eprintln "padding = \(padding)";
+
+    str = manipulate.paddedBy padding str;
+
+    streams = manipulate.duplicated octaves str;
+
+    // forward transform uses the conjugate-transposed kernel, inverse
+    // uses the original
+    kernel = cm.transposed (cm.conjugateTransposed kdata.kernel);
+
+//    eprintln "have kernel";
+
+    fftFunc = fft.forward kdata.fftSize;
+
+    cqblocks =
+        map do octave:
+            frames = map ch.mixedDown //!!! mono for now
+               (framer.frames kdata.fftSize [ Hop kdata.fftHop, Padded false ]
+                   (resample.decimated (pow 2 octave) streams[octave]));
+            map do frame:
+                freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize));
+// eprintln "octave = \(octave), frame = \(vec.list frame)";
+// eprintln "octave = \(octave), freq = \(freq)";
+                cm.product kernel (cm.newComplexColumnVector freq);
+            done frames;
+        done [0..octaves-1];
+
+    // The cqblocks list is a list<list<matrix>>. Each top-level list
+    // corresponds to an octave, from highest to lowest, each having
+    // twice as many elements in its list as the next octave. The
+    // sub-lists are sampled in time with an effective spacing of
+    // fftSize * 2^(octave-1) audio frames, and the matrices are row
+    // vectors with atomsPerFrame * binsPerOctave complex elements.
+    //
+    // ***
+    // 
+    // In a typical constant-Q structure, each (2^(octaves-1) *
+    // fftHop) input frames gives us an output structure conceptually
+    // like this:
+    //
+    // [][][][][][][][]   <- fftHop frames per highest-octave output value
+    // [][][][][][][][]      layered as many times as binsPerOctave (here 2)
+    // [--][--][--][--]   <- fftHop*2 frames for the next lower octave
+    // [--][--][--][--]      etc
+    // [------][------]
+    // [------][------]
+    // [--------------]
+    // [--------------]
+    //
+    // ***
+    //
+    // But the kernel we're using here has more than one temporally
+    // spaced atom; each individual cell is a row vector with
+    // atomsPerFrame * binsPerOctave elements, but that actually
+    // represents a rectangular matrix of result cells with width
+    // atomsPerFrame and height binsPerOctave. The columns of this
+    // matrix (the atoms) then need to be spaced by 2^(octave-1)
+    // relative to those from the highest octave.
+
+    // Reshape each row vector into the appropriate rectangular matrix
+    // and split into single-atom columns
+
+    emptyHops = kdata.firstCentre / kdata.atomSpacing; //!!! int? round?
+//    maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops;
+//    eprintln "maxDrop = \(maxDrop)";
+
+    cqblocks =
+        map do octlist:
+            concatMap do rv:
+                cm.asColumns
+                   (cm.generate do row col:
+                        cm.at rv ((row * kdata.atomsPerFrame) + col) 0
+                    done {
+                        rows = kdata.binsPerOctave,
+                        columns = kdata.atomsPerFrame
+                    })
+            done octlist
+        done cqblocks;
+
+    cqblocks = array (map2 do octlist octave:
+        d = emptyHops * (pow 2 (octaves-octave)) - emptyHops;
+//        eprintln "dropping \(d)";
+        drop d octlist;
+    done cqblocks [1..octaves]);
+
+    assembleBlock bits =
+       (//eprintln "assembleBlock: structure of bits is:";
+        //eprintln (map length bits);
+
+        rows = octaves * kdata.binsPerOctave;
+        columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
+
+        cm.generate do row col:
+
+            // bits structure: [1,2,4,8,...]
+
+            // each elt of bits is a list of the chunks that should
+            // make up this block in that octave (lowest octave first)
+
+            // each chunk has atomsPerFrame * binsPerOctave elts in it
+
+            // row is disposed with 0 at the top, highest octave (in
+            // both pitch and index into bits structure)
+
+            oct = int (row / binsPerOctave);
+            binNo = row % kdata.binsPerOctave;
+
+            chunks = pow 2 oct;
+            colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
+            atomNo = int (col / colsPerAtom);
+            atomOffset = col % colsPerAtom;
+
+            if atomOffset == 0 and atomNo < length bits[oct] then
+                bits[oct][atomNo][binNo];
+            else
+                cplx.zero
+            fi;
+
+        done { rows, columns };
+        );
+
+    assembleBlockSpectrogram bits =
+       (// As assembleBlock, but producing a dense magnitude
+        // spectrogram (rather than a complex output). (todo:
+        // interpolation, smoothing)
+
+        //eprintln "assembleBlockSpectrogram: structure of bits is:";
+        //eprintln (map length bits);
+
+        rows = octaves * kdata.binsPerOctave;
+        columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
+
+        mat.generate do row col:
+
+            oct = int (row / binsPerOctave);
+            binNo = row % kdata.binsPerOctave;
+
+            chunks = pow 2 oct;
+            colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
+            atomNo = int (col / colsPerAtom);
+
+            if atomNo < length bits[oct] then
+                cplx.magnitude bits[oct][atomNo][binNo];
+            else 
+                0
+            fi;
+
+        done { rows, columns };
+        );
+
+    processOctaveLists assembler octs =
+        case octs[0] of
+        block::rest:
+           (toAssemble = array 
+               (map do oct:
+                    n = kdata.atomsPerFrame * pow 2 oct;
+                    if not empty? octs[oct] then
+                        forBlock = array (take n octs[oct]);
+                        octs[oct] := drop n octs[oct];
+                        forBlock
+                    else
+                        array []
+                    fi
+                done (keys octs));
+            assembler toAssemble :. \(processOctaveLists assembler octs));
+         _: []
+        esac;
+
+//eprintln "cqblocks has \(length cqblocks) entries";
+
+    octaveLists = [:];
+
+    cqblocks = array cqblocks;
+    for [1..octaves] do oct:
+        octaveLists[octaves - oct] := cqblocks[oct-1];
+    done;
+/*
+    \() (map2 do octlist octave:
+println "oct \(octaves) - \(octave) = \(octaves - octave)";
+             octaveLists[octaves - octave] := octlist 
+         done cqblocks [1..octaves]);
+*/
+//eprintln "octaveLists keys are: \(keys octaveLists)";
+
+    {
+        kernel = kdata with {
+            binFrequencies = array
+               (concatMap do octave:
+                    map do freq:
+                        freq / (pow 2 octave);
+                    done (reverse (list kdata.binFrequencies))
+                done [0..octaves-1])
+        },
+        sampleRate,
+        octaves,
+        get cqComplex () = processOctaveLists assembleBlock octaveLists,
+        get cqSpectrogram () = processOctaveLists assembleBlockSpectrogram octaveLists,
+    }
+    );
+
+{ cqt }
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/cqtkernel.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,158 @@
+/*
+    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.
+*/
+
+module cqtkernel;
+
+vec = load may.vector;
+complex = load may.complex;
+window = load may.signal.window;
+fft = load may.transform.fft;
+cm = load may.matrix.complex;
+
+{ pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc;
+
+makeKernel { sampleRate, maxFreq, binsPerOctave } =
+   (q = 1;
+    atomHopFactor = 0.25;
+    thresh = 0.0005;
+    minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave));
+    bigQ = q / ((pow 2 (1/binsPerOctave)) - 1);
+
+    maxNK = round(bigQ * sampleRate / minFreq);
+    minNK = round(bigQ * sampleRate /
+                  (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave))));
+
+    atomHop = round(minNK * atomHopFactor);
+    
+    firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop));
+    
+    fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2));
+    
+//    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";
+//    eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)";
+    
+    winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1;
+    
+    lastCentre = firstCentre + (winNr - 1) * atomHop;
+    
+    fftHop = (lastCentre + atomHop) - firstCentre;
+    
+//    eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)";
+    
+    fftFunc = fft.forward fftSize;
+    
+    // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
+    // function. We can't do that here; we need to generate real and imag
+    // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
+    
+    binFrequencies = array [];
+
+    kernels = map do k:
+    
+        nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave))));
+
+        // the cq MATLAB toolbox uses a symmetric window for
+        // blackmanharris -- which is odd because it uses a periodic one
+        // for other types. Oh well
+        win = vec.divideBy nk
+           (vec.sqrt
+               (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk));
+
+        fk = minFreq * (pow 2 ((k-1)/binsPerOctave));
+
+        push binFrequencies fk;
+
+        genKernel f = vec.multiply
+           [win,
+            vec.fromList
+               (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])];
+    
+        reals = genKernel cos;
+        imags = genKernel sin;
+    
+        atomOffset = firstCentre - ceil(nk/2);
+    
+        map do i:
+    
+            shift = vec.zeros (atomOffset + ((i-1) * atomHop));
+    
+            specKernel = fftFunc
+               (complex.complexArray
+                   (vec.concat [shift, reals])
+                   (vec.concat [shift, imags]));
+
+            map do c:
+                if complex.magnitude c <= thresh then complex.zero else c fi
+                done specKernel;
+
+        done [1..winNr];
+    
+    done [1..binsPerOctave];
+    
+    kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels)));
+    
+//    eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))";
+    
+    // Normalisation
+    
+    wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat));
+    wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat));
+
+    subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1));
+    square = cm.product (cm.conjugateTransposed subset) subset;
+
+    diag = complex.magnitudes (cm.getDiagonal 0 square);
+    wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2);
+
+    weight = (fftHop / fftSize) / (vec.mean (vec.abs wK));
+    weight = sqrt(weight);
+
+    kernel = cm.scaled weight kmat;
+
+//    eprintln "weight = \(weight)";
+    
+    {
+        kernel,
+        fftSize,
+        fftHop,
+        binsPerOctave,
+        atomsPerFrame = winNr,
+        atomSpacing = atomHop,
+        firstCentre,
+        maxFrequency = maxFreq,
+        minFrequency = minFreq,
+        binFrequencies,
+        bigQ
+    });
+
+{
+    makeKernel
+}
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/experiment.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,50 @@
+
+program experiment;
+
+af = load may.stream.audiofile;
+plot = load may.plot;
+cm = load may.matrix.complex;
+mat = load may.matrix;
+vec = load may.vector;
+manipulate = load may.stream.manipulate;
+syn = load may.stream.syntheticstream;
+
+{ cqt } = load cqt;
+{ icqt } = load icqt;
+
+//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
+//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
+testStream = af.open "sweep-48000.wav";
+//testStream = af.open "sweep.wav";
+
+eprintln "have test stream";
+
+cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream;
+
+eprintln "bin frequencies: \(cq.kernel.binFrequencies)";
+
+icq = icqt cq;
+
+eprintln "calculated icq...";
+
+outfile = "test-output.wav";
+written = af.write icq outfile;
+
+eprintln "wrote \(written) to \(outfile)";
+
+//for cq.cqSpectrogram do mm:
+//    for (mat.asColumns mm) (println . strJoin "," . vec.list);
+//done;
+
+/*
+bigM = mat.concatHorizontal (map cm.magnitudes cq.output);
+
+eprintln "overall output size = \(mat.size bigM)";
+
+mat.print bigM;
+*/
+//\() (plot.plot [Contour bigM]);
+//\() (plot.plot [Grid bigM]);
+
+()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/icqt.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,145 @@
+/*
+    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.
+*/
+
+module icqt;
+
+cqt = load cqt;
+cm = load may.matrix.complex;
+framer = load may.stream.framer;
+win = load may.signal.window;
+mm = load may.mathmisc;
+vec = load may.vector;
+resamp = load may.stream.resample;
+manip = load may.stream.manipulate;
+mat = load may.matrix;
+
+icqt cq =
+   (kdata = cq.kernel;
+
+    // kdata.kernel is the kernel matrix for a single octave. It has
+    // width kdata.fftSize and height kdata.binsPerOctave *
+    // kdata.atomsPerFrame.
+    //
+    // kdata.fftHop is the overlap between kernel matrices in a single
+    // octave.
+    // 
+    // cq.sampleRate is the output stream sample rate and cq.octaves
+    // the number of octaves.
+    //
+    // cq.cqComplex is the list of complex matrices containing the CQ
+    // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1)
+    // and height kdata.binsPerOctave * cq.octaves.
+
+    bpo = kdata.binsPerOctave;
+    atomsPerFrame = kdata.atomsPerFrame;
+    octaves = cq.octaves;
+
+    // transform a single block, all octaves tall, into an array
+    // (indexed by octave) of lists of individual columns (valid
+    // values for that octave only)
+    decomposeOctaves mat = array
+       (map do oct:
+            inverted = cm.fromRows (reverse (cm.asRows mat));
+            octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1));
+            gap = (mm.pow 2 oct) - 1;
+            pickFrom cols =
+                case cols of
+                c::cs: c::(pickFrom (drop gap cs));
+                _: [];
+                esac;
+            pickFrom (cm.asColumns octMat);
+        done [0..octaves-1]);
+
+    // transform a list of the arrays produced by decomposeOctaves
+    // into a single array (indexed by octave) of lists of the
+    // individual columns
+    flattenOctaves decomposed =
+       (flattenAux acc decomposed = 
+            case decomposed of
+            chunk::rest:
+                flattenAux 
+                   (array
+                       (map do oct:
+                            acc[oct] ++ chunk[oct]
+                            done [0..octaves-1]))
+                    rest;
+            _: acc;
+            esac;
+        flattenAux (array (map \[] [0..octaves-1])) decomposed);
+
+    // given a list of columns, deinterleave and pile up each sequence
+    // of atomsPerFrame columns into a single tall column and return
+    // the resulting list of tall columns
+    pile cols =
+       (pileAux acc cols =
+            if (length cols) >= atomsPerFrame then
+                atoms = take atomsPerFrame cols;
+                juggled = concat (reverse (cm.asRows (cm.fromColumns atoms)));
+                pileAux (juggled :: acc) (drop atomsPerFrame cols);
+            else
+                reverse acc
+            fi;
+        pileAux [] cols);
+
+    octaveColumnLists = 
+        map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex)));
+
+    for octaveColumnLists do l: println "octave column list length: \(length l)" done;
+
+    kernel = cm.transposed kdata.kernel; // right way around for the multiply
+    
+    spectra = 
+        map do l:
+            map do col:
+                res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col));
+                cm.columnSlice res 0 (kdata.fftSize / 2 + 1)
+            done l;
+        done octaveColumnLists;
+
+    eprintln "calculated spectra, now to ifft, overlap-add...";
+
+    rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1];
+    
+    resynthesised =
+        map2 do frames rate:
+            framer.complexStreamed rate kdata.fftSize
+                [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ]
+                frames;
+        done spectra rates;
+
+   eprintln "have streams:";
+   for resynthesised eprintln;
+
+   resampled = map (resamp.resampledTo cq.sampleRate) resynthesised;
+
+   manip.sum resampled;
+);
+
+{icqt}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/nbproject/ide-file-targets.xml	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project basedir=".." name="cqt-IDE">
+    <!-- TODO: edit the following target according to your needs -->
+    <!-- (more info: http://www.netbeans.org/kb/articles/freeform-config.html#debugj2se) -->
+    <property name="maydir" value="${basedir}/../../may"/>
+    <property name="yetidir" value="${basedir}/../../yeti"/>
+    <property name="jardir" value="${maydir}/ext/jar"/>
+        <path id="cp">
+            <fileset dir="${basedir}">
+                <include name="cqt.jar"/>
+            </fileset>
+            <fileset dir="${jardir}">
+                <include name="*.jar"/>
+            </fileset>
+            <fileset dir="${yetidir}">
+                <include name="yeti.jar"/>
+            </fileset>
+            <fileset dir="${maydir}">
+                <include name="may.jar"/>
+            </fileset>
+        </path>
+    <target name="debug-nb">
+        <nbjpdastart addressproperty="jpda.address" name="cqt" transport="dt_socket">
+            <classpath refid="cp"/>
+        </nbjpdastart>
+        <!-- TODO configure the main class for your project here: -->
+        <java classname="plotfile" fork="true">
+            <classpath refid="cp"/>
+            <jvmarg value="-Xdebug"/>
+            <jvmarg value="-Xrunjdwp:transport=dt_socket,address=${jpda.address}"/>
+            <arg value="${basedir}/test.wav"/>
+        </java>
+    </target>
+    <target name="-profile-check">
+        <startprofiler freeform="true"/>
+    </target>
+    <!-- TODO: edit the following target according to your needs -->
+    <!-- (more info: http://netbeans.org/kb/articles/freeform-config.html#profilej2se) -->
+    <target depends="-profile-check" if="profiler.configured" name="profile-nb">
+        <!-- TODO configure the main class for your project here: -->
+        <java classname="plotfile" fork="true">
+            <classpath refid="cp"/>
+            <jvmarg line="${agent.jvmargs}"/>
+            <arg value="${basedir}/test.wav"/>
+        </java>
+    </target>
+</project>

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/nbproject/project.xml	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,53 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project xmlns="http://www.netbeans.org/ns/project/1">
+    <type>org.netbeans.modules.ant.freeform</type>
+    <configuration>
+        <general-data xmlns="http://www.netbeans.org/ns/freeform-project/1">
+            <!-- Do not use Project Properties customizer when editing this file manually. 
+ To prevent the customizer from showing, create nbproject/project.properties file and enter 
+auxiliary.show.customizer=false 
+property there. Adding 
+auxiliary.show.customizer.message=<message>
+ will show your customized message when someone attempts to open the customizer.  -->
+            <name>cqt</name>
+            <properties/>
+            <folders/>
+            <ide-actions>
+                <action name="build">
+                    <target>jar</target>
+                </action>
+                <action name="clean">
+                    <target>clean</target>
+                </action>
+                <action name="rebuild">
+                    <target>clean</target>
+                    <target>jar</target>
+                </action>
+                <action name="debug">
+                    <script>nbproject/ide-file-targets.xml</script>
+                    <target>debug-nb</target>
+                </action>
+                <action name="profile">
+                    <script>nbproject/ide-file-targets.xml</script>
+                    <target>profile-nb</target>
+                </action>
+            </ide-actions>
+            <view>
+                <items>
+                    <source-file>
+                        <location>build.xml</location>
+                    </source-file>
+                </items>
+                <context-menu>
+                    <ide-action name="build"/>
+                    <ide-action name="rebuild"/>
+                    <ide-action name="clean"/>
+                    <ide-action name="debug"/>
+                    <ide-action name="profile"/>
+                </context-menu>
+            </view>
+            <subprojects/>
+        </general-data>
+        <java-data xmlns="http://www.netbeans.org/ns/freeform-project-java/1"/>
+    </configuration>
+</project>

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/plotfile.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,41 @@
+
+program plotfile;
+
+af = load may.stream.audiofile;
+plot = load may.plot;
+cm = load may.matrix.complex;
+mat = load may.matrix;
+
+{ cqt } = load cqt;
+
+minFreq = 27.5;
+maxRateDivisor = 3;
+binsPerOctave = 60;
+
+plotfile f =
+   (testStream = af.open f;
+    eprintln "Opened file stream...";
+    start = System#currentTimeMillis();
+    cq = cqt {
+        maxFreq = testStream.sampleRate/maxRateDivisor, 
+        minFreq, 
+        binsPerOctave
+        } testStream;
+    middle = System#currentTimeMillis();
+    eprintln "Generated kernel \(cm.size cq.kernel.kernel) and primed transform (\(cq.octaves) octaves), took \(middle-start)ms, now calculating...";
+    bigM = mat.concatHorizontal cq.cqSpectrogram;
+    finish = System#currentTimeMillis();
+    eprintln "Done, that part took \(finish-middle)ms, all CQ stuff took \(finish-start)ms";
+    eprintln "Plotting...";
+    \() (plot.plot [Contour bigM]));
+
+usage () =
+   (eprintln "\nUsage: plotfile file.wav";
+    eprintln "\n   Loads audio from file.wav and plots a \(binsPerOctave)bpo Constant-Q spectrogram";
+    eprintln "   from \(minFreq)Hz up to 1/\(maxRateDivisor) of the file samplerate");
+
+case (list _argv) of
+file::[]: plotfile file;
+_: usage ();
+esac;
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/test.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,96 @@
+/*
+    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.
+*/
+
+program test;
+
+{ runTests } = load may.test;
+
+// We want to test:
+// 
+// Kernel design -- check size (number of bins, number of atoms);
+// check an example kernel against known data
+//
+// Time alignment -- feed a dirac train, check that peaks in all bins
+// align
+//
+// Frequency discrimination -- feed a sinusoid, check peaks
+//
+// Latency compensation -- for dirac at 0, check peak can be found at
+// 0 plus the declared latency
+//
+// Signal-noise ratio
+//
+// Specimen output for simple test case
+
+tests = [
+"kernel"    : load test_cqtkernel,
+"frequency" : load test_frequency,
+];
+
+bad = sum (mapHash do name testHash: runTests name testHash done tests);
+
+if (bad > 0) then
+    println "\n** \(bad) test(s) failed!";
+    threadExit 1
+else
+    ()
+fi
+
+
+
+
+/*
+//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
+//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
+testStream = af.open "sweep-48000.wav";
+//testStream = af.open "sweep.wav";
+
+// So the stream is [ 0, 1, 0, -1, 0, 1, 0, -1, ... ] :
+//testStream = manipulate.withDuration 64 (syn.sinusoid 8 2);
+
+testStream = manipulate.withDuration 32 (syn.pulseTrain 8 0.001);
+
+eprintln "have test stream";
+
+cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream;
+
+eprintln "bin frequencies: \(cq.kernel.binFrequencies)";
+
+bigM = mat.concatHorizontal (map cm.magnitudes cq.output);
+
+eprintln "overall output size = \(mat.size bigM)";
+
+mat.print bigM;
+
+//\() (plot.plot [Contour bigM]);
+\() (plot.plot [Grid bigM]);
+*/
+//()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/test_cqtkernel.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,63 @@
+/*
+    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.
+*/
+
+module test_cqtkernel;
+
+cm = load may.matrix.complex;
+mm = load may.mathmisc;
+
+{ compare, compareUsing } = load may.test;
+
+{ makeKernel } = load cqtkernel;
+
+eps = 1e-7;
+
+compareClose = compareUsing do a b: abs (a - b) < eps done;
+
+[
+
+"minimal": \(
+    k = makeKernel { sampleRate = 16, maxFreq = 8, binsPerOctave = 4 };
+    compare k.binsPerOctave 4 and
+        compare (cm.size k.kernel) {
+            rows = k.binsPerOctave * k.atomsPerFrame,
+            columns = k.fftSize
+        } and
+        compareClose k.maxFrequency 8 and
+        compareClose k.minFrequency (4 * (mm.pow 2 (1/4))) and
+        compare k.atomsPerFrame 5 and
+        compare k.fftSize 32 and
+        compare (length k.binFrequencies) k.binsPerOctave and
+        compareClose (head k.binFrequencies) k.minFrequency and
+        compareClose (head (reverse k.binFrequencies)) k.maxFrequency
+),
+
+] is hash<string, () -> boolean>
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/misc/yeti/test_frequency.yeti	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,132 @@
+/*
+    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.
+*/
+
+module test_frequency;
+
+mat = load may.matrix;
+vec = load may.vector;
+win = load may.signal.window;
+mm = load may.mathmisc;
+cm = load may.matrix.complex;
+syn = load may.stream.syntheticstream;
+plot = load may.plot;
+
+{ compare } = load may.test;
+
+{ cqt } = load cqt;
+
+// Test with a single windowed sinusoid, repeating at various frequencies
+
+sinTestStream sampleRate duration signalFreq = // duration is in samples
+   (sin = syn.sinusoid sampleRate signalFreq;
+    chunk = mat.getRow 0 (sin.read duration);
+    syn.precalculatedMono sampleRate (win.windowed win.hann chunk));
+
+// We want to make a CQ transform spanning more than one octave, but
+// not going all the way to fs/2 so we can test it also with
+// frequencies above and below its extents
+
+sampleRate = 100;
+
+// fs/2 = 50 so 10->40 gives us 2 octaves
+cqmin = 10;
+cqmax = 40;
+bpo = 4; // fairly arbitrary
+
+testFreqs = map (* 5) [ 0..10 ];
+duration = sampleRate * 2;
+
+threshold = 0.08;
+
+streamBuilder = sinTestStream sampleRate duration;
+
+binForFreq f =
+    mm.round (bpo * mm.log2 (f / cqmin)) - 1;
+
+report message matrix =
+   (eprintln message;
+    eprintln "matrix is:";
+    mat.eprint matrix);
+//    chart = plot.plot [Grid matrix];
+//    sleep 100;
+//    chart#dispose());
+
+tests = mapIntoHash
+    do f: "freq_\(f)" done
+    do f: \(
+        str = streamBuilder f;
+        cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str;
+        spec = cq.cqSpectrogram;
+        rightSize = all id
+           (map do s:
+            compare (mat.size s) { 
+                rows = cq.kernel.binsPerOctave * cq.octaves,
+                columns = cq.kernel.atomsPerFrame * mm.pow 2 (cq.octaves - 1)
+            }
+            done spec);
+        m = mat.concatHorizontal spec;
+//    println "binFrequencies = \(cq.kernel.binFrequencies)";
+//    println "binForFreq \(f) = \(binForFreq f)";
+        var colno = 0;
+        success = all id
+           (rightSize :: map do c:
+                // The test passes for this column if:
+                //
+                //  * the max bin is the expected one, or
+                //
+                //  * the expected max is out of range entirely (but
+                //  we need to test _something_ in this case --
+                //  what?), or
+                //
+                //  * all bins are below a threshold, or
+                //
+                //  * this is an odd column and the expected max is in
+                //  the lower octave
+                //
+                // We should also check that all values in the lower
+                // octave are zero for odd columns.
+                //
+                expected = binForFreq f;
+                good =
+                   (expected < 0 or expected >= vec.length c) or
+                   ((colno % 2 == 1) and expected < (vec.length c / 2)) or
+                   (vec.max c < threshold) or
+                   (vec.maxindex c == binForFreq f);
+                if (not good) then
+                    report " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)" m;
+                fi;
+                colno := colno + 1;
+                good;
+            done (mat.asColumns m));
+        success;
+    ) done
+    testFreqs;
+
+tests is hash<string, () -> boolean>

Binary file paper/smc2010.pdf has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/CQInverse.cpp	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,411 @@
+/* -*- 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, 50, 0.05);
+
+#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 j = 0; j < m_octaves; ++j) {
+        int factor = pow(2, j);
+        int latency = (j > 0 ? m_upsamplers[j]->getLatency() : 0) / factor;
+        for (int i = 0; i < (latency + m_p.fftSize) / m_p.fftHop; ++i) {
+            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/src/CQKernel.cpp	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,337 @@
+/* -*- 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 "CQKernel.h"
+
+#include "maths/MathUtilities.h"
+#include "dsp/transforms/FFT.h"
+#include "base/Window.h"
+
+#include <cmath>
+#include <cassert>
+#include <vector>
+#include <iostream>
+#include <algorithm>
+
+using std::vector;
+using std::complex;
+using std::cerr;
+using std::endl;
+
+typedef std::complex<double> C;
+
+CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave) :
+    m_fft(0)
+{
+    m_p.sampleRate = sampleRate;
+    m_p.maxFrequency = maxFreq;
+    m_p.binsPerOctave = binsPerOctave;
+    generateKernel();
+}
+
+CQKernel::~CQKernel()
+{
+    delete m_fft;
+}
+
+void
+CQKernel::generateKernel()
+{
+    double q = 1;
+    double atomHopFactor = 0.25;
+    double thresh = 0.0005;
+
+    double bpo = m_p.binsPerOctave;
+
+    m_p.minFrequency = (m_p.maxFrequency / 2) * pow(2, 1.0/bpo);
+    m_p.Q = q / (pow(2, 1.0/bpo) - 1.0);
+
+    double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency);
+    double minNK = round
+        (m_p.Q * m_p.sampleRate /
+         (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo)));
+
+    if (minNK == 0 || maxNK == 0) {
+        // most likely pathological parameters of some sort
+        cerr << "WARNING: CQKernel::generateKernel: minNK or maxNK is zero (minNK == " << minNK << ", maxNK == " << maxNK << "), not generating a kernel" << endl;
+        m_p.atomSpacing = 0;
+        m_p.firstCentre = 0;
+        m_p.fftSize = 0;
+        m_p.atomsPerFrame = 0;
+        m_p.lastCentre = 0;
+        m_p.fftHop = 0;
+        return;
+    }
+
+    m_p.atomSpacing = round(minNK * atomHopFactor);
+    m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing);
+    m_p.fftSize = MathUtilities::nextPowerOfTwo
+        (m_p.firstCentre + ceil(maxNK / 2.0));
+
+    m_p.atomsPerFrame = floor
+        (1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing);
+
+    cerr << "atomsPerFrame = " << m_p.atomsPerFrame << " (atomHopFactor = " << atomHopFactor << ", atomSpacing = " << m_p.atomSpacing << ", fftSize = " << m_p.fftSize << ", maxNK = " << maxNK << ", firstCentre = " << m_p.firstCentre << ")" << endl;
+
+    m_p.lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing;
+
+    m_p.fftHop = (m_p.lastCentre + m_p.atomSpacing) - m_p.firstCentre;
+
+    cerr << "fftHop = " << m_p.fftHop << endl;
+
+    m_fft = new FFT(m_p.fftSize);
+
+    for (int k = 1; k <= m_p.binsPerOctave; ++k) {
+        
+        int nk = round(m_p.Q * m_p.sampleRate /
+                       (m_p.minFrequency * pow(2, ((k-1.0) / bpo))));
+
+        // The MATLAB version uses a symmetric window, but our windows
+        // are periodic. A symmetric window of size N is a periodic
+        // one of size N-1 with the first element stuck on the end
+        Window<double> w(BlackmanHarrisWindow, nk-1);
+        vector<double> win = w.getWindowData();
+        win.push_back(win[0]);
+
+        for (int i = 0; i < (int)win.size(); ++i) {
+            win[i] = sqrt(win[i]) / nk;
+        }
+
+        double fk = m_p.minFrequency * pow(2, ((k-1.0) / bpo));
+
+        vector<double> reals, imags;
+        
+        for (int i = 0; i < nk; ++i) {
+            double arg = (2.0 * M_PI * fk * i) / m_p.sampleRate;
+            reals.push_back(win[i] * cos(arg));
+            imags.push_back(win[i] * sin(arg));
+        }
+
+        int atomOffset = m_p.firstCentre - int(ceil(nk/2.0));
+
+        for (int i = 0; i < m_p.atomsPerFrame; ++i) {
+
+            int shift = atomOffset + (i * m_p.atomSpacing);
+
+            vector<double> rin(m_p.fftSize, 0.0);
+            vector<double> iin(m_p.fftSize, 0.0);
+
+            for (int j = 0; j < nk; ++j) {
+                rin[j + shift] = reals[j];
+                iin[j + shift] = imags[j];
+            }
+
+            vector<double> rout(m_p.fftSize, 0.0);
+            vector<double> iout(m_p.fftSize, 0.0);
+
+            m_fft->process(false,
+                           rin.data(), iin.data(),
+                           rout.data(), iout.data());
+
+            // Keep this dense for the moment (until after
+            // normalisation calculations)
+
+            vector<C> row;
+
+            for (int j = 0; j < m_p.fftSize; ++j) {
+                if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) {
+                    row.push_back(C(0, 0));
+                } else {
+                    row.push_back(C(rout[j] / m_p.fftSize,
+                                    iout[j] / m_p.fftSize));
+                }
+            }
+
+            m_kernel.origin.push_back(0);
+            m_kernel.data.push_back(row);
+        }
+    }
+
+    assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+
+    // print density as diagnostic
+
+    int nnz = 0;
+    for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
+        for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) {
+            if (m_kernel.data[i][j] != C(0, 0)) {
+                ++nnz;
+            }
+        }
+    }
+
+    cerr << "size = " << m_kernel.data.size() << "*" << m_kernel.data[0].size() << " (fft size = " << m_p.fftSize << ")" << endl;
+
+    assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+    assert((int)m_kernel.data[0].size() == m_p.fftSize);
+
+    cerr << "density = " << double(nnz) / double(m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize) << " (" << nnz << " of " << m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize << ")" << endl;
+
+    finaliseKernel();
+}
+
+static bool ccomparator(C &c1, C &c2)
+{
+    return abs(c1) < abs(c2);
+}
+
+static int maxidx(vector<C> &v)
+{
+    return std::max_element(v.begin(), v.end(), ccomparator) - v.begin();
+}
+
+void
+CQKernel::finaliseKernel()
+{
+    // calculate weight for normalisation
+
+    int wx1 = maxidx(m_kernel.data[0]);
+    int wx2 = maxidx(m_kernel.data[m_kernel.data.size()-1]);
+
+    vector<vector<C> > subset(m_kernel.data.size());
+    for (int j = wx1; j <= wx2; ++j) {
+        for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
+            subset[i].push_back(m_kernel.data[i][j]);
+        }
+    }
+
+    int nrows = subset.size();
+    int ncols = subset[0].size();
+    vector<vector<C> > square(ncols); // conjugate transpose of subset * subset
+
+    for (int i = 0; i < nrows; ++i) {
+        assert((int)subset[i].size() == ncols);
+    }
+
+    for (int j = 0; j < ncols; ++j) {
+        for (int i = 0; i < ncols; ++i) {
+            C v(0, 0);
+            for (int k = 0; k < nrows; ++k) {
+                v += subset[k][i] * conj(subset[k][j]);
+            }
+            square[i].push_back(v);
+        }
+    }
+
+    vector<double> wK;
+    double q = 1; //!!! duplicated from constructor
+    for (int i = round(1.0/q); i < ncols - round(1.0/q) - 2; ++i) {
+        wK.push_back(abs(square[i][i]));
+    }
+
+    double weight = double(m_p.fftHop) / m_p.fftSize;
+    weight /= MathUtilities::mean(wK.data(), wK.size());
+    weight = sqrt(weight);
+    
+    cerr << "weight = " << weight << endl;
+
+    // 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;
+
+    for (int i = 0; i < (int)m_kernel.data.size(); ++i) {
+
+        sk.origin.push_back(0);
+        sk.data.push_back(vector<C>());
+
+        int lastNZ = 0;
+        for (int j = (int)m_kernel.data[i].size()-1; j >= 0; --j) {
+            if (abs(m_kernel.data[i][j]) != 0.0) {
+                lastNZ = j;
+                break;
+            }
+        }
+
+        bool haveNZ = false;
+        for (int j = 0; j <= lastNZ; ++j) {
+            if (haveNZ || abs(m_kernel.data[i][j]) != 0.0) {
+                if (!haveNZ) sk.origin[i] = j;
+                haveNZ = true;
+                sk.data[i].push_back(conj(m_kernel.data[i][j]) * weight);
+            }
+        }
+    }
+
+    m_kernel = sk;
+}
+
+vector<C>
+CQKernel::processForward(const vector<C> &cv)
+{
+    // straightforward matrix multiply (taking into account m_kernel's
+    // slightly-sparse representation)
+
+    if (m_kernel.data.empty()) return vector<C>();
+
+    int nrows = m_p.binsPerOctave * m_p.atomsPerFrame;
+
+    vector<C> rv(nrows, C());
+
+    for (int i = 0; i < nrows; ++i) {
+        int len = m_kernel.data[i].size();
+        for (int j = 0; j < len; ++j) {
+            rv[i] += cv[j + m_kernel.origin[i]] * m_kernel.data[i][j];
+        }
+    }
+
+    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.
+
+    if (m_kernel.data.empty()) return vector<C>();
+
+    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;
+}
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/CQSpectrogram.cpp	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,261 @@
+/* -*- 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;
+	    }
+	}
+    }
+
+//    cerr << "fetchLinear: firstFullHeight = " << firstFullHeight << ", secondFullHeight = " << secondFullHeight << endl;
+
+    if (firstFullHeight < 0) {
+	if (insist) {
+            return fetchHold(true);
+	} 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) {
+            return fetchHold(true);
+	} 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/src/ConstantQ.cpp	Thu May 15 12:04:00 2014 +0100
@@ -0,0 +1,352 @@
+/* -*- 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 "ConstantQ.h"
+
+#include "CQKernel.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
+
+ConstantQ::ConstantQ(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();
+}
+
+ConstantQ::~ConstantQ()
+{
+    delete m_fft;
+    for (int i = 0; i < (int)m_decimators.size(); ++i) {
+        delete m_decimators[i];
+    }
+    delete m_kernel;
+}
+
+double
+ConstantQ::getMinFrequency() const
+{
+    return m_p.minFrequency / pow(2.0, m_octaves - 1);
+}
+
+double
+ConstantQ::getBinFrequency(int bin) const
+{
+    return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave()));
+}
+
+void
+ConstantQ::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_decimators.push_back(0);
+
+    for (int i = 1; i < m_octaves; ++i) {
+
+        int factor = pow(2, i);
+
+        Resampler *r = new Resampler
+            (sourceRate, sourceRate / factor, 50, 0.05);
+
+#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.
+        // 
+        // Our decimator uses a linear phase filter, but being causal
+        // it is not zero phase: it has a latency that depends on the
+        // decimation factor. Those latencies have been calculated
+        // per-octave and are available to us in the latencies
+        // array. Left to its own devices, the first input sample will
+        // appear at output sample 0 in the highest octave (where no
+        // decimation is needed), sample number latencies[1] in the
+        // next octave down, latencies[2] in the next one, etc. We get
+        // to apply some artificial per-octave latency after the
+        // decimator in the processing chain, in order to compensate
+        // for the differing latencies associated with different
+        // decimation factors. How much should we insert?
+        //
+        // The outputs of the decimators are at different rates (in
+        // terms of the relation between clock time and samples) and
+        // we want them aligned in terms of time. So, for example, a
+        // latency of 10 samples with a decimation factor of 2 is
+        // equivalent to a latency of 20 with no decimation -- they
+        // both result in the first output sample happening at the
+        // same equivalent time in milliseconds.
+	// 
+	// So here we record the latency added by the decimator, in
+	// terms of the sample rate of the undecimated signal. Then we
+	// use that to compensate in a moment, when we've discovered
+	// what the longest latency across all octaves is.
+
+        latencies.push_back(r->getLatency() * factor);
+        m_decimators.push_back(r);
+    }
+
+    m_bigBlockSize = m_p.fftSize * pow(2, m_octaves - 1);
+
+    // Now add in the extra padding and compensate for hops that must
+    // be dropped in order to align the atom centres across
+    // octaves. Again this is a bit trickier because we are doing it
+    // at input rather than output and so must work in per-octave
+    // sample rates rather than output blocks
+
+    int emptyHops = m_p.firstCentre / m_p.atomSpacing;
+
+    vector<int> drops;
+    for (int i = 0; i < m_octaves; ++i) {
+	int factor = pow(2, i);
+	int dropHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops;
+	int drop = ((dropHops * m_p.fftHop) * factor) / m_p.atomsPerFrame;
+	drops.push_back(drop);
+    }
+
+    int maxLatPlusDrop = 0;
+    for (int i = 0; i < m_octaves; ++i) {
+	int latPlusDrop = latencies[i] + drops[i];
+	if (latPlusDrop > maxLatPlusDrop) maxLatPlusDrop = latPlusDrop;
+    }
+
+    int totalLatency = maxLatPlusDrop;
+
+    int lat0 = totalLatency - latencies[0] - drops[0];
+    totalLatency = ceil(double(lat0 / m_p.fftHop) * m_p.fftHop)
+	+ latencies[0] + drops[0];
+
+    // 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);
+
+#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) {
+
+	double factor = pow(2, i);
+
+	// Calculate the difference between the total latency applied
+	// across all octaves, and the existing latency due to the
+	// decimator for this octave, and then convert it back into
+	// the sample rate appropriate for the output latency of this
+	// decimator -- including one additional big block of padding
+	// (as in the reference).
+
+	double octaveLatency =
+	    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
+            (RealSequence(int(round(octaveLatency)), 0.0));
+    }
+
+    m_fft = new FFTReal(m_p.fftSize);
+}
+
+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) {
+        RealSequence dec = m_decimators[i]->process(td.data(), td.size());
+        m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end());
+    }
+
+    ComplexBlock out;
+
+    while (true) {
+
+	// We could have quite different remaining sample counts in
+	// different octaves, because (apart from the predictable
+	// added counts for decimator output on each block) we also
+	// have variable additional latency per octave
+	bool enough = true;
+	for (int i = 0; i < m_octaves; ++i) {
+	    int required = m_p.fftSize * pow(2, m_octaves - i - 1);
+	    if ((int)m_buffers[i].size() < required) {
+		enough = false;
+	    }
+	}
+	if (!enough) break;
+
+        int base = out.size();
+        int totalColumns = pow(2, m_octaves - 1) * m_p.atomsPerFrame;
+        for (int i = 0; i < totalColumns; ++i) {
+            out.push_back(ComplexColumn());
+        }
+
+        for (int octave = 0; octave < m_octaves; ++octave) {
+
+            int blocksThisOctave = pow(2, (m_octaves - octave - 1));
+
+            for (int b = 0; b < blocksThisOctave; ++b) {
+                ComplexBlock block = processOctaveBlock(octave);
+                
+                for (int j = 0; j < m_p.atomsPerFrame; ++j) {
+
+                    int target = base +
+			    (b * (totalColumns / blocksThisOctave) + 
+			     (j * ((totalColumns / blocksThisOctave) /
+				   m_p.atomsPerFrame)));
+
+                    while (int(out[target].size()) < 
+                           m_p.binsPerOctave * (octave + 1)) {
+                        out[target].push_back(Complex());
+                    }
+                    
+                    for (int i = 0; i < m_p.binsPerOctave; ++i) {
+                        out[target][m_p.binsPerOctave * octave + i] = 
+                            block[j][m_p.binsPerOctave - i - 1];
+                    }
+                }
+            }
+        }
+    }
+
+    return out;
+}
+
+ConstantQ::ComplexBlock
+ConstantQ::getRemainingOutput()
+{
+    // Same as padding added at start, though rounded up
+    int pad = ceil(double(m_outputLatency) / m_bigBlockSize) * m_bigBlockSize;
+    RealSequence zeros(pad, 0.0);
+    return process(zeros);
+}
+
+ConstantQ::ComplexBlock
+ConstantQ::processOctaveBlock(int octave)
+{
+    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());
+
+    RealSequence shifted;
+    shifted.insert(shifted.end(), 
+                   m_buffers[octave].begin() + m_p.fftHop,
+                   m_buffers[octave].end());
+    m_buffers[octave] = shifted;
+
+    ComplexSequence cv;
+    for (int i = 0; i < m_p.fftSize; ++i) {
+        cv.push_back(Complex(ro[i], io[i]));
+    }
+
+    ComplexSequence cqrowvec = m_kernel->processForward(cv);
+
+    // Reform into a column matrix
+    ComplexBlock cqblock;
+    for (int j = 0; j < m_p.atomsPerFrame; ++j) {
+        cqblock.push_back(ComplexColumn());
+        for (int i = 0; i < m_p.binsPerOctave; ++i) {
+            cqblock[j].push_back(cqrowvec[i * m_p.atomsPerFrame + j]);
+        }
+    }
+
+    return cqblock;
+}
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/processfile.cpp	Thu May 15 12:04:00 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/src/test.cpp	Thu May 15 12:04:00 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");
+	}
+    }
+}
+

--- a/yeti/Makefile	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-
-MAY_DIR	:= ../../may
-MAY_YC	:= $(MAY_DIR)/bin/yc
-
-test:
-	$(MAY_YC) ./test.yeti

--- a/yeti/build.xml	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,74 +0,0 @@
-<project name="cqt" default="test" basedir=".">
-
-  <property name="maydir" value="${basedir}/../../may"/>
-  <property name="yetidir" value="${basedir}/../../yeti"/>
-
-  <property name="jardir" value="${maydir}/ext/jar"/>
-
-  <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/jzy3d-swt-0.9.1.jar:${jardir}/jzy3d-api-0.9.1.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/>
-
-  <condition property="archtag" value="linux32">
-    <os family="unix" arch="i386"/>
-  </condition>
-  <condition property="archtag" value="linux64">
-    <os family="unix" arch="amd64"/>
-  </condition>
-  <condition property="archtag" value="win32">
-    <os family="windows" arch="x86"/>
-  </condition>
-  <condition property="archtag" value="win64">
-    <os family="windows" arch="amd64"/>
-  </condition>
-  <condition property="archtag" value="osx">
-    <os family="mac"/>
-  </condition>
-
-  <target name="taskdef">
-    <taskdef name="yetic" classname="yeti.lang.compiler.YetiTask"
-	     classpath="${yetidir}/yeti.jar:${maydir}/may.jar:${extjars}" />
-  </target>
-
-  <target name="prepare">
-    <mkdir dir="${basedir}/classes"/>
-  </target>
-
-  <target name="yeticlasses" depends="taskdef">
-    <yetic srcdir="${basedir}" 
-	   destdir="${basedir}/classes"
-	   includes="**/*.yeti"
-	   preload="yeti/lang/std:yeti/lang/io"/>
-  </target>
-
-  <target name="classes" depends="prepare,yeticlasses"/>
-
-  <target name="jar" depends="classes,taskdef">
-    <jar jarfile="${basedir}/cqt.jar">
-      <fileset dir="${basedir}/classes" 
-	       includes="**/*.class" 
-	       excludes="**/test*.class"/>
-    </jar>
-  </target>
-
-  <target name="testjar" depends="classes,taskdef">
-    <jar jarfile="${basedir}/test.jar">
-      <fileset dir="${basedir}/classes" 
-	       includes="**/test*.class"/>
-    </jar>
-  </target>
-
-  <target name="test" depends="jar,testjar,taskdef">
-    <java classpath="${basedir}/test.jar:${basedir}/cqt.jar:${maydir}/may.jar:${yetidir}/yeti.jar:${extjars}"
-	  classname="test"
-	  fork="true" failonerror="true">
-      <sysproperty key="java.library.path" path="${maydir}/ext/native/${archtag}"/>
-    </java>
-  </target>
-
-  <target name="clean">
-    <delete dir="${basedir}/classes"/>
-  </target>
-
-  <target name="rebuild" depends="clean,jar"/>
-
-</project>
-

--- a/yeti/cqt.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,260 +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.
-*/
-
-module cqt;
-
-cqtkernel = load cqtkernel;
-resample = load may.stream.resample;
-manipulate = load may.stream.manipulate;
-mat = load may.matrix;
-cm = load may.matrix.complex;
-framer = load may.stream.framer;
-cplx = load may.complex;
-fft = load may.transform.fft;
-vec = load may.vector;
-ch = load may.stream.channels;
-
-{ pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc;
-
-cqt { maxFreq, minFreq, binsPerOctave } str =
-   (sampleRate = str.sampleRate;
-    octaves = ceil (log2 (maxFreq / minFreq));
-//    actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave));
-
-    kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave };
-
-//    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave), fftSize = \(kdata.fftSize), hop = \(kdata.fftHop)";
-
-//    eprintln "atomsPerFrame = \(kdata.atomsPerFrame)";
-
-    padding = (kdata.fftSize * (pow 2 (octaves-1)));
-
-//    eprintln "padding = \(padding)";
-
-    str = manipulate.paddedBy padding str;
-
-    streams = manipulate.duplicated octaves str;
-
-    // forward transform uses the conjugate-transposed kernel, inverse
-    // uses the original
-    kernel = cm.transposed (cm.conjugateTransposed kdata.kernel);
-
-//    eprintln "have kernel";
-
-    fftFunc = fft.forward kdata.fftSize;
-
-    cqblocks =
-        map do octave:
-            frames = map ch.mixedDown //!!! mono for now
-               (framer.frames kdata.fftSize [ Hop kdata.fftHop, Padded false ]
-                   (resample.decimated (pow 2 octave) streams[octave]));
-            map do frame:
-                freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize));
-// eprintln "octave = \(octave), frame = \(vec.list frame)";
-// eprintln "octave = \(octave), freq = \(freq)";
-                cm.product kernel (cm.newComplexColumnVector freq);
-            done frames;
-        done [0..octaves-1];
-
-    // The cqblocks list is a list<list<matrix>>. Each top-level list
-    // corresponds to an octave, from highest to lowest, each having
-    // twice as many elements in its list as the next octave. The
-    // sub-lists are sampled in time with an effective spacing of
-    // fftSize * 2^(octave-1) audio frames, and the matrices are row
-    // vectors with atomsPerFrame * binsPerOctave complex elements.
-    //
-    // ***
-    // 
-    // In a typical constant-Q structure, each (2^(octaves-1) *
-    // fftHop) input frames gives us an output structure conceptually
-    // like this:
-    //
-    // [][][][][][][][]   <- fftHop frames per highest-octave output value
-    // [][][][][][][][]      layered as many times as binsPerOctave (here 2)
-    // [--][--][--][--]   <- fftHop*2 frames for the next lower octave
-    // [--][--][--][--]      etc
-    // [------][------]
-    // [------][------]
-    // [--------------]
-    // [--------------]
-    //
-    // ***
-    //
-    // But the kernel we're using here has more than one temporally
-    // spaced atom; each individual cell is a row vector with
-    // atomsPerFrame * binsPerOctave elements, but that actually
-    // represents a rectangular matrix of result cells with width
-    // atomsPerFrame and height binsPerOctave. The columns of this
-    // matrix (the atoms) then need to be spaced by 2^(octave-1)
-    // relative to those from the highest octave.
-
-    // Reshape each row vector into the appropriate rectangular matrix
-    // and split into single-atom columns
-
-    emptyHops = kdata.firstCentre / kdata.atomSpacing; //!!! int? round?
-//    maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops;
-//    eprintln "maxDrop = \(maxDrop)";
-
-    cqblocks =
-        map do octlist:
-            concatMap do rv:
-                cm.asColumns
-                   (cm.generate do row col:
-                        cm.at rv ((row * kdata.atomsPerFrame) + col) 0
-                    done {
-                        rows = kdata.binsPerOctave,
-                        columns = kdata.atomsPerFrame
-                    })
-            done octlist
-        done cqblocks;
-
-    cqblocks = array (map2 do octlist octave:
-        d = emptyHops * (pow 2 (octaves-octave)) - emptyHops;
-//        eprintln "dropping \(d)";
-        drop d octlist;
-    done cqblocks [1..octaves]);
-
-    assembleBlock bits =
-       (//eprintln "assembleBlock: structure of bits is:";
-        //eprintln (map length bits);
-
-        rows = octaves * kdata.binsPerOctave;
-        columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
-
-        cm.generate do row col:
-
-            // bits structure: [1,2,4,8,...]
-
-            // each elt of bits is a list of the chunks that should
-            // make up this block in that octave (lowest octave first)
-
-            // each chunk has atomsPerFrame * binsPerOctave elts in it
-
-            // row is disposed with 0 at the top, highest octave (in
-            // both pitch and index into bits structure)
-
-            oct = int (row / binsPerOctave);
-            binNo = row % kdata.binsPerOctave;
-
-            chunks = pow 2 oct;
-            colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
-            atomNo = int (col / colsPerAtom);
-            atomOffset = col % colsPerAtom;
-
-            if atomOffset == 0 and atomNo < length bits[oct] then
-                bits[oct][atomNo][binNo];
-            else
-                cplx.zero
-            fi;
-
-        done { rows, columns };
-        );
-
-    assembleBlockSpectrogram bits =
-       (// As assembleBlock, but producing a dense magnitude
-        // spectrogram (rather than a complex output). (todo:
-        // interpolation, smoothing)
-
-        //eprintln "assembleBlockSpectrogram: structure of bits is:";
-        //eprintln (map length bits);
-
-        rows = octaves * kdata.binsPerOctave;
-        columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
-
-        mat.generate do row col:
-
-            oct = int (row / binsPerOctave);
-            binNo = row % kdata.binsPerOctave;
-
-            chunks = pow 2 oct;
-            colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
-            atomNo = int (col / colsPerAtom);
-
-            if atomNo < length bits[oct] then
-                cplx.magnitude bits[oct][atomNo][binNo];
-            else 
-                0
-            fi;
-
-        done { rows, columns };
-        );
-
-    processOctaveLists assembler octs =
-        case octs[0] of
-        block::rest:
-           (toAssemble = array 
-               (map do oct:
-                    n = kdata.atomsPerFrame * pow 2 oct;
-                    if not empty? octs[oct] then
-                        forBlock = array (take n octs[oct]);
-                        octs[oct] := drop n octs[oct];
-                        forBlock
-                    else
-                        array []
-                    fi
-                done (keys octs));
-            assembler toAssemble :. \(processOctaveLists assembler octs));
-         _: []
-        esac;
-
-//eprintln "cqblocks has \(length cqblocks) entries";
-
-    octaveLists = [:];
-
-    cqblocks = array cqblocks;
-    for [1..octaves] do oct:
-        octaveLists[octaves - oct] := cqblocks[oct-1];
-    done;
-/*
-    \() (map2 do octlist octave:
-println "oct \(octaves) - \(octave) = \(octaves - octave)";
-             octaveLists[octaves - octave] := octlist 
-         done cqblocks [1..octaves]);
-*/
-//eprintln "octaveLists keys are: \(keys octaveLists)";
-
-    {
-        kernel = kdata with {
-            binFrequencies = array
-               (concatMap do octave:
-                    map do freq:
-                        freq / (pow 2 octave);
-                    done (reverse (list kdata.binFrequencies))
-                done [0..octaves-1])
-        },
-        sampleRate,
-        octaves,
-        get cqComplex () = processOctaveLists assembleBlock octaveLists,
-        get cqSpectrogram () = processOctaveLists assembleBlockSpectrogram octaveLists,
-    }
-    );
-
-{ cqt }
-

--- a/yeti/cqtkernel.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,158 +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.
-*/
-
-module cqtkernel;
-
-vec = load may.vector;
-complex = load may.complex;
-window = load may.signal.window;
-fft = load may.transform.fft;
-cm = load may.matrix.complex;
-
-{ pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc;
-
-makeKernel { sampleRate, maxFreq, binsPerOctave } =
-   (q = 1;
-    atomHopFactor = 0.25;
-    thresh = 0.0005;
-    minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave));
-    bigQ = q / ((pow 2 (1/binsPerOctave)) - 1);
-
-    maxNK = round(bigQ * sampleRate / minFreq);
-    minNK = round(bigQ * sampleRate /
-                  (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave))));
-
-    atomHop = round(minNK * atomHopFactor);
-    
-    firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop));
-    
-    fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2));
-    
-//    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";
-//    eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)";
-    
-    winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1;
-    
-    lastCentre = firstCentre + (winNr - 1) * atomHop;
-    
-    fftHop = (lastCentre + atomHop) - firstCentre;
-    
-//    eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)";
-    
-    fftFunc = fft.forward fftSize;
-    
-    // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
-    // function. We can't do that here; we need to generate real and imag
-    // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
-    
-    binFrequencies = array [];
-
-    kernels = map do k:
-    
-        nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave))));
-
-        // the cq MATLAB toolbox uses a symmetric window for
-        // blackmanharris -- which is odd because it uses a periodic one
-        // for other types. Oh well
-        win = vec.divideBy nk
-           (vec.sqrt
-               (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk));
-
-        fk = minFreq * (pow 2 ((k-1)/binsPerOctave));
-
-        push binFrequencies fk;
-
-        genKernel f = vec.multiply
-           [win,
-            vec.fromList
-               (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])];
-    
-        reals = genKernel cos;
-        imags = genKernel sin;
-    
-        atomOffset = firstCentre - ceil(nk/2);
-    
-        map do i:
-    
-            shift = vec.zeros (atomOffset + ((i-1) * atomHop));
-    
-            specKernel = fftFunc
-               (complex.complexArray
-                   (vec.concat [shift, reals])
-                   (vec.concat [shift, imags]));
-
-            map do c:
-                if complex.magnitude c <= thresh then complex.zero else c fi
-                done specKernel;
-
-        done [1..winNr];
-    
-    done [1..binsPerOctave];
-    
-    kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels)));
-    
-//    eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))";
-    
-    // Normalisation
-    
-    wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat));
-    wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat));
-
-    subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1));
-    square = cm.product (cm.conjugateTransposed subset) subset;
-
-    diag = complex.magnitudes (cm.getDiagonal 0 square);
-    wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2);
-
-    weight = (fftHop / fftSize) / (vec.mean (vec.abs wK));
-    weight = sqrt(weight);
-
-    kernel = cm.scaled weight kmat;
-
-//    eprintln "weight = \(weight)";
-    
-    {
-        kernel,
-        fftSize,
-        fftHop,
-        binsPerOctave,
-        atomsPerFrame = winNr,
-        atomSpacing = atomHop,
-        firstCentre,
-        maxFrequency = maxFreq,
-        minFrequency = minFreq,
-        binFrequencies,
-        bigQ
-    });
-
-{
-    makeKernel
-}
-

--- a/yeti/experiment.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,50 +0,0 @@
-
-program experiment;
-
-af = load may.stream.audiofile;
-plot = load may.plot;
-cm = load may.matrix.complex;
-mat = load may.matrix;
-vec = load may.vector;
-manipulate = load may.stream.manipulate;
-syn = load may.stream.syntheticstream;
-
-{ cqt } = load cqt;
-{ icqt } = load icqt;
-
-//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
-//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
-testStream = af.open "sweep-48000.wav";
-//testStream = af.open "sweep.wav";
-
-eprintln "have test stream";
-
-cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream;
-
-eprintln "bin frequencies: \(cq.kernel.binFrequencies)";
-
-icq = icqt cq;
-
-eprintln "calculated icq...";
-
-outfile = "test-output.wav";
-written = af.write icq outfile;
-
-eprintln "wrote \(written) to \(outfile)";
-
-//for cq.cqSpectrogram do mm:
-//    for (mat.asColumns mm) (println . strJoin "," . vec.list);
-//done;
-
-/*
-bigM = mat.concatHorizontal (map cm.magnitudes cq.output);
-
-eprintln "overall output size = \(mat.size bigM)";
-
-mat.print bigM;
-*/
-//\() (plot.plot [Contour bigM]);
-//\() (plot.plot [Grid bigM]);
-
-()
-

--- a/yeti/icqt.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,145 +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.
-*/
-
-module icqt;
-
-cqt = load cqt;
-cm = load may.matrix.complex;
-framer = load may.stream.framer;
-win = load may.signal.window;
-mm = load may.mathmisc;
-vec = load may.vector;
-resamp = load may.stream.resample;
-manip = load may.stream.manipulate;
-mat = load may.matrix;
-
-icqt cq =
-   (kdata = cq.kernel;
-
-    // kdata.kernel is the kernel matrix for a single octave. It has
-    // width kdata.fftSize and height kdata.binsPerOctave *
-    // kdata.atomsPerFrame.
-    //
-    // kdata.fftHop is the overlap between kernel matrices in a single
-    // octave.
-    // 
-    // cq.sampleRate is the output stream sample rate and cq.octaves
-    // the number of octaves.
-    //
-    // cq.cqComplex is the list of complex matrices containing the CQ
-    // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1)
-    // and height kdata.binsPerOctave * cq.octaves.
-
-    bpo = kdata.binsPerOctave;
-    atomsPerFrame = kdata.atomsPerFrame;
-    octaves = cq.octaves;
-
-    // transform a single block, all octaves tall, into an array
-    // (indexed by octave) of lists of individual columns (valid
-    // values for that octave only)
-    decomposeOctaves mat = array
-       (map do oct:
-            inverted = cm.fromRows (reverse (cm.asRows mat));
-            octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1));
-            gap = (mm.pow 2 oct) - 1;
-            pickFrom cols =
-                case cols of
-                c::cs: c::(pickFrom (drop gap cs));
-                _: [];
-                esac;
-            pickFrom (cm.asColumns octMat);
-        done [0..octaves-1]);
-
-    // transform a list of the arrays produced by decomposeOctaves
-    // into a single array (indexed by octave) of lists of the
-    // individual columns
-    flattenOctaves decomposed =
-       (flattenAux acc decomposed = 
-            case decomposed of
-            chunk::rest:
-                flattenAux 
-                   (array
-                       (map do oct:
-                            acc[oct] ++ chunk[oct]
-                            done [0..octaves-1]))
-                    rest;
-            _: acc;
-            esac;
-        flattenAux (array (map \[] [0..octaves-1])) decomposed);
-
-    // given a list of columns, deinterleave and pile up each sequence
-    // of atomsPerFrame columns into a single tall column and return
-    // the resulting list of tall columns
-    pile cols =
-       (pileAux acc cols =
-            if (length cols) >= atomsPerFrame then
-                atoms = take atomsPerFrame cols;
-                juggled = concat (reverse (cm.asRows (cm.fromColumns atoms)));
-                pileAux (juggled :: acc) (drop atomsPerFrame cols);
-            else
-                reverse acc
-            fi;
-        pileAux [] cols);
-
-    octaveColumnLists = 
-        map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex)));
-
-    for octaveColumnLists do l: println "octave column list length: \(length l)" done;
-
-    kernel = cm.transposed kdata.kernel; // right way around for the multiply
-    
-    spectra = 
-        map do l:
-            map do col:
-                res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col));
-                cm.columnSlice res 0 (kdata.fftSize / 2 + 1)
-            done l;
-        done octaveColumnLists;
-
-    eprintln "calculated spectra, now to ifft, overlap-add...";
-
-    rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1];
-    
-    resynthesised =
-        map2 do frames rate:
-            framer.complexStreamed rate kdata.fftSize
-                [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ]
-                frames;
-        done spectra rates;
-
-   eprintln "have streams:";
-   for resynthesised eprintln;
-
-   resampled = map (resamp.resampledTo cq.sampleRate) resynthesised;
-
-   manip.sum resampled;
-);
-
-{icqt}

--- a/yeti/nbproject/ide-file-targets.xml	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project basedir=".." name="cqt-IDE">
-    <!-- TODO: edit the following target according to your needs -->
-    <!-- (more info: http://www.netbeans.org/kb/articles/freeform-config.html#debugj2se) -->
-    <property name="maydir" value="${basedir}/../../may"/>
-    <property name="yetidir" value="${basedir}/../../yeti"/>
-    <property name="jardir" value="${maydir}/ext/jar"/>
-        <path id="cp">
-            <fileset dir="${basedir}">
-                <include name="cqt.jar"/>
-            </fileset>
-            <fileset dir="${jardir}">
-                <include name="*.jar"/>
-            </fileset>
-            <fileset dir="${yetidir}">
-                <include name="yeti.jar"/>
-            </fileset>
-            <fileset dir="${maydir}">
-                <include name="may.jar"/>
-            </fileset>
-        </path>
-    <target name="debug-nb">
-        <nbjpdastart addressproperty="jpda.address" name="cqt" transport="dt_socket">
-            <classpath refid="cp"/>
-        </nbjpdastart>
-        <!-- TODO configure the main class for your project here: -->
-        <java classname="plotfile" fork="true">
-            <classpath refid="cp"/>
-            <jvmarg value="-Xdebug"/>
-            <jvmarg value="-Xrunjdwp:transport=dt_socket,address=${jpda.address}"/>
-            <arg value="${basedir}/test.wav"/>
-        </java>
-    </target>
-    <target name="-profile-check">
-        <startprofiler freeform="true"/>
-    </target>
-    <!-- TODO: edit the following target according to your needs -->
-    <!-- (more info: http://netbeans.org/kb/articles/freeform-config.html#profilej2se) -->
-    <target depends="-profile-check" if="profiler.configured" name="profile-nb">
-        <!-- TODO configure the main class for your project here: -->
-        <java classname="plotfile" fork="true">
-            <classpath refid="cp"/>
-            <jvmarg line="${agent.jvmargs}"/>
-            <arg value="${basedir}/test.wav"/>
-        </java>
-    </target>
-</project>

--- a/yeti/nbproject/project.xml	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,53 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project xmlns="http://www.netbeans.org/ns/project/1">
-    <type>org.netbeans.modules.ant.freeform</type>
-    <configuration>
-        <general-data xmlns="http://www.netbeans.org/ns/freeform-project/1">
-            <!-- Do not use Project Properties customizer when editing this file manually. 
- To prevent the customizer from showing, create nbproject/project.properties file and enter 
-auxiliary.show.customizer=false 
-property there. Adding 
-auxiliary.show.customizer.message=<message>
- will show your customized message when someone attempts to open the customizer.  -->
-            <name>cqt</name>
-            <properties/>
-            <folders/>
-            <ide-actions>
-                <action name="build">
-                    <target>jar</target>
-                </action>
-                <action name="clean">
-                    <target>clean</target>
-                </action>
-                <action name="rebuild">
-                    <target>clean</target>
-                    <target>jar</target>
-                </action>
-                <action name="debug">
-                    <script>nbproject/ide-file-targets.xml</script>
-                    <target>debug-nb</target>
-                </action>
-                <action name="profile">
-                    <script>nbproject/ide-file-targets.xml</script>
-                    <target>profile-nb</target>
-                </action>
-            </ide-actions>
-            <view>
-                <items>
-                    <source-file>
-                        <location>build.xml</location>
-                    </source-file>
-                </items>
-                <context-menu>
-                    <ide-action name="build"/>
-                    <ide-action name="rebuild"/>
-                    <ide-action name="clean"/>
-                    <ide-action name="debug"/>
-                    <ide-action name="profile"/>
-                </context-menu>
-            </view>
-            <subprojects/>
-        </general-data>
-        <java-data xmlns="http://www.netbeans.org/ns/freeform-project-java/1"/>
-    </configuration>
-</project>

--- a/yeti/plotfile.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-
-program plotfile;
-
-af = load may.stream.audiofile;
-plot = load may.plot;
-cm = load may.matrix.complex;
-mat = load may.matrix;
-
-{ cqt } = load cqt;
-
-minFreq = 27.5;
-maxRateDivisor = 3;
-binsPerOctave = 60;
-
-plotfile f =
-   (testStream = af.open f;
-    eprintln "Opened file stream...";
-    start = System#currentTimeMillis();
-    cq = cqt {
-        maxFreq = testStream.sampleRate/maxRateDivisor, 
-        minFreq, 
-        binsPerOctave
-        } testStream;
-    middle = System#currentTimeMillis();
-    eprintln "Generated kernel \(cm.size cq.kernel.kernel) and primed transform (\(cq.octaves) octaves), took \(middle-start)ms, now calculating...";
-    bigM = mat.concatHorizontal cq.cqSpectrogram;
-    finish = System#currentTimeMillis();
-    eprintln "Done, that part took \(finish-middle)ms, all CQ stuff took \(finish-start)ms";
-    eprintln "Plotting...";
-    \() (plot.plot [Contour bigM]));
-
-usage () =
-   (eprintln "\nUsage: plotfile file.wav";
-    eprintln "\n   Loads audio from file.wav and plots a \(binsPerOctave)bpo Constant-Q spectrogram";
-    eprintln "   from \(minFreq)Hz up to 1/\(maxRateDivisor) of the file samplerate");
-
-case (list _argv) of
-file::[]: plotfile file;
-_: usage ();
-esac;
-

--- a/yeti/test.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,96 +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.
-*/
-
-program test;
-
-{ runTests } = load may.test;
-
-// We want to test:
-// 
-// Kernel design -- check size (number of bins, number of atoms);
-// check an example kernel against known data
-//
-// Time alignment -- feed a dirac train, check that peaks in all bins
-// align
-//
-// Frequency discrimination -- feed a sinusoid, check peaks
-//
-// Latency compensation -- for dirac at 0, check peak can be found at
-// 0 plus the declared latency
-//
-// Signal-noise ratio
-//
-// Specimen output for simple test case
-
-tests = [
-"kernel"    : load test_cqtkernel,
-"frequency" : load test_frequency,
-];
-
-bad = sum (mapHash do name testHash: runTests name testHash done tests);
-
-if (bad > 0) then
-    println "\n** \(bad) test(s) failed!";
-    threadExit 1
-else
-    ()
-fi
-
-
-
-
-/*
-//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
-//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
-testStream = af.open "sweep-48000.wav";
-//testStream = af.open "sweep.wav";
-
-// So the stream is [ 0, 1, 0, -1, 0, 1, 0, -1, ... ] :
-//testStream = manipulate.withDuration 64 (syn.sinusoid 8 2);
-
-testStream = manipulate.withDuration 32 (syn.pulseTrain 8 0.001);
-
-eprintln "have test stream";
-
-cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream;
-
-eprintln "bin frequencies: \(cq.kernel.binFrequencies)";
-
-bigM = mat.concatHorizontal (map cm.magnitudes cq.output);
-
-eprintln "overall output size = \(mat.size bigM)";
-
-mat.print bigM;
-
-//\() (plot.plot [Contour bigM]);
-\() (plot.plot [Grid bigM]);
-*/
-//()
-

--- a/yeti/test_cqtkernel.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,63 +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.
-*/
-
-module test_cqtkernel;
-
-cm = load may.matrix.complex;
-mm = load may.mathmisc;
-
-{ compare, compareUsing } = load may.test;
-
-{ makeKernel } = load cqtkernel;
-
-eps = 1e-7;
-
-compareClose = compareUsing do a b: abs (a - b) < eps done;
-
-[
-
-"minimal": \(
-    k = makeKernel { sampleRate = 16, maxFreq = 8, binsPerOctave = 4 };
-    compare k.binsPerOctave 4 and
-        compare (cm.size k.kernel) {
-            rows = k.binsPerOctave * k.atomsPerFrame,
-            columns = k.fftSize
-        } and
-        compareClose k.maxFrequency 8 and
-        compareClose k.minFrequency (4 * (mm.pow 2 (1/4))) and
-        compare k.atomsPerFrame 5 and
-        compare k.fftSize 32 and
-        compare (length k.binFrequencies) k.binsPerOctave and
-        compareClose (head k.binFrequencies) k.minFrequency and
-        compareClose (head (reverse k.binFrequencies)) k.maxFrequency
-),
-
-] is hash<string, () -> boolean>
-

--- a/yeti/test_frequency.yeti	Thu May 15 11:59:11 2014 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,132 +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.
-*/
-
-module test_frequency;
-
-mat = load may.matrix;
-vec = load may.vector;
-win = load may.signal.window;
-mm = load may.mathmisc;
-cm = load may.matrix.complex;
-syn = load may.stream.syntheticstream;
-plot = load may.plot;
-
-{ compare } = load may.test;
-
-{ cqt } = load cqt;
-
-// Test with a single windowed sinusoid, repeating at various frequencies
-
-sinTestStream sampleRate duration signalFreq = // duration is in samples
-   (sin = syn.sinusoid sampleRate signalFreq;
-    chunk = mat.getRow 0 (sin.read duration);
-    syn.precalculatedMono sampleRate (win.windowed win.hann chunk));
-
-// We want to make a CQ transform spanning more than one octave, but
-// not going all the way to fs/2 so we can test it also with
-// frequencies above and below its extents
-
-sampleRate = 100;
-
-// fs/2 = 50 so 10->40 gives us 2 octaves
-cqmin = 10;
-cqmax = 40;
-bpo = 4; // fairly arbitrary
-
-testFreqs = map (* 5) [ 0..10 ];
-duration = sampleRate * 2;
-
-threshold = 0.08;
-
-streamBuilder = sinTestStream sampleRate duration;
-
-binForFreq f =
-    mm.round (bpo * mm.log2 (f / cqmin)) - 1;
-
-report message matrix =
-   (eprintln message;
-    eprintln "matrix is:";
-    mat.eprint matrix);
-//    chart = plot.plot [Grid matrix];
-//    sleep 100;
-//    chart#dispose());
-
-tests = mapIntoHash
-    do f: "freq_\(f)" done
-    do f: \(
-        str = streamBuilder f;
-        cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str;
-        spec = cq.cqSpectrogram;
-        rightSize = all id
-           (map do s:
-            compare (mat.size s) { 
-                rows = cq.kernel.binsPerOctave * cq.octaves,
-                columns = cq.kernel.atomsPerFrame * mm.pow 2 (cq.octaves - 1)
-            }
-            done spec);
-        m = mat.concatHorizontal spec;
-//    println "binFrequencies = \(cq.kernel.binFrequencies)";
-//    println "binForFreq \(f) = \(binForFreq f)";
-        var colno = 0;
-        success = all id
-           (rightSize :: map do c:
-                // The test passes for this column if:
-                //
-                //  * the max bin is the expected one, or
-                //
-                //  * the expected max is out of range entirely (but
-                //  we need to test _something_ in this case --
-                //  what?), or
-                //
-                //  * all bins are below a threshold, or
-                //
-                //  * this is an odd column and the expected max is in
-                //  the lower octave
-                //
-                // We should also check that all values in the lower
-                // octave are zero for odd columns.
-                //
-                expected = binForFreq f;
-                good =
-                   (expected < 0 or expected >= vec.length c) or
-                   ((colno % 2 == 1) and expected < (vec.length c / 2)) or
-                   (vec.max c < threshold) or
-                   (vec.maxindex c == binForFreq f);
-                if (not good) then
-                    report " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)" m;
-                fi;
-                colno := colno + 1;
-                good;
-            done (mat.asColumns m));
-        success;
-    ) done
-    testFreqs;
-
-tests is hash<string, () -> boolean>