Start adding Loris files

1 files changed, 317 insertions, 0 deletions

diff --git a/src/loris/AssociateBandwidth.C b/src/loris/AssociateBandwidth.C
new file mode 100644
index 0000000..e72c52f
--- /dev/null
+++ b/src/loris/AssociateBandwidth.C
@@ -0,0 +1,317 @@
+/*
+ * 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
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY, without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *
+ *
+ * AssociateBandwidth.C
+ *
+ * Implementation of a class representing a policy for associating noise
+ * (bandwidth) energy with reassigned spectral peaks to be used in
+ * Partial formation.
+ *
+ * Kelly Fitz, 20 Jan 2000
+ * loris@cerlsoundgroup.org
+ *
+ * http://www.cerlsoundgroup.org/Loris/
+ *
+ */
+
+#if HAVE_CONFIG_H
+	#include "config.h"
+#endif
+
+#include "AssociateBandwidth.h"
+#include "Breakpoint.h"
+#include "BreakpointUtils.h"
+#include "LorisExceptions.h"
+#include "Notifier.h"
+#include "SpectralPeaks.h"
+
+#include <algorithm>
+#include <cmath>
+
+using namespace Loris;
+
+// ---------------------------------------------------------------------------
+//	AssociateBandwidth constructor
+// ---------------------------------------------------------------------------
+//	Association regions are centered on all integer bin frequencies, regionWidth
+//	is the total width (in Hz) of the overlapping bandwidth association regions, 
+//	the region centers are spaced at half this width.
+//
+AssociateBandwidth::AssociateBandwidth( double regionWidth, double srate ) :
+	_regionRate( 0 )
+{
+	if ( ! (regionWidth>0) )
+		Throw( InvalidArgument, "The regionWidth must be greater than 0 Hz." );
+	if ( ! (srate>0) )
+		Throw( InvalidArgument, "The sample rate must be greater than 0 Hz." );
+		
+		
+	_weights.resize( int(srate/regionWidth) );
+	_surplus.resize( int(srate/regionWidth) );
+	_regionRate = 2./regionWidth;
+}	
+
+// ---------------------------------------------------------------------------
+//	AssociateBandwidth destructor
+// ---------------------------------------------------------------------------
+//
+AssociateBandwidth::~AssociateBandwidth( void )
+{
+}
+
+// ---------------------------------------------------------------------------
+//	binFrequency
+// ---------------------------------------------------------------------------
+//	Compute the warped fractional bin/region frequency corresponding to 
+//	freqHz. (_regionRate is the number of regions per hertz.)
+//
+//	Once, we used bark frequency scale warping here, but there seems to be
+//	no reason to do so. The best results seem to be indistinguishable from
+// 	plain 'ol 1k bins, and some results are much worse.
+//	
+static double binFrequency( double freqHz, double regionRate )
+{
+//#define Use_Barks
+#ifndef Use_Barks
+	return freqHz * regionRate;
+#else
+	//	Compute Bark frequency from Hertz.
+	//	Got this formula for Bark frequency from Sciarraba's thesis.
+	//	Ignore region rate when using barks
+	double tmp = std::atan( ( 0.001 / 7.5 ) * freqHz );
+	return  13. * std::atan( 0.76 * 0.001 * freqHz ) + 3.5 * ( tmp * tmp );
+#endif
+	
+}
+
+// ---------------------------------------------------------------------------
+//	findRegionBelow
+// ---------------------------------------------------------------------------
+//	Return the index of the last region having center frequency less than
+//	or equal to freq, or -1 if no region is low enough. 
+//
+//	Note: the zeroeth region is centered at bin frequency 1 and tapers
+//	to zero at bin frequency 0! (when booger is 1.)
+//
+static int findRegionBelow( double binfreq, unsigned int howManyBins )
+{
+	const double booger = 0.;
+	if ( binfreq < booger ) 
+	{
+		return -1;
+	}
+	else 
+	{
+		return int( std::min( std::floor(binfreq - booger), howManyBins - 1. ) );
+	}
+}
+
+// ---------------------------------------------------------------------------
+//	computeAlpha
+// ---------------------------------------------------------------------------
+//	binfreq is a warped, fractional bin frequency, and bin frequencies 
+//	are integers. Return the relative contribution of a component at
+//	binfreq to the bins (bw association regions) below and above
+//	binfreq. 
+//
+static double computeAlpha( double binfreq, unsigned int howManyBins )
+{
+	//	everything above the center of the highest
+	//	bin is lumped into that bin; i.e it does
+	//	not taper off at higher frequencies:	
+	if ( binfreq > howManyBins ) 
+	{
+		return 0.;
+	}
+	else 
+	{
+		return binfreq - std::floor( binfreq );
+	}
+}
+
+// ---------------------------------------------------------------------------
+//	distribute
+// ---------------------------------------------------------------------------
+//
+static void distribute( double fractionalBin, double x, std::vector<double> & regions )
+{
+	//	contribute x to two regions having center
+	//	frequencies less and greater than freqHz:
+	int posBelow = findRegionBelow( fractionalBin, regions.size() );
+	int posAbove = posBelow + 1;
+	
+	double alpha = computeAlpha( fractionalBin, regions.size() );
+	
+	if ( posAbove < regions.size() )
+		regions[posAbove] += alpha * x;
+	
+	if ( posBelow >= 0 )
+		regions[posBelow] += (1. - alpha) * x;
+}
+
+// ---------------------------------------------------------------------------
+//	computeNoiseEnergy
+// ---------------------------------------------------------------------------
+//	Return the noise energy to be associated with a component at freqHz.
+//	_surplus contains the surplus spectral energy in each region, which is,
+//	by defintion, non-negative.
+//
+double 
+AssociateBandwidth::computeNoiseEnergy( double freq, double amp )
+{
+	//	don't mess with negative frequencies:
+	if ( freq < 0. )
+		return 0.;
+	
+	//	compute the fractional bin frequency 
+	//	corresponding to freqHz:
+	double bin = binFrequency( freq, _regionRate );
+	
+	//	contribute x to two regions having center
+	//	frequencies less and greater than freqHz:
+	int posBelow = findRegionBelow( bin, _surplus.size() );
+	int posAbove = posBelow + 1;
+
+	double alpha = computeAlpha( bin, _surplus.size() );
+
+	double noise = 0.;
+	//	Have to check for alpha == 0, because 
+	//	the weights will be zero (see computeAlpha()):
+	//	(ignore lowest regions)
+	const int LowestRegion = 2;
+	/*
+	if ( posAbove < _surplus.size() && alpha != 0. && posAbove >= LowestRegion )
+		noise += _surplus[posAbove] * alpha / _weights[posAbove];
+	
+	if ( posBelow >= LowestRegion )
+		noise += _surplus[posBelow] * (1. - alpha) / _weights[posBelow];
+	*/
+	//	new idea, weight Partials by amplitude:	
+	if ( posAbove < _surplus.size() && alpha != 0. && posAbove >= LowestRegion )
+		noise += _surplus[posAbove] * alpha * amp / _weights[posAbove];
+	
+	if ( posBelow >= LowestRegion )
+		noise += _surplus[posBelow] * (1. - alpha) * amp / _weights[posBelow];
+		
+	return noise;
+}
+
+// ---------------------------------------------------------------------------
+//	accumulateSinusoid
+// ---------------------------------------------------------------------------
+//	Accumulate sinusoidal energy at frequency f and amplitude a.
+//	The amplitude isn't used for anything.
+//	
+void
+AssociateBandwidth::accumulateSinusoid( double freq, double amp )
+{
+	//	distribute weight at the peak frequency,
+	//	don't mess with negative frequencies:
+	if ( freq > 0. )
+	{
+		//distribute( binFrequency( freq, _regionRate ), 1., _weights );
+		//	new idea: weight Partials by amplitude:
+		distribute( binFrequency( freq, _regionRate ), amp, _weights );
+	}
+}
+
+// ---------------------------------------------------------------------------
+//	accumulateNoise
+// ---------------------------------------------------------------------------
+//	Accumulate a rejected spectral peak as surplus (noise) energy.
+//
+void
+AssociateBandwidth::accumulateNoise( double freq, double amp )
+{
+	//	compute energy contribution and distribute 
+	//	at frequency f, don't mess with negative 
+	//	frequencies:
+	if ( freq > 0. )
+    {
+		distribute( binFrequency( freq, _regionRate ), amp * amp, _surplus  );
+    }
+}
+
+// ---------------------------------------------------------------------------
+//	associate
+// ---------------------------------------------------------------------------
+//	Associate bandwidth with a single SpectralPeak.
+//
+void 
+AssociateBandwidth::associate( SpectralPeak & pk )
+{		
+    pk.setBandwidth(0);
+    pk.addNoiseEnergy( computeNoiseEnergy( pk.frequency(), pk.amplitude() ) );
+}
+
+// ---------------------------------------------------------------------------
+//	reset
+// ---------------------------------------------------------------------------
+//	This is called after each distribution of bandwidth energy.
+//
+void
+AssociateBandwidth::reset( void )
+{
+	std::fill( _weights.begin(), _weights.end(), 0. );
+	std::fill( _surplus.begin(), _surplus.end(), 0. );
+}
+
+// ---------------------------------------------------------------------------
+//	associate
+// ---------------------------------------------------------------------------
+//	Perform bandwidth association on a collection of reassigned spectral peaks
+//	or ridges. The range [begin, rejected) spans the Peaks selected to form
+//	Partials. The range [rejected, end) spans the Peaks that were found in
+//	the reassigned spectrum, but rejected as too weak or too close (in 
+//	frequency) to another stronger Peak. 
+//
+void 
+AssociateBandwidth::associateBandwidth( Peaks::iterator begin, 		//	beginning of Peaks
+										Peaks::iterator rejected, 	//	first rejected Peak
+										Peaks::iterator end )		//	end of Peaks
+{		
+	if ( begin == rejected )
+		return;
+		
+	//	accumulate retained Breakpoints as sinusoids, 
+	for ( Peaks::iterator it = begin; it != rejected; ++it )
+	{
+		accumulateSinusoid( it->frequency(), it->amplitude() );
+	}
+	
+	//	accumulate rejected breakpoints as noise:
+	for ( Peaks::iterator it = rejected; it != end; ++it )
+	{
+		accumulateNoise( it->frequency(), it->amplitude() );
+	}
+
+	//	associate bandwidth with each retained Breakpoint:
+	for ( Peaks::iterator it = begin; it != rejected; ++it )
+	{
+        //  sets bandwidth to zero, then calls addNoiseEnergy()
+		associate( *it );
+	}
+	
+	//	reset after association, yuk:
+	reset();
+
+}