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@135: 
c@135:     This program is free software; you can redistribute it and/or
c@135:     modify it under the terms of the GNU General Public License as
c@135:     published by the Free Software Foundation; either version 2 of the
c@135:     License, or (at your option) any later version.  See the file
c@135:     COPYING included with this distribution for more information.
c@92: */
c@92: 
c@92: #include "AdaptiveSpectrogram.h"
c@92: 
c@92: #include <cstdlib>
c@133: #include <cstdio>
c@92: #include <cstring>
c@114: #include <cfloat>
c@92: 
c@92: #include <iostream>
c@92: 
c@92: #include <dsp/transforms/FFT.h>
c@156: #include <dsp/rateconversion/Decimator.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@114:     m_n(2),
c@114:     m_coarse(false),
c@109:     m_threaded(true),
c@156:     m_decFactor(1),
c@178:     m_buffer(0),
c@156:     m_buflen(0),
c@178:     m_decimator(0),
c@178:     m_threadsInUse(false)
c@92: {
c@92: }
c@92: 
c@92: AdaptiveSpectrogram::~AdaptiveSpectrogram()
c@92: {
c@178:     for (int i = 0; i < int(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@156: 
c@156:     delete[] m_buffer;
c@156:     delete m_decimator;
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@156:     return 2;
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@156:     return getPreferredBlockSize();
c@92: }
c@92: 
c@92: size_t
c@92: AdaptiveSpectrogram::getPreferredBlockSize() const
c@92: {
c@156:     // the /2 at the end is for 50% overlap (we handle framing ourselves)
c@156:     return ((2 << m_w) << m_n) * m_decFactor / 2;
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@156:     if (blockSize != getPreferredBlockSize()) {
c@156:         std::cerr << "AdaptiveSpectrogram::initialise: Block size " << blockSize << " does not match required block size of " << getPreferredBlockSize() << std::endl;
c@156:         return false;
c@156:     }
c@156:     if (stepSize != getPreferredStepSize()) {
c@156:         std::cerr << "AdaptiveSpectrogram::initialise: Step size " << stepSize << " does not match required step size of " << getPreferredStepSize() << std::endl;
c@156:         return false;
c@156:     }
c@156: 
c@156:     if (m_decFactor > 1) {
c@156:         m_decimator = new Decimator(blockSize, m_decFactor);
c@156:     }
c@156: 
c@156:     m_buflen = (blockSize * 2) / m_decFactor; // *2 for 50% overlap
c@156:     m_buffer = new float[m_buflen];
c@156: 
c@156:     reset();
c@156:     
c@92:     return true;
c@92: }
c@92: 
c@92: void
c@92: AdaptiveSpectrogram::reset()
c@92: {
c@156:     if (m_decimator) m_decimator->resetFilter();
c@156:     for (int i = 0; i < m_buflen; ++i) m_buffer[i] = 0.f;
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@114:     desc.description = "Number of consecutive powers of two in the range to be used as spectrogram resolutions, starting with the minimum resolution specified";
c@92:     desc.unit = "";
c@114:     desc.minValue = 2;
c@92:     desc.maxValue = 10;
c@114:     desc.defaultValue = 3;
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@156:     desc2.identifier = "dec";
c@156:     desc2.name = "Decimation factor";
c@156:     desc2.description = "Factor to down-sample by, increasing speed but lowering maximum frequency";
c@156:     desc2.unit = "";
c@156:     desc2.minValue = 0;
c@156:     desc2.maxValue = 3;
c@156:     desc2.defaultValue = 0;
c@156:     desc2.isQuantized = true;
c@156:     desc2.quantizeStep = 1;
c@156:     desc2.valueNames.clear();
c@156:     desc2.valueNames.push_back("No decimation");
c@156:     desc2.valueNames.push_back("2");
c@156:     desc2.valueNames.push_back("4");
c@156:     desc2.valueNames.push_back("8");
c@156:     list.push_back(desc2);
c@156: 
c@109:     ParameterDescriptor desc3;
c@114:     desc3.identifier = "coarse";
c@114:     desc3.name = "Omit alternate resolutions";
c@114:     desc3.description = "Generate a coarser spectrogram faster by excluding every alternate resolution (first and last resolution are always retained)";
c@114:     desc3.unit = "";
c@114:     desc3.minValue = 0;
c@114:     desc3.maxValue = 1;
c@114:     desc3.defaultValue = 0;
c@114:     desc3.isQuantized = true;
c@114:     desc3.quantizeStep = 1;
c@114:     list.push_back(desc3);
c@114: 
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@114:     else if (id == "coarse") return (m_coarse ? 1 : 0);
c@156:     else if (id == "dec") {
c@156:         int f = m_decFactor;
c@156:         int p = 0;
c@156:         while (f > 1) {
c@156:             f >>= 1;
c@156:             p += 1;
c@156:         }
c@156:         return p;
c@156:     }
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@114:     } else if (id == "coarse") {
c@114:         m_coarse = (value > 0.5);
c@156:     } else if (id == "dec") {
c@156:         int p = lrintf(value);
c@156:         if (p >= 0 && p <= 3) m_decFactor = 1 << p;
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@156:     d.binCount = getPreferredBlockSize() / (m_decFactor * 2);
c@92:     d.hasKnownExtents = false;
c@92:     d.isQuantized = false;
c@92:     d.sampleType = OutputDescriptor::FixedSampleRate;
c@156:     d.sampleRate = m_inputSampleRate / (m_decFactor * ((2 << m_w) / 2));
c@92:     d.hasDuration = false;
c@112:     char name[20];
c@178:     for (int i = 0; i < int(d.binCount); ++i) {
c@156:         float freq = (m_inputSampleRate / (m_decFactor * (d.binCount * 2)) * (i + 1)); // no DC bin
c@157:         sprintf(name, "%.1f Hz", 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@178: AdaptiveSpectrogram::process(const float *const *inputBuffers, RealTime)
c@100: {
c@156:     // framing: shift and write the new data to right half
c@156:     for (int i = 0; i < m_buflen/2; ++i) {
c@156:         m_buffer[i] = m_buffer[m_buflen/2 + i];
c@156:     }
c@156: 
c@156:     if (m_decFactor == 1) {
c@156:         for (int i = 0; i < m_buflen/2; ++i) {
c@156:             m_buffer[m_buflen/2 + i] = inputBuffers[0][i];
c@156:         }
c@156:     } else {
c@156:         m_decimator->process(inputBuffers[0], m_buffer + m_buflen/2);
c@156:     }
c@156: 
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@114: 
c@114:         if (!isResolutionWanted(s, w/2)) {
c@114:             w *= 2;
c@114:             ++index;
c@114:             continue;
c@114:         }
c@114: 
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@156:             m_fftThreads[w]->startCalculation(m_buffer, s, index, maxwid);
c@109:         } else {
c@156:             m_fftThreads[w]->setParameters(m_buffer, 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@114:         index = 0;
c@109:         while (w <= maxwid) {
c@114:             if (!isResolutionWanted(s, w/2)) {
c@114:                 w *= 2;
c@114:                 ++index;
c@114:                 continue;
c@114:             }
c@109:             m_fftThreads[w]->await();
c@109:             w *= 2;
c@114:             ++index;
c@109:         }
c@105:     }
c@102: 
c@109:     m_threadsInUse = false;
c@104: 
c@114: //    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@114:     int cutwid = maxwid/2;
c@114:     Cutting *cutting = cut(s, cutwid, 0, 0, cutwid, 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@114:     assemble(s, cutting, rmat, 0, 0, maxwid/minwid, cutwid);
c@100: 
c@110:     cutting->erase();
c@100: 
c@178:     for (int i = 0; i < int(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@114:                                 Cutting **top, Cutting **bottom,
c@114:                                 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@114:         if (top) m_cutThreads[0]->cut(s, res, x, y + h/2, h/2);
c@114:         if (bottom) m_cutThreads[1]->cut(s, res, x, y, h/2);
c@104: 
c@114:         if (left) m_cutThreads[2]->cut(s, res/2, 2 * x, y/2, h/2);
c@114:         if (right) m_cutThreads[3]->cut(s, res/2, 2 * x + 1, y/2, h/2);
c@114: 
c@114:         if (top) *top = m_cutThreads[0]->get();
c@114:         if (bottom) *bottom = m_cutThreads[1]->get();
c@114:         if (left) *left = m_cutThreads[2]->get();
c@114:         if (right) *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@114:         if (top) *top = cut(s, res, x, y + h/2, h/2, allocator);
c@114:         if (bottom) *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@114:         if (left) *left = cut(s, res/2, 2 * x, y/2, h/2, allocator);
c@114:         if (right) *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@114:         if (!isResolutionWanted(s, res)) {
c@100: 
c@114:             Cutting *left = 0, *right = 0;
c@114:             getSubCuts(s, res, x, y, h, 0, 0, &left, &right, allocator);
c@114:                 
c@114:             double hcost = left->cost + right->cost;
c@101:             double henergy = left->value + right->value;
c@114:             hcost = normalize(hcost, henergy);
c@114:                 
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@100: 
c@114:         } else if (h == 2 && !isResolutionWanted(s, res/2)) {
c@100: 
c@114:             Cutting *top = 0, *bottom = 0;
c@114:             getSubCuts(s, res, x, y, h, &top, &bottom, 0, 0, allocator);
c@114:                 
c@114:             double vcost = top->cost + bottom->cost;
c@114:             double venergy = top->value + bottom->value;
c@114:             vcost = normalize(vcost, venergy);
c@114:                 
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@114: 
c@114:         } else {
c@114: 
c@114:             Cutting *top = 0, *bottom = 0, *left = 0, *right = 0;
c@114:             getSubCuts(s, res, x, y, h, &top, &bottom, &left, &right, allocator);
c@114: 
c@114:             double vcost = top->cost + bottom->cost;
c@114:             double venergy = top->value + bottom->value;
c@114:             vcost = normalize(vcost, venergy);
c@114:             
c@114:             double hcost = left->cost + right->cost;
c@114:             double henergy = left->value + right->value;
c@114:             hcost = normalize(hcost, henergy);
c@114:             
c@114:             if (vcost > hcost) {
c@114:                 cutting->cut = Cutting::Horizontal;
c@114:                 cutting->first = left;
c@114:                 cutting->second = right;
c@114:                 cutting->cost = hcost;
c@114:                 cutting->value = left->value + right->value;
c@114:                 top->erase();
c@114:                 bottom->erase();
c@114:                 return cutting;
c@114:             } else {
c@114:                 cutting->cut = Cutting::Vertical;
c@114:                 cutting->first = top;
c@114:                 cutting->second = bottom;
c@114:                 cutting->cost = vcost;
c@114:                 cutting->value = top->value + bottom->value;
c@114:                 left->erase();
c@114:                 right->erase();
c@114:                 return cutting;
c@114:             }
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@114:         for (int r = res; r > s.minres; r >>= 1) ++n;
c@100:         const Spectrogram *spectrogram = s.spectrograms[n];
c@100:         cutting->cost = cost(*spectrogram, x, y);
c@100:         cutting->value = value(*spectrogram, x, y);
c@114:     }
c@100: 
c@114:     return cutting;
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@114:                 rmat[x+i][y+j] = 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: