From c0415ce3b3b41a03ec9184ccfd1064786cf72442 Mon Sep 17 00:00:00 2001 From: Vesa Date: Thu, 18 Dec 2014 22:36:39 +0200 Subject: [PATCH] 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. --- include/Delay.h | 363 ++++++++++++++++++++++++++++++++++++++++ include/interpolation.h | 12 +- 2 files changed, 369 insertions(+), 6 deletions(-) create mode 100644 include/Delay.h diff --git a/include/Delay.h b/include/Delay.h new file mode 100644 index 000000000..1c51f1ba8 --- /dev/null +++ b/include/Delay.h @@ -0,0 +1,363 @@ +/* + * Delay.h - Delay effect objects to use as building blocks in DSP + * + * Copyright (c) 2014 Vesa Kivimäki + * Copyright (c) 2006-2014 Tobias Doerffel + * + * 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 +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( 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 +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( 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 +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( ceil( delay1 ) ); + m_fraction1 = 1.0 - ( delay1 - floor( delay1 ) ); + + m_delay2 = static_cast( 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 +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( 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 diff --git a/include/interpolation.h b/include/interpolation.h index cbe274d42..693897587 100644 --- a/include/interpolation.h +++ b/include/interpolation.h @@ -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 ); }