Changes to interpolation (explicitly use fma more), some new DSP building blocks

New file: Delay.h - contains some simple delay effects for use in DSP - perhaps for designing reverbs or similar. All are in double precision because why not.
This commit is contained in:
Vesa
2014-12-18 22:36:39 +02:00
parent 8402bcb525
commit c0415ce3b3
2 changed files with 369 additions and 6 deletions

363
include/Delay.h Normal file
View File

@@ -0,0 +1,363 @@
/*
* Delay.h - Delay effect objects to use as building blocks in DSP
*
* Copyright (c) 2014 Vesa Kivimäki <contact/dot/diizy/at/nbl/dot/fi>
* Copyright (c) 2006-2014 Tobias Doerffel <tobydox/at/users.sourceforge.net>
*
* This file is part of LMMS - http://lmms.io
*
* 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 (see COPYING); if not, write to the
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301 USA.
*
*/
#ifndef DELAY_H
#define DELAY_H
#include "lmms_basics.h"
#include "lmms_math.h"
#include "interpolation.h"
#include "MemoryManager.h"
// brief usage
// Classes:
// CombFeedback: a feedback comb filter - basically a simple delay line, makes a comb shape in the freq response
// CombFeedfwd: a feed-forward comb filter - an "inverted" comb filter, can be combined with CombFeedback to create a net allpass if negative gain is used
// CombFeedbackDualtap: same as CombFeedback but takes two delay values
// AllpassDelay: an allpass delay - combines feedback and feed-forward - has flat frequency response
// all classes are templated with channel count, any arbitrary channel count can be used for each fx
// Methods (for all classes):
// setDelay sets delay amount in frames. It's up to you to make this samplerate-agnostic.
// Fractions are allowed - linear interpolation is used to deal with them
// CombFeedbackDualTap is a special case: it requires 2 delay times
// setMaxDelay (re)sets the maximum allowed delay, in frames
// NOTE: for performance reasons, there's no bounds checking at setDelay, so make sure you set maxDelay >= delay!
// clearHistory clears the delay buffer
// setGain sets the feedback/feed-forward gain, in linear amplitude, negative values are allowed
// 1.0 is full feedback/feed-forward, -1.0 is full negative feedback/feed-forward
// update runs the fx for one frame - takes as arguments input and number of channel to run, returns output
template<ch_cnt_t CHANNELS>
class CombFeedback
{
public:
typedef double frame[CHANNELS];
CombFeedback( int maxDelay ) :
m_size( maxDelay ),
m_position( 0 ),
m_feedBack( 0.0 ),
m_delay( 0 ),
m_fraction( 0.0 )
{
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
virtual ~CombFeedback()
{
MM_FREE( m_buffer );
}
inline void setMaxDelay( int maxDelay )
{
if( maxDelay > m_size )
{
MM_FREE( m_buffer );
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
m_size = maxDelay;
m_position %= m_size;
}
inline void clearHistory()
{
memset( m_buffer, 0, sizeof( frame ) * m_size );
}
inline void setDelay( double delay )
{
m_delay = static_cast<int>( ceil( delay ) );
m_fraction = 1.0 - ( delay - floor( delay ) );
}
inline void setGain( double gain )
{
m_gain = gain;
}
inline double update( double in, ch_cnt_t ch )
{
int readPos = m_position - m_delay;
if( readPos < 0 ) { readPos += m_size; }
const double y = linearInterpolate( m_buffer[readPos][ch], m_buffer[( readPos + 1 ) % m_size][ch], m_fraction );
++m_position %= m_size;
m_buffer[m_position][ch] = in + m_gain * y;
return y;
}
private:
frame * m_buffer;
int m_size;
int m_position;
double m_gain;
int m_delay;
double m_fraction;
};
template<ch_cnt_t CHANNELS>
class CombFeedfwd
{
typedef double frame[CHANNELS];
CombFeedfwd( int maxDelay ) :
m_size( maxDelay ),
m_position( 0 ),
m_feedBack( 0.0 ),
m_delay( 0 ),
m_fraction( 0.0 )
{
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
virtual ~CombFeedfwd()
{
MM_FREE( m_buffer );
}
inline void setMaxDelay( int maxDelay )
{
if( maxDelay > m_size )
{
MM_FREE( m_buffer );
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
m_size = maxDelay;
m_position %= m_size;
}
inline void clearHistory()
{
memset( m_buffer, 0, sizeof( frame ) * m_size );
}
inline void setDelay( double delay )
{
m_delay = static_cast<int>( ceil( delay ) );
m_fraction = 1.0 - ( delay - floor( delay ) );
}
inline void setGain( double gain )
{
m_gain = gain;
}
inline double update( double in, ch_cnt_t ch )
{
int readPos = m_position - m_delay;
if( readPos < 0 ) { readPos += m_size; }
const double y = linearInterpolate( m_buffer[readPos][ch], m_buffer[( readPos + 1 ) % m_size][ch], m_fraction ) + in * m_gain;
++m_position %= m_size;
m_buffer[m_position][ch] = in;
return y;
}
private:
frame * m_buffer;
int m_size;
int m_position;
double m_gain;
int m_delay;
double m_fraction;
};
template<ch_cnt_t CHANNELS>
class CombFeedbackDualtap
{
typedef double frame[CHANNELS];
CombFeedbackDualtap( int maxDelay ) :
m_size( maxDelay ),
m_position( 0 ),
m_feedBack( 0.0 ),
m_delay( 0 ),
m_fraction( 0.0 )
{
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
virtual ~CombFeedbackDualtap()
{
MM_FREE( m_buffer );
}
inline void setMaxDelay( int maxDelay )
{
if( maxDelay > m_size )
{
MM_FREE( m_buffer );
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
m_size = maxDelay;
m_position %= m_size;
}
inline void clearHistory()
{
memset( m_buffer, 0, sizeof( frame ) * m_size );
}
inline void setDelays( double delay1, double delay2 )
{
m_delay1 = static_cast<int>( ceil( delay1 ) );
m_fraction1 = 1.0 - ( delay1 - floor( delay1 ) );
m_delay2 = static_cast<int>( ceil( delay2 ) );
m_fraction2 = 2.0 - ( delay2 - floor( delay2 ) );
}
inline void setGain( double gain )
{
m_gain = gain;
}
inline double update( double in, ch_cnt_t ch )
{
int readPos1 = m_position - m_delay1;
if( readPos1 < 0 ) { readPos1 += m_size; }
int readPos2 = m_position - m_delay2;
if( readPos2 < 0 ) { readPos2 += m_size; }
const double y = linearInterpolate( m_buffer[readPos1][ch], m_buffer[( readPos1 + 1 ) % m_size][ch], m_fraction1 ) +
linearInterpolate( m_buffer[readPos2][ch], m_buffer[( readPos2 + 1 ) % m_size][ch], m_fraction2 );
++m_position %= m_size;
m_buffer[m_position][ch] = in + m_gain * y;
return y;
}
private:
frame * m_buffer;
int m_size;
int m_position;
double m_gain;
int m_delay1;
int m_delay2;
double m_fraction1;
double m_fraction2;
};
template<ch_cnt_t CHANNELS>
class AllpassDelay
{
public:
typedef double frame[CHANNELS];
AllpassDelay( int maxDelay ) :
m_size( maxDelay ),
m_position( 0 ),
m_feedBack( 0.0 ),
m_delay( 0 ),
m_fraction( 0.0 )
{
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
virtual ~AllpassDelay()
{
MM_FREE( m_buffer );
}
inline void setMaxDelay( int maxDelay )
{
if( maxDelay > m_size )
{
MM_FREE( m_buffer );
m_buffer = MM_ALLOC( frame, maxDelay );
memset( m_buffer, 0, sizeof( frame ) * maxDelay );
}
m_size = maxDelay;
m_position %= m_size;
}
inline void clearHistory()
{
memset( m_buffer, 0, sizeof( frame ) * m_size );
}
inline void setDelay( double delay )
{
m_delay = static_cast<int>( ceil( delay ) );
m_fraction = 1.0 - ( delay - floor( delay ) );
}
inline void setGain( double gain )
{
m_gain = gain;
}
inline double update( double in, ch_cnt_t ch )
{
int readPos = m_position - m_delay;
if( readPos < 0 ) { readPos += m_size; }
const double y = linearInterpolate( m_buffer[readPos][ch], m_buffer[( readPos + 1 ) % m_size][ch], m_fraction ) + in * -m_gain;
const double x = in + m_gain * y;
++m_position %= m_size;
m_buffer[m_position][ch] = x;
return y;
}
private:
frame * m_buffer;
int m_size;
int m_position;
double m_gain;
int m_delay;
double m_fraction;
};
// convenience typedefs for stereo effects
typedef CombFeedback<2> StereoCombFeedback;
typedef CombFeedfwd<2> StereoCombFeedfwd;
typedef CombFeedbackDualtap<2> StereoCombFeedbackDualtap;
typedef AllpassDelay<2> StereoAllpassDelay;
#endif

View File

@@ -71,9 +71,9 @@ inline float cubicInterpolate( float v0, float v1, float v2, float v3, float x )
float frcu = frsq*v0;
float t1 = v3 + 3*v1;
return( v1 + 0.5f * frcu + x * ( v2 - frcu * ( 1.0f/6.0f ) -
t1 * ( 1.0f/6.0f ) - v0 / 3.0f ) + frsq * x * ( t1 *
( 1.0f/6.0f ) - 0.5f * v2 ) + frsq * ( 0.5f * v2 - v1 ) );
return( v1 + fastFmaf( 0.5f, frcu, x ) * ( v2 - frcu * ( 1.0f/6.0f ) -
fastFmaf( t1, ( 1.0f/6.0f ), -v0 ) * ( 1.0f/3.0f ) ) + frsq * x * ( t1 *
( 1.0f/6.0f ) - 0.5f * v2 ) + frsq * fastFmaf( 0.5f, v2, -v1 ) );
}
@@ -102,7 +102,7 @@ inline float optimalInterpolate( float v0, float v1, float x )
const float c2 = even * -0.004541102062639801;
const float c3 = odd * -1.57015627178718420;
return ( ( c3*z + c2 ) * z + c1 ) * z + c0;
return fastFmaf( fastFmaf( fastFmaf( c3, z, c2 ), z, c1 ), z, c0 );
}
@@ -119,7 +119,7 @@ inline float optimal4pInterpolate( float v0, float v1, float v2, float v3, float
const float c2 = even1 * -0.246185007019907091 + even2 * 0.24614027139700284;
const float c3 = odd1 * -0.36030925263849456 + odd2 * 0.10174985775982505;
return ( ( c3*z + c2 ) * z + c1 ) * z + c0;
return fastFmaf( fastFmaf( fastFmaf( c3, z, c2 ), z, c1 ), z, c0 );
}
@@ -130,7 +130,7 @@ inline float lagrangeInterpolate( float v0, float v1, float v2, float v3, float
const float c1 = v2 - v0 * ( 1.0f / 3.0f ) - v1 * 0.5f - v3 * ( 1.0f / 6.0f );
const float c2 = 0.5f * (v0 + v2) - v1;
const float c3 = ( 1.0f/6.0f ) * ( v3 - v0 ) + 0.5f * ( v1 - v2 );
return ( ( c3*x + c2 ) * x + c1 ) * x + c0;
return fastFmaf( fastFmaf( fastFmaf( c3, x, c2 ), x, c1 ), x, c0 );
}