# HG changeset patch
# User Chris Cannam
# Date 1380908792 -3600
# Node ID 263181813eec4bbd6c901baf241ef2e65924f3e4
# Parent  a37635bbb2c1aa9aba827d45f10808c0bce97bf7
Add Kaiser window (and some tests for it)

diff -r a37635bbb2c1 -r 263181813eec base/KaiserWindow.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/KaiserWindow.cpp	Fri Oct 04 18:46:32 2013 +0100
@@ -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.
+*/
+
+#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);
+    for (int i = 0; i < m_length; ++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);
+    }
+}
diff -r a37635bbb2c1 -r 263181813eec base/KaiserWindow.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/KaiserWindow.h	Fri Oct 04 18:46:32 2013 +0100
@@ -0,0 +1,100 @@
+/* -*- 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>
+
+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 -r a37635bbb2c1 -r 263181813eec qm-dsp.pro
--- a/qm-dsp.pro	Fri Oct 04 18:34:30 2013 +0100
+++ b/qm-dsp.pro	Fri Oct 04 18:46:32 2013 +0100
@@ -44,6 +44,7 @@
 # Input
 HEADERS += base/Pitch.h \
            base/Window.h \
+           base/KaiserWindow.h \
            dsp/chromagram/Chromagram.h \
            dsp/chromagram/ConstantQ.h \
            dsp/keydetection/GetKeyMode.h \
@@ -82,6 +83,7 @@
            thread/BlockAllocator.h \
            thread/Thread.h
 SOURCES += base/Pitch.cpp \
+           base/KaiserWindow.cpp \
            dsp/chromagram/Chromagram.cpp \
            dsp/chromagram/ConstantQ.cpp \
            dsp/keydetection/GetKeyMode.cpp \
diff -r a37635bbb2c1 -r 263181813eec tests/TestWindow.cpp
--- a/tests/TestWindow.cpp	Fri Oct 04 18:34:30 2013 +0100
+++ b/tests/TestWindow.cpp	Fri Oct 04 18:46:32 2013 +0100
@@ -1,6 +1,7 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
 
 #include "base/Window.h"
+#include "base/KaiserWindow.h"
 
 #include <iostream>
 
@@ -117,6 +118,40 @@
     };
     testWindow<10>(BlackmanWindow, e10);
 }
+
+BOOST_AUTO_TEST_CASE(kaiser_4_10)
+{
+    double d[10];
+    for (int i = 0; i < 10; ++i) d[i] = 1.0;
+    double e[] = { 
+        0.0885, 0.2943, 0.5644, 0.8216, 0.9789,
+        0.9789, 0.8216, 0.5644, 0.2943, 0.0885
+    };
+    KaiserWindow::Parameters p;
+    p.length = 10;
+    p.beta = 4;
+    KaiserWindow k(p);
+    k.cut(d);
+    COMPARE_ARRAY(d, e);
+}
     
+BOOST_AUTO_TEST_CASE(kaiser_2p5_11)
+{
+    double d[11];
+    for (int i = 0; i < 11; ++i) d[i] = 1.0;
+    double e[] = { 
+        0.3040, 0.5005, 0.6929, 0.8546, 0.9622,
+        1.0000, 0.9622, 0.8546, 0.6929, 0.5005, 0.3040
+    };
+    KaiserWindow::Parameters p;
+    p.length = 11;
+    p.beta = 2.5;
+    KaiserWindow k(p);
+    k.cut(d);
+    COMPARE_ARRAY(d, e);
+}
+
+//!!! todo: tests for kaiser with attenuation and bandwidth parameters
+        
 BOOST_AUTO_TEST_SUITE_END()