Merge branch 'master' into migrate-to-QJSEngine

48 files changed, 1898 insertions, 626 deletions

diff --git a/lib/qm-dsp/CONTRIBUTING.md b/lib/qm-dsp/CONTRIBUTING.md
new file mode 100644
index 0000000000..6f6229e610
--- /dev/null
+++ b/lib/qm-dsp/CONTRIBUTING.md
@@ -0,0 +1,79 @@
+
+Contributing
+============
+
+The qm-dsp library is maintained in a Github repository at
+https://github.com/c4dm/qm-dsp.
+
+
+Reporting bugs
+--------------
+
+Please use Github issues for bug reports. Try to make them as specific
+as possible. For example, describe an input that triggers some
+particular behaviour, and tell us how that behaviour differs from what
+you expected.
+
+If your bug can be reproduced by processing an audio file using one of
+the QM Vamp Plugins (https://github.com/c4dm/qm-vamp-plugins), which
+are built using this library, that might be a good way to illustrate
+the problem.
+
+
+Pull requests
+-------------
+
+We're happy to see pull requests, and can pull them directly in some
+circumstances.
+
+ * Please make sure your change compiles without warnings and passes
+   the existing tests.
+
+ * Please follow the code style guidelines (see below).
+
+ * Please make it as easy as possible to verify the behaviour of the
+   pull request, for example by adding a new test in the `tests`
+   directory. This library has only limited test coverage, but we
+   would like to expand it, and prefer not to make changes unless they
+   are backed by tests.
+
+ * Please provide your changes under terms which permit Queen Mary
+   University of London to relicense the code for commercial
+   purposes. The qm-dsp library as a whole is provided under the GPL,
+   but QM also make commercial licences available separately, and
+   cannot accept any pull request whose copyright status would prevent
+   that. In practice, this means any non-trivial change not
+   originating from QM must be explicitly licensed using a BSD-like
+   licence text, either in the source file itself or in an
+   accompanying file. See `thread/BlockAllocator.h` for an example of
+   typical language.
+
+Please note also that fixes which change the behaviour of the existing
+QM Vamp Plugins will need particularly close scrutiny - these are
+reasonably widely used and, even where they have defects, changes may
+cause problems for users and will at least need to be documented with
+the plugins. For this reason it may take some time for such changes to
+be reviewed or integrated.
+
+
+Code style
+----------
+
+ * C++ code must compile with the C++98 standard, except for the unit
+   tests which are C++14
+
+ * Classes are named `LikeThis` - functions, methods, and local
+   variables `likeThis` - class member variables `m_likeThis`
+
+ * Indentation is four spaces at a time (no tabs)
+
+ * The opening brace for a block goes at the end of the line, except
+   at the start of a function or class definition where it gets a line
+   of its own
+
+ * Please use braces around any conditional or loop block that
+   occupies its own line
+
+Some of the older code in this library does not follow these
+guidelines - usually this means the code needs to be updated.
+
diff --git a/lib/qm-dsp/README.md b/lib/qm-dsp/README.md
new file mode 100644
index 0000000000..20bcb60f04
--- /dev/null
+++ b/lib/qm-dsp/README.md
@@ -0,0 +1,59 @@
+
+QM-DSP library
+==============
+
+This is a C++ library of functions for Digital Signal Processing and
+Music Informatics purposes developed in the [Centre for Digital
+Music](http://c4dm.eecs.qmul.ac.uk) at Queen Mary, University of
+London.
+
+It is used by [QM Vamp Plugins](http://isophonics.net/QMVampPlugins)
+amongst other things.
+
+Despite the assertive name "qm-dsp", it is not "the official QM DSP
+library", just one library for DSP that happens to have been written
+at QM. It got this name because nothing else was using it at the time.
+
+
+Compiling the library
+---------------------
+
+ - Linux: `make -f build/linux/Makefile.linux64`
+
+ - Mac: `make -f build/osx/Makefile.osx`
+
+ - Windows (MSVC): Use the project file `build/msvc/QMDSP.vcxproj`
+
+To build and run unit tests as well, add the `test` target to your
+Make invocation, e.g. `make -f build/linux/Makefile.linux64
+test`. Tests require the Boost library.
+
+
+Licence
+-------
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at
+your option) any later version.  See the file COPYING included with
+this distribution for more information.
+
+This code is Copyright (c) 2006-2019 Queen Mary, University of London,
+with the following exceptions:
+
+ - `ext/kissfft` - Copyright (c) 2003-2010 Mark Borgerding
+
+ - `maths/pca/pca.c` - Fionn Murtagh, from StatLib, used with permission
+
+ - `maths/Polyfit.h` - by Allen Miller, David J Taylor and others;
+also for Delphi in the the JEDI Math Library, under the Mozilla Public
+License
+
+ - `thread/BlockAllocator.h` - derived from FSB Allocator by Juha
+Nieminen, under a BSD-style license
+
+See individual files for further authorship details.
+
+If you wish to use this code in a proprietary application or product
+for which the terms of the GPL are not appropriate, please contact QM
+Innovation https://www.qminnovation.co.uk/ for licensing terms.
diff --git a/lib/qm-dsp/README.txt b/lib/qm-dsp/README.txt
deleted file mode 100644
index 6927c7a1b8..0000000000
--- a/lib/qm-dsp/README.txt
+++ /dev/null
@@ -1,35 +0,0 @@
-
-
-QM-DSP library
-==============
-
-This is a C++ library of functions for DSP and Music Informatics
-purposes developed at Queen Mary, University of London.
-It is used by the QM Vamp Plugins (q.v.) amongst other things.
-
-This code is Copyright (c) 2006-2015 Queen Mary, University of London,
-with the following exceptions:
-
-ext/kissfft -- Copyright (c) 2003-2010 Mark Borgerding
-
-maths/pca.c -- Fionn Murtagh, from StatLib; with permission
-
-maths/Polyfit.h -- Allen Miller, David J Taylor and others; also for
-Delphi in the the JEDI Math Library, under the Mozilla Public License
-
-thread/BlockAllocator.h -- derived from FSB Allocator by Juha Nieminen,
-under BSD-style license
-
-See individual files for further authorship details.
-
-
-License
-=======
-
-This program is free software; you can redistribute it and/or
-modify it under the terms of the GNU General Public License as
-published by the Free Software Foundation; either version 2 of the
-License, or (at your option) any later version.  See the file
-COPYING included with this distribution for more information.
-
-
diff --git a/lib/qm-dsp/base/Restrict.h b/lib/qm-dsp/base/Restrict.h
new file mode 100644
index 0000000000..723a8b3b67
--- /dev/null
+++ b/lib/qm-dsp/base/Restrict.h
@@ -0,0 +1,17 @@
+

+#ifndef QM_DSP_RESTRICT_H

+#define QM_DSP_RESTRICT_H

+

+#ifdef _MSC_VER

+#define QM_R__ __restrict

+#endif

+

+#ifdef __GNUC__

+#define QM_R__ __restrict__

+#endif

+

+#ifndef QM_R__

+#define QM_R__

+#endif

+

+#endif

diff --git a/lib/qm-dsp/dsp/chromagram/Chromagram.cpp b/lib/qm-dsp/dsp/chromagram/Chromagram.cpp
index a8597a5ddd..9f40b1d2f1 100644
--- a/lib/qm-dsp/dsp/chromagram/Chromagram.cpp
+++ b/lib/qm-dsp/dsp/chromagram/Chromagram.cpp
@@ -33,8 +33,9 @@ int Chromagram::initialise( ChromaConfig Config )
     m_BPO  = Config.BPO;		// bins per octave
     m_normalise = Config.normalise;     // if frame normalisation is required
 
-    // No. of constant Q bins
-    m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0));	
+    // Extend range to a full octave
+    double octaves = log(m_FMax / m_FMin) / log(2.0);
+    m_FMax = m_FMin * pow(2.0, ceil(octaves));
 
     // Create array for chroma result
     m_chromadata = new double[ m_BPO ];
@@ -44,7 +45,7 @@ int Chromagram::initialise( ChromaConfig Config )
 
     // Populate CQ config structure with parameters
     // inherited from the Chroma config
-    ConstantQConfig.FS	 = Config.FS;
+    ConstantQConfig.FS = Config.FS;
     ConstantQConfig.min = m_FMin;
     ConstantQConfig.max = m_FMax;
     ConstantQConfig.BPO = m_BPO;
@@ -53,6 +54,9 @@ int Chromagram::initialise( ChromaConfig Config )
     // Initialise ConstantQ operator
     m_ConstantQ = new ConstantQ( ConstantQConfig );
 
+    // No. of constant Q bins
+    m_uK = m_ConstantQ->getK();
+
     // Initialise working arrays
     m_frameSize = m_ConstantQ->getfftlength();
     m_hopSize = m_ConstantQ->gethop();
@@ -112,7 +116,7 @@ void Chromagram::unityNormalise(double *src)
 
     MathUtilities::getFrameMinMax( src, m_BPO, & min, &max );
 
-    for( unsigned int i = 0; i < m_BPO; i++ )
+    for (int i = 0; i < m_BPO; i++)
     {
 	val = src[ i ] / max;
 
@@ -121,7 +125,7 @@ void Chromagram::unityNormalise(double *src)
 }
 
 
-double* Chromagram::process( const double *data )
+double *Chromagram::process(const double *data)
 {
     if (!m_skGenerated) {
         // Generate CQ Kernel 
@@ -139,12 +143,20 @@ double* Chromagram::process( const double *data )
     }
     m_window->cut(m_windowbuf);
 
+    // The frequency-domain version expects pre-fftshifted input - so
+    // we must do the same here
+    for (int i = 0; i < m_frameSize/2; ++i) {
+        double tmp = m_windowbuf[i];
+        m_windowbuf[i] = m_windowbuf[i + m_frameSize/2];
+        m_windowbuf[i + m_frameSize/2] = tmp;
+    }
+
     m_FFT->forward(m_windowbuf, m_FFTRe, m_FFTIm);
 
     return process(m_FFTRe, m_FFTIm);
 }
 
-double* Chromagram::process( const double *real, const double *imag )
+double *Chromagram::process(const double *real, const double *imag)
 {
     if (!m_skGenerated) {
         // Generate CQ Kernel 
@@ -153,19 +165,17 @@ double* Chromagram::process( const double *real, const double *imag )
     }
 
     // initialise chromadata to 0
-    for (unsigned i = 0; i < m_BPO; i++) m_chromadata[i] = 0;
+    for (int i = 0; i < m_BPO; i++) m_chromadata[i] = 0;
 
-    double cmax = 0.0;
-    double cval = 0;
     // Calculate ConstantQ frame
     m_ConstantQ->process( real, imag, m_CQRe, m_CQIm );
 	
     // add each octave of cq data into Chromagram
-    const unsigned octaves = (int)floor(double( m_uK/m_BPO))-1;
-    for (unsigned octave = 0; octave <= octaves; octave++) 
+    const int octaves = m_uK / m_BPO;
+    for (int octave = 0; octave < octaves; octave++) 
     {
-	unsigned firstBin = octave*m_BPO;
-	for (unsigned i = 0; i < m_BPO; i++) 
+	int firstBin = octave*m_BPO;
+	for (int i = 0; i < m_BPO; i++) 
 	{
 	    m_chromadata[i] += kabs( m_CQRe[ firstBin + i ], m_CQIm[ firstBin + i ]);
 	}
diff --git a/lib/qm-dsp/dsp/chromagram/Chromagram.h b/lib/qm-dsp/dsp/chromagram/Chromagram.h
index bd928f5d48..ca68ee9071 100644
--- a/lib/qm-dsp/dsp/chromagram/Chromagram.h
+++ b/lib/qm-dsp/dsp/chromagram/Chromagram.h
@@ -20,11 +20,11 @@
 #include "base/Window.h"
 #include "ConstantQ.h"
 
-struct ChromaConfig{
-    unsigned int FS;
+struct ChromaConfig {
+    double FS;
     double min;
     double max;
-    unsigned int BPO;
+    int BPO;
     double CQThresh;
     MathUtilities::NormaliseType normalise;
 };
@@ -35,19 +35,44 @@ class Chromagram
 public:	
     Chromagram( ChromaConfig Config );
     ~Chromagram();
-	
-    double* process( const double *data ); // time domain
-    double* process( const double *real, const double *imag ); // frequency domain
-    void unityNormalise( double* src );
+
+    /**
+     * Process a time-domain input signal of length equal to
+     * getFrameSize().
+     * 
+     * The returned buffer contains the chromagram values indexed by
+     * bin, with the number of values corresponding to the BPO field
+     * in the ChromaConfig supplied at construction. It is owned by
+     * the Chromagram object and is reused from one process call to
+     * the next.
+     */
+    double *process(const double *data);
+    
+    /**
+     * Process a frequency-domain input signal generated from a
+     * time-domain signal of length equal to getFrameSize() that has
+     * been windowed and "fftshifted" to place the zero index in the
+     * centre of the frame. The real and imag buffers must each
+     * contain the full getFrameSize() frequency bins.
+     * 
+     * The returned buffer contains the chromagram values indexed by
+     * bin, with the number of values corresponding to the BPO field
+     * in the ChromaConfig supplied at construction. It is owned by
+     * the Chromagram object and is reused from one process call to
+     * the next.
+     */
+    double *process(const double *real, const double *imag);
+    
+    void unityNormalise(double* src);
 
     // Complex arithmetic
     double kabs( double real, double imag );
 	
     // Results
-    unsigned int getK() { return m_uK;}
-    unsigned int getFrameSize() { return m_frameSize; }
-    unsigned int getHopSize()   { return m_hopSize; }
-
+    int getK() { return m_uK;}
+    int getFrameSize() { return m_frameSize; }
+    int getHopSize()   { return m_hopSize; }
+    
 private:
     int initialise( ChromaConfig Config );
     int deInitialise();
@@ -58,13 +83,13 @@ private:
     double* m_chromadata;
     double m_FMin;
     double m_FMax;
-    unsigned int m_BPO;
-    unsigned int m_uK;
+    int m_BPO;
+    int m_uK;
 
     MathUtilities::NormaliseType m_normalise;
 
-    unsigned int m_frameSize;
-    unsigned int m_hopSize;
+    int m_frameSize;
+    int m_hopSize;
 
     FFTReal* m_FFT;
     ConstantQ* m_ConstantQ;
diff --git a/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp b/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp
index b764235fcd..a5b1c12105 100644
--- a/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp
+++ b/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp
@@ -117,28 +117,29 @@ void ConstantQ::sparsekernel()
 
     FFT m_FFT(m_FFTLength);
 	
-    for (unsigned k = m_uK; k--; ) 
-    {
-        for (unsigned u=0; u < m_FFTLength; u++) 
-        {
+    for (unsigned k = m_uK; k--;) {
+        for (unsigned u=0; u < m_FFTLength; u++) {
             hammingWindowRe[u] = 0;
             hammingWindowIm[u] = 0;
         }
         
-	// Computing a hamming window
-	const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
+        const double samplesPerCycle =
+                m_FS / (m_FMin * pow(2, (double)k / (double)m_BPO));
+
+        // Computing a hamming window
+        const unsigned hammingLength = (int) ceil(
+                m_dQ * samplesPerCycle);
 
         unsigned origin = m_FFTLength/2 - hammingLength/2;
 
-	for (unsigned i=0; i<hammingLength; i++) 
-	{
-	    const double angle = 2*PI*m_dQ*i/hammingLength;
-	    const double real = cos(angle);
-	    const double imag = sin(angle);
-	    const double absol = hamming(hammingLength, i)/hammingLength;
-	    hammingWindowRe[ origin + i ] = absol*real;
-	    hammingWindowIm[ origin + i ] = absol*imag;
-	}
+        for (unsigned i=0; i<hammingLength; i++) {
+            const double angle = 2*PI*i/samplesPerCycle;
+            const double real = cos(angle);
+            const double imag = sin(angle);
+            const double absol = hamming(hammingLength, i)/hammingLength;
+            hammingWindowRe[ origin + i ] = absol*real;
+            hammingWindowIm[ origin + i ] = absol*imag;
+        }
 
         for (unsigned i = 0; i < m_FFTLength/2; ++i) {
             double temp = hammingWindowRe[i];
@@ -149,24 +150,24 @@ void ConstantQ::sparsekernel()
             hammingWindowIm[i + m_FFTLength/2] = temp;
         }
     
-	//do fft of hammingWindow
-	m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
-
-		
-	for (unsigned j=0; j<( m_FFTLength ); j++) 
-	{
-	    // perform thresholding
-	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
-	    if (squaredBin <= squareThreshold) continue;
-		
-	    // Insert non-zero position indexes, doubled because they are floats
-	    sk->is.push_back(j);
-	    sk->js.push_back(k);
-
-	    // take conjugate, normalise and add to array sparkernel
-	    sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
-	    sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
-	}
+        //do fft of hammingWindow
+        m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
+
+
+        for (unsigned j=0; j<( m_FFTLength ); j++) {
+            // perform thresholding
+            const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
+            if (squaredBin <= squareThreshold) {
+                continue;
+            }
+            // Insert non-zero position indexes
+            sk->is.push_back(j);
+            sk->js.push_back(k);
+
+            // take conjugate, normalise and add to array sparkernel
+            sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
+            sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
+        }
 
     }
 
@@ -256,10 +257,9 @@ double* ConstantQ::process( const double* fftdata )
 
     SparseKernel *sk = m_sparseKernel;
 
-    for (unsigned row=0; row<2*m_uK; row++) 
-    {
-	m_CQdata[ row ] = 0;
-	m_CQdata[ row+1 ] = 0;
+    for (unsigned row=0; row<2*m_uK; row++) {
+        m_CQdata[ row ] = 0;
+        m_CQdata[ row+1 ] = 0;
     }
     const unsigned *fftbin = &(sk->is[0]);
     const unsigned *cqbin  = &(sk->js[0]);
@@ -267,17 +267,19 @@ double* ConstantQ::process( const double* fftdata )
     const double   *imag   = &(sk->imag[0]);
     const unsigned int sparseCells = sk->real.size();
 	
-    for (unsigned i = 0; i<sparseCells; i++)
-    {
-	const unsigned row = cqbin[i];
-	const unsigned col = fftbin[i];
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
-	const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
-	// add the multiplication
-	m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
-	m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
+    for (unsigned i = 0; i<sparseCells; i++) {
+        const unsigned row = cqbin[i];
+        const unsigned col = fftbin[i];
+        if (col == 0) {
+            continue;
+        }
+        const double & r1  = real[i];
+        const double & i1  = imag[i];
+        const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
+        const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
+        // add the multiplication
+        m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
+        m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
     }
 
     return m_CQdata;
@@ -300,7 +302,7 @@ void ConstantQ::initialise( CQConfig Config )
     // work out length of fft required for this constant Q Filter bank
     m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
 
-    m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32
+    m_hop = m_FFTLength/8;
 
 //    std::cerr << "ConstantQ::initialise: -> fft length = " << m_FFTLength << ", hop = " << m_hop << std::endl;
 
@@ -324,10 +326,9 @@ void ConstantQ::process(const double *FFTRe, const double* FFTIm,
 
     SparseKernel *sk = m_sparseKernel;
 
-    for (unsigned row=0; row<m_uK; row++) 
-    {
-	CQRe[ row ] = 0;
-	CQIm[ row ] = 0;
+    for (unsigned row=0; row<m_uK; row++) {
+        CQRe[ row ] = 0;
+        CQIm[ row ] = 0;
     }
 
     const unsigned *fftbin = &(sk->is[0]);
@@ -336,16 +337,18 @@ void ConstantQ::process(const double *FFTRe, const double* FFTIm,
     const double   *imag   = &(sk->imag[0]);
     const unsigned int sparseCells = sk->real.size();
 	
-    for (unsigned i = 0; i<sparseCells; i++)
-    {
-	const unsigned row = cqbin[i];
-	const unsigned col = fftbin[i];
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = FFTRe[ m_FFTLength - col - 1 ];
-	const double & i2  = FFTIm[ m_FFTLength - col - 1 ];
-	// add the multiplication
-	CQRe[ row ] += (r1*r2 - i1*i2);
-	CQIm[ row ] += (r1*i2 + i1*r2);
+    for (unsigned i = 0; i<sparseCells; i++) {
+        const unsigned row = cqbin[i];
+        const unsigned col = fftbin[i];
+        if (col == 0) {
+            continue;
+        }
+        const double & r1  = real[i];
+        const double & i1  = imag[i];
+        const double & r2  = FFTRe[ m_FFTLength - col ];
+        const double & i2  = FFTIm[ m_FFTLength - col ];
+        // add the multiplication
+        CQRe[ row ] += (r1*r2 - i1*i2);
+        CQIm[ row ] += (r1*i2 + i1*r2);
     }
 }
diff --git a/lib/qm-dsp/dsp/chromagram/ConstantQ.h b/lib/qm-dsp/dsp/chromagram/ConstantQ.h
index c06f60a9d0..7507a1f964 100644
--- a/lib/qm-dsp/dsp/chromagram/ConstantQ.h
+++ b/lib/qm-dsp/dsp/chromagram/ConstantQ.h
@@ -20,8 +20,8 @@
 #include "maths/MathAliases.h"
 #include "maths/MathUtilities.h"
 
-struct CQConfig{
-    unsigned int FS;   // samplerate
+struct CQConfig {
+    double FS;         // samplerate
     double min;        // minimum frequency
     double max;        // maximum frequency
     unsigned int BPO;  // bins per octave
@@ -58,7 +58,7 @@ private:
     void deInitialise();
 	
     double* m_CQdata;
-    unsigned int m_FS;
+