c@92: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@92: 
c@92: /*
c@92:     QM Vamp Plugin Set
c@92: 
c@92:     Centre for Digital Music, Queen Mary, University of London.
c@92:     All rights reserved.
c@92: */
c@92: 
c@92: #include "AdaptiveSpectrogram.h"
c@92: 
c@92: #include <cstdlib>
c@92: #include <cstring>
c@92: 
c@92: #include <iostream>
c@92: 
c@92: #include <dsp/transforms/FFT.h>
c@92: 
c@92: using std::string;
c@92: using std::vector;
c@92: using std::cerr;
c@92: using std::endl;
c@92: 
c@92: using Vamp::RealTime;
c@92: 
c@99: //#define DEBUG_VERBOSE 1
c@99: 
c@92: AdaptiveSpectrogram::AdaptiveSpectrogram(float inputSampleRate) :
c@92:     Plugin(inputSampleRate),
c@104:     m_w(8),
c@104:     m_n(3),
c@109:     m_threaded(true),
c@109:     m_threadsInUse(false)
c@92: {
c@92: }
c@92: 
c@92: AdaptiveSpectrogram::~AdaptiveSpectrogram()
c@92: {
c@104:     for (int i = 0; i < m_cutThreads.size(); ++i) {
c@104:         delete m_cutThreads[i];
c@104:     }
c@104:     m_cutThreads.clear();
c@105: 
c@110:     for (FFTMap::iterator i = m_fftThreads.begin();
c@110:          i != m_fftThreads.end(); ++i) {
c@106:         delete i->second;
c@105:     }
c@105:     m_fftThreads.clear();
c@92: }
c@92: 
c@92: string
c@92: AdaptiveSpectrogram::getIdentifier() const
c@92: {
c@93:     return "qm-adaptivespectrogram";
c@92: }
c@92: 
c@92: string
c@92: AdaptiveSpectrogram::getName() const
c@92: {
c@92:     return "Adaptive Spectrogram";
c@92: }
c@92: 
c@92: string
c@92: AdaptiveSpectrogram::getDescription() const
c@92: {
c@92:     return "Produce an adaptive spectrogram by adaptive selection from spectrograms at multiple resolutions";
c@92: }
c@92: 
c@92: string
c@92: AdaptiveSpectrogram::getMaker() const
c@92: {
c@92:     return "Queen Mary, University of London";
c@92: }
c@92: 
c@92: int
c@92: AdaptiveSpectrogram::getPluginVersion() const
c@92: {
c@92:     return 1;
c@92: }
c@92: 
c@92: string
c@92: AdaptiveSpectrogram::getCopyright() const
c@92: {
c@92:     return "Plugin by Wen Xue and Chris Cannam.  Copyright (c) 2009 Wen Xue and QMUL - All Rights Reserved";
c@92: }
c@92: 
c@92: size_t
c@92: AdaptiveSpectrogram::getPreferredStepSize() const
c@92: {
c@92:     return ((2 << m_w) << m_n) / 2;
c@92: }
c@92: 
c@92: size_t
c@92: AdaptiveSpectrogram::getPreferredBlockSize() const
c@92: {
c@92:     return (2 << m_w) << m_n;
c@92: }
c@92: 
c@92: bool
c@92: AdaptiveSpectrogram::initialise(size_t channels, size_t stepSize, size_t blockSize)
c@92: {
c@92:     if (channels < getMinChannelCount() ||
c@92: 	channels > getMaxChannelCount()) return false;
c@92: 
c@92:     return true;
c@92: }
c@92: 
c@92: void
c@92: AdaptiveSpectrogram::reset()
c@92: {
c@92: 
c@92: }
c@92: 
c@92: AdaptiveSpectrogram::ParameterList
c@92: AdaptiveSpectrogram::getParameterDescriptors() const
c@92: {
c@92:     ParameterList list;
c@92: 
c@92:     ParameterDescriptor desc;
c@92:     desc.identifier = "n";
c@92:     desc.name = "Number of resolutions";
c@92:     desc.description = "Number of consecutive powers of two to use as spectrogram resolutions, starting with the minimum resolution specified";
c@92:     desc.unit = "";
c@92:     desc.minValue = 1;
c@92:     desc.maxValue = 10;
c@104:     desc.defaultValue = 4;
c@92:     desc.isQuantized = true;
c@92:     desc.quantizeStep = 1;
c@92:     list.push_back(desc);
c@92: 
c@92:     ParameterDescriptor desc2;
c@92:     desc2.identifier = "w";
c@92:     desc2.name = "Smallest resolution";
c@92:     desc2.description = "Smallest of the consecutive powers of two to use as spectrogram resolutions";
c@92:     desc2.unit = "";
c@92:     desc2.minValue = 1;
c@92:     desc2.maxValue = 14;
c@104:     desc2.defaultValue = 9;
c@92:     desc2.isQuantized = true;
c@92:     desc2.quantizeStep = 1;
c@92:     // I am so lazy
c@92:     desc2.valueNames.push_back("2");
c@92:     desc2.valueNames.push_back("4");
c@92:     desc2.valueNames.push_back("8");
c@92:     desc2.valueNames.push_back("16");
c@92:     desc2.valueNames.push_back("32");
c@92:     desc2.valueNames.push_back("64");
c@92:     desc2.valueNames.push_back("128");
c@92:     desc2.valueNames.push_back("256");
c@92:     desc2.valueNames.push_back("512");
c@92:     desc2.valueNames.push_back("1024");
c@92:     desc2.valueNames.push_back("2048");
c@92:     desc2.valueNames.push_back("4096");
c@92:     desc2.valueNames.push_back("8192");
c@92:     desc2.valueNames.push_back("16384");
c@92:     list.push_back(desc2);
c@92: 
c@109:     ParameterDescriptor desc3;
c@109:     desc3.identifier = "threaded";
c@109:     desc3.name = "Multi-threaded processing";
c@110:     desc3.description = "Perform calculations using several threads in parallel";
c@109:     desc3.unit = "";
c@109:     desc3.minValue = 0;
c@109:     desc3.maxValue = 1;
c@109:     desc3.defaultValue = 1;
c@109:     desc3.isQuantized = true;
c@109:     desc3.quantizeStep = 1;
c@109:     list.push_back(desc3);
c@109: 
c@92:     return list;
c@92: }
c@92: 
c@92: float
c@92: AdaptiveSpectrogram::getParameter(std::string id) const
c@92: {
c@92:     if (id == "n") return m_n+1;
c@92:     else if (id == "w") return m_w+1;
c@109:     else if (id == "threaded") return (m_threaded ? 1 : 0);
c@92:     return 0.f;
c@92: }
c@92: 
c@92: void
c@92: AdaptiveSpectrogram::setParameter(std::string id, float value)
c@92: {
c@92:     if (id == "n") {
c@92:         int n = lrintf(value);
c@92:         if (n >= 1 && n <= 10) m_n = n-1;
c@92:     } else if (id == "w") {
c@92:         int w = lrintf(value);
c@92:         if (w >= 1 && w <= 14) m_w = w-1;
c@109:     } else if (id == "threaded") {
c@109:         m_threaded = (value > 0.5);
c@109:     }
c@92: }
c@92: 
c@92: AdaptiveSpectrogram::OutputList
c@92: AdaptiveSpectrogram::getOutputDescriptors() const
c@92: {
c@92:     OutputList list;
c@92: 
c@92:     OutputDescriptor d;
c@92:     d.identifier = "output";
c@92:     d.name = "Output";
c@92:     d.description = "The output of the plugin";
c@92:     d.unit = "";
c@92:     d.hasFixedBinCount = true;
c@92:     d.binCount = ((2 << m_w) << m_n) / 2;
c@92:     d.hasKnownExtents = false;
c@92:     d.isQuantized = false;
c@92:     d.sampleType = OutputDescriptor::FixedSampleRate;
c@92:     d.sampleRate = m_inputSampleRate / ((2 << m_w) / 2);
c@92:     d.hasDuration = false;
c@112:     char name[20];
c@112:     for (int i = 0; i < d.binCount; ++i) {
c@112:         float freq = (m_inputSampleRate / d.binCount) * (i + 1); // no DC bin
c@112:         sprintf(name, "%d Hz", int(freq));
c@112:         d.binNames.push_back(name);
c@112:     }
c@92:     list.push_back(d);
c@92: 
c@92:     return list;
c@92: }
c@92: 
c@92: AdaptiveSpectrogram::FeatureSet
c@92: AdaptiveSpectrogram::getRemainingFeatures()
c@92: {
c@92:     FeatureSet fs;
c@92:     return fs;
c@92: }
c@92: 
c@100: AdaptiveSpectrogram::FeatureSet
c@100: AdaptiveSpectrogram::process(const float *const *inputBuffers, RealTime ts)
c@100: {
c@100:     FeatureSet fs;
c@100: 
c@100:     int minwid = (2 << m_w), maxwid = ((2 << m_w) << m_n);
c@100: 
c@101: #ifdef DEBUG_VERBOSE
c@100:     cerr << "widths from " << minwid << " to " << maxwid << " ("
c@100:          << minwid/2 << " to " << maxwid/2 << " in real parts)" << endl;
c@101: #endif
c@100: 
c@100:     Spectrograms s(minwid/2, maxwid/2, 1);
c@100: 
c@100:     int w = minwid;
c@100:     int index = 0;
c@100: 
c@100:     while (w <= maxwid) {
c@106:         if (m_fftThreads.find(w) == m_fftThreads.end()) {
c@106:             m_fftThreads[w] = new FFTThread(w);
c@106:         }
c@109:         if (m_threaded) {
c@109:             m_fftThreads[w]->startCalculation(inputBuffers[0], s, index, maxwid);
c@109:         } else {
c@109:             m_fftThreads[w]->setParameters(inputBuffers[0], s, index, maxwid);
c@109:             m_fftThreads[w]->performTask();
c@109:         }
c@100:         w *= 2;
c@100:         ++index;
c@100:     }
c@100: 
c@109:     if (m_threaded) {
c@109:         w = minwid;
c@109:         while (w <= maxwid) {
c@109:             m_fftThreads[w]->await();
c@109:             w *= 2;
c@109:         }
c@105:     }
c@102: 
c@109:     m_threadsInUse = false;
c@104: 
c@111:     std::cerr << "maxwid/2 = " << maxwid/2 << ", minwid/2 = " << minwid/2 << ", n+1 = " << m_n+1 << ", 2^(n+1) = " << (2<<m_n) << std::endl;
c@110: 
c@110:     Cutting *cutting = cut(s, maxwid/2, 0, 0, maxwid/2, 0);
c@100: 
c@101: #ifdef DEBUG_VERBOSE
c@100:     printCutting(cutting, "  ");
c@101: #endif
c@100: 
c@100:     vector<vector<float> > rmat(maxwid/minwid);
c@100:     for (int i = 0; i < maxwid/minwid; ++i) {
c@100:         rmat[i] = vector<float>(maxwid/2);
c@100:     }
c@100:     
c@100:     assemble(s, cutting, rmat, 0, 0, maxwid/minwid, maxwid/2);
c@100: 
c@110:     cutting->erase();
c@100: 
c@100:     for (int i = 0; i < rmat.size(); ++i) {
c@100:         Feature f;
c@100:         f.hasTimestamp = false;
c@100:         f.values = rmat[i];
c@100:         fs[0].push_back(f);
c@100:     }
c@100: 
c@104: //    std::cerr << "process returning!\n" << std::endl;
c@104: 
c@100:     return fs;
c@100: }
c@100: 
c@100: void
c@104: AdaptiveSpectrogram::printCutting(Cutting *c, string pfx) const
c@100: {
c@100:     if (c->first) {
c@100:         if (c->cut == Cutting::Horizontal) {
c@100:             cerr << pfx << "H" << endl;
c@100:         } else if (c->cut == Cutting::Vertical) {
c@100:             cerr << pfx << "V" << endl;
c@100:         }
c@100:         printCutting(c->first, pfx + "  ");
c@100:         printCutting(c->second, pfx + "  ");
c@100:     } else {
c@100:         cerr << pfx << "* " << c->value << endl;
c@100:     }
c@100: }
c@100: 
c@104: void
c@104: AdaptiveSpectrogram::getSubCuts(const Spectrograms &s,
c@104:                                 int res,
c@104:                                 int x, int y, int h,
c@104:                                 Cutting *&top, Cutting *&bottom,
c@113:                                 Cutting *&left, Cutting *&right,
c@113:                                 BlockAllocator *allocator) const
c@104: {
c@109:     if (m_threaded && !m_threadsInUse) {
c@104: 
c@109:         m_threadsInUse = true;
c@104: 
c@104:         if (m_cutThreads.empty()) {
c@104:             for (int i = 0; i < 4; ++i) {
c@104:                 CutThread *t = new CutThread(this);
c@104:                 m_cutThreads.push_back(t);
c@104:             }
c@104:         }
c@104: 
c@109:         // Cut threads 0 and 1 calculate the top and bottom halves;
c@110:         // threads 2 and 3 calculate left and right.  See notes in
c@110:         // unthreaded code below for more information.
c@104: 
c@110:         m_cutThreads[0]->cut(s, res, x, y + h/2, h/2);        // top
c@110:         m_cutThreads[1]->cut(s, res, x, y, h/2);              // bottom
c@110:         m_cutThreads[2]->cut(s, res/2, 2 * x, y/2, h/2);      // left
c@110:         m_cutThreads[3]->cut(s, res/2, 2 * x + 1, y/2, h/2);  // right
c@104: 
c@104:         top    = m_cutThreads[0]->get();
c@104:         bottom = m_cutThreads[1]->get();
c@104:         left   = m_cutThreads[2]->get();
c@104:         right  = m_cutThreads[3]->get();
c@104: 
c@104:     } else {
c@104: 
c@110:         // Unthreaded version
c@104: 
c@104:         // The "vertical" division is a top/bottom split.
c@104:         // Splitting this way keeps us in the same resolution,
c@104:         // but with two vertical subregions of height h/2.
c@104: 
c@113:         top    = cut(s, res, x, y + h/2, h/2, allocator);
c@113:         bottom = cut(s, res, x, y, h/2, allocator);
c@104: 
c@104:         // The "horizontal" division is a left/right split.  Splitting
c@104:         // this way places us in resolution res/2, which has lower
c@104:         // vertical resolution but higher horizontal resolution.  We
c@104:         // need to double x accordingly.
c@104:         
c@113:         left   = cut(s, res/2, 2 * x, y/2, h/2, allocator);
c@113:         right  = cut(s, res/2, 2 * x + 1, y/2, h/2, allocator);
c@104:     }
c@104: }
c@104: 
c@100: AdaptiveSpectrogram::Cutting *
c@100: AdaptiveSpectrogram::cut(const Spectrograms &s,
c@100:                          int res,
c@110:                          int x, int y, int h,
c@110:                          BlockAllocator *allocator) const
c@100: {
c@100: //    cerr << "res = " << res << ", x = " << x << ", y = " << y << ", h = " << h << endl;
c@100: 
c@110:     Cutting *cutting;
c@110:     if (allocator) {
c@110:         cutting = (Cutting *)(allocator->allocate());
c@110:         cutting->allocator = allocator;
c@110:     } else {
c@110:         cutting = new Cutting;
c@110:         cutting->allocator = 0;
c@110:     }
c@110: 
c@100:     if (h > 1 && res > s.minres) {
c@100: 
c@104:         Cutting *top = 0, *bottom = 0, *left = 0, *right = 0;
c@113:         getSubCuts(s, res, x, y, h, top, bottom, left, right, allocator);
c@100: 
c@100:         double vcost = top->cost + bottom->cost;
c@100:         double hcost = left->cost + right->cost;
c@100: 
c@101:         bool normalize = true;
c@101: 
c@101:         if (normalize) {
c@101: 
c@101:             double venergy = top->value + bottom->value;
c@101:             vcost = (vcost + (venergy * log(venergy))) / venergy;
c@101: 
c@101:             double henergy = left->value + right->value;
c@101:             hcost = (hcost + (henergy * log(henergy))) / henergy;
c@101:         }            
c@101: 
c@100:         if (vcost > hcost) {
c@100: 
c@100:             // cut horizontally (left/right)
c@100:             cutting->cut = Cutting::Horizontal;
c@100:             cutting->first = left;
c@100:             cutting->second = right;
c@100:             cutting->cost = hcost;
c@111:             cutting->value = left->value + right->value;
c@110:             top->erase();
c@110:             bottom->erase();
c@100:             return cutting;
c@100: 
c@100:         } else {
c@100: 
c@110:             // cut vertically (top/bottom)
c@100:             cutting->cut = Cutting::Vertical;
c@100:             cutting->first = top;
c@100:             cutting->second = bottom;
c@100:             cutting->cost = vcost;
c@111:             cutting->value = top->value + bottom->value;
c@110:             left->erase();
c@110:             right->erase();
c@100:             return cutting;
c@100:         }
c@100: 
c@100:     } else {
c@100: 
c@100:         // no cuts possible from this level
c@100: 
c@100:         cutting->cut = Cutting::Finished;
c@100:         cutting->first = 0;
c@100:         cutting->second = 0;
c@100: 
c@100:         int n = 0;
c@100:         for (int r = res; r > s.minres; r /= 2) ++n;
c@100:         const Spectrogram *spectrogram = s.spectrograms[n];
c@100: 
c@100:         cutting->cost = cost(*spectrogram, x, y);
c@100:         cutting->value = value(*spectrogram, x, y);
c@100: 
c@100: //        cerr << "cost for this cell: " << cutting->cost << endl;
c@100: 
c@100:         return cutting;
c@100:     }
c@100: }
c@100: 
c@100: void
c@100: AdaptiveSpectrogram::assemble(const Spectrograms &s,
c@100:                               const Cutting *cutting,
c@100:                               vector<vector<float> > &rmat,
c@104:                               int x, int y, int w, int h) const
c@100: {
c@100:     switch (cutting->cut) {
c@100: 
c@100:     case Cutting::Finished:
c@100:         for (int i = 0; i < w; ++i) {
c@100:             for (int j = 0; j < h; ++j) {
c@112:                 rmat[x+i][y+j] = cutting->value * cutting->value;
c@100:             }
c@100:         }
c@100:         return;
c@100: 
c@100:     case Cutting::Horizontal:
c@100:         assemble(s, cutting->first, rmat, x, y, w/2, h);
c@100:         assemble(s, cutting->second, rmat, x+w/2, y, w/2, h);
c@100:         break;
c@100:         
c@100:     case Cutting::Vertical:
c@100:         assemble(s, cutting->first, rmat, x, y+h/2, w, h/2);
c@100:         assemble(s, cutting->second, rmat, x, y, w, h/2);
c@100:         break;
c@100:     }        
c@100: }
c@100: