# HG changeset patch
# User Chris Cannam
# Date 1444208135 -3600
# Node ID 857ca50ca25f075edfa4cc083ea7a1f932d657a4
# Parent  e4a57215ddee51b617a8fb3557d11ff3ddb82bac
Add DCT

diff -r e4a57215ddee -r 857ca50ca25f .hgignore
--- a/.hgignore	Mon Sep 28 12:33:17 2015 +0100
+++ b/.hgignore	Wed Oct 07 09:55:35 2015 +0100
@@ -5,5 +5,7 @@
 *.a
 doc/html/
 tests/test-*
+tests/test.log
 *.rej
+ext/uncertain/
 
diff -r e4a57215ddee -r 857ca50ca25f build/general/Makefile.inc
--- a/build/general/Makefile.inc	Mon Sep 28 12:33:17 2015 +0100
+++ b/build/general/Makefile.inc	Wed Oct 07 09:55:35 2015 +0100
@@ -38,6 +38,7 @@
            dsp/tonal/ChangeDetectionFunction.h \
            dsp/tonal/TCSgram.h \
            dsp/tonal/TonalEstimator.h \
+           dsp/transforms/DCT.h \
            dsp/transforms/FFT.h \
            dsp/wavelet/Wavelet.h \
            hmm/hmm.h \
@@ -83,6 +84,7 @@
            dsp/tonal/ChangeDetectionFunction.cpp \
            dsp/tonal/TCSgram.cpp \
            dsp/tonal/TonalEstimator.cpp \
+           dsp/transforms/DCT.cpp \
            dsp/transforms/FFT.cpp \
            dsp/wavelet/Wavelet.cpp \
            hmm/hmm.c \
diff -r e4a57215ddee -r 857ca50ca25f dsp/transforms/DCT.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/transforms/DCT.cpp	Wed Oct 07 09:55:35 2015 +0100
@@ -0,0 +1,91 @@
+/* -*- 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 file by Chris Cannam.
+
+    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 "DCT.h"
+
+#include <cmath>
+
+DCT::DCT(int n) :
+    m_n(n),
+    m_scaled(n, 0.0),
+    m_time2(n*4, 0.0),
+    m_freq2r(n*4, 0.0),
+    m_freq2i(n*4, 0.0),
+    m_fft(n*4)
+{
+    m_scale = m_n * sqrt(2.0 / m_n);
+}
+
+DCT::~DCT()
+{
+}
+
+void
+DCT::forward(const double *in, double *out)
+{
+    for (int i = 0; i < m_n; ++i) {
+	m_time2[i*2 + 1] = in[i];
+	m_time2[m_n*4 - i*2 - 1] = in[i];
+    }
+
+    m_fft.forward(m_time2.data(), m_freq2r.data(), m_freq2i.data());
+
+    for (int i = 0; i < m_n; ++i) {
+	out[i] = m_freq2r[i];
+    }
+}
+
+void
+DCT::forwardUnitary(const double *in, double *out)
+{
+    forward(in, out);
+    for (int i = 0; i < m_n; ++i) {
+	out[i] /= m_scale;
+    }
+    out[0] /= sqrt(2.0);
+}
+
+void
+DCT::inverse(const double *in, double *out)
+{
+    for (int i = 0; i < m_n; ++i) {
+	m_freq2r[i] = in[i];
+    }
+    for (int i = 0; i < m_n; ++i) {
+	m_freq2r[m_n*2 - i] = -in[i];
+    }
+    m_freq2r[m_n] = 0.0;
+
+    for (int i = 0; i <= m_n*2; ++i) {
+	m_freq2i[i] = 0.0;
+    }
+    
+    m_fft.inverse(m_freq2r.data(), m_freq2i.data(), m_time2.data());
+
+    for (int i = 0; i < m_n; ++i) {
+	out[i] = m_time2[i*2 + 1];
+    }
+}
+
+void
+DCT::inverseUnitary(const double *in, double *out)
+{
+    for (int i = 0; i < m_n; ++i) {
+	m_scaled[i] = in[i] * m_scale;
+    }
+    m_scaled[0] *= sqrt(2.0);
+    inverse(m_scaled.data(), out);
+}
+
diff -r e4a57215ddee -r 857ca50ca25f dsp/transforms/DCT.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/transforms/DCT.h	Wed Oct 07 09:55:35 2015 +0100
@@ -0,0 +1,85 @@
+/* -*- 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 file by Chris Cannam.
+
+    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 DCT_H
+#define DCT_H
+
+#include "FFT.h"
+
+#include <vector>
+
+class DCT
+{
+public:
+    /**
+     * Construct a DCT object to calculate the Discrete Cosine
+     * Transform given input of size n samples. The transform is
+     * implemented using an FFT of size 4n, for simplicity.
+     */
+    DCT(int n);
+
+    ~DCT();
+
+    /**
+     * Carry out a type-II DCT of size n, where n is the value
+     * provided to the constructor above.
+     *
+     * The in and out pointers must point to (enough space for) n
+     * values.
+     */
+    void forward(const double *in, double *out);
+
+    /**
+     * Carry out a type-II unitary DCT of size n, where n is the value
+     * provided to the constructor above. This is a scaled version of
+     * the type-II DCT corresponding to a transform with an orthogonal
+     * matrix. This is the transform implemented by the dct() function
+     * in MATLAB.
+     *
+     * The in and out pointers must point to (enough space for) n
+     * values.
+     */
+    void forwardUnitary(const double *in, double *out);
+
+    /**
+     * Carry out a type-III (inverse) DCT of size n, where n is the
+     * value provided to the constructor above.
+     *
+     * The in and out pointers must point to (enough space for) n
+     * values.
+     */
+    void inverse(const double *in, double *out);
+
+    /**
+     * Carry out a type-III (inverse) unitary DCT of size n, where n
+     * is the value provided to the constructor above. This is the
+     * inverse of forwardUnitary().
+     *
+     * The in and out pointers must point to (enough space for) n
+     * values.
+     */
+    void inverseUnitary(const double *in, double *out);
+
+private:
+    int m_n;
+    double m_scale;
+    std::vector<double> m_scaled;
+    std::vector<double> m_time2;
+    std::vector<double> m_freq2r;
+    std::vector<double> m_freq2i;
+    FFTReal m_fft;
+};
+
+#endif
diff -r e4a57215ddee -r 857ca50ca25f dsp/transforms/FFT.h
--- a/dsp/transforms/FFT.h	Mon Sep 28 12:33:17 2015 +0100
+++ b/dsp/transforms/FFT.h	Wed Oct 07 09:55:35 2015 +0100
@@ -4,6 +4,12 @@
     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 FFT_H
diff -r e4a57215ddee -r 857ca50ca25f tests/Makefile
--- a/tests/Makefile	Mon Sep 28 12:33:17 2015 +0100
+++ b/tests/Makefile	Wed Oct 07 09:55:35 2015 +0100
@@ -1,11 +1,11 @@
 
 CFLAGS := -I.. $(CFLAGS)
-CXXFLAGS := -I.. -Wall -Wextra -g $(CXXFLAGS)
+CXXFLAGS := -I.. -Wall -Wextra -std=c++11 -g $(CXXFLAGS)
 
 LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework -lpthread
 LIBS := ../libqm-dsp.a
 
-TESTS	:= test-mathutilities test-window test-fft test-pvoc test-resampler test-medianfilter
+TESTS	:= test-mathutilities test-window test-fft test-dct test-pvoc test-resampler test-medianfilter
 
 VG	:= valgrind -q
 
@@ -24,6 +24,9 @@
 test-fft: TestFFT.o $(LIBS)
 	$(CXX) -o $@ $^ $(LDFLAGS)
 
+test-dct: TestDCT.o $(LIBS)
+	$(CXX) -o $@ $^ $(LDFLAGS)
+
 test-pvoc: TestPhaseVocoder.o $(LIBS)
 	$(CXX) -o $@ $^ $(LDFLAGS)
 
@@ -34,6 +37,7 @@
 TestMedianFilter.o: $(LIBS)
 TestWindow.o: $(LIBS)
 TestFFT.o: $(LIBS)
+TestDCT.o: $(LIBS)
 TestPhaseVocoder.o: $(LIBS)
 TestResampler.o: $(LIBS)
 
diff -r e4a57215ddee -r 857ca50ca25f tests/TestDCT.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/TestDCT.cpp	Wed Oct 07 09:55:35 2015 +0100
@@ -0,0 +1,196 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "dsp/transforms/DCT.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+#include <stdexcept>
+
+using namespace std;
+
+BOOST_AUTO_TEST_SUITE(TestDCT)
+
+BOOST_AUTO_TEST_CASE(forwardArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double out[] = { 999, 999, 999, 999, 999, 999 };
+    DCT(4).forward(in, out+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(out[0], 999.0);
+    BOOST_CHECK_EQUAL(out[5], 999.0);
+    delete[] in;
+}
+
+BOOST_AUTO_TEST_CASE(u_forwardArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double out[] = { 999, 999, 999, 999, 999, 999 };
+    DCT(4).forwardUnitary(in, out+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(out[0], 999.0);
+    BOOST_CHECK_EQUAL(out[5], 999.0);
+    delete[] in;
+}
+
+BOOST_AUTO_TEST_CASE(inverseArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double out[] = { 999, 999, 999, 999, 999, 999 };
+    DCT(4).inverse(in, out+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(out[0], 999.0);
+    BOOST_CHECK_EQUAL(out[5], 999.0);
+    delete[] in;
+
+}
+
+BOOST_AUTO_TEST_CASE(u_inverseArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell if
+    // they haven't been written; and allocate the inputs on the heap
+    // so that, if running under valgrind, we get warnings about
+    // overruns
+    double *in = new double[4];
+    in[0] = 1;
+    in[1] = 1;
+    in[2] = -1;
+    in[3] = -1;
+    double out[] = { 999, 999, 999, 999, 999, 999 };
+    DCT(4).inverseUnitary(in, out+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(out[0], 999.0);
+    BOOST_CHECK_EQUAL(out[5], 999.0);
+    delete[] in;
+}
+
+BOOST_AUTO_TEST_CASE(nonU4)
+{
+    int n = 4;
+    vector<double> in { 1, 2, 3, 5 };
+    vector<double> out(n);
+    vector<double> inv(n);
+    vector<double> expected { 22.0, -8.1564, 1.4142, -1.2137 };
+    
+    DCT d(n);
+
+    // our expected values are very approximate
+    double thresh = 1e-4;
+
+    d.forward(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+
+    d.inverse(out.data(), inv.data());
+
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(in[i] - inv[i], thresh);
+    }
+    
+    // do it again, just in case some internal state was modified in inverse
+    
+    d.forward(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(u4)
+{
+    int n = 4;
+    vector<double> in { 1, 2, 3, 5 };
+    vector<double> out(n);
+    vector<double> inv(n);
+    vector<double> expected { 5.5, -2.8837, 0.5, -0.4291 };
+    
+    DCT d(n);
+
+    // our expected values are very approximate
+    double thresh = 1e-4;
+
+    d.forwardUnitary(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+
+    d.inverseUnitary(out.data(), inv.data());
+
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(in[i] - inv[i], thresh);
+    }
+    
+    // do it again, just in case some internal state was modified in inverse
+    
+    d.forwardUnitary(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(u5)
+{
+    int n = 5;
+    vector<double> in { 1, 2, 3, 5, 6 };
+    vector<double> out(n);
+    vector<double> inv(n);
+    vector<double> expected { 7.6026, -4.1227, 0.3162, -0.0542, -0.3162 };
+    
+    DCT d(n);
+
+    // our expected values are very approximate
+    double thresh = 1e-4;
+
+    d.forwardUnitary(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+
+    d.inverseUnitary(out.data(), inv.data());
+
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(in[i] - inv[i], thresh);
+    }
+    
+    // do it again, just in case some internal state was modified in inverse
+    
+    d.forwardUnitary(in.data(), out.data());
+    
+    for (int i = 0; i < n; ++i) {
+	BOOST_CHECK_SMALL(expected[i] - out[i], thresh);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+