c@178: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@178: 
c@178: /*
c@178:   QM Vamp Plugin Set
c@178: 
c@178:   Centre for Digital Music, Queen Mary, University of London.
c@178:   This file copyright 2009 Thomas Wilmering.
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@178:     COPYING included with this distribution for more information.
c@178: */
c@178: 
c@178: #include "DWT.h"
c@178: 
c@178: #include <cmath>
c@178: 
c@178: using std::string;
c@178: using std::vector;
c@178: using std::cerr;
c@178: using std::endl;
c@178: 
c@178: DWT::DWT(float inputSampleRate) :
c@178:     Plugin(inputSampleRate),
c@178:     m_stepSize(0),
c@178:     m_blockSize(0)
c@178: {
c@178:     m_scales = 10;
c@178:     m_flength = 0;
c@178:     m_wavelet = Wavelet::Haar;
c@178:     m_threshold = 0;
c@178:     m_absolute = 0;
c@178: }
c@178: 
c@178: DWT::~DWT()
c@178: {
c@178: }
c@178: 
c@178: string
c@178: DWT::getIdentifier() const
c@178: {
c@178:     return "qm-dwt";
c@178: }
c@178: 
c@178: string
c@178: DWT::getName() const
c@178: {
c@178:     return "Discrete Wavelet Transform";
c@178: }
c@178: 
c@178: string
c@178: DWT::getDescription() const
c@178: {
c@178:     return "Visualisation by scalogram";
c@178: }
c@178: 
c@178: string
c@178: DWT::getMaker() const
c@178: {
c@178:     return "Queen Mary, University of London";
c@178: }
c@178: 
c@178: int
c@178: DWT::getPluginVersion() const
c@178: {
c@178:     return 1;
c@178: }
c@178: 
c@178: string
c@178: DWT::getCopyright() const
c@178: {
c@178:     return "Plugin by Thomas Wilmering.  Copyright (c) 2009 Thomas Wilmering and QMUL - All Rights Reserved";
c@178: }
c@178: 
c@178: size_t 
c@178: DWT::getPreferredBlockSize() const 
c@178: { 
c@178:     size_t s = (1 << m_scales);
c@178:     while (s < 1024) s *= 2;
c@178:     return s;
c@178: } 
c@178: 
c@178: size_t 
c@178: DWT::getPreferredStepSize() const 
c@178: { 
c@178:     return 0;  
c@178: } 
c@178: 
c@178: bool
c@178: DWT::initialise(size_t channels, size_t stepSize, size_t blockSize)
c@178: {
c@178:     if (channels < getMinChannelCount() ||
c@178:         channels > getMaxChannelCount()) return false;
c@210: 
c@210:     if (m_inputSampleRate > 1000000) { // somewhat arbitrarily
c@210:         std::cerr << "DWT::initialise: ERROR: Maximum sample rate exceeded" << std::endl;
c@210:         return false;
c@210:     }
c@178: 	
c@178:     if ((1U << m_scales) > blockSize) {
c@178:         std::cerr << "DWT::initialise: ERROR: Block size must be at least 2^scales (specified block size " << blockSize << " < " << (1 << m_scales) << ")" << std::endl;
c@178:         return false;
c@178:     }
c@178: 
c@178:     m_stepSize = stepSize;
c@178:     m_blockSize = blockSize;
c@178: 
c@178:     Wavelet::createDecompositionFilters(m_wavelet, m_lpd, m_hpd);
c@178: 
c@178:     m_flength = m_lpd.size(); // or m_hpd.size()
c@178: 	
c@178:     m_samplePass.resize(m_scales);				// resize buffer for samples to pass to next block
c@178: 	
c@178:     for (int i=0; i<m_scales; ++i) {
c@178:         m_samplePass[i].resize(m_flength-2, 0.0);
c@178:     }
c@178: 
c@178:     return true;	
c@178: }
c@178: 
c@178: void
c@178: DWT::reset()
c@178: {
c@178:     m_samplePass.clear();
c@178: 
c@178:     m_samplePass.resize(m_scales);
c@178: 	
c@178:     for (int i=0; i<m_scales; ++i) {
c@178:         m_samplePass[i].resize(m_flength-2, 0.0);
c@178:     }
c@178: }
c@178: 
c@178: DWT::OutputList
c@178: DWT::getOutputDescriptors() const
c@178: {
c@178:     OutputList list;
c@178: 	
c@178:     OutputDescriptor sg;
c@178:     sg.identifier = "wcoeff";
c@178:     sg.name = "Wavelet Coefficients";
c@178:     sg.description = "Wavelet coefficients";
c@178:     sg.unit = "";
c@178:     sg.hasFixedBinCount = true;              // depends on block size
c@178:     sg.binCount = m_scales; 	// number of scales
c@178:     sg.hasKnownExtents = false;
c@178:     sg.isQuantized = false;
c@178:     sg.sampleType = OutputDescriptor::FixedSampleRate; 
c@178:     sg.sampleRate = .5 * m_inputSampleRate;
c@178: 	
c@178:     list.push_back(sg);
c@178: 	
c@178:     return list;
c@178: }
c@178: 
c@178: 
c@178: DWT::ParameterList
c@178: DWT::getParameterDescriptors() const
c@178: {
c@178:     ParameterList list;
c@178:  
c@178:     ParameterDescriptor d;
c@178:     d.identifier = "scales";
c@178:     d.name = "Scales";
c@178:     d.description = "Scale depth";
c@178:     d.unit = "";
c@178:     d.minValue = 1.0f;
c@178:     d.maxValue = 16.0f;
c@178:     d.defaultValue = 10.0f;
c@178:     d.isQuantized = true;
c@178:     d.quantizeStep = 1.0f;
c@178:     list.push_back(d);
c@178:  
c@178:     d.identifier = "wavelet";
c@178:     d.name = "Wavelet";
c@178:     d.description = "Wavelet type to use";
c@178:     d.unit = "";
c@178:     d.minValue = 0.f;
c@178:     d.maxValue = int(Wavelet::LastType);
c@178:     d.defaultValue = int(Wavelet::Haar);
c@178:     d.isQuantized = true;
c@178:     d.quantizeStep = 1.0f;
c@178: 
c@178:     for (int i = 0; i <= int(Wavelet::LastType); ++i) {
c@178:         d.valueNames.push_back(Wavelet::getWaveletName(Wavelet::Type(i)));
c@178:     }
c@178:     list.push_back(d);
c@178:     d.valueNames.clear();
c@178: 
c@178:     d.identifier = "threshold";
c@178:     d.name = "Threshold";
c@178:     d.description = "Wavelet coefficient threshold";
c@178:     d.unit = "";
c@178:     d.minValue = 0.0f;
c@178:     d.maxValue = 0.01f;
c@178:     d.defaultValue = 0.0f;
c@178:     d.isQuantized = false;
c@178:     list.push_back(d);
c@178: 
c@178:     d.identifier = "absolute";
c@178:     d.name = "Absolute values";
c@178:     d.description = "Return absolute values";
c@178:     d.unit = "";
c@178:     d.minValue = 0.0f;
c@178:     d.maxValue = 1.00f;
c@178:     d.defaultValue = 0.0f;
c@178:     d.isQuantized = true;
c@178:     d.quantizeStep = 1.0f;
c@178:     list.push_back(d);
c@178:  
c@178:     return list;
c@178: }
c@178:  
c@178: void DWT::setParameter(std::string paramid, float newval)
c@178: {
c@178:     if (paramid == "scales") {
c@178:         m_scales = newval;
c@178:     }
c@178:     else if (paramid == "wavelet") {
c@178:         m_wavelet = (Wavelet::Type)(int(newval + 0.1));
c@178:     }
c@178:     else if (paramid == "threshold") {
c@178:         m_threshold = newval;
c@178:     }
c@178:     else if (paramid == "absolute") {
c@178:         m_absolute = newval;
c@178:     }
c@178: }
c@178:  
c@178: float DWT::getParameter(std::string paramid) const
c@178: {
c@178:     if (paramid == "scales") {
c@178:         return m_scales;
c@178:     } 
c@178:     else if (paramid == "wavelet") {
c@178:         return int(m_wavelet);
c@178:     }
c@178:     else if (paramid == "threshold") {
c@178:         return m_threshold;
c@178:     } 
c@178:     else if (paramid == "absolute") {
c@178:         return m_absolute;
c@178:     } 
c@178:  
c@178:     return 0.0f; 
c@178: }
c@178:  
c@178:  
c@178: DWT::FeatureSet
c@178: DWT::process(const float *const *inputBuffers,
c@178:              Vamp::RealTime)
c@178: {
c@178:     FeatureSet fs;
c@178: 	
c@178:     if (m_blockSize == 0) {
c@178:         cerr << "ERROR: DWT::process: Not initialised" << endl;
c@178:         return fs;
c@178:     } 
c@178: 	
c@178:     int s = m_scales;
c@178:     int b = m_blockSize;
c@178:     int b_init = b;
c@178: 
c@178:     if ((1 << s) > b) b = 1 << s;												// correct blocksize if smaller than 2^(max scale)
c@178: 
c@178: //--------------------------------------------------------------------------------------------------	
c@178: 
c@178:     float tempDet;
c@178:     float aTempDet;
c@178:     int outloc;
c@178:     int halfblocksize = int(.5 * b);
c@178:     int fbufloc;
c@178:     int fbufloc2;    
c@178: 		
c@178:     vector< vector<float> > wCoefficients(m_scales);						// result
c@178:     vector<float> tempAprx(halfblocksize,0.0);								// approximation
c@178:     vector<float> fbuf(b+m_flength-2,0.0);									// input buffer
c@178: 
c@178:     for (int n=m_flength-2; n<b+m_flength-2; n++)								// copy input buffer to dwt input
c@178:         fbuf[n] = inputBuffers[0][n-m_flength+2];		
c@178: 	
c@178:     for (int scale=0; scale<m_scales; ++scale)								// do for each scale
c@178:     {	
c@178:         for (int n=0; n<m_flength-2; ++n)										// get samples from previous block
c@178:             fbuf[n]  = m_samplePass[scale][n];
c@178: 		
c@178: 		
c@178:         if ((m_flength-2)<b)													// pass samples to next block
c@178:             for (int n=0; n<m_flength-2; ++n)					
c@178:                 m_samplePass[scale][n] = fbuf[b+n];
c@178:         else {
c@178:             for (int n=0; n<b; ++n)										// if number of samples to pass > blocksize
c@178:                 m_samplePass[scale].push_back(fbuf[m_flength-2+n]);
c@178:             m_samplePass[scale].erase (m_samplePass[scale].begin(),m_samplePass[scale].begin()+b);
c@178:         }
c@178: 					
c@178:         for (int n=0; n<halfblocksize; ++n)	{								// do for every other sample of the input buffer
c@178:             tempDet = 0;
c@178:             fbufloc = 2*n+m_flength-1;
c@178:             for (int m=0; m<m_flength; ++m)	{								// Convolve the sample with filter coefficients
c@178:                 fbufloc2 = fbufloc - m;
c@178:                 tempAprx[n] += fbuf[fbufloc2] * m_lpd[m];						// approximation
c@178:                 tempDet += fbuf[fbufloc2] * m_hpd[m];							// detail
c@178:             }
c@178: 		
c@178:             aTempDet = fabs(tempDet);
c@178:             if (m_absolute == 1) tempDet = aTempDet;
c@178: 			
c@178: 			
c@178:             if (aTempDet < m_threshold) tempDet = 0;							// simple hard thresholding, same for each scale
c@178:             wCoefficients[scale].push_back(tempDet);
c@178:         }
c@178: 																				
c@178:         if (scale+1<m_scales) {												// prepare variables for next scale
c@178:             b = b >> 1;														// the approximation in tmpfwd is stored as 
c@178:             halfblocksize = halfblocksize >> 1;								// input for next level
c@178: 	
c@178:             for (int n=m_flength-2; n<b+m_flength-2; n++)						// copy approximation to dwt input
c@178:                 fbuf[n] = tempAprx[n-m_flength+2];
c@178: 			
c@178:             //vector<float>(b+m_flength-2).swap(fbuf);
c@178:             vector<float>(halfblocksize).swap(tempAprx);					// set new size with zeros
c@178:         }
c@178:     }
c@178: 	
c@178: 	
c@178: //-----------------------------------------------------------------------------------------
c@178: 	
c@178:     halfblocksize = int(.5 * b_init);
c@178: 	
c@178:     for (int m = 0; m<halfblocksize; m++) {
c@178: 		
c@178:         Feature feature;
c@178:         feature.hasTimestamp = false;
c@178: 		
c@178:         for (int j = 0; j < s; j++) {
c@178:             outloc = m / (1 << j);									// This one pushes a single result bin 
c@178:             // onto the top of a feature column
c@178:             feature.values.push_back(wCoefficients[j][outloc]);				// each coefficient on higher scales need 
c@178:         }																	// to be copied multiple times to feature columns
c@178:         fs[0].push_back(feature);
c@178:     }
c@178:     return fs;
c@178: }
c@178: 
c@178: 
c@178: 
c@178: DWT::FeatureSet
c@178: DWT::getRemainingFeatures()
c@178: {
c@178:     int s = m_scales;
c@178: 
c@178:     FeatureSet fs;
c@178: 	
c@178: /*	
c@178: 	int b = 1;							
c@178: 	while (b<((m_flength-1) * (1 << s))) {				//set blocksize to tail length
c@178:         b= (b << 1);
c@178: 	}
c@178: 	int b_init = b;
c@178: 	
c@178: */
c@178:     int b = m_blockSize;
c@178:     int b_init = b;
c@178:     int tailIterations = int(((m_flength-1) * (1 << s)) / b) + 1;   // number of iterations for tail 
c@178: 
c@178: 	
c@178:     for(int m=0; m<tailIterations; ++m)
c@178:     {
c@178: 
c@178: 	b = b_init;
c@178: 	
c@178: 	//-------------------------------------------------------------------------------------------	
c@178: 	float tempDet;
c@178:         float aTempDet;
c@178: 	int outloc;
c@178: 	int halfblocksize = int(.5 * b);
c@178: 	int fbufloc;
c@178: 	int fbufloc2;    
c@178:         int len = m_flength;
c@178: 	
c@178: 	vector< vector<float> > wCoefficients(m_scales);						// result
c@178: 	vector<float> tempAprx(halfblocksize,0.0);								// approximation
c@178: 	vector<float> fbuf(b+len-2,0.0);									// input buffer
c@178: 	
c@178: 	//for (int n=len-2; n<b+len-2; n++)								// copy input buffer to dwt input
c@178: 	//	fbuf[n] = 0; //inputBuffers[0][n-len+2];		
c@178: 	
c@178: 	for (int scale=0; scale<m_scales; ++scale)								// do for each scale
c@178: 	{	
c@178:             for (int n=0; n<len-2; ++n)										// get samples from previous block
c@178:                 fbuf[n]  = m_samplePass[scale][n];
c@178: 		
c@178: 		
c@178:             if ((len-2)<b)													// pass samples to next block
c@178:                 for (int n=0; n<len-2; ++n)					
c@178:                     m_samplePass[scale][n] = fbuf[b+n];
c@178:             else {
c@178:                 for (int n=0; n<b; ++n)										// if number of samples to pass > blocksize
c@178:                     m_samplePass[scale].push_back(fbuf[len-2+n]);
c@178:                 m_samplePass[scale].erase (m_samplePass[scale].begin(),m_samplePass[scale].begin()+b);
c@178:             }
c@178: 		
c@178:             for (int n=0; n<halfblocksize; ++n)	{								// do for every other sample of the input buffer
c@178:                 tempDet = 0;
c@178:                 fbufloc = 2*n+len-1;
c@178:                 for (int m=0; m<len; ++m)	{								// Convolve the sample with filter coefficients
c@178:                     fbufloc2 = fbufloc - m;
c@178:                     tempAprx[n] += fbuf[fbufloc2] * m_lpd[m];						// approximation
c@178:                     tempDet += fbuf[fbufloc2] * m_hpd[m];							// detail
c@178:                 }
c@178: 
c@178:                 aTempDet = fabs(tempDet);
c@178:                 if (m_absolute == 1) tempDet = aTempDet;
c@178:                 if (aTempDet < m_threshold) tempDet = 0;							// simple hard thresholding, same for each scale
c@178:                 wCoefficients[scale].push_back(tempDet);
c@178:             }
c@178: 		
c@178: 	    if (scale+1<m_scales) {												// prepare variables for next scale
c@178:                 b = b >> 1;														// the approximation in tmpfwd is stored as 
c@178:                 halfblocksize = halfblocksize >> 1;								// input for next level
c@178: 			
c@178:                 for (int n=len-2; n<b+len-2; n++)						// copy approximation to dwt input
c@178:                     fbuf[n] = tempAprx[n-len+2];
c@178: 			
c@178:                 //vector<float>(b+len-2).swap(fbuf);
c@178:                 vector<float>(halfblocksize).swap(tempAprx);					// set new size with zeros
c@178:             }
c@178: 	
c@178: 	}
c@178: 	
c@178: //-----------------------------------------------------------------------------------------
c@178: 	
c@178: 	halfblocksize = int(.5 * b_init + 0.1);
c@178: 	
c@178: 	for (int m = 0; m<halfblocksize; m++) {
c@178: 		
c@178:             Feature feature;
c@178:             feature.hasTimestamp = false;
c@178: 		
c@178:             for (int j = 0; j < s; j++) {
c@178:                 outloc = m / (1 << j);									// This one pushes a single result bin 
c@178:                 // onto the top of a feature column
c@178:                 feature.values.push_back(wCoefficients[j][outloc]);				// each coefficient on higher scales need 
c@178:             }																	// to be copied multiple times to feature columns
c@178:             fs[0].push_back(feature);
c@178:         }
c@178:     }
c@178:     return fs;
c@178: 
c@178: }
c@178: