Mercurial > hg > qm-dsp

changeset 124:263181813eec

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Kaiser window (and some tests for it)
author Chris Cannam
date Fri, 04 Oct 2013 18:46:32 +0100
parents a37635bbb2c1
children 5351b5e9ad9f
files base/KaiserWindow.cpp base/KaiserWindow.h qm-dsp.pro tests/TestWindow.cpp
diffstat 4 files changed, 198 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /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);
+    }
+}
--- /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
--- 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 \
--- 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()