# HG changeset patch
# User Chris Cannam
# Date 1318413352 -3600
# Node ID 887c629502a9135f2d0f86a66b78810c622badca
# Parent f1252b6a7cf5a4a94be6a46747d3a75110347ee9
refactor: pull method implementations into .cpp files
diff -r f1252b6a7cf5 -r 887c629502a9 Agent.cpp
--- a/Agent.cpp Fri Oct 07 14:07:42 2011 +0100
+++ b/Agent.cpp Wed Oct 12 10:55:52 2011 +0100
@@ -31,6 +31,30 @@
double Agent::expiryTime = 0.0;
double Agent::decayFactor = 0.0;
+void Agent::accept(Event e, double err, int beats) {
+ beatTime = e.time;
+ events.push_back(e);
+ if (fabs(initialBeatInterval - beatInterval -
+ err / correctionFactor) < MAX_CHANGE * initialBeatInterval)
+ beatInterval += err / correctionFactor;// Adjust tempo
+ beatCount += beats;
+ double conFactor = 1.0 - CONF_FACTOR * err /
+ (err>0? postMargin: -preMargin);
+ if (decayFactor > 0) {
+ double memFactor = 1. - 1. / threshold((double)beatCount,1,decayFactor);
+ phaseScore = memFactor * phaseScore +
+ (1.0 - memFactor) * conFactor * e.salience;
+ } else
+ phaseScore += conFactor * e.salience;
+
+#ifdef DEBUG_BEATROOT
+ std::cerr << "Ag#" << idNumber << ": " << beatInterval << std::endl;
+ std::cerr << " Beat" << beatCount << " Time=" << beatTime
+ << " Score=" << tempoScore << ":P" << phaseScore << ":"
+ << topScoreTime << std::endl;
+#endif
+} // accept()
+
bool Agent::considerAsBeat(Event e, AgentList &a) {
double err;
if (beatTime < 0) { // first event
diff -r f1252b6a7cf5 -r 887c629502a9 Agent.h
--- a/Agent.h Fri Oct 07 14:07:42 2011 +0100
+++ b/Agent.h Wed Oct 12 10:55:52 2011 +0100
@@ -157,28 +157,7 @@
* @param err The difference between the predicted and actual beat times.
* @param beats The number of beats since the last beat that matched an Event.
*/
- void accept(Event e, double err, int beats) {
- beatTime = e.time;
- events.push_back(e);
- if (fabs(initialBeatInterval - beatInterval -
- err / correctionFactor) < MAX_CHANGE * initialBeatInterval)
- beatInterval += err / correctionFactor;// Adjust tempo
- beatCount += beats;
- double conFactor = 1.0 - CONF_FACTOR * err /
- (err>0? postMargin: -preMargin);
- if (decayFactor > 0) {
- double memFactor = 1. - 1. / threshold((double)beatCount,1,decayFactor);
- phaseScore = memFactor * phaseScore +
- (1.0 - memFactor) * conFactor * e.salience;
- } else
- phaseScore += conFactor * e.salience;
-#ifdef DEBUG_BEATROOT
- std::cerr << "Ag#" << idNumber << ": " << beatInterval << std::endl;
- std::cerr << " Beat" << beatCount << " Time=" << beatTime
- << " Score=" << tempoScore << ":P" << phaseScore << ":"
- << topScoreTime << std::endl;
-#endif
- } // accept()
+ void accept(Event e, double err, int beats);
/** The given Event is tested for a possible beat time. The following situations can occur:
* 1) The Agent has no beats yet; the Event is accepted as the first beat.
diff -r f1252b6a7cf5 -r 887c629502a9 AgentList.cpp
--- a/AgentList.cpp Fri Oct 07 14:07:42 2011 +0100
+++ b/AgentList.cpp Wed Oct 12 10:55:52 2011 +0100
@@ -19,3 +19,117 @@
const double AgentList::DEFAULT_BI = 0.02;
const double AgentList::DEFAULT_BT = 0.04;
+void AgentList::removeDuplicates()
+{
+ sort();
+ for (iterator itr = begin(); itr != end(); ++itr) {
+ if (itr->phaseScore < 0.0) // already flagged for deletion
+ continue;
+ iterator itr2 = itr;
+ for (++itr2; itr2 != end(); ++itr2) {
+ if (itr2->beatInterval - itr->beatInterval > DEFAULT_BI)
+ break;
+ if (fabs(itr->beatTime - itr2->beatTime) > DEFAULT_BT)
+ continue;
+ if (itr->phaseScore < itr2->phaseScore) {
+ itr->phaseScore = -1.0; // flag for deletion
+ if (itr2->topScoreTime < itr->topScoreTime)
+ itr2->topScoreTime = itr->topScoreTime;
+ break;
+ } else {
+ itr2->phaseScore = -1.0; // flag for deletion
+ if (itr->topScoreTime < itr2->topScoreTime)
+ itr->topScoreTime = itr2->topScoreTime;
+ }
+ }
+ }
+ int removed = 0;
+ for (iterator itr = begin(); itr != end(); ) {
+ if (itr->phaseScore < 0.0) {
+ ++removed;
+ list.erase(itr);
+ } else {
+ ++itr;
+ }
+ }
+#ifdef DEBUG_BEATROOT
+ if (removed > 0) {
+ std::cerr << "removeDuplicates: removed " << removed << ", have "
+ << list.size() << " agent(s) remaining" << std::endl;
+ }
+ int n = 0;
+ for (Container::iterator i = list.begin(); i != list.end(); ++i) {
+ std::cerr << "agent " << n++ << ": time " << i->beatTime << std::endl;
+ }
+#endif
+} // removeDuplicates()
+
+
+void AgentList::beatTrack(EventList el, double stop)
+{
+ EventList::iterator ei = el.begin();
+ bool phaseGiven = !empty() && (begin()->beatTime >= 0); // if given for one, assume given for others
+ while (ei != el.end()) {
+ Event ev = *ei;
+ ++ei;
+ if ((stop > 0) && (ev.time > stop))
+ break;
+ bool created = phaseGiven;
+ double prevBeatInterval = -1.0;
+ // cc: Duplicate our list of agents, and scan through the
+ // copy. This means we can safely add agents to our own
+ // list while scanning without disrupting our scan. Each
+ // agent needs to be re-added to our own list explicitly
+ // (since it is modified by e.g. considerAsBeat)
+ Container currentAgents = list;
+ list.clear();
+ for (Container::iterator ai = currentAgents.begin();
+ ai != currentAgents.end(); ++ai) {
+ Agent currentAgent = *ai;
+ if (currentAgent.beatInterval != prevBeatInterval) {
+ if ((prevBeatInterval>=0) && !created && (ev.time<5.0)) {
+#ifdef DEBUG_BEATROOT
+ std::cerr << "Creating a new agent" << std::endl;
+#endif
+ // Create new agent with different phase
+ Agent newAgent(prevBeatInterval);
+ // This may add another agent to our list as well
+ newAgent.considerAsBeat(ev, *this);
+ add(newAgent);
+ }
+ prevBeatInterval = currentAgent.beatInterval;
+ created = phaseGiven;
+ }
+ if (currentAgent.considerAsBeat(ev, *this))
+ created = true;
+ add(currentAgent);
+ } // loop for each agent
+ removeDuplicates();
+ } // loop for each event
+} // beatTrack()
+
+Agent *AgentList::bestAgent()
+{
+ double best = -1.0;
+ Agent *bestAg = 0;
+ for (iterator itr = begin(); itr != end(); ++itr) {
+ if (itr->events.empty()) continue;
+ double startTime = itr->events.begin()->time;
+ double conf = (itr->phaseScore + itr->tempoScore) /
+ (useAverageSalience? (double)itr->beatCount: 1.0);
+ if (conf > best) {
+ bestAg = &(*itr);
+ best = conf;
+ }
+ }
+#ifdef DEBUG_BEATROOT
+ if (bestAg) {
+ std::cerr << "Best agent: Ag#" << bestAg->idNumber << std::endl;
+ std::cerr << " Av-salience = " << best << std::endl;
+ } else {
+ std::cerr << "No surviving agent - beat tracking failed" << std::endl;
+ }
+#endif
+ return bestAg;
+} // bestAgent()
+
diff -r f1252b6a7cf5 -r 887c629502a9 AgentList.h
--- a/AgentList.h Fri Oct 07 14:07:42 2011 +0100
+++ b/AgentList.h Wed Oct 12 10:55:52 2011 +0100
@@ -93,49 +93,7 @@
* A duplicate is defined by the tempo and phase thresholds
* thresholdBI and thresholdBT respectively.
*/
- void removeDuplicates() {
- sort();
- for (iterator itr = begin(); itr != end(); ++itr) {
- if (itr->phaseScore < 0.0) // already flagged for deletion
- continue;
- iterator itr2 = itr;
- for (++itr2; itr2 != end(); ++itr2) {
- if (itr2->beatInterval - itr->beatInterval > DEFAULT_BI)
- break;
- if (fabs(itr->beatTime - itr2->beatTime) > DEFAULT_BT)
- continue;
- if (itr->phaseScore < itr2->phaseScore) {
- itr->phaseScore = -1.0; // flag for deletion
- if (itr2->topScoreTime < itr->topScoreTime)
- itr2->topScoreTime = itr->topScoreTime;
- break;
- } else {
- itr2->phaseScore = -1.0; // flag for deletion
- if (itr->topScoreTime < itr2->topScoreTime)
- itr->topScoreTime = itr2->topScoreTime;
- }
- }
- }
- int removed = 0;
- for (iterator itr = begin(); itr != end(); ) {
- if (itr->phaseScore < 0.0) {
- ++removed;
- list.erase(itr);
- } else {
- ++itr;
- }
- }
-#ifdef DEBUG_BEATROOT
- if (removed > 0) {
- std::cerr << "removeDuplicates: removed " << removed << ", have "
- << list.size() << " agent(s) remaining" << std::endl;
- }
- int n = 0;
- for (Container::iterator i = list.begin(); i != list.end(); ++i) {
- std::cerr << "agent " << n++ << ": time " << i->beatTime << std::endl;
- }
-#endif
- } // removeDuplicates()
+ void removeDuplicates();
public:
/** Perform beat tracking on a list of events (onsets).
@@ -149,75 +107,12 @@
* @param el The list of onsets (or events or peaks) to beat track.
* @param stop Do not find beats after stop seconds.
*/
- void beatTrack(EventList el, double stop) {
- EventList::iterator ei = el.begin();
- bool phaseGiven = !empty() && (begin()->beatTime >= 0); // if given for one, assume given for others
- while (ei != el.end()) {
- Event ev = *ei;
- ++ei;
- if ((stop > 0) && (ev.time > stop))
- break;
- bool created = phaseGiven;
- double prevBeatInterval = -1.0;
- // cc: Duplicate our list of agents, and scan through the
- // copy. This means we can safely add agents to our own
- // list while scanning without disrupting our scan. Each
- // agent needs to be re-added to our own list explicitly
- // (since it is modified by e.g. considerAsBeat)
- Container currentAgents = list;
- list.clear();
- for (Container::iterator ai = currentAgents.begin();
- ai != currentAgents.end(); ++ai) {
- Agent currentAgent = *ai;
- if (currentAgent.beatInterval != prevBeatInterval) {
- if ((prevBeatInterval>=0) && !created && (ev.time<5.0)) {
-#ifdef DEBUG_BEATROOT
- std::cerr << "Creating a new agent" << std::endl;
-#endif
- // Create new agent with different phase
- Agent newAgent(prevBeatInterval);
- // This may add another agent to our list as well
- newAgent.considerAsBeat(ev, *this);
- add(newAgent);
- }
- prevBeatInterval = currentAgent.beatInterval;
- created = phaseGiven;
- }
- if (currentAgent.considerAsBeat(ev, *this))
- created = true;
- add(currentAgent);
- } // loop for each agent
- removeDuplicates();
- } // loop for each event
- } // beatTrack()
+ void beatTrack(EventList el, double stop);
/** Finds the Agent with the highest score in the list, or NULL if beat tracking has failed.
* @return The Agent with the highest score
*/
- Agent *bestAgent() {
- double best = -1.0;
- Agent *bestAg = 0;
- for (iterator itr = begin(); itr != end(); ++itr) {
- if (itr->events.empty()) continue;
- double startTime = itr->events.begin()->time;
- double conf = (itr->phaseScore + itr->tempoScore) /
- (useAverageSalience? (double)itr->beatCount: 1.0);
- if (conf > best) {
- bestAg = &(*itr);
- best = conf;
- }
- }
-#ifdef DEBUG_BEATROOT
- if (bestAg) {
- std::cerr << "Best agent: Ag#" << bestAg->idNumber << std::endl;
- std::cerr << " Av-salience = " << best << std::endl;
- } else {
- std::cerr << "No surviving agent - beat tracking failed" << std::endl;
- }
-#endif
- return bestAg;
- } // bestAgent()
-
+ Agent *bestAgent();
}; // class AgentList
diff -r f1252b6a7cf5 -r 887c629502a9 BeatRootProcessor.cpp
--- a/BeatRootProcessor.cpp Fri Oct 07 14:07:42 2011 +0100
+++ b/BeatRootProcessor.cpp Wed Oct 12 10:55:52 2011 +0100
@@ -18,3 +18,56 @@
bool
BeatRootProcessor::silent = true;
+void BeatRootProcessor::processFrame(const float *const *inputBuffers) {
+ double flux = 0;
+ for (int i = 0; i <= fftSize/2; i++) {
+ double mag = sqrt(inputBuffers[0][i*2] * inputBuffers[0][i*2] +
+ inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1]);
+ if (mag > prevFrame[i]) flux += mag - prevFrame[i];
+ prevFrame[i] = mag;
+ }
+
+ spectralFlux.push_back(flux);
+
+} // processFrame()
+
+EventList BeatRootProcessor::beatTrack() {
+
+#ifdef DEBUG_BEATROOT
+ std::cerr << "Spectral flux:" << std::endl;
+ for (int i = 0; i < spectralFlux.size(); ++i) {
+ if ((i % 8) == 0) std::cerr << "\n";
+ std::cerr << spectralFlux[i] << " ";
+ }
+#endif
+
+ double hop = hopTime;
+ Peaks::normalise(spectralFlux);
+ vector peaks = Peaks::findPeaks(spectralFlux, (int)lrint(0.06 / hop), 0.35, 0.84, true);
+ onsets.clear();
+ onsets.resize(peaks.size(), 0);
+ vector::iterator it = peaks.begin();
+ onsetList.clear();
+ double minSalience = Peaks::min(spectralFlux);
+ for (int i = 0; i < onsets.size(); i++) {
+ int index = *it;
+ ++it;
+ onsets[i] = index * hop;
+ Event e = BeatTracker::newBeat(onsets[i], 0);
+// if (debug)
+// System.err.printf("Onset: %8.3f %8.3f %8.3f\n",
+// onsets[i], energy[index], slope[index]);
+// e.salience = slope[index]; // or combination of energy + slope??
+ // Note that salience must be non-negative or the beat tracking system fails!
+ e.salience = spectralFlux[index] - minSalience;
+ onsetList.push_back(e);
+ }
+
+#ifdef DEBUG_BEATROOT
+ std::cerr << "Onsets: " << onsetList.size() << std::endl;
+#endif
+
+ return BeatTracker::beatTrack(onsetList);
+
+} // processFile()
+
diff -r f1252b6a7cf5 -r 887c629502a9 BeatRootProcessor.h
--- a/BeatRootProcessor.h Fri Oct 07 14:07:42 2011 +0100
+++ b/BeatRootProcessor.h Wed Oct 12 10:55:52 2011 +0100
@@ -108,57 +108,11 @@
* then computing the spectral flux then (optionally) normalising
* and calculating onsets.
*/
- void processFrame(const float *const *inputBuffers) {
- double flux = 0;
- for (int i = 0; i <= fftSize/2; i++) {
- double mag = sqrt(inputBuffers[0][i*2] * inputBuffers[0][i*2] +
- inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1]);
- if (mag > prevFrame[i]) flux += mag - prevFrame[i];
- prevFrame[i] = mag;
- }
-
- spectralFlux.push_back(flux);
-
- } // processFrame()
+ void processFrame(const float *const *inputBuffers);
/** Tracks beats once all frames have been processed by processFrame
*/
- EventList beatTrack() {
-
- for (int i = 0; i < spectralFlux.size(); ++i) {
- if ((i % 8) == 0) std::cerr << "\n";
- std::cerr << spectralFlux[i] << " ";
- }
-
- double hop = hopTime;
- Peaks::normalise(spectralFlux);
- vector peaks = Peaks::findPeaks(spectralFlux, (int)lrint(0.06 / hop), 0.35, 0.84, true);
- onsets.clear();
- onsets.resize(peaks.size(), 0);
- vector::iterator it = peaks.begin();
- onsetList.clear();
- double minSalience = Peaks::min(spectralFlux);
- for (int i = 0; i < onsets.size(); i++) {
- int index = *it;
- ++it;
- onsets[i] = index * hop;
- Event e = BeatTracker::newBeat(onsets[i], 0);
-// if (debug)
-// System.err.printf("Onset: %8.3f %8.3f %8.3f\n",
-// onsets[i], energy[index], slope[index]);
-// e.salience = slope[index]; // or combination of energy + slope??
- // Note that salience must be non-negative or the beat tracking system fails!
- e.salience = spectralFlux[index] - minSalience;
- onsetList.push_back(e);
- }
-
-#ifdef DEBUG_BEATROOT
- std::cerr << "Onsets: " << onsetList.size() << std::endl;
-#endif
-
- return BeatTracker::beatTrack(onsetList);
-
- } // processFile()
+ EventList beatTrack();
protected:
/** Allocates memory for arrays, based on parameter settings */
diff -r f1252b6a7cf5 -r 887c629502a9 BeatTracker.h
--- a/BeatTracker.h Fri Oct 07 14:07:42 2011 +0100
+++ b/BeatTracker.h Wed Oct 12 10:55:52 2011 +0100
@@ -65,64 +65,7 @@
* @param beats The initial beats which are given, if any
* @return The list of beats, or an empty list if beat tracking fails
*/
- static EventList beatTrack(EventList events, EventList beats) {
- AgentList agents;
- int count = 0;
- double beatTime = -1;
- if (!beats.empty()) {
- count = beats.size() - 1;
- EventList::iterator itr = beats.end();
- --itr;
- beatTime = itr->time;
- }
- if (count > 0) { // tempo given by mean of initial beats
- double ioi = (beatTime - beats.begin()->time) / count;
- agents.push_back(Agent(ioi));
- } else // tempo not given; use tempo induction
- agents = Induction::beatInduction(events);
- if (!beats.empty())
- for (AgentList::iterator itr = agents.begin(); itr != agents.end();
- ++itr) {
- itr->beatTime = beatTime;
- itr->beatCount = count;
- itr->events = beats;
- }
- agents.beatTrack(events, -1);
- Agent *best = agents.bestAgent();
- if (best) {
- best->fillBeats(beatTime);
- return best->events;
- }
- return EventList();
- } // beatTrack()/1
-
- /** Finds the mean tempo (as inter-beat interval) from an array of beat times
- * @param d An array of beat times
- * @return The average inter-beat interval
- */
- static double getAverageIBI(vector d) {
- if (d.size() < 2)
- return -1.0;
- return (d[d.size() - 1] - d[0]) / (d.size() - 1);
- } // getAverageIBI()
-
- /** Finds the median tempo (as inter-beat interval) from an array of beat times
- * @param d An array of beat times
- * @return The median inter-beat interval
- */
- static double getMedianIBI(vector d) {
- if (d.size() < 2)
- return -1.0;
- vector ibi;
- ibi.resize(d.size()-1);
- for (int i = 1; i < d.size(); i++)
- ibi[i-1] = d[i] - d[i-1];
- std::sort(ibi.begin(), ibi.end());
- if (ibi.size() % 2 == 0)
- return (ibi[ibi.size() / 2] + ibi[ibi.size() / 2 - 1]) / 2;
- else
- return ibi[ibi.size() / 2];
- } // getAverageIBI()
+ static EventList beatTrack(EventList events, EventList beats);
// Various get and set methods
diff -r f1252b6a7cf5 -r 887c629502a9 Induction.cpp
--- a/Induction.cpp Fri Oct 07 14:07:42 2011 +0100
+++ b/Induction.cpp Wed Oct 12 10:55:52 2011 +0100
@@ -22,3 +22,173 @@
double Induction::maxIBI = 1.0;
int Induction::topN = 10;
+
+AgentList Induction::beatInduction(EventList events) {
+ int i, j, b, bestCount;
+ bool submult;
+ int intervals = 0; // number of interval clusters
+ vector bestn;// count of high-scoring clusters
+ bestn.resize(topN);
+
+ double ratio, err;
+ int degree;
+ int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth);
+ vector clusterMean;
+ clusterMean.resize(maxClusterCount);
+ vector clusterSize;
+ clusterSize.resize(maxClusterCount);
+ vector clusterScore;
+ clusterScore.resize(maxClusterCount);
+
+ EventList::iterator ptr1, ptr2;
+ Event e1, e2;
+ ptr1 = events.begin();
+ while (ptr1 != events.end()) {
+ e1 = *ptr1;
+ ++ptr1;
+ ptr2 = events.begin();
+ e2 = *ptr2;
+ ++ptr2;
+ while (e2 != e1 && ptr2 != events.end()) {
+ e2 = *ptr2;
+ ++ptr2;
+ }
+ while (ptr2 != events.end()) {
+ e2 = *ptr2;
+ ++ptr2;
+ double ioi = e2.time - e1.time;
+ if (ioi < minIOI) // skip short intervals
+ continue;
+ if (ioi > maxIOI) // ioi too long
+ break;
+ for (b = 0; b < intervals; b++) // assign to nearest cluster
+ if (fabs(clusterMean[b] - ioi) < clusterWidth) {
+ if ((b < intervals - 1) && (
+ fabs(clusterMean[b+1] - ioi) <
+ fabs(clusterMean[b] - ioi)))
+ b++; // next cluster is closer
+ clusterMean[b] = (clusterMean[b] * clusterSize[b] +ioi)/
+ (clusterSize[b] + 1);
+ clusterSize[b]++;
+ break;
+ }
+ if (b == intervals) { // no suitable cluster; create new one
+ if (intervals == maxClusterCount) {
+// System.err.println("Warning: Too many clusters");
+ continue; // ignore this IOI
+ }
+ intervals++;
+ for ( ; (b>0) && (clusterMean[b-1] > ioi); b--) {
+ clusterMean[b] = clusterMean[b-1];
+ clusterSize[b] = clusterSize[b-1];
+ }
+ clusterMean[b] = ioi;
+ clusterSize[b] = 1;
+ }
+ }
+ }
+ for (b = 0; b < intervals; b++) // merge similar intervals
+ // TODO: they are now in order, so don't need the 2nd loop
+ // TODO: check BOTH sides before averaging or upper gps don't work
+ for (i = b+1; i < intervals; i++)
+ if (fabs(clusterMean[b] - clusterMean[i]) < clusterWidth) {
+ clusterMean[b] = (clusterMean[b] * clusterSize[b] +
+ clusterMean[i] * clusterSize[i]) /
+ (clusterSize[b] + clusterSize[i]);
+ clusterSize[b] = clusterSize[b] + clusterSize[i];
+ --intervals;
+ for (j = i+1; j <= intervals; j++) {
+ clusterMean[j-1] = clusterMean[j];
+ clusterSize[j-1] = clusterSize[j];
+ }
+ }
+ if (intervals == 0)
+ return AgentList();
+ for (b = 0; b < intervals; b++)
+ clusterScore[b] = 10 * clusterSize[b];
+ bestn[0] = 0;
+ bestCount = 1;
+ for (b = 0; b < intervals; b++)
+ for (i = 0; i <= bestCount; i++)
+ if ((i < topN) && ((i == bestCount) ||
+ (clusterScore[b] > clusterScore[bestn[i]]))){
+ if (bestCount < topN)
+ bestCount++;
+ for (j = bestCount - 1; j > i; j--)
+ bestn[j] = bestn[j-1];
+ bestn[i] = b;
+ break;
+ }
+ for (b = 0; b < intervals; b++) // score intervals
+ for (i = b+1; i < intervals; i++) {
+ ratio = clusterMean[b] / clusterMean[i];
+ submult = ratio < 1;
+ if (submult)
+ degree = (int) nearbyint(1/ratio);
+ else
+ degree = (int) nearbyint(ratio);
+ if ((degree >= 2) && (degree <= 8)) {
+ if (submult)
+ err = fabs(clusterMean[b]*degree - clusterMean[i]);
+ else
+ err = fabs(clusterMean[b] - clusterMean[i]*degree);
+ if (err < (submult? clusterWidth : clusterWidth * degree)) {
+ if (degree >= 5)
+ degree = 1;
+ else
+ degree = 6 - degree;
+ clusterScore[b] += degree * clusterSize[i];
+ clusterScore[i] += degree * clusterSize[b];
+ }
+ }
+ }
+
+ AgentList a;
+ for (int index = 0; index < bestCount; index++) {
+ b = bestn[index];
+ // Adjust it, using the size of super- and sub-intervals
+ double newSum = clusterMean[b] * clusterScore[b];
+ int newCount = clusterSize[b];
+ int newWeight = clusterScore[b];
+ for (i = 0; i < intervals; i++) {
+ if (i == b)
+ continue;
+ ratio = clusterMean[b] / clusterMean[i];
+ if (ratio < 1) {
+ degree = (int) nearbyint(1 / ratio);
+ if ((degree >= 2) && (degree <= 8)) {
+ err = fabs(clusterMean[b]*degree - clusterMean[i]);
+ if (err < clusterWidth) {
+ newSum += clusterMean[i] / degree * clusterScore[i];
+ newCount += clusterSize[i];
+ newWeight += clusterScore[i];
+ }
+ }
+ } else {
+ degree = (int) nearbyint(ratio);
+ if ((degree >= 2) && (degree <= 8)) {
+ err = fabs(clusterMean[b] - degree*clusterMean[i]);
+ if (err < clusterWidth * degree) {
+ newSum += clusterMean[i] * degree * clusterScore[i];
+ newCount += clusterSize[i];
+ newWeight += clusterScore[i];
+ }
+ }
+ }
+ }
+ double beat = newSum / newWeight;
+ // Scale within range ... hope the grouping isn't ternary :(
+ while (beat < minIBI) // Maximum speed
+ beat *= 2.0;
+ while (beat > maxIBI) // Minimum speed
+ beat /= 2.0;
+ if (beat >= minIBI) {
+ a.push_back(Agent(beat));
+ }
+ }
+#ifdef DEBUG_BEATROOT
+ std::cerr << "Induction complete, returning " << a.size() << " agent(s)" << std::endl;
+#endif
+ return a;
+} // beatInduction()
+
diff -r f1252b6a7cf5 -r 887c629502a9 Induction.h
--- a/Induction.h Fri Oct 07 14:07:42 2011 +0100
+++ b/Induction.h Wed Oct 12 10:55:52 2011 +0100
@@ -68,174 +68,7 @@
* @return A list of beat tracking agents, where each is initialised with one
* of the top tempo hypotheses but no beats
*/
- static AgentList beatInduction(EventList events) {
- int i, j, b, bestCount;
- bool submult;
- int intervals = 0; // number of interval clusters
- vector bestn;// count of high-scoring clusters
- bestn.resize(topN);
-
- double ratio, err;
- int degree;
- int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth);
- vector clusterMean;
- clusterMean.resize(maxClusterCount);
- vector clusterSize;
- clusterSize.resize(maxClusterCount);
- vector clusterScore;
- clusterScore.resize(maxClusterCount);
-
- EventList::iterator ptr1, ptr2;
- Event e1, e2;
- ptr1 = events.begin();
- while (ptr1 != events.end()) {
- e1 = *ptr1;
- ++ptr1;
- ptr2 = events.begin();
- e2 = *ptr2;
- ++ptr2;
- while (e2 != e1 && ptr2 != events.end()) {
- e2 = *ptr2;
- ++ptr2;
- }
- while (ptr2 != events.end()) {
- e2 = *ptr2;
- ++ptr2;
- double ioi = e2.time - e1.time;
- if (ioi < minIOI) // skip short intervals
- continue;
- if (ioi > maxIOI) // ioi too long
- break;
- for (b = 0; b < intervals; b++) // assign to nearest cluster
- if (fabs(clusterMean[b] - ioi) < clusterWidth) {
- if ((b < intervals - 1) && (
- fabs(clusterMean[b+1] - ioi) <
- fabs(clusterMean[b] - ioi)))
- b++; // next cluster is closer
- clusterMean[b] = (clusterMean[b] * clusterSize[b] +ioi)/
- (clusterSize[b] + 1);
- clusterSize[b]++;
- break;
- }
- if (b == intervals) { // no suitable cluster; create new one
- if (intervals == maxClusterCount) {
-// System.err.println("Warning: Too many clusters");
- continue; // ignore this IOI
- }
- intervals++;
- for ( ; (b>0) && (clusterMean[b-1] > ioi); b--) {
- clusterMean[b] = clusterMean[b-1];
- clusterSize[b] = clusterSize[b-1];
- }
- clusterMean[b] = ioi;
- clusterSize[b] = 1;
- }
- }
- }
- for (b = 0; b < intervals; b++) // merge similar intervals
- // TODO: they are now in order, so don't need the 2nd loop
- // TODO: check BOTH sides before averaging or upper gps don't work
- for (i = b+1; i < intervals; i++)
- if (fabs(clusterMean[b] - clusterMean[i]) < clusterWidth) {
- clusterMean[b] = (clusterMean[b] * clusterSize[b] +
- clusterMean[i] * clusterSize[i]) /
- (clusterSize[b] + clusterSize[i]);
- clusterSize[b] = clusterSize[b] + clusterSize[i];
- --intervals;
- for (j = i+1; j <= intervals; j++) {
- clusterMean[j-1] = clusterMean[j];
- clusterSize[j-1] = clusterSize[j];
- }
- }
- if (intervals == 0)
- return AgentList();
- for (b = 0; b < intervals; b++)
- clusterScore[b] = 10 * clusterSize[b];
- bestn[0] = 0;
- bestCount = 1;
- for (b = 0; b < intervals; b++)
- for (i = 0; i <= bestCount; i++)
- if ((i < topN) && ((i == bestCount) ||
- (clusterScore[b] > clusterScore[bestn[i]]))){
- if (bestCount < topN)
- bestCount++;
- for (j = bestCount - 1; j > i; j--)
- bestn[j] = bestn[j-1];
- bestn[i] = b;
- break;
- }
- for (b = 0; b < intervals; b++) // score intervals
- for (i = b+1; i < intervals; i++) {
- ratio = clusterMean[b] / clusterMean[i];
- submult = ratio < 1;
- if (submult)
- degree = (int) nearbyint(1/ratio);
- else
- degree = (int) nearbyint(ratio);
- if ((degree >= 2) && (degree <= 8)) {
- if (submult)
- err = fabs(clusterMean[b]*degree - clusterMean[i]);
- else
- err = fabs(clusterMean[b] - clusterMean[i]*degree);
- if (err < (submult? clusterWidth : clusterWidth * degree)) {
- if (degree >= 5)
- degree = 1;
- else
- degree = 6 - degree;
- clusterScore[b] += degree * clusterSize[i];
- clusterScore[i] += degree * clusterSize[b];
- }
- }
- }
-
- AgentList a;
- for (int index = 0; index < bestCount; index++) {
- b = bestn[index];
- // Adjust it, using the size of super- and sub-intervals
- double newSum = clusterMean[b] * clusterScore[b];
- int newCount = clusterSize[b];
- int newWeight = clusterScore[b];
- for (i = 0; i < intervals; i++) {
- if (i == b)
- continue;
- ratio = clusterMean[b] / clusterMean[i];
- if (ratio < 1) {
- degree = (int) nearbyint(1 / ratio);
- if ((degree >= 2) && (degree <= 8)) {
- err = fabs(clusterMean[b]*degree - clusterMean[i]);
- if (err < clusterWidth) {
- newSum += clusterMean[i] / degree * clusterScore[i];
- newCount += clusterSize[i];
- newWeight += clusterScore[i];
- }
- }
- } else {
- degree = (int) nearbyint(ratio);
- if ((degree >= 2) && (degree <= 8)) {
- err = fabs(clusterMean[b] - degree*clusterMean[i]);
- if (err < clusterWidth * degree) {
- newSum += clusterMean[i] * degree * clusterScore[i];
- newCount += clusterSize[i];
- newWeight += clusterScore[i];
- }
- }
- }
- }
- double beat = newSum / newWeight;
- // Scale within range ... hope the grouping isn't ternary :(
- while (beat < minIBI) // Maximum speed
- beat *= 2.0;
- while (beat > maxIBI) // Minimum speed
- beat /= 2.0;
- if (beat >= minIBI) {
- a.push_back(Agent(beat));
- }
- }
-#ifdef DEBUG_BEATROOT
- std::cerr << "Induction complete, returning " << a.size() << " agent(s)" << std::endl;
-#endif
- return a;
- } // beatInduction()
+ static AgentList beatInduction(EventList events);
protected:
/** For variable cluster widths in newInduction().
diff -r f1252b6a7cf5 -r 887c629502a9 Makefile
--- a/Makefile Fri Oct 07 14:07:42 2011 +0100
+++ b/Makefile Wed Oct 12 10:55:52 2011 +0100
@@ -1,7 +1,8 @@
-CXXFLAGS := -g -DDEBUG_BEATROOT
+#CXXFLAGS := -fPIC -g -DDEBUG_BEATROOT
+CXXFLAGS := -fPIC -O3
-OBJECTS := BeatRootProcessor.o BeatRootVampPlugin.o Peaks.o Agent.o AgentList.o Induction.o
+OBJECTS := BeatRootProcessor.o BeatRootVampPlugin.o Peaks.o Agent.o AgentList.o Induction.o BeatTracker.o
HEADERS := Agent.h AgentList.h BeatRootProcessor.h BeatRootVampPlugin.h BeatTracker.h Event.h Induction.h Peaks.h
@@ -19,5 +20,12 @@
BeatRootProcessor.o: Agent.h AgentList.h Induction.h
BeatRootVampPlugin.o: BeatRootVampPlugin.h BeatRootProcessor.h Peaks.h
BeatRootVampPlugin.o: Event.h BeatTracker.h Agent.h AgentList.h Induction.h
+BeatTracker.o: BeatTracker.h Event.h Agent.h AgentList.h Induction.h
Induction.o: Induction.h Agent.h Event.h AgentList.h
Peaks.o: Peaks.h
+Agent.o: Event.h
+AgentList.o: Agent.h Event.h
+BeatRootProcessor.o: Peaks.h Event.h BeatTracker.h Agent.h AgentList.h
+BeatRootProcessor.o: Induction.h
+BeatTracker.o: Event.h Agent.h AgentList.h Induction.h
+Induction.o: Agent.h Event.h AgentList.h
diff -r f1252b6a7cf5 -r 887c629502a9 Peaks.cpp
--- a/Peaks.cpp Fri Oct 07 14:07:42 2011 +0100
+++ b/Peaks.cpp Wed Oct 12 10:55:52 2011 +0100
@@ -18,3 +18,163 @@
int Peaks::pre = 3;
int Peaks::post = 1;
+
+int Peaks::findPeaks(const vector &data, vector peaks, int width) {
+ int peakCount = 0;
+ int maxp = 0;
+ int mid = 0;
+ int end = data.size();
+ while (mid < end) {
+ int i = mid - width;
+ if (i < 0)
+ i = 0;
+ int stop = mid + width + 1;
+ if (stop > data.size())
+ stop = data.size();
+ maxp = i;
+ for (i++; i < stop; i++)
+ if (data[i] > data[maxp])
+ maxp = i;
+ if (maxp == mid) {
+ int j;
+ for (j = peakCount; j > 0; j--) {
+ if (data[maxp] <= data[peaks[j-1]])
+ break;
+ else if (j < peaks.size())
+ peaks[j] = peaks[j-1];
+ }
+ if (j != peaks.size())
+ peaks[j] = maxp;
+ if (peakCount != peaks.size())
+ peakCount++;
+ }
+ mid++;
+ }
+ return peakCount;
+} // findPeaks()
+
+vector Peaks::findPeaks(const vector &data, int width,
+ double threshold, double decayRate, bool isRelative) {
+ vector peaks;
+ int maxp = 0;
+ int mid = 0;
+ int end = data.size();
+ double av = data[0];
+ while (mid < end) {
+ av = decayRate * av + (1 - decayRate) * data[mid];
+ if (av < data[mid])
+ av = data[mid];
+ int i = mid - width;
+ if (i < 0)
+ i = 0;
+ int stop = mid + width + 1;
+ if (stop > data.size())
+ stop = data.size();
+ maxp = i;
+ for (i++; i < stop; i++)
+ if (data[i] > data[maxp])
+ maxp = i;
+ if (maxp == mid) {
+ if (overThreshold(data, maxp, width, threshold, isRelative,av)){
+ peaks.push_back(maxp);
+ }
+ }
+ mid++;
+ }
+ return peaks;
+} // findPeaks()
+
+double Peaks::expDecayWithHold(double av, double decayRate,
+ const vector &data, int start, int stop) {
+ while (start < stop) {
+ av = decayRate * av + (1 - decayRate) * data[start];
+ if (av < data[start])
+ av = data[start];
+ start++;
+ }
+ return av;
+} // expDecayWithHold()
+
+bool Peaks::overThreshold(const vector &data, int index, int width,
+ double threshold, bool isRelative,
+ double av) {
+ if (data[index] < av)
+ return false;
+ if (isRelative) {
+ int iStart = index - pre * width;
+ if (iStart < 0)
+ iStart = 0;
+ int iStop = index + post * width;
+ if (iStop > data.size())
+ iStop = data.size();
+ double sum = 0;
+ int count = iStop - iStart;
+ while (iStart < iStop)
+ sum += data[iStart++];
+ return (data[index] > sum / count + threshold);
+ } else
+ return (data[index] > threshold);
+} // overThreshold()
+
+void Peaks::normalise(vector &data) {
+ double sx = 0;
+ double sxx = 0;
+ for (int i = 0; i < data.size(); i++) {
+ sx += data[i];
+ sxx += data[i] * data[i];
+ }
+ double mean = sx / data.size();
+ double sd = sqrt((sxx - sx * mean) / data.size());
+ if (sd == 0)
+ sd = 1; // all data[i] == mean -> 0; avoids div by 0
+ for (int i = 0; i < data.size(); i++) {
+ data[i] = (data[i] - mean) / sd;
+ }
+} // normalise()
+
+/** Uses an n-point linear regression to estimate the slope of data.
+ * @param data input data
+ * @param hop spacing of data points
+ * @param n length of linear regression
+ * @param slope output data
+ */
+void Peaks::getSlope(const vector &data, double hop, int n,
+ vector &slope) {
+ int i = 0, j = 0;
+ double t;
+ double sx = 0, sxx = 0, sy = 0, sxy = 0;
+ for ( ; i < n; i++) {
+ t = i * hop;
+ sx += t;
+ sxx += t * t;
+ sy += data[i];
+ sxy += t * data[i];
+ }
+ double delta = n * sxx - sx * sx;
+ for ( ; j < n / 2; j++)
+ slope[j] = (n * sxy - sx * sy) / delta;
+ for ( ; j < data.size() - (n + 1) / 2; j++, i++) {
+ slope[j] = (n * sxy - sx * sy) / delta;
+ sy += data[i] - data[i - n];
+ sxy += hop * (n * data[i] - sy);
+ }
+ for ( ; j < data.size(); j++)
+ slope[j] = (n * sxy - sx * sy) / delta;
+} // getSlope()
+
+int Peaks::imin(const vector &arr) {
+ int i = 0;
+ for (int j = 1; j < arr.size(); j++)
+ if (arr[j] < arr[i])
+ i = j;
+ return i;
+} // imin()
+
+int Peaks::imax(const vector &arr) {
+ int i = 0;
+ for (int j = 1; j < arr.size(); j++)
+ if (arr[j] > arr[i])
+ i = j;
+ return i;
+} // imax()
+
diff -r f1252b6a7cf5 -r 887c629502a9 Peaks.h
--- a/Peaks.h Fri Oct 07 14:07:42 2011 +0100
+++ b/Peaks.h Wed Oct 12 10:55:52 2011 +0100
@@ -33,39 +33,7 @@
* @param peaks list of peak indexes
* @param width minimum distance between peaks
*/
- static int findPeaks(const vector &data, vector peaks, int width) {
- int peakCount = 0;
- int maxp = 0;
- int mid = 0;
- int end = data.size();
- while (mid < end) {
- int i = mid - width;
- if (i < 0)
- i = 0;
- int stop = mid + width + 1;
- if (stop > data.size())
- stop = data.size();
- maxp = i;
- for (i++; i < stop; i++)
- if (data[i] > data[maxp])
- maxp = i;
- if (maxp == mid) {
- int j;
- for (j = peakCount; j > 0; j--) {
- if (data[maxp] <= data[peaks[j-1]])
- break;
- else if (j < peaks.size())
- peaks[j] = peaks[j-1];
- }
- if (j != peaks.size())
- peaks[j] = maxp;
- if (peakCount != peaks.size())
- peakCount++;
- }
- mid++;
- }
- return peakCount;
- } // findPeaks()
+ static int findPeaks(const vector &data, vector peaks, int width);
/** General peak picking method for finding local maxima in an array
* @param data input data
@@ -87,83 +55,16 @@
* @return list of peak indexes
*/
static vector findPeaks(const vector &data, int width,
- double threshold, double decayRate, bool isRelative) {
- vector peaks;
- int maxp = 0;
- int mid = 0;
- int end = data.size();
- double av = data[0];
- while (mid < end) {
- av = decayRate * av + (1 - decayRate) * data[mid];
- if (av < data[mid])
- av = data[mid];
- int i = mid - width;
- if (i < 0)
- i = 0;
- int stop = mid + width + 1;
- if (stop > data.size())
- stop = data.size();
- maxp = i;
- for (i++; i < stop; i++)
- if (data[i] > data[maxp])
- maxp = i;
- if (maxp == mid) {
- if (overThreshold(data, maxp, width, threshold, isRelative,av)){
- peaks.push_back(maxp);
- }
- }
- mid++;
- }
- return peaks;
- } // findPeaks()
+ double threshold, double decayRate, bool isRelative);
static double expDecayWithHold(double av, double decayRate,
- const vector &data, int start, int stop) {
- while (start < stop) {
- av = decayRate * av + (1 - decayRate) * data[start];
- if (av < data[start])
- av = data[start];
- start++;
- }
- return av;
- } // expDecayWithHold()
+ const vector &data, int start, int stop);
static bool overThreshold(const vector &data, int index, int width,
- double threshold, bool isRelative,
- double av) {
- if (data[index] < av)
- return false;
- if (isRelative) {
- int iStart = index - pre * width;
- if (iStart < 0)
- iStart = 0;
- int iStop = index + post * width;
- if (iStop > data.size())
- iStop = data.size();
- double sum = 0;
- int count = iStop - iStart;
- while (iStart < iStop)
- sum += data[iStart++];
- return (data[index] > sum / count + threshold);
- } else
- return (data[index] > threshold);
- } // overThreshold()
+ double threshold, bool isRelative,
+ double av);
- static void normalise(vector &data) {
- double sx = 0;
- double sxx = 0;
- for (int i = 0; i < data.size(); i++) {
- sx += data[i];
- sxx += data[i] * data[i];
- }
- double mean = sx / data.size();
- double sd = sqrt((sxx - sx * mean) / data.size());
- if (sd == 0)
- sd = 1; // all data[i] == mean -> 0; avoids div by 0
- for (int i = 0; i < data.size(); i++) {
- data[i] = (data[i] - mean) / sd;
- }
- } // normalise()
+ static void normalise(vector &data);
/** Uses an n-point linear regression to estimate the slope of data.
* @param data input data
@@ -172,48 +73,13 @@
* @param slope output data
*/
static void getSlope(const vector &data, double hop, int n,
- vector &slope) {
- int i = 0, j = 0;
- double t;
- double sx = 0, sxx = 0, sy = 0, sxy = 0;
- for ( ; i < n; i++) {
- t = i * hop;
- sx += t;
- sxx += t * t;
- sy += data[i];
- sxy += t * data[i];
- }
- double delta = n * sxx - sx * sx;
- for ( ; j < n / 2; j++)
- slope[j] = (n * sxy - sx * sy) / delta;
- for ( ; j < data.size() - (n + 1) / 2; j++, i++) {
- slope[j] = (n * sxy - sx * sy) / delta;
- sy += data[i] - data[i - n];
- sxy += hop * (n * data[i] - sy);
- }
- for ( ; j < data.size(); j++)
- slope[j] = (n * sxy - sx * sy) / delta;
- } // getSlope()
+ vector &slope);
static double min(const vector &arr) { return arr[imin(arr)]; }
-
static double max(const vector &arr) { return arr[imax(arr)]; }
- static int imin(const vector &arr) {
- int i = 0;
- for (int j = 1; j < arr.size(); j++)
- if (arr[j] < arr[i])
- i = j;
- return i;
- } // imin()
-
- static int imax(const vector &arr) {
- int i = 0;
- for (int j = 1; j < arr.size(); j++)
- if (arr[j] > arr[i])
- i = j;
- return i;
- } // imax()
+ static int imin(const vector &arr);
+ static int imax(const vector &arr);
}; // class Peaks