Chris@7: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@7: 
Chris@7: /*
Chris@7:   Vamp feature extraction plugin for the BeatRoot beat tracker.
Chris@7: 
Chris@7:   Centre for Digital Music, Queen Mary, University of London.
Chris@7:   This file copyright 2011 Simon Dixon, Chris Cannam and QMUL.
Chris@7:     
Chris@7:   This program is free software; you can redistribute it and/or
Chris@7:   modify it under the terms of the GNU General Public License as
Chris@7:   published by the Free Software Foundation; either version 2 of the
Chris@7:   License, or (at your option) any later version.  See the file
Chris@7:   COPYING included with this distribution for more information.
Chris@7: */
Chris@7: 
Chris@7: #ifndef _INDUCTION_H_
Chris@7: #define _INDUCTION_H_
Chris@7: 
Chris@7: #include "Agent.h"
Chris@7: #include "Event.h"
Chris@7: 
Chris@7: #include <vector>
Chris@7: 
Chris@7: using std::vector;
Chris@7: 
Chris@7: /** Performs tempo induction by finding clusters of similar
Chris@7:  *  inter-onset intervals (IOIs), ranking them according to the number
Chris@7:  *  of intervals and relationships between them, and returning a set
Chris@7:  *  of tempo hypotheses for initialising the beat tracking agents.
Chris@7:  */
Chris@7: class Induction
Chris@7: {
Chris@7: public:
Chris@7:     /** The maximum difference in IOIs which are in the same cluster */ 
Chris@7:     static double clusterWidth;
Chris@7:     
Chris@7:     /** The minimum IOI for inclusion in a cluster */
Chris@7:     static double minIOI;
Chris@7: 	
Chris@7:     /** The maximum IOI for inclusion in a cluster */
Chris@7:     static double maxIOI;
Chris@7: 	
Chris@7:     /** The minimum inter-beat interval (IBI), i.e. the maximum tempo
Chris@7:      *  hypothesis that can be returned.
Chris@7:      *  0.30 seconds == 200 BPM
Chris@7:      *  0.25 seconds == 240 BPM
Chris@7:      */
Chris@7:     static double minIBI; 
Chris@7: 
Chris@7:     /** The maximum inter-beat interval (IBI), i.e. the minimum tempo
Chris@7:      *  hypothesis that can be returned.
Chris@7:      *  1.00 seconds ==  60 BPM
Chris@7:      *  0.75 seconds ==  80 BPM
Chris@7:      *  0.60 seconds == 100 BPM
Chris@7:      */
Chris@7:     static double maxIBI;	//  60BPM	// was 0.75 =>  80
Chris@7: 	
Chris@7:     /** The maximum number of tempo hypotheses to return */
Chris@7:     static int topN;
Chris@7: 	
Chris@7:     /** Performs tempo induction (see JNMR 2001 paper by Simon Dixon for details). 
Chris@7:      *  @param events The onsets (or other events) from which the tempo is induced
Chris@7:      *  @return A list of beat tracking agents, where each is initialised with one
Chris@7:      *          of the top tempo hypotheses but no beats
Chris@7:      */
Chris@7:     static AgentList beatInduction(EventList events) {
Chris@7: 	int i, j, b, bestCount;
Chris@7: 	bool submult;
Chris@7: 	int intervals = 0;			// number of interval clusters
Chris@7: 	vector<int> bestn;// count of high-scoring clusters
Chris@7: 	bestn.resize(topN);
Chris@7: 
Chris@7: 	double ratio, err;
Chris@7: 	int degree;
Chris@7: 	int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth);
Chris@7: 	vector<double> clusterMean;
Chris@7: 	clusterMean.resize(maxClusterCount);
Chris@7: 	vector<int> clusterSize;
Chris@7: 	clusterSize.resize(maxClusterCount);
Chris@7: 	vector<int> clusterScore;
Chris@7: 	clusterScore.resize(maxClusterCount);
Chris@7: 		
Chris@7: 	EventList::iterator ptr1, ptr2;
Chris@7: 	Event e1, e2;
Chris@7: 	ptr1 = events.begin();
Chris@7: 	while (ptr1 != events.end()) {
Chris@7: 	    e1 = *ptr1;
Chris@7: 	    ++ptr1;
Chris@7: 	    ptr2 = events.begin();
Chris@7: 	    e2 = *ptr2;
Chris@7: 	    ++ptr2;
Chris@7: 	    while (e2 != e1 && ptr2 != events.end()) {
Chris@7: 		e2 = *ptr2;
Chris@7: 		++ptr2;
Chris@7: 	    }
Chris@7: 	    while (ptr2 != events.end()) {
Chris@7: 		e2 = *ptr2;
Chris@7: 		++ptr2;
Chris@7: 		double ioi = e2.time - e1.time;
Chris@7: 		if (ioi < minIOI)		// skip short intervals
Chris@7: 		    continue;
Chris@7: 		if (ioi > maxIOI)		// ioi too long
Chris@7: 		    break;
Chris@7: 		for (b = 0; b < intervals; b++)		// assign to nearest cluster
Chris@7: 		    if (fabs(clusterMean[b] - ioi) < clusterWidth) {
Chris@7: 			if ((b < intervals - 1) && (
Chris@7: 				fabs(clusterMean[b+1] - ioi) <
Chris@7: 				fabs(clusterMean[b] - ioi)))
Chris@7: 			    b++;		// next cluster is closer
Chris@7: 			clusterMean[b] = (clusterMean[b] * clusterSize[b] +ioi)/
Chris@7: 			    (clusterSize[b] + 1);
Chris@7: 			clusterSize[b]++;
Chris@7: 			break;
Chris@7: 		    }
Chris@7: 		if (b == intervals) {	// no suitable cluster; create new one
Chris@7: 		    if (intervals == maxClusterCount) {
Chris@7: //			System.err.println("Warning: Too many clusters");
Chris@7: 			continue; // ignore this IOI
Chris@7: 		    }
Chris@7: 		    intervals++;
Chris@7: 		    for ( ; (b>0) && (clusterMean[b-1] > ioi); b--) {
Chris@7: 			clusterMean[b] = clusterMean[b-1];
Chris@7: 			clusterSize[b] = clusterSize[b-1];
Chris@7: 		    }
Chris@7: 		    clusterMean[b] = ioi;
Chris@7: 		    clusterSize[b] = 1;
Chris@7: 		}
Chris@7: 	    }
Chris@7: 	}
Chris@7: 	for (b = 0; b < intervals; b++)	// merge similar intervals
Chris@7: 	    // TODO: they are now in order, so don't need the 2nd loop
Chris@7: 	    // TODO: check BOTH sides before averaging or upper gps don't work
Chris@7: 	    for (i = b+1; i < intervals; i++)
Chris@7: 		if (fabs(clusterMean[b] - clusterMean[i]) < clusterWidth) {
Chris@7: 		    clusterMean[b] = (clusterMean[b] * clusterSize[b] +
Chris@7: 				      clusterMean[i] * clusterSize[i]) /
Chris@7: 			(clusterSize[b] + clusterSize[i]);
Chris@7: 		    clusterSize[b] = clusterSize[b] + clusterSize[i];
Chris@7: 		    --intervals;
Chris@7: 		    for (j = i+1; j <= intervals; j++) {
Chris@7: 			clusterMean[j-1] = clusterMean[j];
Chris@7: 			clusterSize[j-1] = clusterSize[j];
Chris@7: 		    }
Chris@7: 		}
Chris@7: 	if (intervals == 0)
Chris@7: 	    return AgentList();
Chris@7: 	for (b = 0; b < intervals; b++)
Chris@7: 	    clusterScore[b] = 10 * clusterSize[b];
Chris@7: 	bestn[0] = 0;
Chris@7: 	bestCount = 1;
Chris@7: 	for (b = 0; b < intervals; b++)
Chris@7: 	    for (i = 0; i <= bestCount; i++)
Chris@7: 		if ((i < topN) && ((i == bestCount) ||
Chris@7: 				   (clusterScore[b] > clusterScore[bestn[i]]))){
Chris@7: 		    if (bestCount < topN)
Chris@7: 			bestCount++;
Chris@7: 		    for (j = bestCount - 1; j > i; j--)
Chris@7: 			bestn[j] = bestn[j-1];
Chris@7: 		    bestn[i] = b;
Chris@7: 		    break;
Chris@7: 		}
Chris@7: 	for (b = 0; b < intervals; b++)	// score intervals
Chris@7: 	    for (i = b+1; i < intervals; i++) {
Chris@7: 		ratio = clusterMean[b] / clusterMean[i];
Chris@7: 		submult = ratio < 1;
Chris@7: 		if (submult)
Chris@7: 		    degree = (int) nearbyint(1/ratio);
Chris@7: 		else
Chris@7: 		    degree = (int) nearbyint(ratio);
Chris@7: 		if ((degree >= 2) && (degree <= 8)) {
Chris@7: 		    if (submult)
Chris@7: 			err = fabs(clusterMean[b]*degree - clusterMean[i]);
Chris@7: 		    else
Chris@7: 			err = fabs(clusterMean[b] - clusterMean[i]*degree);
Chris@7: 		    if (err < (submult? clusterWidth : clusterWidth * degree)) {
Chris@7: 			if (degree >= 5)
Chris@7: 			    degree = 1;
Chris@7: 			else
Chris@7: 			    degree = 6 - degree;
Chris@7: 			clusterScore[b] += degree * clusterSize[i];
Chris@7: 			clusterScore[i] += degree * clusterSize[b];
Chris@7: 		    }
Chris@7: 		}
Chris@7: 	    }
Chris@7: 
Chris@7: 	AgentList a;
Chris@7: 	for (int index = 0; index < bestCount; index++) {
Chris@7: 	    b = bestn[index];
Chris@7: 	    // Adjust it, using the size of super- and sub-intervals
Chris@7: 	    double newSum = clusterMean[b] * clusterScore[b];
Chris@7: 	    int newCount = clusterSize[b];
Chris@7: 	    int newWeight = clusterScore[b];
Chris@7: 	    for (i = 0; i < intervals; i++) {
Chris@7: 		if (i == b)
Chris@7: 		    continue;
Chris@7: 		ratio = clusterMean[b] / clusterMean[i];
Chris@7: 		if (ratio < 1) {
Chris@7: 		    degree = (int) nearbyint(1 / ratio);
Chris@7: 		    if ((degree >= 2) && (degree <= 8)) {
Chris@7: 			err = fabs(clusterMean[b]*degree - clusterMean[i]);
Chris@7: 			if (err < clusterWidth) {
Chris@7: 			    newSum += clusterMean[i] / degree * clusterScore[i];
Chris@7: 			    newCount += clusterSize[i];
Chris@7: 			    newWeight += clusterScore[i];
Chris@7: 			}
Chris@7: 		    }
Chris@7: 		} else {
Chris@7: 		    degree = (int) nearbyint(ratio);
Chris@7: 		    if ((degree >= 2) && (degree <= 8)) {
Chris@7: 			err = fabs(clusterMean[b] - degree*clusterMean[i]);
Chris@7: 			if (err < clusterWidth * degree) {
Chris@7: 			    newSum += clusterMean[i] * degree * clusterScore[i];
Chris@7: 			    newCount += clusterSize[i];
Chris@7: 			    newWeight += clusterScore[i];
Chris@7: 			}
Chris@7: 		    }
Chris@7: 		}
Chris@7: 	    }
Chris@7: 	    double beat = newSum / newWeight;
Chris@7: 	    // Scale within range ... hope the grouping isn't ternary :(
Chris@7: 	    while (beat < minIBI)		// Maximum speed
Chris@7: 		beat *= 2.0;
Chris@7: 	    while (beat > maxIBI)		// Minimum speed
Chris@7: 		beat /= 2.0;
Chris@7: 	    if (beat >= minIBI) {
Chris@7: 		a.push_back(Agent(beat));
Chris@7: 	    }
Chris@7: 	}
Chris@7: 	return a;
Chris@7:     } // beatInduction()
Chris@7: 
Chris@7: protected:
Chris@7:     /** For variable cluster widths in newInduction().
Chris@7:      * @param low The lowest IOI allowed in the cluster
Chris@7:      * @return The highest IOI allowed in the cluster
Chris@7:      */
Chris@7:     static int top(int low) {
Chris@7: 	return low + 25; // low/10;
Chris@7:     } // top()
Chris@7: 
Chris@7: }; // class Induction
Chris@7: 
Chris@7: #endif