Mercurial > hg > tipic

changeset 18:c785eaaeac40

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add DCT implementation
author Chris Cannam
date Thu, 24 Sep 2015 12:45:43 +0100
parents a8ca8a90257c
children 51eb8b1a1910
files Makefile.inc src/DCT.cpp src/DCT.h src/PitchFilterbank.cpp src/test-dct.cpp src/test-filter.cpp
diffstat 6 files changed, 263 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile.inc	Thu Aug 20 17:22:25 2015 +0100
+++ b/Makefile.inc	Thu Sep 24 12:45:43 2015 +0100
@@ -24,38 +24,59 @@
 
 PUBLIC_HEADERS	:=
 
-LIB_HEADERS	:= $(SRC_DIR)/delays.h $(SRC_DIR)/filter-a.h $(SRC_DIR)/filter-b.h $(SRC_DIR)/Filter.h $(SRC_DIR)/PitchFilterbank.h
-LIB_SOURCES	:= $(SRC_DIR)/Filter.cpp $(SRC_DIR)/PitchFilterbank.cpp
+LIB_HEADERS	:= $(SRC_DIR)/delays.h $(SRC_DIR)/filter-a.h $(SRC_DIR)/filter-b.h $(SRC_DIR)/Filter.h $(SRC_DIR)/PitchFilterbank.h $(SRC_DIR)/DCT.h
+LIB_SOURCES	:= $(SRC_DIR)/Filter.cpp $(SRC_DIR)/PitchFilterbank.cpp $(SRC_DIR)/DCT.cpp
 LIB_OBJECTS	:= $(LIB_SOURCES:.cpp=.o)
 LIB_OBJECTS	:= $(LIB_OBJECTS:.c=.o)
 
 PLUGIN_HEADERS  := $(SRC_DIR)/TipicVampPlugin.h
 PLUGIN_SOURCES  := $(SRC_DIR)/TipicVampPlugin.cpp $(SRC_DIR)/libmain.cpp
+PLUGIN_OBJECTS	:= $(PLUGIN_SOURCES:.cpp=.o)
+PLUGIN_OBJECTS	:= $(PLUGIN_OBJECTS:.c=.o)
 
 BQVEC_HEADERS	:= $(BQVEC_DIR)/Allocators.h $(BQVEC_DIR)/Restrict.h $(BQVEC_DIR)/VectorOps.h
 BQVEC_SOURCES	:= $(BQVEC_DIR)/src/Allocators.cpp
+BQVEC_OBJECTS	:= $(BQVEC_SOURCES:.cpp=.o)
+BQVEC_OBJECTS	:= $(BQVEC_OBJECTS:.c=.o)
+
+TEST_SOURCES	:= $(SRC_DIR)/test-filter.cpp $(SRC_DIR)/test-dct.cpp
+TEST_OBJECTS	:= $(TEST_SOURCES:.cpp=.o)
+TEST_OBJECTS	:= $(TEST_OBJECTS:.c=.o)
+
+OTHER_OBJECTS	:= $(CQ_DIR)/src/dsp/FFT.o $(CQ_DIR)/src/ext/kissfft/kiss_fft.o $(CQ_DIR)/src/ext/kissfft/tools/kiss_fftr.o
 
 HEADERS	     := $(PUBLIC_HEADERS) $(LIB_HEADERS) $(PLUGIN_HEADERS) $(BQVEC_HEADERS)
-SOURCES	     := $(PUBLIC_SOURCES) $(LIB_SOURCES) $(PLUGIN_SOURCES) $(BQVEC_SOURCES)
+SOURCES	     := $(PUBLIC_SOURCES) $(LIB_SOURCES) $(PLUGIN_SOURCES) $(BQVEC_SOURCES) $(TEST_SOURCES)
 OBJECTS	     := $(SOURCES:.cpp=.o)
 OBJECTS	     := $(OBJECTS:.c=.o)
 
 LIBS	:= $(CQ_DIR)/libcq.a $(VAMPSDK_DIR)/libvamp-sdk.a
 
-all: constant-q-cpp $(LIBRARY) $(PLUGIN)
+all: constant-q-cpp $(LIBRARY) $(PLUGIN) tests
 
 .PHONY: constant-q-cpp
 constant-q-cpp: 
 	$(MAKE) -C $@ -f Makefile$(MAKEFILE_EXT) libcq.a
 
-$(PLUGIN):	$(OBJECTS) $(LIBS)
+$(PLUGIN):	$(PLUGIN_OBJECTS) $(LIB_OBJECTS) $(BQVEC_OBJECTS) $(LIBS)
 	$(CXX) -o $@ $^ $(LIBS) $(PLUGIN_LDFLAGS)
 
-$(LIBRARY):    $(LIB_OBJECTS)
+$(LIBRARY):    $(LIB_OBJECTS) $(BQVEC_OBJECTS) $(OTHER_OBJECTS)
 	$(RM) -f $@
 	$(AR) cr $@ $^
 	$(RANLIB) $@
 
+.PHONY:	tests
+tests:	test-dct test-filter
+	./test-dct
+	./test-filter
+
+test-dct:	  $(TEST_OBJECTS) $(LIBRARY)
+	$(CXX) -o $@ src/test-dct.o $(LIBRARY)
+
+test-filter:	  $(TEST_OBJECTS) $(LIBRARY)
+	$(CXX) -o $@ src/test-filter.o $(LIBRARY)
+
 clean:		
 	rm -f $(OBJECTS)
 	$(MAKE) -C constant-q-cpp -f Makefile$(MAKEFILE_EXT) clean
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DCT.cpp	Thu Sep 24 12:45:43 2015 +0100
@@ -0,0 +1,78 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#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);
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DCT.h	Thu Sep 24 12:45:43 2015 +0100
@@ -0,0 +1,72 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#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
--- a/src/PitchFilterbank.cpp	Thu Aug 20 17:22:25 2015 +0100
+++ b/src/PitchFilterbank.cpp	Thu Sep 24 12:45:43 2015 +0100
@@ -31,6 +31,13 @@
 	//!!! todo: tuning frequency adjustment
 	// * resample input by a small amount
 	// * adjust output block timings by a small amount
+
+	//!!! nb "we use forward-backward filtering such that the
+	// resulting output signal has precisely zero phase distortion
+	// and a magnitude modified by the square of the filter’s
+	// magnitude response" -- we are not doing forward-backward
+	// here & so need to adapt magnitudes in compensation to match
+	// original
 	
 	m_resamplers[882] = new Resampler(sampleRate, 882);
 	m_resamplers[4410] = new Resampler(sampleRate, 4410);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/test-dct.cpp	Thu Sep 24 12:45:43 2015 +0100
@@ -0,0 +1,72 @@
+
+#include "DCT.h"
+
+#include <iostream>
+#include <cmath>
+
+using namespace std;
+
+void check(string context, vector<double> got, vector<double> expected, bool &good)
+{
+//    cerr << "Checking " << context << "..." << endl;
+    double thresh = 1e-4;
+    for (int i = 0; i < int(got.size()); ++i) {
+	if (fabs(got[i] - expected[i]) > thresh) {
+	    cerr << "ERROR: " << context << "[" << i << "] (" << got[i]
+		 << ") differs from expected " << expected[i] << endl;
+	    good = false;
+	}
+    }
+}
+
+int main(int argc, char **argv)
+{
+    bool good = true;
+    
+    vector<double> in4 { 1, 2, 3, 5 };
+    vector<double> out4(4), inv4(4);
+    vector<double> expected4nu { 22.0, -8.1564, 1.4142, -1.2137 };
+    vector<double> expected4u { 5.5, -2.8837, 0.5, -0.4291 };
+
+    DCT d4(4);
+
+    d4.forward(in4.data(), out4.data());
+    check("out4", out4, expected4nu, good);
+
+    d4.inverse(out4.data(), inv4.data());
+    check("inverse4", inv4, in4, good);
+
+    d4.forwardUnitary(in4.data(), out4.data());
+    check("out4u", out4, expected4u, good);
+
+    d4.inverseUnitary(out4.data(), inv4.data());
+    check("inverse4u", inv4, in4, good);
+
+    // do it again, just in case some internal state was modified in inverse
+    
+    d4.forwardUnitary(in4.data(), out4.data());
+    check("out4u", out4, expected4u, good);
+
+    d4.inverseUnitary(out4.data(), inv4.data());
+    check("inverse4u", inv4, in4, good);
+    
+    vector<double> in5 { 1, 2, 3, 5, 6 };
+    vector<double> out5(5), inv5(5);
+    vector<double> expected5u { 7.6026, -4.1227, 0.3162, -0.0542, -0.3162 };
+
+    DCT d5(5);
+
+    d5.forwardUnitary(in5.data(), out5.data());
+    check("out5u", out5, expected5u, good);
+
+    d5.inverseUnitary(out5.data(), inv5.data());
+    check("inverse5u", inv5, in5, good);
+    
+    if (good) {
+	cerr << "Success" << endl;
+	return 0;
+    } else {
+	return 1;
+    }
+}
+
--- a/src/test-filter.cpp	Thu Aug 20 17:22:25 2015 +0100
+++ b/src/test-filter.cpp	Thu Sep 24 12:45:43 2015 +0100
@@ -31,6 +31,12 @@
 	    good = false;
 	}
     }
-    if (good) cerr << "SUCCESS" << endl;
+
+    if (good) {
+	cerr << "Success" << endl;
+	return 0;
+    } else {
+	return 1;
+    }
 }