/*
* 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();
}