# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1381424021 -3600
# Node ID 5d9489187abdccee14504119e9f1cbfa910f9e29
# Parent  a067c2eeb13c3049e4b940ae012fcf510f1e7929
Add sinc

diff -r a067c2eeb13c -r 5d9489187abd base/SincWindow.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/SincWindow.cpp	Thu Oct 10 17:53:41 2013 +0100
@@ -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 -r a067c2eeb13c -r 5d9489187abd base/SincWindow.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/SincWindow.h	Thu Oct 10 17:53:41 2013 +0100
@@ -0,0 +1,57 @@
+/* -*- 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>
+
+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 -r a067c2eeb13c -r 5d9489187abd qm-dsp.pro
--- a/qm-dsp.pro	Tue Oct 08 17:23:17 2013 +0100
+++ b/qm-dsp.pro	Thu Oct 10 17:53:41 2013 +0100
@@ -45,6 +45,7 @@
 HEADERS += base/Pitch.h \
            base/Window.h \
            base/KaiserWindow.h \
+           base/SincWindow.h \
            dsp/chromagram/Chromagram.h \
            dsp/chromagram/ConstantQ.h \
            dsp/keydetection/GetKeyMode.h \
@@ -84,6 +85,7 @@
            thread/Thread.h
 SOURCES += base/Pitch.cpp \
            base/KaiserWindow.cpp \
+           base/SincWindow.cpp \
            dsp/chromagram/Chromagram.cpp \
            dsp/chromagram/ConstantQ.cpp \
            dsp/keydetection/GetKeyMode.cpp \
diff -r a067c2eeb13c -r 5d9489187abd tests/TestWindow.cpp
--- a/tests/TestWindow.cpp	Tue Oct 08 17:23:17 2013 +0100
+++ b/tests/TestWindow.cpp	Thu Oct 10 17:53:41 2013 +0100
@@ -2,6 +2,7 @@
 
 #include "base/Window.h"
 #include "base/KaiserWindow.h"
+#include "base/SincWindow.h"
 
 #include <iostream>
 
@@ -152,6 +153,38 @@
 }
 
 //!!! todo: tests for kaiser with attenuation and bandwidth parameters
+
+template <int N>
+void testSinc(double p, const double expected[N])
+{
+    double d[N];
+    for (int i = 0; i < N; ++i) d[i] = 1.0;
+    SincWindow w(N, p);
+    w.cut(d);
+    COMPARE_ARRAY(d, expected);
+
+    double d0[N], d1[N];
+    for (int i = 0; i < N; ++i) d0[i] = 0.5 + (1.0 / (N * 2)) * (i + 1);
+    w.cut(d0, d1);
+    for (int i = 0; i < N; ++i) {
+	BOOST_CHECK_SMALL(d1[i] - d0[i] * expected[i], 1e-4);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(sinc)
+{
+    double e1[] = { 0, 0, 1, 0, 0 };
+    testSinc<5>(1, e1);
+
+    double e2[] = { 0, 0, 1, 0, 0 };
+    testSinc<5>(2, e2);
+
+    double e3[] = { -0.2122, 0.0, 0.6366, 1.0, 0.6366, 0, -0.2122 };
+    testSinc<7>(4, e3);
+
+    double e4[] = { -0.2122, 0, 0.6366, 1, 0.6366, 0 };
+    testSinc<6>(4, e4);
+}
         
 BOOST_AUTO_TEST_SUITE_END()