Mercurial > hg > beatroot-vamp
changeset 15:887c629502a9
refactor: pull method implementations into .cpp files
| author | Chris Cannam |
|---|---|
| date | Wed, 12 Oct 2011 10:55:52 +0100 |
| parents | f1252b6a7cf5 |
| children | 33d0b18b2509 |
| files | Agent.cpp Agent.h AgentList.cpp AgentList.h BeatRootProcessor.cpp BeatRootProcessor.h BeatTracker.h Induction.cpp Induction.h Makefile Peaks.cpp Peaks.h |
| diffstat | 12 files changed, 548 insertions(+), 549 deletions(-) [+] |
line wrap: on
line diff
--- 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
--- 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.
--- 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() +
--- 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 <code>stop</code> 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
--- 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<int> peaks = Peaks::findPeaks(spectralFlux, (int)lrint(0.06 / hop), 0.35, 0.84, true); + onsets.clear(); + onsets.resize(peaks.size(), 0); + vector<int>::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() +
--- 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<int> peaks = Peaks::findPeaks(spectralFlux, (int)lrint(0.06 / hop), 0.35, 0.84, true); - onsets.clear(); - onsets.resize(peaks.size(), 0); - vector<int>::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 */
--- 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<double> 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<double> d) { - if (d.size() < 2) - return -1.0; - vector<double> 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
--- 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<int> bestn;// count of high-scoring clusters + bestn.resize(topN); + + double ratio, err; + int degree; + int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth); + vector<double> clusterMean; + clusterMean.resize(maxClusterCount); + vector<int> clusterSize; + clusterSize.resize(maxClusterCount); + vector<int> 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() +
--- 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<int> bestn;// count of high-scoring clusters - bestn.resize(topN); - - double ratio, err; - int degree; - int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth); - vector<double> clusterMean; - clusterMean.resize(maxClusterCount); - vector<int> clusterSize; - clusterSize.resize(maxClusterCount); - vector<int> 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().
--- 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
--- 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<double> &data, vector<int> 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<int> Peaks::findPeaks(const vector<double> &data, int width, + double threshold, double decayRate, bool isRelative) { + vector<int> 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<double> &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<double> &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<double> &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<double> &data, double hop, int n, + vector<double> &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<double> &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<double> &arr) { + int i = 0; + for (int j = 1; j < arr.size(); j++) + if (arr[j] > arr[i]) + i = j; + return i; +} // imax() +
--- 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<double> &data, vector<int> 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<double> &data, vector<int> 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<int> findPeaks(const vector<double> &data, int width, - double threshold, double decayRate, bool isRelative) { - vector<int> 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<double> &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<double> &data, int start, int stop); static bool overThreshold(const vector<double> &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<double> &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<double> &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<double> &data, double hop, int n, - vector<double> &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<double> &slope); static double min(const vector<double> &arr) { return arr[imin(arr)]; } - static double max(const vector<double> &arr) { return arr[imax(arr)]; } - static int imin(const vector<double> &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<double> &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<double> &arr); + static int imax(const vector<double> &arr); }; // class Peaks