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