diff options
Diffstat (limited to 'src/loris/AssociateBandwidth.C')
| -rw-r--r-- | src/loris/AssociateBandwidth.C | 317 |
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(); + +} |