changeset 51:0997774f5fdc

Test cepstrum

--- a/CepstralPitchTracker.cpp	Tue Sep 11 20:42:52 2012 +0100
+++ b/CepstralPitchTracker.cpp	Tue Sep 11 21:55:05 2012 +0100
@@ -23,6 +23,7 @@
 */
 
 #include "CepstralPitchTracker.h"
+#include "Cepstrum.h"
 #include "MeanFilter.h"
 #include "PeakInterpolator.h"
 
@@ -263,43 +264,14 @@
 {
     FeatureSet fs;
 
-    int bs = m_blockSize;
-    int hs = m_blockSize/2 + 1;
-
-    double *rawcep = new double[bs];
-    double *io = new double[bs];
-    double *logmag = new double[bs];
-
-    // The "inverse symmetric" method. Seems to be the most reliable
-        
-    double magmean = 0.0;
-
-    for (int i = 0; i < hs; ++i) {
-
-	double power =
-	    inputBuffers[0][i*2  ] * inputBuffers[0][i*2  ] +
-	    inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1];
-	double mag = sqrt(power);
-
-        magmean += mag;
-
-	double lm = log(mag + 0.00000001);
-	
-	logmag[i] = lm;
-	if (i > 0) logmag[bs - i] = lm;
-    }
-
-    magmean /= hs;
-    double threshold = 0.1; // for magmean
-    
-    Vamp::FFT::inverse(bs, logmag, 0, rawcep, io);
-    
-    delete[] logmag;
-    delete[] io;
+    double *rawcep = new double[m_blockSize];
+    double magmean = Cepstrum(m_blockSize).process(inputBuffers[0], rawcep);
 
     int n = m_bins;
     double *data = new double[n];
-    MeanFilter(m_vflen).filterSubsequence(rawcep, data, m_blockSize, n, m_binFrom);
+    MeanFilter(m_vflen).filterSubsequence
+        (rawcep, data, m_blockSize, n, m_binFrom);
+
     delete[] rawcep;
 
     double maxval = 0.0;
@@ -332,6 +304,8 @@
     double peakfreq = m_inputSampleRate / (cimax + m_binFrom);
 
     double confidence = 0.0;
+    double threshold = 0.1; // for magmean
+
     if (nextPeakVal != 0.0) {
         confidence = (maxval - nextPeakVal) * 10.0;
         if (magmean < threshold) confidence = 0.0;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Cepstrum.h	Tue Sep 11 21:55:05 2012 +0100
@@ -0,0 +1,105 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#ifndef _CEPSTRUM_H_
+#define _CEPSTRUM_H_
+
+#include "vamp-sdk/FFT.h"
+#include <cmath>
+
+#include <iostream>
+#include <exception>
+
+class Cepstrum
+{
+public:
+    /**
+     * Construct a cepstrum converter based on an n-point FFT.
+     */
+    Cepstrum(int n) : m_n(n) {
+	if (n & (n-1)) {
+	    throw "N must be a power of two";
+	}
+    }
+    ~Cepstrum() { }
+    
+    /**
+     * Convert the given frequency-domain data to the cepstral domain.
+     *
+     * The input must be in the format used for FrequencyDomain data
+     * in the Vamp SDK: n/2+1 consecutive pairs of real and imaginary
+     * component floats corresponding to bins 0..(n/2) of the FFT
+     * output. Thus, n+2 values in total.
+     *
+     * The output consists of the raw cepstrum of length n.
+     *
+     * The cepstrum is calculated as the inverse FFT of a
+     * synthetically symmetrical base-10 log magnitude spectrum.
+     *
+     * Returns the mean magnitude of the input spectrum.
+     */
+    double process(const float *in, double *out) {
+
+	int hs = m_n/2 + 1;
+	double *io = new double[m_n];
+	double *logmag = new double[m_n];
+	double epsilon = 1e-10;
+
+	double magmean = 0.0;
+
+	for (int i = 0; i < hs; ++i) {
+
+	    double power = in[i*2] * in[i*2] + in[i*2+1] * in[i*2+1];
+	    double mag = sqrt(power);
+	    magmean += mag;
+
+	    logmag[i] = log10(mag + epsilon);
+
+	    if (i > 0) {
+		// make the log magnitude spectrum symmetrical
+		logmag[m_n - i] = logmag[i];
+	    }
+	}
+	
+	magmean /= hs;
+	/*
+	std::cerr << "logmags:" << std::endl;
+	for (int i = 0; i < m_n; ++i) {
+	    std::cerr << logmag[i] << " ";
+	}
+	std::cerr << std::endl;
+	*/
+	Vamp::FFT::inverse(m_n, logmag, 0, out, io);
+    
+	delete[] logmag;
+	delete[] io;
+
+	return magmean;
+    }
+
+private:
+    int m_n;
+};
+
+#endif

--- a/Makefile.inc	Tue Sep 11 20:42:52 2012 +0100
+++ b/Makefile.inc	Tue Sep 11 21:55:05 2012 +0100
@@ -24,11 +24,12 @@
 
 PLUGIN_MAIN := libmain.cpp
 
-TESTS := test/test-notehypothesis \
-         test/test-meanfilter \
+TESTS := test/test-meanfilter \
          test/test-fft \
-         test/test-peakinterpolator
-
+	 test/test-cepstrum \
+         test/test-peakinterpolator \
+	 test/test-notehypothesis
+         
 OBJECTS := $(SOURCES:.cpp=.o)
 OBJECTS := $(OBJECTS:.c=.o)
 
@@ -46,6 +47,9 @@
 test/test-meanfilter: test/TestMeanFilter.o $(OBJECTS)
 	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
 
+test/test-cepstrum: test/TestCepstrum.o $(OBJECTS)
+	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
+
 test/test-fft: test/TestFFT.o $(OBJECTS)
 	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
 
@@ -58,11 +62,18 @@
 distclean:	clean
 		rm -f $(PLUGIN) $(TESTS)
 
+depend:
+		makedepend -Y -fMakefile.inc *.cpp test/*.cpp *.h test/*.h
+
 # DO NOT DELETE
 
-CepstralPitchTracker.o: CepstralPitchTracker.h NoteHypothesis.h
+CepstralPitchTracker.o: CepstralPitchTracker.h NoteHypothesis.h Cepstrum.h
+CepstralPitchTracker.o: MeanFilter.h PeakInterpolator.h
 libmain.o: CepstralPitchTracker.h NoteHypothesis.h
 NoteHypothesis.o: NoteHypothesis.h
 PeakInterpolator.o: PeakInterpolator.h
+test/TestCepstrum.o: Cepstrum.h
+test/TestMeanFilter.o: MeanFilter.h
 test/TestNoteHypothesis.o: NoteHypothesis.h
 test/TestPeakInterpolator.o: PeakInterpolator.h
+CepstralPitchTracker.o: NoteHypothesis.h

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestCepstrum.cpp	Tue Sep 11 21:55:05 2012 +0100
@@ -0,0 +1,89 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#include "Cepstrum.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestCepstrum)
+
+BOOST_AUTO_TEST_CASE(cosine)
+{
+    // input format is re0, im0, re1, im1...
+    float in[] = { 0,0, 10,0, 0,0 };
+    double out[4];
+    double mm = Cepstrum(4).process(in, out);
+    BOOST_CHECK_SMALL(out[0] - (-4.5), 1e-10);
+    BOOST_CHECK_EQUAL(out[1], 0);
+    BOOST_CHECK_SMALL(out[2] - (-5.5), 1e-10);
+    BOOST_CHECK_EQUAL(out[3], 0);
+    BOOST_CHECK_EQUAL(mm, 10.0/3.0);
+}
+
+BOOST_AUTO_TEST_CASE(symmetry)
+{
+    // Cepstrum output bins 1..n-1 are symmetric about bin n/2
+    float in[] = { 1,2,3,4,5,6,7,8,9,10 };
+    double out[8];
+    double mm = Cepstrum(8).process(in, out);
+    BOOST_CHECK_SMALL(out[1] - out[7], 1e-14);
+    BOOST_CHECK_SMALL(out[2] - out[6], 1e-14);
+    BOOST_CHECK_SMALL(out[3] - out[5], 1e-14);
+}
+
+BOOST_AUTO_TEST_CASE(oneHarmonic)
+{
+    // input format is re0, im0, re1, im1, re2, im2
+    // freq for bin i is i * samplerate / n
+    // freqs:        0  sr/n  2sr/n  3sr/n  4sr/n  5sr/n  6sr/n  7sr/n  sr/2
+    float in[] = { 0,0,  0,0,  10,0,   0,0,  10,0,   0,0,  10,0,   0,0,  0,0 };
+    double out[16];
+    double mm = Cepstrum(16).process(in, out);
+    BOOST_CHECK_EQUAL(mm, 30.0/9.0);
+    // peak is at 8
+    BOOST_CHECK(out[8] > 0);
+    // odd bins are all zero
+    BOOST_CHECK_EQUAL(out[1], 0);
+    BOOST_CHECK_EQUAL(out[3], 0);
+    BOOST_CHECK_EQUAL(out[5], 0);
+    BOOST_CHECK_EQUAL(out[7], 0);
+    BOOST_CHECK_EQUAL(out[9], 0);
+    BOOST_CHECK_EQUAL(out[11], 0);
+    BOOST_CHECK_EQUAL(out[13], 0);
+    BOOST_CHECK_EQUAL(out[15], 0);
+    // the rest are negative
+    BOOST_CHECK(out[0] < 0);
+    BOOST_CHECK(out[2] < 0);
+    BOOST_CHECK(out[4] < 0);
+    BOOST_CHECK(out[6] < 0);
+    BOOST_CHECK(out[10] < 0);
+    BOOST_CHECK(out[12] < 0);
+    BOOST_CHECK(out[14] < 0);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- a/test/TestFFT.cpp	Tue Sep 11 20:42:52 2012 +0100
+++ b/test/TestFFT.cpp	Tue Sep 11 21:55:05 2012 +0100
@@ -112,24 +112,6 @@
     COMPARE_ARRAY(back, in);
 }
 
-BOOST_AUTO_TEST_CASE(sineCosineDC)
-{
-    // Sine and cosine mixed
-    double in[] = { 0.5, 1.0, -0.5, -1.0 };
-    double re[4], im[4];
-    Vamp::FFT::forward(4, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 0.0);
-    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
-    BOOST_CHECK_EQUAL(re[2], 0.0);
-    BOOST_CHECK_EQUAL(im[0], 0.0);
-    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
-    BOOST_CHECK_EQUAL(im[2], 0.0);
-    double back[4];
-    double backim[4];
-    Vamp::FFT::inverse(4, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-
 BOOST_AUTO_TEST_CASE(nyquist)
 {
     double in[] = { 1, -1, 1, -1 };