[partial_tracking] Bug fix: Add custom implementation of

Loris PartialBuilder::buildPartials that works in 
real-time.

5 files changed, 284 insertions, 94 deletions

diff --git a/src/loris/PartialBuilder.C b/src/loris/PartialBuilder.C
index f1ce3be..87fa0e1 100644
--- a/src/loris/PartialBuilder.C
+++ b/src/loris/PartialBuilder.C
@@ -1,6 +1,6 @@
 /*
- * This is the Loris C++ Class Library, implementing analysis, 
- * manipulation, and synthesis of digitized sounds using the Reassigned 
+ * This is the Loris C++ Class Library, implementing analysis,
+ * manipulation, and synthesis of digitized sounds using the Reassigned
  * Bandwidth-Enhanced Additive Sound Model.
  *
  * Loris is Copyright (c) 1999-2010 by Kelly Fitz and Lippold Haken
@@ -26,12 +26,12 @@
  * extracted from a reassigned time-frequency spectrum to form ridges
  * and construct Partials.
  *
- * This strategy attemps to follow a mFreqWarping frequency envelope when 
+ * This strategy attemps to follow a mFreqWarping frequency envelope when
  * forming Partials, by prewarping all peak frequencies according to the
- * (inverse of) frequency mFreqWarping envelope. At the end of the analysis, 
+ * (inverse of) frequency mFreqWarping envelope. At the end of the analysis,
  * Partial frequencies need to be un-warped by calling fixPartialFrequencies().
  *
- * The first attempt was the same as the basic partial formation strategy, 
+ * The first attempt was the same as the basic partial formation strategy,
  * but for purposes of matching, peak frequencies are scaled by the ratio
  * of the mFreqWarping envelope's value at the previous frame to its value
  * at the current frame. This was not adequate, didn't store enough history
@@ -82,6 +82,7 @@ PartialBuilder::PartialBuilder( double drift ) :
 	mFreqWarping( new BreakpointEnvelope(1.0) ),
 	mFreqDrift( drift )
 {
+    reset();
 }
 
 // ---------------------------------------------------------------------------
@@ -99,6 +100,7 @@ PartialBuilder::PartialBuilder( double drift, const Envelope & env ) :
 	mFreqWarping( env.clone() ),
 	mFreqDrift( drift )
 {
+    reset();
 }
 
 // --- local helpers for Partial building ---
@@ -117,20 +119,26 @@ static inline double end_frequency( const Partial & partial )
 //	freq_distance
 // ---------------------------------------------------------------------------
 //	Helper function, used in formPartials().
-//	Returns the (positive) frequency distance between a Breakpoint 
+//	Returns the (positive) frequency distance between a Breakpoint
 //	and the last Breakpoint in a Partial.
 //
-inline double 
+inline double
 PartialBuilder::freq_distance( const Partial & partial, const SpectralPeak & pk )
 {
     double normBpFreq = pk.frequency() / mFreqWarping->valueAt( pk.time() );
-    
-    double normPartialEndFreq = 
+
+    double normPartialEndFreq =
         partial.last().frequency() / mFreqWarping->valueAt( partial.endTime() );
-    
+
 	return std::fabs( normPartialEndFreq - normBpFreq );
 }
 
+inline double
+PartialBuilder::freq_distance( const SpectralPeak & pk1, const SpectralPeak & pk2 )
+{
+    return std::fabs( pk1.frequency() - pk2.frequency() );
+}
+
 // ---------------------------------------------------------------------------
 //	better_match
 // ---------------------------------------------------------------------------
@@ -148,18 +156,55 @@ bool PartialBuilder::better_match( const Partial & part, const SpectralPeak & pk
                                    const SpectralPeak & pk2 )
 {
 	Assert( part.numBreakpoints() > 0 );
-	
+
 	return freq_distance( part, pk1 ) < freq_distance( part, pk2 );
-}	                                   
-                                   
-bool PartialBuilder::better_match( const Partial & part1, 
+}
+
+bool PartialBuilder::better_match( const Partial & part1,
                                    const Partial & part2, const SpectralPeak & pk )
 {
 	Assert( part1.numBreakpoints() > 0 );
 	Assert( part2.numBreakpoints() > 0 );
-	
+
 	return freq_distance( part1, pk ) < freq_distance( part2, pk );
-}	
+}
+
+bool PartialBuilder::better_match( const SpectralPeak & pk1,
+                                   const SpectralPeak & pk2,
+                                   const SpectralPeak & pk3 )
+{
+	return freq_distance( pk1, pk2 ) < freq_distance( pk2, pk3 );
+}
+
+// ---------------------------------------------------------------------------
+//	getNextActive
+// ---------------------------------------------------------------------------
+int PartialBuilder::getNextActive(int start)
+{
+    for(int i = start + 1; i < mCurrentPartials.size(); i++)
+    {
+        if(mActivePartials[i] && !mMatchedPartials[i])
+        {
+            return i;
+        }
+    }
+    return mCurrentPartials.size();
+}
+
+// ---------------------------------------------------------------------------
+//	getNextInActive
+// ---------------------------------------------------------------------------
+int PartialBuilder::getNextInactive(int start)
+{
+    for(int i = start + 1; i < mCurrentPartials.size(); i++)
+    {
+        if(!mActivePartials[i] && !mMatchedPartials[i])
+        {
+            return i;
+        }
+    }
+    return mCurrentPartials.size();
+}
 
 // --- Partial building members ---
 
@@ -171,31 +216,30 @@ bool PartialBuilder::better_match( const Partial & part1,
 //	be used to extend eliglble Partials spawn new Partials.
 //
 //	This is similar to the basic MQ partial formation strategy, except that
-//	before matching, all frequencies are normalized by the value of the 
+//	before matching, all frequencies are normalized by the value of the
 //	warping envelope at the time of the current frame. This means that
-//	the frequency envelopes of all the Partials are warped, and need to 
+//	the frequency envelopes of all the Partials are warped, and need to
 //	be un-normalized by calling finishBuilding at the end of the building
 //  process.
 //
-void 
+void
 PartialBuilder::buildPartials( Peaks & peaks, double frameTime )
-{	
+{
 	mNewlyEligible.clear();
-	
+
 	unsigned int matchCount = 0;	//	for debugging
-		
+
 	//	frequency-sort the spectral peaks:
 	//	(the eligible partials are always sorted by
 	//	increasing frequency if we always sort the
 	//	peaks this way)
 	std::sort( peaks.begin(), peaks.end(), SpectralPeak::sort_increasing_freq );
-	
+
 	PartialPtrs::iterator eligible = mEligiblePartials.begin();
-	for ( Peaks::iterator bpIter = peaks.begin(); bpIter != peaks.end(); ++bpIter ) 
+	for ( Peaks::iterator bpIter = peaks.begin(); bpIter != peaks.end(); ++bpIter )
 	{
-		//const Breakpoint & bp = bpIter->breakpoint;
 		const double peakTime = frameTime + bpIter->time();
-		
+
 		// 	find the Partial that is nearest in frequency to the Peak:
 		PartialPtrs::iterator nextEligible = eligible;
 		if ( eligible != mEligiblePartials.end() &&
@@ -208,39 +252,29 @@ PartialBuilder::buildPartials( Peaks & peaks, double frameTime )
 				++nextEligible;
 				++eligible;
 			}
-			
+
 			if ( nextEligible != mEligiblePartials.end() &&
 				 better_match( **nextEligible, **eligible, *bpIter ) )
 			{
 				eligible = nextEligible;
 			}
 		}
-		
+
 		// 	INVARIANT:
 		//
 		//	eligible is the position of the nearest (in frequency)
 		//	eligible Partial (pointer) or it is mEligiblePartials.end().
 		//
-		//	nextEligible is the eligible Partial with frequency 
-		//	greater than bp, or it is mEligiblePartials.end().  
-              
-#if defined(Debug_Loris) && Debug_Loris
-        /*
-		if ( nextEligible != mEligiblePartials.end() )
-		{
-			debugger << matchFrequency << "( " << end_frequency( **eligible )
-					 << ", " << end_frequency( **nextEligible ) << ")" << endl;
-		}
-        */
-#endif
-								
+		//	nextEligible is the eligible Partial with frequency
+		//	greater than bp, or it is mEligiblePartials.end().
+
+
 		//	create a new Partial if there is no eligible Partial,
-		//	or the frequency difference to the eligible Partial is 
-		//	too great, or the next peak is a better match for the 
+		//	or the frequency difference to the eligible Partial is
+		//	too great, or the next peak is a better match for the
 		//	eligible Partial, otherwise add this peak to the eligible
 		//	Partial:
-		Peaks::iterator nextPeak = //Peaks::iterator( bpIter ); ++nextPeak;
-								   ++Peaks::iterator( bpIter ); //  some compilers choke on this?
+		Peaks::iterator nextPeak = ++Peaks::iterator( bpIter );
 
         //  decide whether this match should be made:
         //  - can only make the match if eligible is not the end of the list
@@ -249,71 +283,208 @@ PartialBuilder::buildPartials( Peaks & peaks, double frameTime )
         bool makeMatch = false;
         if ( eligible != mEligiblePartials.end() )
         {
-            bool matchIsGood =  mFreqDrift > 
+            bool matchIsGood =  mFreqDrift >
                 std::fabs( end_frequency( **eligible ) - bpIter->frequency() );
             if ( matchIsGood )
             {
                 bool nextIsBetter = ( nextPeak != peaks.end() &&
-                                      better_match( **eligible, *nextPeak, *bpIter ) ); 
+                                      better_match( **eligible, *nextPeak, *bpIter ) );
                 if ( ! nextIsBetter )
                 {
                     makeMatch = true;
                 }
-            }                       
+            }
         }
-        
+
         Breakpoint bp = bpIter->createBreakpoint();
-        
+
         if ( makeMatch )
         {
             //  invariant:
             //  if makeMatch is true, then eligible is the position of a valid Partial
             (*eligible)->insert( peakTime, bp );
 			mNewlyEligible.push_back( *eligible );
-			
+
 			++matchCount;
         }
         else
         {
             Partial p;
             p.insert( peakTime, bp );
-            mCollectedPartials.push_back( p );
-            mNewlyEligible.push_back( & mCollectedPartials.back() );
+            /* mCollectedPartials.push_back( p ); */
+            /* mNewlyEligible.push_back( & mCollectedPartials.back() ); */
+            mNewlyEligible.push_back( & p );
         }
-        
+
 		//	update eligible, nextEligible is the eligible Partial
 		//	with frequency greater than bp, or it is mEligiblePartials.end():
 		eligible = nextEligible;
-	}			 
-	 	
+	}
+
 	mEligiblePartials = mNewlyEligible;
-	
-    /*
-	debugger << "PartialBuilder::buildPartials: matched " << matchCount << endl;
-	debugger << "PartialBuilder::buildPartials: " << mNewlyEligible.size() << " newly eligible partials" << endl;
-    */
+}
+
+void
+PartialBuilder::buildPartials( Peaks & peaks )
+{
+	unsigned int matchCount = 0;
+	unsigned int Np = mCurrentPartials.size();
+
+    for(int i = 0; i < mCurrentPartials.size(); i++)
+    {
+        mMatchedPartials[i] = false;
+    }
+
+    for(int p = 0; p < peaks.size(); p++)
+	{
+		// 	find the Partial that is nearest in frequency to the Peak
+        int eligible = getNextActive(-1);
+		int nextEligible = eligible;
+
+		if ( eligible < Np )
+		{
+            nextEligible = getNextActive(nextEligible);
+			while ( nextEligible < Np )
+			{
+                if ( better_match( mCurrentPartials[nextEligible],
+                                   peaks[p],
+                                   mCurrentPartials[eligible] ) )
+                {
+                    eligible = nextEligible;
+                }
+
+                nextEligible = getNextActive(nextEligible);
+			}
+		}
+
+		//	create a new Partial if there is no eligible Partial,
+		//	or the frequency difference to the eligible Partial is
+		//	too great, or the next peak is a better match for the
+		//	eligible Partial, otherwise add this peak to the eligible
+		//	Partial:
+		int nextPeak = p + 1;
+
+        //  decide whether this match should be made:
+        //  - can only make the match if eligible is not the end of the list
+        //  - the match is only good if it is close enough in frequency
+        //  - even if the match is good, only match if the next one is not better
+        bool makeMatch = false;
+        if ( eligible < Np )
+        {
+            bool matchIsGood = mFreqDrift >
+                std::fabs( mCurrentPartials[eligible].frequency() -
+                           peaks[p].frequency() );
+            if ( matchIsGood )
+            {
+                bool nextIsBetter = ( nextPeak < peaks.size() &&
+                                      better_match( peaks[nextPeak],
+                                                    mCurrentPartials[eligible],
+                                                    peaks[p] ) );
+                if ( ! nextIsBetter )
+                {
+                    makeMatch = true;
+                }
+            }
+        }
+
+        if ( makeMatch )
+        {
+            //  if makeMatch is true then eligible is the position of a
+            //  valid Partial
+            mCurrentPartials[eligible] = peaks[p];
+            mActivePartials[eligible] = true;
+            mMatchedPartials[eligible] = true;
+            ++matchCount;
+        }
+        else
+        {
+            // save to first inactive partial (if any)
+            int inactive = getNextInactive(-1);
+            if(inactive < Np)
+            {
+                mCurrentPartials[inactive] = peaks[p];
+                mActivePartials[inactive] = true;
+                mMatchedPartials[inactive] = true;
+            }
+        }
+	}
+
+    // kill inactive partials
+    for(int i = 0; i < mCurrentPartials.size(); i++)
+    {
+        if(!mMatchedPartials[i])
+        {
+            mCurrentPartials[i].setAmplitude(0.f);
+            mActivePartials[i] = false;
+        }
+    }
 }
 
 // ---------------------------------------------------------------------------
 //	finishBuilding
 // ---------------------------------------------------------------------------
-//	Un-do the frequency warping performed in buildPartials, and return 
+//	Un-do the frequency warping performed in buildPartials, and return
 //	the Partials that were built. After calling finishBuilding, the
 //  builder is returned to its initial state, and ready to build another
-//  set of Partials. Partials are returned by appending them to the 
+//  set of Partials. Partials are returned by appending them to the
 //  supplied PartialList.
 //
 void
 PartialBuilder::finishBuilding( PartialList & product )
-{	
+{
     //  append the collected Partials to the product list:
 	product.splice( product.end(), mCollectedPartials );
-    
+
     //  reset the builder state:
     mEligiblePartials.clear();
     mNewlyEligible.clear();
 }
 
+// ---------------------------------------------------------------------------
+//	getPartials
+// ---------------------------------------------------------------------------
+// Return partials by appending them to the supplised PartialList
+void
+PartialBuilder::getPartials( PartialList & product )
+{
+    //  append the collected Partials to the product list:
+	product.splice( product.end(), mCollectedPartials );
+}
+
+Peaks &
+PartialBuilder::getPartials()
+{
+    return mCurrentPartials;
+}
+
+// ---------------------------------------------------------------------------
+//	maxPartials
+// ---------------------------------------------------------------------------
+// Change the maximum number of partials per frame
+void
+PartialBuilder::maxPartials(int max)
+{
+    mCurrentPartials.resize(max);
+    mActivePartials.resize(max);
+    mMatchedPartials.resize(max);
+}
 
+// ---------------------------------------------------------------------------
+//	reset
+// ---------------------------------------------------------------------------
+// Reset the current partial list
+void
+PartialBuilder::reset()
+{
+    for(int i = 0; i < mCurrentPartials.size(); i++)
+    {
+        mCurrentPartials[i].setAmplitude(0.f);
+        mCurrentPartials[i].setFrequency(0.f);
+        mCurrentPartials[i].setPhase(0.f);
+        mCurrentPartials[i].setBandwidth(0.f);
+        mActivePartials[i] = false;
+        mMatchedPartials[i] = false;
+    }
+}
 
 }	//	end of namespace Loris
diff --git a/src/loris/PartialBuilder.h b/src/loris/PartialBuilder.h
index 2592ffc..0459010 100644
--- a/src/loris/PartialBuilder.h
+++ b/src/loris/PartialBuilder.h
@@ -94,6 +94,7 @@ public:
     //	be un-normalized by calling finishBuilding at the end of the building
     //  process.
 	void buildPartials( Peaks & peaks, double frameTime );
+	void buildPartials( Peaks & peaks );
 
     //  finishBuilding
     //
@@ -104,16 +105,38 @@ public:
     //  supplied PartialList.
 	void finishBuilding( PartialList & product );
 
+    // getPartials
+    //
+    // Return partials
+	void getPartials( PartialList & product );
+	Peaks & getPartials();
+
+    // maxPartials
+    //
+    // Change the maximum number of partials per frame
+    void maxPartials(int max);
+
+    // reset
+    //
+    // Reset the current partial list
+    void reset();
+
 private:
 
 // --- auxiliary member functions ---
 
     double freq_distance( const Partial & partial, const SpectralPeak & pk );
+    double freq_distance( const SpectralPeak & pk1, const SpectralPeak & pk2 );
 
     bool better_match( const Partial & part, const SpectralPeak & pk1,
 		   	           const SpectralPeak & pk2 );
     bool better_match( const Partial & part1, 
 	  		           const Partial & part2, const SpectralPeak & pk );                       
+    bool better_match( const SpectralPeak & pk1, const SpectralPeak & pk2,
+		   	           const SpectralPeak & pk3 );
+
+    int getNextActive(int start);
+    int getNextInactive(int start);
                        
                        
 // --- collected partials ---
@@ -122,8 +145,11 @@ private:
 
 // --- builder state variables ---
 		
+    std::vector<bool> mActivePartials;
+    std::vector<bool> mMatchedPartials;
 	PartialPtrs mEligiblePartials;
     PartialPtrs mNewlyEligible;                 // 	keep track of eligible partials here
+    Peaks mCurrentPartials;
 
 // --- parameters ---
     	
diff --git a/src/loris/SpectralPeaks.h b/src/loris/SpectralPeaks.h
index eeb9b7e..bdf3479 100644
--- a/src/loris/SpectralPeaks.h
+++ b/src/loris/SpectralPeaks.h
@@ -69,11 +69,14 @@ public:
     
     double amplitude( void ) const { return m_breakpoint.amplitude(); }
     double frequency( void ) const { return m_breakpoint.frequency(); }
+    double phase( void ) const { return m_breakpoint.phase(); }
     double bandwidth( void ) const { return m_breakpoint.bandwidth(); }
     
     //  --- mutation ---
     
     void setAmplitude( double x ) { m_breakpoint.setAmplitude(x); }
+    void setFrequency( double x ) { m_breakpoint.setFrequency(x); }
+    void setPhase( double x ) { m_breakpoint.setPhase(x); }
     void setBandwidth( double x ) { m_breakpoint.setBandwidth(x); }
     
     //  this REALLY shouldn't be in here...
diff --git a/src/simpl/partial_tracking.cpp b/src/simpl/partial_tracking.cpp
index 42df76b..99a8f34 100644
--- a/src/simpl/partial_tracking.cpp
+++ b/src/simpl/partial_tracking.cpp
@@ -308,11 +308,12 @@ Peaks SndObjPartialTracking::update_partials(Frame* frame) {
 // ---------------------------------------------------------------------------
 // LorisPartialTracking
 // ---------------------------------------------------------------------------
-SimplLorisPTAnalyzer::SimplLorisPTAnalyzer() :
+SimplLorisPTAnalyzer::SimplLorisPTAnalyzer(int max_partials) :
     Loris::Analyzer(50, 100),
     _env(1.0) {
     buildFundamentalEnv(false);
-    _partial_builder = new Loris::PartialBuilder(m_freqDrift, _env);
+    _partial_builder = new Loris::PartialBuilder(m_freqDrift);
+    _partial_builder->maxPartials(max_partials);
 }
 
 SimplLorisPTAnalyzer::~SimplLorisPTAnalyzer() {
@@ -322,20 +323,16 @@ SimplLorisPTAnalyzer::~SimplLorisPTAnalyzer() {
 void SimplLorisPTAnalyzer::analyze() {
     m_ampEnvBuilder->reset();
     m_f0Builder->reset();
-    m_partials.clear();
 
-    //  estimate the amplitude in this frame:
+    // estimate the amplitude in this frame
     m_ampEnvBuilder->build(peaks, 0);
 
-    //  collect amplitudes and frequencies and try to
-    //  estimate the fundamental
+    // estimate the fundamental frequency
     m_f0Builder->build(peaks, 0);
 
-    //  form Partials from the extracted Breakpoints:
-    _partial_builder->buildPartials(peaks, 0);
-
-    //  unwarp the Partial frequency envelopes:
-    _partial_builder->finishBuilding(m_partials);
+    // form partials
+    _partial_builder->buildPartials(peaks);
+    partials = _partial_builder->getPartials();
 }
 
 // ---------------------------------------------------------------------------
@@ -355,7 +352,7 @@ void LorisPartialTracking::reset() {
     if(_analyzer) {
         delete _analyzer;
     }
-    _analyzer = new SimplLorisPTAnalyzer();
+    _analyzer = new SimplLorisPTAnalyzer(_max_partials);
 }
 
 void LorisPartialTracking::max_partials(int new_max_partials) {
@@ -377,25 +374,18 @@ Peaks LorisPartialTracking::update_partials(Frame* frame) {
                                                  frame->peak(i)->amplitude,
                                                  frame->peak(i)->bandwidth,
                                                  frame->peak(i)->phase);
-        _analyzer->peaks.push_back(Loris::SpectralPeak(1, bp));
+        _analyzer->peaks.push_back(Loris::SpectralPeak(0, bp));
     }
 
     _analyzer->analyze();
-    _partials = _analyzer->partials();
-    int num_partials = _partials.size();
-
-    for(Loris::PartialListIterator i = _partials.begin(); i != _partials.end(); ++i) {
-        Peak* p = new Peak();
-        p->amplitude =  i->amplitudeAt(1);
-        p->frequency =  i->frequencyAt(1);
-        p->phase =  i->phaseAt(1);
-        p->bandwidth =  i->bandwidthAt(1);
-        peaks.push_back(p);
-        frame->add_partial(p);
-    }
+    int num_partials = _analyzer->partials.size();
 
-    for(int i = num_partials; i < _max_partials; i++) {
+    for(int i = 0; i < num_partials; i++) {
         Peak* p = new Peak();
+        p->amplitude =  _analyzer->partials[i].amplitude();
+        p->frequency =  _analyzer->partials[i].frequency();
+        p->phase =  _analyzer->partials[i].phase();
+        p->bandwidth =  _analyzer->partials[i].bandwidth();
         peaks.push_back(p);
         frame->add_partial(p);
     }
diff --git a/src/simpl/partial_tracking.h b/src/simpl/partial_tracking.h
index b9f2d28..4f8d60e 100644
--- a/src/simpl/partial_tracking.h
+++ b/src/simpl/partial_tracking.h
@@ -112,9 +112,10 @@ class SimplLorisPTAnalyzer : public Loris::Analyzer {
         Loris::PartialBuilder* _partial_builder;
 
     public:
-        SimplLorisPTAnalyzer();
+        SimplLorisPTAnalyzer(int max_partials);
         ~SimplLorisPTAnalyzer();
         Loris::Peaks peaks;
+        Loris::Peaks partials;
         void analyze();
 };
 
@@ -122,12 +123,11 @@ class SimplLorisPTAnalyzer : public Loris::Analyzer {
 class LorisPartialTracking : public PartialTracking {
     private:
         SimplLorisPTAnalyzer* _analyzer;
-        Loris::PartialList _partials;
-        void reset();
 
     public:
         LorisPartialTracking();
         ~LorisPartialTracking();
+        void reset();
         void max_partials(int new_max_partials);
         Peaks update_partials(Frame* frame);
 };