c@131: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@131: 
c@131: #include "cq/CQSpectrogram.h"
c@131: 
c@131: #include "dsp/Window.h"
c@131: 
c@131: #include <cmath>
c@131: #include <vector>
c@132: #include <iostream>
c@131: 
c@131: using std::vector;
c@132: using std::cerr;
c@132: using std::endl;
c@131: 
c@131: #define BOOST_TEST_DYN_LINK
c@131: #define BOOST_TEST_MAIN
c@131: 
c@131: #include <boost/test/unit_test.hpp>
c@131: 
c@131: BOOST_AUTO_TEST_SUITE(TestCQFrequency)
c@131: 
c@131: // The principle here is to feed a single windowed sinusoid into a
c@131: // small CQ transform and check that the output has its peak bin at
c@133: // the correct frequency. 
c@131: 
c@131: // Set up fs/2 = 50, frequency range 10 -> 40 i.e. 2 octaves, fixed
c@131: // duration of 2 seconds
c@131: static const double sampleRate = 100;
c@139: static const double cqmin = 11.8921;
c@131: static const double cqmax = 40;
c@131: static const double bpo = 4;
c@131: static const int duration = sampleRate * 2;
c@131: 
c@132: // Threshold below which to ignore a column completely
c@131: static const double threshold = 0.08;
c@131: 
c@132: int
c@132: binForFrequency(double freq)
c@132: {
c@139:     int bin = (bpo * 2) - round(bpo * log2(freq / cqmin)) - 1;
c@132:     return bin;
c@132: }
c@132: 
c@131: void
c@135: checkCQFreqColumn(int i, vector<double> column,
c@135:                   double freq, CQSpectrogram::Interpolation interp)
c@131: {
c@132:     double maxval = 0.0;
c@132:     int maxidx = -1;
c@132:     int height = column.size();
c@139: 
c@139:     int nonZeroHeight = ((i % 2 == 1) ? height/2 : height);
c@139: 
c@139:     for (int j = 0; j < nonZeroHeight; ++j) {
c@132:         if (j == 0 || column[j] > maxval) {
c@132:             maxval = column[j];
c@132:             maxidx = j;
c@132:         }
c@132:     }
c@135: 
c@132:     int expected = binForFrequency(freq);
c@132:     if (maxval < threshold) {
c@132:         return; // ignore these columns at start and end
c@135:     } else if (expected < nonZeroHeight && maxidx != expected) {
c@139:         cerr << "ERROR: In column " << i << " with interpolation " << interp
c@139:              << ", maximum value for frequency " << freq
c@139:              << "\n       found at index " << maxidx
c@139:              << " (expected index " << expected << ")" << endl;
c@135:         cerr << "column contains: ";
c@135:         for (int j = 0; j < height; ++j) {
c@135:             cerr << column[j] << " ";
c@133:         }
c@135:         cerr << endl;
c@132:         BOOST_CHECK_EQUAL(maxidx, expected);
c@132:     }
c@131: }
c@131: 
c@131: void
c@144: testCQFrequencyWith(CQParameters params,
c@144:                     CQSpectrogram::Interpolation interp,
c@144:                     double freq)
c@144: {
c@144:     CQSpectrogram cq(params, interp);
c@144: 
c@144:     BOOST_CHECK_EQUAL(cq.getBinsPerOctave(), bpo);
c@144:     BOOST_CHECK_EQUAL(cq.getOctaves(), 2);
c@144:     BOOST_CHECK_CLOSE(cq.getBinFrequency(0), 40, 1e-10);
c@144:     BOOST_CHECK_CLOSE(cq.getBinFrequency(4), 20, 1e-10);
c@144:     BOOST_CHECK_CLOSE(cq.getBinFrequency(7), cqmin, 1e-3);
c@144:     
c@144:     vector<double> input;
c@144:     for (int i = 0; i < duration; ++i) {
c@144:         input.push_back(sin((i * 2 * M_PI * freq) / sampleRate));
c@144:     }
c@144:     Window<double>(HanningWindow, duration).cut(input.data());
c@144:     
c@144:     CQSpectrogram::RealBlock output = cq.process(input);
c@144:     CQSpectrogram::RealBlock rest = cq.getRemainingOutput();
c@144:     output.insert(output.end(), rest.begin(), rest.end());
c@144:     
c@144:     BOOST_CHECK_EQUAL(output[0].size(), 
c@144:                       cq.getBinsPerOctave() * cq.getOctaves());
c@144:     
c@144:     for (int i = 0; i < int(output.size()); ++i) {
c@144:         checkCQFreqColumn(i, output[i], freq, interp);
c@144:     }
c@144: }
c@144: 
c@144: void
c@131: testCQFrequency(double freq)
c@131: {
c@135:     vector<CQSpectrogram::Interpolation> interpolationTypes;
c@135:     interpolationTypes.push_back(CQSpectrogram::InterpolateZeros);
c@135:     interpolationTypes.push_back(CQSpectrogram::InterpolateHold);
c@135:     interpolationTypes.push_back(CQSpectrogram::InterpolateLinear);
c@131: 
c@135:     for (int k = 0; k < int(interpolationTypes.size()); ++k) {
c@135:         CQSpectrogram::Interpolation interp = interpolationTypes[k];
c@135:         CQParameters params(sampleRate, cqmin, cqmax, bpo);
c@144:         testCQFrequencyWith(params, interp, freq);
c@132:     }
c@131: }
c@131: 
c@133: BOOST_AUTO_TEST_CASE(freq_11) { testCQFrequency(11); }
c@140: BOOST_AUTO_TEST_CASE(freq_17) { testCQFrequency(17); }
c@140: BOOST_AUTO_TEST_CASE(freq_24) { testCQFrequency(24); }
c@139: BOOST_AUTO_TEST_CASE(freq_27) { testCQFrequency(27); }
c@139: BOOST_AUTO_TEST_CASE(freq_33) { testCQFrequency(33); }
c@131: BOOST_AUTO_TEST_CASE(freq_40) { testCQFrequency(40); }
c@131: 
c@131: BOOST_AUTO_TEST_SUITE_END()
c@131: