/*
* 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
*
*
* Channelizer.C
*
* Implementation of class Channelizer.
*
* Kelly Fitz, 21 July 2000
* loris@cerlsoundgroup.org
*
* http://www.cerlsoundgroup.org/Loris/
*
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include "Channelizer.h"
#include "Envelope.h"
#include "LinearEnvelope.h"
#include "Partial.h"
#include "PartialList.h"
#include "Notifier.h"
#include <cmath>
// begin namespace
namespace Loris {
// ---------------------------------------------------------------------------
// Channelizer constructor from reference envelope
// ---------------------------------------------------------------------------
//! Construct a new Channelizer using the specified reference
//! Envelope to represent the a numbered channel. If the sound
//! being channelized is known to have detuned harmonics, a
//! stretching factor can be specified (defaults to 0 for no
//! stretching). The stretching factor can be computed using
//! the static member computeStretchFactor.
//!
//! \param refChanFreq is an Envelope representing the center frequency
//! of a channel.
//! \param refChanLabel is the corresponding channel number (i.e. 1
//! if refChanFreq is the lowest-frequency channel, and all
//! other channels are harmonics of refChanFreq, or 2 if
//! refChanFreq tracks the second harmonic, etc.).
//! \param stretchFactor is a stretching factor to account for detuned
//! harmonics, default is 0.
//!
//! \throw InvalidArgument if refChanLabel is not positive.
//! \throw InvalidArgument if stretchFactor is negative.
//
Channelizer::Channelizer( const Envelope & refChanFreq, int refChanLabel, double stretchFactor ) :
_refChannelFreq( refChanFreq.clone() ),
_refChannelLabel( refChanLabel ),
_stretchFactor( stretchFactor ),
_ampWeighting( 0 )
{
if ( refChanLabel <= 0 )
{
Throw( InvalidArgument, "Channelizer reference label must be positive." );
}
if ( stretchFactor < 0. )
{
Throw( InvalidArgument, "Channelizer stretch factor must be non-negative." );
}
}
// ---------------------------------------------------------------------------
// Channelizer constructor from constant reference frequency
// ---------------------------------------------------------------------------
//! Construct a new Channelizer having a constant reference frequency.
//! The specified frequency is the center frequency of the lowest-frequency
//! channel (for a harmonic sound, the channel containing the fundamental
//! Partial.
//!
//! \param refFreq is the reference frequency (in Hz) corresponding
//! to the first frequency channel.
//! \param stretchFactor is a stretching factor to account for detuned
//! harmonics, default is 0.
//!
//! \throw InvalidArgument is the reference frequency is not positive
Channelizer::Channelizer( double refFreq, double stretchFactor ) :
_refChannelFreq( new LinearEnvelope( refFreq ) ),
_refChannelLabel( 1 ),
_stretchFactor( stretchFactor ),
_ampWeighting( 0 )
{
if ( refFreq <= 0 )
{
Throw( InvalidArgument, "Channelizer reference frequency must be positive." );
}
if ( stretchFactor < 0. )
{
Throw( InvalidArgument, "Channelizer stretch factor must be non-negative." );
}
}
// ---------------------------------------------------------------------------
// Channelizer copy constructor
// ---------------------------------------------------------------------------
//! Construct a new Channelizer that is an exact copy of another.
//! The copy represents the same set of frequency channels, constructed
//! from the same reference Envelope and channel number.
//!
//! \param other is the Channelizer to copy
//
Channelizer::Channelizer( const Channelizer & other ) :
_refChannelFreq( other._refChannelFreq->clone() ),
_refChannelLabel( other._refChannelLabel ),
_stretchFactor( other._stretchFactor ),
_ampWeighting( other._ampWeighting )
{
}
// ---------------------------------------------------------------------------
// Channelizer assignment
// ---------------------------------------------------------------------------
//! Assignment operator: make this Channelizer an exact copy of another.
//! This Channelizer is made to represent the same set of frequency channels,
//! constructed from the same reference Envelope and channel number as @a rhs.
//!
//! \param rhs is the Channelizer to copy
//
Channelizer &
Channelizer::operator=( const Channelizer & rhs )
{
if ( &rhs != this )
{
_refChannelFreq.reset( rhs._refChannelFreq->clone() );
_refChannelLabel = rhs._refChannelLabel;
_stretchFactor = rhs._stretchFactor;
_ampWeighting = rhs._ampWeighting;
}
return *this;
}
// ---------------------------------------------------------------------------
// Channelizer destructor
// ---------------------------------------------------------------------------
//! Destroy this Channelizer.
//
Channelizer::~Channelizer( void )
{
}
// ---------------------------------------------------------------------------
// amplitudeWeighting
// ---------------------------------------------------------------------------
//! Return the exponent applied to amplitude before weighting
//! the instantaneous estimate of the frequency channel number
//! for a Partial. zero (default) for no weighting, 1 for linear
//! amplitude weighting, 2 for power weighting, etc.
//! Amplitude weighting is a bad idea for many sounds, particularly
//! those with transients, for which it may emphasize the part of
//! the Partial having the least reliable frequency estimate.
//
double Channelizer::amplitudeWeighting( void ) const
{
return _ampWeighting;
}
// ---------------------------------------------------------------------------
// setAmplitudeWeighting
// ---------------------------------------------------------------------------
//! Set the exponent applied to amplitude before weighting
//! the instantaneous estimate of the frequency channel number
//! for a Partial. zero (default) for no weighting, 1 for linear
//! amplitude weighting, 2 for power weighting, etc.
//! Amplitude weighting is a bad idea for many sounds, particularly
//! those with transients, for which it may emphasize the part of
//! the Partial having the least reliable frequency estimate.
//
void Channelizer::setAmplitudeWeighting( double x )
{
_ampWeighting = x;
}
// ---------------------------------------------------------------------------
// stretchFactor
// ---------------------------------------------------------------------------
//! Return the stretching factor used to account for detuned
//! harmonics, as in a piano tone. Normally set to 0 for
//! in-tune harmonics.
//!
//! The stretching factor is a small positive number for
//! heavy vibrating strings (as in pianos) for which the
//! mass of the string significantly affects the frequency
//! of the vibrating modes. See Martin Keane, "Understanding
//! the complex nature of the piano tone", 2004, for a discussion
//! and the source of the mode frequency stretching algorithms
//! implemented here.
//
double Channelizer::stretchFactor( void ) const
{
return _stretchFactor;
}
// ---------------------------------------------------------------------------
// setStretchFactor
// ---------------------------------------------------------------------------
//! Set the stretching factor used to account for detuned
//! harmonics, as in a piano tone. Normally set to 0 for
//! in-tune harmonics. The stretching factor for massy
//! vibrating strings (like pianos) can be computed from
//! the physical characteristics of the string, or using
//! computeStretchFactor().
//!
//! The stretching factor is a small positive number for
//! heavy vibrating strings (as in pianos) for which the
//! mass of the string significantly affects the frequency
//! of the vibrating modes. See Martin Keane, "Understanding
//! the complex nature of the piano tone", 2004, for a discussion
//! and the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! \throw InvalidArgument if stretch is negative.
//
void Channelizer::setStretchFactor( double stretch )
{
if ( stretch < 0. )
{
Throw( InvalidArgument, "Channelizer stretch factor must be non-negative." );
}
_stretchFactor = stretch;
}
// ---------------------------------------------------------------------------
// setStretchFactor
// ---------------------------------------------------------------------------
//! Set the stretching factor used to account for (consistently)
//! detuned harmonics, as in a piano tone, from a pair of
//! mode (harmonic) frequencies and numbers.
//!
//! The stretching factor is a small positive number for
//! heavy vibrating strings (as in pianos) for which the
//! mass of the string significantly affects the frequency
//! of the vibrating modes. See Martin Keane, "Understanding
//! the complex nature of the piano tone", 2004, for a discussion
//! and the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! The stretching factor is computed using computeStretchFactor,
//! but only a valid stretch factor will ever be assigned. If an
//! invalid (negative) stretching factor is computed for the
//! specified frequencies and mode numbers, the stretch factor
//! will be set to zero.
//!
//! \param fm is the frequency of the Mth stretched harmonic
//! \param m is the harmonic number of the harmonic whose frequnecy is fm
//! \param fn is the frequency of the Nth stretched harmonic
//! \param n is the harmonic number of the harmonic whose frequnecy is fn
void Channelizer::setStretchFactor( double fm, int m, double fn, int n )
{
const double B = computeStretchFactor( fm, m, fn, n );
if ( 0 < B )
{
_stretchFactor = B;
}
else
{
_stretchFactor = 0;
}
}
// ---------------------------------------------------------------------------
// computeStretchFactor (STATIC class member)
// ---------------------------------------------------------------------------
//! Static member to compute the stretch factor for a sound having
//! (consistently) detuned harmonics, like piano tones.
//!
//! The stretching factor is a small positive number for
//! heavy vibrating strings (as in pianos) for which the
//! mass of the string significantly affects the frequency
//! of the vibrating modes. See Martin Keane, "Understanding
//! the complex nature of the piano tone", 2004, for a discussion
//! and the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! The value returned by this function MAY NOT be a valid stretch
//! factor. If this function returns a negative stretch factor,
//! then the specified pair of frequencies and mode numbers cannot
//! be used to estimate the effects of string mass on mode frequency
//! (because the negative stretch factor implies a physical
//! impossibility, like negative mass or negative length).
//!
//! \param fm is the frequency of the Mth stretched harmonic
//! \param m is the harmonic number of the harmonic whose frequnecy is fm
//! \param fn is the frequency of the Nth stretched harmonic
//! \param n is the harmonic number of the harmonic whose frequnecy is fn
//! \returns the stretching factor, usually a very small positive
//! floating point number, or 0 for pefectly tuned harmonics
//! (that is, if fn = n*f1).
//
double
Channelizer::computeStretchFactor( double fm, int m, double fn, int n )
{
if ( fm <= 0. || fn <= 0. )
{
Throw( InvalidArgument, "Channelizer stretched harmonic frequencies must be positive." );
}
if ( m <= 0 || n <= 0 )
{
Throw( InvalidArgument, "Channelizer stretched harmonic numbers must be positive." );
}
// K is a factor that depends on the frequencies
// of the two stretched harmonics, equal to 1.0 for
// perfectly tuned (not stretched) harmonics
const double K = (m*fn) / (n*fm);
const double num = 1. - (K*K);
const double denom = (K*K*m*m) - (n*n);
/*
OLD and wrong I think
double num = (fn*fn) - (n*n*fref*fref);
double denom = (n*n*n*n)*(fref*fref);
*/
return num / denom;
}
// ---------------------------------------------------------------------------
// computeStretchFactor (STATIC class member)
// ---------------------------------------------------------------------------
//! Static member to compute the stretch factor for a sound having
//! (consistently) detuned harmonics, like piano tones. Legacy version
//! that assumes the first argument corresponds to the first partial.
//!
//! \param f1 is the frequency of the lowest numbered (1) partial.
//! \param fn is the frequency of the Nth stretched harmonic
//! \param n is the harmonic number of the harmonic whose frequnecy is fn
//! \returns the stretching factor, usually a very small positive
//! floating point number, or 0 for pefectly tuned harmonics
//! (that is, for harmonic frequencies fn = n*f1).
//
double
Channelizer::computeStretchFactor( double f1, double fn, double n )
{
return computeStretchFactor( f1, 1, fn, int(n + 0.5) );
}
// ---------------------------------------------------------------------------
// referenceFrequencyAt
// ---------------------------------------------------------------------------
//! Compute the reference frequency at the specified time. For non-stretched
//! harmonics, this is simply the ratio of the reference envelope evaluated
//! at that time to the reference channel number, and is the center frequecy
//! for the lowest channel. For stretched harmonics, the reference frequency
//! is NOT equal to the center frequency of any of the channels, and is also
//! a function of the stretch factor.
//
double Channelizer::referenceFrequencyAt( double time ) const
{
const double N = _refChannelLabel;
double fref = _refChannelFreq->valueAt( time ) / N;
if ( 0 != _stretchFactor )
{
double divisor = std::sqrt( 1.0 + ( _stretchFactor*N*N) );
fref = fref / divisor;
}
return fref;
}
// ---------------------------------------------------------------------------
// computeFractionalChannelNumber
// ---------------------------------------------------------------------------
//! Compute the (fractional) channel number estimate for a Partial having a
//! given frequency at a specified time. For ordinary harmonics, this
//! is simply the ratio of the specified frequency to the reference
//! frequency at the specified time. For stretched harmonics (as in
//! a piano), the stretching factor is used to compute the frequency
//! of the corresponding modes of a massy string. See Martin Keane,
//! "Understanding the complex nature of the piano tone", 2004, for
//! the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! The fractional channel number is used internally to determine
//! a best estimate for the channel number (label) for a Partial
//! having time-varying frequency.
//!
//! \param time is the time (in seconds) at which to evalute
//! the reference envelope
//! \param frequency is the frequency (in Hz) for wihch the channel
//! number is to be determined
//! \return the fractional channel number corresponding to the specified
//! frequency and time
//
double
Channelizer::computeFractionalChannelNumber( double time, double frequency ) const
{
double refFreq = referenceFrequencyAt( time );
if ( 0 == _stretchFactor )
{
return frequency / refFreq;
}
/*
const double frefsqrd = fref*fref;
double num = sqrt( (frefsqrd*frefsqrd) + (4*stretch*frefsqrd*fn*fn) ) - (frefsqrd);
double denom = 2*stretch*frefsqrd;
return sqrt( num / denom );
*/
// else:
// avoid squaring big numbers... two sqrts kind of sucks too.
const double rB = 1. / _stretchFactor; // reciprocal of B, the stretch factor
const double fratio = frequency / refFreq;
return std::sqrt( std::sqrt( (.25 * rB * rB) + (fratio * fratio * rB) ) - (.5 * rB) );
}
// ---------------------------------------------------------------------------
// computeChannelNumber
// ---------------------------------------------------------------------------
//! Compute the (fractional) channel number estimate for a Partial having a
//! given frequency at a specified time. For ordinary harmonics, this
//! is simply the ratio of the specified frequency to the reference
//! frequency at the specified time. For stretched harmonics (as in
//! a piano), the stretching factor is used to compute the frequency
//! of the corresponding modes of a massy string. See Martin Keane,
//! "Understanding the complex nature of the piano tone", 2004, for
//! the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! \param time is the time (in seconds) at which to evalute
//! the reference envelope
//! \param frequency is the frequency (in Hz) for wihch the channel
//! number is to be determined
//! \return the channel number corresponding to the specified
//! frequency and time
//
int
Channelizer::computeChannelNumber( double time, double frequency ) const
{
return int( computeFractionalChannelNumber( time, frequency ) + 0.5 );
}
// ---------------------------------------------------------------------------
// channelFrequencyAt
// ---------------------------------------------------------------------------
//! Compute the center frequency of one a channel at the specified
//! time. For non-stretched harmonics, this is simply the value
//! of the reference envelope scaled by the ratio of the specified
//! channel number to the reference channel number. For stretched
//! harmonics, the channel center frequency is computed using the
//! stretch factor. See Martin Keane, "Understanding
//! the complex nature of the piano tone", 2004, for a discussion
//! and the source of the mode frequency stretching algorithms
//! implemented here.
//!
//! \param time is the time (in seconds) at which to evalute
//! the reference envelope
//! \param channel is the frequency channel (or harmonic, or vibrational
//! mode) number whose frequency is to be determined
//! \return the center frequency in Hz of the specified frequency channel
//! at the specified time
//
double
Channelizer::channelFrequencyAt( double time, int channel ) const
{
const double fref = referenceFrequencyAt( time );
double fn = channel * fref;
if ( 0 != _stretchFactor )
{
const double scale = std::sqrt( 1.0 + (_stretchFactor*channel*channel) );
fn = fn * scale;
}
return fn;
}
// ---------------------------------------------------------------------------
// channelize (one Partial)
// ---------------------------------------------------------------------------
//! Label a Partial with the number of the frequency channel corresponding to
//! the average frequency over all the Partial's Breakpoints.
//!
//! \param partial is the Partial to label.
//
void
Channelizer::channelize( Partial & partial ) const
{
using std::pow;
debugger << "channelizing Partial with " << partial.numBreakpoints() << " Breakpoints" << endl;
// compute an amplitude-weighted average channel
// label for each Partial:
//double ampsum = 0.;
double weightedlabel = 0.;
Partial::const_iterator bp;
for ( bp = partial.begin(); bp != partial.end(); ++bp )
{
double f = bp.breakpoint().frequency();
double t = bp.time();
double weight = 1;
if ( 0 != _ampWeighting )
{
// This used to be an amplitude-weighted avg, but for many sounds,
// particularly those for which the weighted avg would be very
// different from the simple avg, the amplitude-weighted avg
// emphasized the part of the sound in which the frequency estimates
// are least reliable (e.g. a piano tone). The unweighted
// average should give more intuitive results in most cases.
// use sinusoidal amplitude:
double a = bp.breakpoint().amplitude() * std::sqrt( 1. - bp.breakpoint().bandwidth() );
weight = pow( a, _ampWeighting );
}
weightedlabel += weight * computeFractionalChannelNumber( t, f );
}
int label = 0;
if ( 0 < partial.numBreakpoints() ) // should always be the case
{
label = (int)((weightedlabel / partial.numBreakpoints()) + 0.5);
}
Assert( label >= 0 );
// assign label, and remember it, but
// only if it is a valid (positive)
// distillation label:
partial.setLabel( label );
}
// -- simplified interface --
// ---------------------------------------------------------------------------
// channelize (static simplified interface)
// ---------------------------------------------------------------------------
//! Static member that constructs an instance and applies
//! it to a PartialList (simplified interface).
//!
//! Construct a Channelizer using the specified Envelope
//! and reference label, and use it to channelize a
//! sequence of Partials.
//!
//! \param partials is the sequence of Partials to
//! channelize.
//! \param refChanFreq is an Envelope representing the center frequency
//! of a channel.
//! \param refChanLabel is the corresponding channel number (i.e. 1
//! if refChanFreq is the lowest-frequency channel, and all
//! other channels are harmonics of refChanFreq, or 2 if
//! refChanFreq tracks the second harmonic, etc.).
//! \throw InvalidArgument if refChanLabel is not positive.
//
void
Channelizer::channelize( PartialList & partials,
const Envelope & refChanFreq, int refChanLabel )
{
Channelizer instance( refChanFreq, refChanLabel );
for ( PartialList::iterator it = partials.begin(); it != partials.end(); ++it )
{
instance.channelize( *it );
}
}
} // end of namespace Loris