update qm-dsp library

7 files changed, 351 insertions, 41 deletions

diff --git a/libs/qm-dsp/base/KaiserWindow.cpp b/libs/qm-dsp/base/KaiserWindow.cpp
new file mode 100644
index 0000000000..4fe838e090
--- /dev/null
+++ b/libs/qm-dsp/base/KaiserWindow.cpp
@@ -0,0 +1,65 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+ 
+    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.
+*/
+
+#include "KaiserWindow.h"
+
+#include "maths/MathUtilities.h"
+
+KaiserWindow::Parameters
+KaiserWindow::parametersForTransitionWidth(double attenuation,
+					   double transition)
+{
+    Parameters p;
+    p.length = 1 + (attenuation > 21.0 ?
+		    ceil((attenuation - 7.95) / (2.285 * transition)) :
+		    ceil(5.79 / transition));
+    p.beta = (attenuation > 50.0 ? 
+	      0.1102 * (attenuation - 8.7) :
+	      attenuation > 21.0 ? 
+	      0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) :
+	      0);
+    return p;
+}
+
+static double besselTerm(double x, int i)
+{
+    if (i == 0) {
+	return 1;
+    } else {
+	double f = MathUtilities::factorial(i);
+	return pow(x/2, i*2) / (f*f);
+    }
+}
+
+static double bessel0(double x)
+{
+    double b = 0.0;
+    for (int i = 0; i < 20; ++i) {
+	b += besselTerm(x, i);
+    }
+    return b;
+}
+
+void
+KaiserWindow::init()
+{
+    double denominator = bessel0(m_beta);
+    bool even = (m_length % 2 == 0);
+    for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) {
+	double k = double(2*i) / double(m_length-1) - 1.0;
+	m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator);
+    }
+    for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) {
+        m_window.push_back(m_window[int(m_length/2) - i - 1]);
+    }
+}
diff --git a/libs/qm-dsp/base/KaiserWindow.h b/libs/qm-dsp/base/KaiserWindow.h
new file mode 100644
index 0000000000..f16a4b6c16
--- /dev/null
+++ b/libs/qm-dsp/base/KaiserWindow.h
@@ -0,0 +1,105 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+ 
+    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.
+*/
+
+#ifndef KAISER_WINDOW_H
+#define KAISER_WINDOW_H
+
+#include <vector>
+#include <cmath>
+
+/**
+ * Kaiser window: A windower whose bandwidth and sidelobe height
+ * (signal-noise ratio) can be specified. These parameters are traded
+ * off against the window length.
+ */
+class KaiserWindow
+{
+public:
+    struct Parameters {
+	int length;
+	double beta;
+    };
+
+    /**
+     * Construct a Kaiser windower with the given length and beta
+     * parameter.
+     */
+    KaiserWindow(Parameters p) : m_length(p.length), m_beta(p.beta) { init(); }
+
+    /**
+     * Construct a Kaiser windower with the given attenuation in dB
+     * and transition width in samples.
+     */
+    static KaiserWindow byTransitionWidth(double attenuation,
+					  double transition) {
+	return KaiserWindow
+	    (parametersForTransitionWidth(attenuation, transition));
+    }
+
+    /**
+     * Construct a Kaiser windower with the given attenuation in dB
+     * and transition bandwidth in Hz for the given samplerate.
+     */
+    static KaiserWindow byBandwidth(double attenuation,
+				    double bandwidth,
+				    double samplerate) {
+	return KaiserWindow
+	    (parametersForBandwidth(attenuation, bandwidth, samplerate));
+    }
+
+    /**
+     * Obtain the parameters necessary for a Kaiser window of the
+     * given attenuation in dB and transition width in samples.
+     */
+    static Parameters parametersForTransitionWidth(double attenuation,
+						   double transition);
+
+    /**
+     * Obtain the parameters necessary for a Kaiser window of the
+     * given attenuation in dB and transition bandwidth in Hz for the
+     * given samplerate.
+     */
+    static Parameters parametersForBandwidth(double attenuation,
+					     double bandwidth,
+					     double samplerate) {
+	return parametersForTransitionWidth
+	    (attenuation, (bandwidth * 2 * M_PI) / samplerate);
+    } 
+
+    int getLength() const {
+	return m_length;
+    }
+
+    const double *getWindow() const { 
+	return m_window.data();
+    }
+
+    void cut(double *src) const { 
+	cut(src, src); 
+    }
+
+    void cut(const double *src, double *dst) const {
+	for (int i = 0; i < m_length; ++i) {
+	    dst[i] = src[i] * m_window[i];
+	}
+    }
+
+private:
+    int m_length;
+    double m_beta;
+    std::vector<double> m_window;
+
+    void init();
+};
+
+#endif
diff --git a/libs/qm-dsp/base/Pitch.cpp b/libs/qm-dsp/base/Pitch.cpp
index e86d6792b7..9ba8b4b3f7 100644
--- a/libs/qm-dsp/base/Pitch.cpp
+++ b/libs/qm-dsp/base/Pitch.cpp
@@ -39,7 +39,7 @@ Pitch::getPitchForFrequency(float frequency,
 	midiPitch = midiPitch + 1;
 	centsOffset = -(100.0 - centsOffset);
     }
-
+    
     if (centsOffsetReturn) *centsOffsetReturn = centsOffset;
     return midiPitch;
 }
diff --git a/libs/qm-dsp/base/Pitch.h b/libs/qm-dsp/base/Pitch.h
index 7b1d6c96aa..5b55ecc4a3 100644
--- a/libs/qm-dsp/base/Pitch.h
+++ b/libs/qm-dsp/base/Pitch.h
@@ -15,6 +15,10 @@
 #ifndef _PITCH_H_
 #define _PITCH_H_
 
+/**
+ * Convert between musical pitch (i.e. MIDI pitch number) and
+ * fundamental frequency.
+ */
 class Pitch
 {
 public:
diff --git a/libs/qm-dsp/base/SincWindow.cpp b/libs/qm-dsp/base/SincWindow.cpp
new file mode 100644
index 0000000000..5d31a631e3
--- /dev/null
+++ b/libs/qm-dsp/base/SincWindow.cpp
@@ -0,0 +1,45 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+ 
+    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.
+*/
+
+#include "SincWindow.h"
+
+#include <cmath>
+
+void
+SincWindow::init()
+{
+    if (m_length < 1) {
+	return;
+    } else if (m_length < 2) {
+	m_window.push_back(1);
+	return;
+    } else {
+
+	int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2);
+	int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2);
+	double m = 2 * M_PI / m_p;
+
+	for (int i = 0; i < n0; ++i) {
+	    double x = ((m_length / 2) - i) * m;
+	    m_window.push_back(sin(x) / x);
+	}
+
+	m_window.push_back(1.0);
+
+	for (int i = 1; i < n1; ++i) {
+	    double x = i * m;
+	    m_window.push_back(sin(x) / x);
+	}
+    }
+}
+
diff --git a/libs/qm-dsp/base/SincWindow.h b/libs/qm-dsp/base/SincWindow.h
new file mode 100644
index 0000000000..bb35d90c20
--- /dev/null
+++ b/libs/qm-dsp/base/SincWindow.h
@@ -0,0 +1,61 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+ 
+    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.
+*/
+
+#ifndef SINC_WINDOW_H
+#define SINC_WINDOW_H
+
+#include <vector>
+
+/**
+ * A window containing values of the sinc function, i.e. sin(x)/x with
+ * sinc(0) == 1, with x == 0 at the centre.
+ */
+class SincWindow
+{
+public:
+    /**
+     * Construct a windower of the given length, containing the values
+     * of sinc(x) with x=0 in the middle, i.e. at sample (length-1)/2
+     * for odd or (length/2)+1 for even length, such that the distance
+     * from -pi to pi (the nearest zero crossings either side of the
+     * peak) is p samples.
+     */
+    SincWindow(int length, double p) : m_length(length), m_p(p) { init(); }
+
+    int getLength() const {
+	return m_length;
+    }
+
+    const double *getWindow() const { 
+	return m_window.data();
+    }
+
+    void cut(double *src) const { 
+	cut(src, src); 
+    }
+
+    void cut(const double *src, double *dst) const {
+	for (int i = 0; i < m_length; ++i) {
+	    dst[i] = src[i] * m_window[i];
+	}
+    }
+
+private:
+    int m_length;
+    double m_p;
+    std::vector<double> m_window;
+
+    void init();
+};
+
+#endif
diff --git a/libs/qm-dsp/base/Window.h b/libs/qm-dsp/base/Window.h
index 02ca426d66..93909115c0 100644
--- a/libs/qm-dsp/base/Window.h
+++ b/libs/qm-dsp/base/Window.h
@@ -18,6 +18,7 @@
 #include <cmath>
 #include <iostream>
 #include <map>
+#include <vector>
 
 enum WindowType {
     RectangularWindow,
@@ -25,18 +26,28 @@ enum WindowType {
     HammingWindow,
     HanningWindow,
     BlackmanWindow,
-    GaussianWindow,
-    ParzenWindow
+    BlackmanHarrisWindow,
+
+    FirstWindow = RectangularWindow,
+    LastWindow = BlackmanHarrisWindow
 };
 
+/**
+ * Various shaped windows for sample frame conditioning, including
+ * cosine windows (Hann etc) and triangular and rectangular windows.
+ */
 template <typename T>
 class Window
 {
 public:
     /**
-     * Construct a windower of the given type.
+     * Construct a windower of the given type and size. 
+     *
+     * Note that the cosine windows are periodic by design, rather
+     * than symmetrical. (A window of size N is equivalent to a
+     * symmetrical window of size N+1 with the final element missing.)
      */
-    Window(WindowType type, size_t size) : m_type(type), m_size(size) { encache(); }
+    Window(WindowType type, int size) : m_type(type), m_size(size) { encache(); }
     Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); }
     Window &operator=(const Window &w) {
 	if (&w == this) return *this;
@@ -46,79 +57,98 @@ public:
 	return *this;
     }
     virtual ~Window() { delete[] m_cache; }
-
+    
     void cut(T *src) const { cut(src, src); }
     void cut(const T *src, T *dst) const {
-	for (size_t i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
+	for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
     }
 
     WindowType getType() const { return m_type; }
-    size_t getSize() const { return m_size; }
+    int getSize() const { return m_size; }
+
+    std::vector<T> getWindowData() const {
+        std::vector<T> d;
+        for (int i = 0; i < m_size; ++i) {
+            d.push_back(m_cache[i]);
+        }
+        return d;
+    }
 
 protected:
     WindowType m_type;
-    size_t m_size;
+    int m_size;
     T *m_cache;
-
+    
     void encache();
 };
 
 template <typename T>
 void Window<T>::encache()
 {
-    size_t n = m_size;
+    int n = m_size;
     T *mult = new T[n];
-    size_t i;
+    int i;
     for (i = 0; i < n; ++i) mult[i] = 1.0;
 
     switch (m_type) {
-
+		
     case RectangularWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * 0.5;
+        for (i = 0; i < n; ++i) {
+            mult[i] = mult[i] * 0.5;
 	}
 	break;
-
+	    
     case BartlettWindow:
-	for (i = 0; i < n/2; ++i) {
-	    mult[i] = mult[i] * (i / T(n/2));
-	    mult[i + n/2] = mult[i + n/2] * (1.0 - (i / T(n/2)));
+        if (n == 2) {
+            mult[0] = mult[1] = 0; // "matlab compatible"
+        } else if (n == 3) {
+            mult[0] = 0;
+            mult[1] = mult[2] = 2./3.;
+        } else if (n > 3) {
+            for (i = 0; i < n/2; ++i) {
+                mult[i] = mult[i] * (i / T(n/2));
+                mult[i + n - n/2] = mult[i + n - n/2] * (1.0 - (i / T(n/2)));
+            }
 	}
 	break;
-
+	    
     case HammingWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n));
+        if (n > 1) {
+            for (i = 0; i < n; ++i) {
+                mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n));
+            }
 	}
 	break;
-
+	    
     case HanningWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n));
+        if (n > 1) {
+            for (i = 0; i < n; ++i) {
+                mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n));
+            }
 	}
 	break;
-
+	    
     case BlackmanWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n)
-				 + 0.08 * cos(4 * M_PI * i / n));
+        if (n > 1) {
+            for (i = 0; i < n; ++i) {
+                mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n)
+                                     + 0.08 * cos(4 * M_PI * i / n));
+            }
 	}
 	break;
-
-    case GaussianWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * exp((-1.0 / (n*n)) * ((T(2*i) - n) *
-						      (T(2*i) - n)));
-	}
-	break;
-
-    case ParzenWindow:
-	for (i = 0; i < n; ++i) {
-	    mult[i] = mult[i] * (1.0 - fabs((T(2*i) - n) / T(n + 1)));
+	    
+    case BlackmanHarrisWindow:
+        if (n > 1) {
+            for (i = 0; i < n; ++i) {
+                mult[i] = mult[i] * (0.35875
+                                     - 0.48829 * cos(2 * M_PI * i / n)
+                                     + 0.14128 * cos(4 * M_PI * i / n)
+                                     - 0.01168 * cos(6 * M_PI * i / n));
+            }
 	}
 	break;
     }
-
+	   
     m_cache = mult;
 }