Commit 6c9da095 authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files

Adding single-in multi-out variable delay line (unfinished)

parent 62c3c542
## ITADSP
ITADSP is a C++ digital signal processing library for acoustics. It is a component from [ITACoreLibs](https://git.rwth-aachen.de/ita/ITACoreLibs), a collection of C++ libraries for virtual acoustics.
ITADSP is a C++ digital signal processing library for acoustics.
It is a component of [ITACoreLibs](https://git.rwth-aachen.de/ita/ITACoreLibs), a collection of C++ libraries for virtual acoustics.
### License
See [LICENSE](LICENSE.md) file.
Copyright 2015-2017 Institute of Technical Acoustics, RWTH Aachen University
Licensed under the Apache License, Version 2.0 (the "License");
you may not use files of this project except in compliance with the License.
You may obtain a copy of the License at
<http://www.apache.org/licenses/LICENSE-2.0>
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
### Quick build guide
Follow instructions from Wiki pages of [ITABase](https://git.rwth-aachen.de/ita/ITABase/wikis/home) project.
It is recommended to clone and follow the build guide of the parent project [ITACoreLibs](https://git.rwth-aachen.de/ita/ITACoreLibs/wikis/home), which includes this project as a submodule.
......@@ -276,88 +276,4 @@ private:
};
//! Abstract interface for interpolation algorithms
/**
* For equidistant sampling.
*
*/
class IITAVDLInterpolationRoutine
{
public:
//! Human-readable name
virtual std::string GetName() const = 0;
//! Anzahl an Stützwerten, die linksseitig (Anfang) bzw. rechtseitig (Ende) benötigt werden
/**
* Auf dem Überlappungsbereich wird die Interpolation angewendet, auch wenn sie
* auf diesen keine gültigen Daten liefern kann. Erst ab dem Offset iLeft und bis
* zum Ende abzüglich iRight werden in der Interpolationsroutine korrekte Daten
* abgelegt.
*
* \param iLeft Anzahl Samples zur linken Seite (call-by-reference)
* \param iRight Anzahl Samples zur rechten Seite (call-by-reference)
*/
virtual void GetOverlapSamples( int& iLeft, int& iRight ) const = 0;
//! Interpolation/Resampling
/**
* An input buffer with a given length will be interpolated starting at a given input offset.
* As many samples will be generated, as defined for the output buffer sample length. In- and
* output bufers are samples equidistant. An overlapping area is required in the source buffer
* to ensure a sufficient covererage of base samples (depending on algorithm). Those left-side
* and right-side samples can be retrieved via the function GetOverlapSamples().
*
* \param pInput In buffer pointer
* \param iInputLength Length of input samples
* \param iInputStartOffset Offset on input buffer (usually = left-side overlap samples)
* \param pOutput Out buffer pointer
* \param iOutputLength Length of requested output sample
* \param iOutputOffset Offset in output buffer
*
* \return True, if interpolation algorithm could be executed
*
*/
virtual bool Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset = 0 ) const = 0;
};
//! Linear interpolation (fast)
class CITAVDLLinearInterpolation : public IITAVDLInterpolationRoutine
{
public:
inline CITAVDLLinearInterpolation() {};
inline ~CITAVDLLinearInterpolation() {};
inline std::string GetName() const { return "Linear Interpolation"; };
inline void GetOverlapSamples( int& iLeft, int& iRight ) const { iLeft = 1; iRight = 0; };
bool Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset = 0 ) const;
};
//! Cubic-Spline interpolation (efficient)
class CITAVDLCubicSplineInterpolation : public IITAVDLInterpolationRoutine
{
public:
inline CITAVDLCubicSplineInterpolation() {};
inline ~CITAVDLCubicSplineInterpolation() {};
inline std::string GetName() const { return "Cubic Spline Interpolation"; };
inline void GetOverlapSamples( int& iLeft, int& iRight ) const { iLeft = 2; iRight = 2; };
bool Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset = 0 ) const;
};
// Window size (entire samples of window, forced to be uneven)
#define WINDOWED_SINC_INTERPOLATION_WINDOW_SIZE 13
//! Windowed Sinc-Interpolation (costly)
class CITAVDLWindowedSincInterpolation : public IITAVDLInterpolationRoutine
{
public:
inline CITAVDLWindowedSincInterpolation() : m_iWindowSize( WINDOWED_SINC_INTERPOLATION_WINDOW_SIZE ) { m_iWindowSize = m_iWindowSize + m_iWindowSize % 1; };
inline ~CITAVDLWindowedSincInterpolation() {};
inline std::string GetName() const { return "Windowed Sinc-Interpolation"; };
inline void GetOverlapSamples( int& iLeft, int& iRight ) const { iLeft = m_iWindowSize / 2; iRight = m_iWindowSize / 2; };
bool Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset = 0 ) const;
private:
int m_iWindowSize;
};
#endif // IW_ITA_SIMO_VARIABLE_DELAY_LINE
......@@ -3,12 +3,12 @@
#include <ITAFastMath.h>
#include <ITANumericUtils.h>
#include <ITASampleBuffer.h>
#include <ITAInterpolation.h>
#include <spline.h>
#include <cassert>
#include <cmath>
#include "ITAConstants.h"
// Interpolation resampling factor range (otherwise crossfading is used)
#define MIN_RESAMPLING_FACTOR 0.0f // [0, inf) ... aber sinnvoll z.B. 0, 0.1, 0.5 (< 1)
......@@ -396,300 +396,3 @@ void CITASIMOVariableDelayLine::ReadBlock( const int iCursorID, ITASampleBuffer*
double t = m_swProcess.stop();
return;
}
// Calculate cubic spline set (second derivatives, ypp) for equidistant data
void spline_cubic_set_equidistant( const int n, const float* y, float* ypp );
// Evaluate cubic spline interpolation point using only base values (y) and second derivatives (ypp) for equidistant data
float spline_cubic_val_equidistant( const int n, const float fX, const float* y, const float* ypp );
bool CITAVDLLinearInterpolation::Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset /*=0*/ ) const
{
assert( pOutput->length() >= iOutputLength + iOutputOffset );
// Interpolation ist nur mglich, wenn der Eingabe-Offset auf dem Eingabepuffer
// grer oder gleich 1 ist, da sonst keine Sttzwerte vorhanden sind.
assert( iInputStartOffset > 0 );
if( iInputStartOffset < 1 )
return false;
// Samples, die aufgeholt (>0) bzw. abgebremst (<0) werden mssen
int iDeltaDelay = ( iInputLength - iInputStartOffset ) - iOutputLength;
// Resamplingfaktor
float r = 1 + iDeltaDelay / ( float ) iOutputLength;
// Grenzen prfen, in denen Resampling betrieben werden kann.
// Dies kann schon mal auftreten, wenn eine Quelle manuell versetzt wird
// und sich dadurch eine zu hohe Relativgeschwindigkeit ergibt.
if( ( r > 1.5 ) || ( r < 0.5 ) )
return false;
// TODO klren ob diese resampling grenzen berhaupt gelten, theor. gibt es keine einschrnkung hier!
// --= Resampling =--
if( r < 1 ) // Oversampling
{
// ber Ausgabesamples (Laufvariable i ber den Ausgangspuffer)
// laufen und interpolierte Werte eintragen, dabei evtl.
// Samples des berlappungspuffers am Anfang einbeziehen
for( int i = 0; i < iOutputLength; i++ )
{
// Position in der Eingabe berechnen
float x_input = ( i + 1 )*r;
// Die Position des letzten Samples muss bereinstimmen mit dem letzten
// Eingabesample, dann wurde die Zeit erfolgreich eingeholt bzw.
// abgebremst.
if( i == iOutputLength - 1 )
{
assert( iOutputLength + iDeltaDelay == iInputLength - iInputStartOffset );
}
// Linkes/rechtes Nachbarsample in der Eingabe
int a = ( int ) floor( x_input );
int b = ( int ) ceil( x_input );
// Relative Position zwischen den Nachbarsamples in der Eingabe
float frac = x_input - ( float ) a;
// Die Linke Seite liegt immer innerhalb des Ausgangs
assert( ( a >= 0 ) && ( a < iInputLength ) );
// Sample genau getroffen (z.B. a=b=0)
if( a == b )
{
( *pOutput )[ i + iOutputOffset ] = ( *pInput )[ a ];
}
else
{
// Zwischenwert anhand linkem und rechtem Sttzwert ermitteln
assert( ( b >= 1 ) && ( b < iInputLength ) );
// Steigung (leicht lesbare Implementierung, wird vom Compiler wegoptimiert)
float left_val = ( *pInput )[ a ];
float right_val = ( *pInput )[ b ];
float m = right_val - left_val;
float val = left_val + m*frac;
( *pOutput )[ i + iOutputOffset ] = val;
}
}
}
else if( r > 1 ) // Undersampling
{
// ber Ausgabesamples (Laufvariable i ber den Ausgangspuffer)
// laufen und interpolierte Werte eintragen, dabei evtl.
// Samples des berlappungspuffers am Anfang einbeziehen
for( int i = 0; i < iOutputLength; i++ )
{
// Position in der Eingabe berechnen
float x_input = ( i + 1 )*r;
if( i == iOutputLength - 1 )
{
assert( x_input == ( float ) iInputLength - iInputStartOffset );
}
// Linkes/rechtes Nachbarsample in der Eingabe
int a = ( int ) floor( x_input );
int b = ( int ) ceil( x_input );
// Relative Position zwischen den Nachbarsamples in der Eingabe
float frac = x_input - ( float ) a;
// Die Linke Seite liegt immer innerhalb des Ausgangs
assert( ( a >= 0 ) && ( a < ( iInputLength ) ) );
// Sample genau getroffen (z.B. a=b=0)
if( a == b )
{
( *pOutput )[ i + iOutputOffset ] = ( *pInput )[ a ];
}
else
{
// Zwischenwert anhand linkem und rechtem Sttzwert ermitteln
assert( ( b >= 1 ) && ( b < iInputLength ) );
// Steigung
float m = ( *pInput )[ b ] - ( *pInput )[ a ];
( *pOutput )[ i + iOutputOffset ] = ( *pInput )[ a ] + m*frac;
}
}
}
else // Kein Resampling; TODO evtl immer resamplen, dieser Fall sollte vor Aufruf der Interpolation bereits abgefangen werden
{
// Nur Kopieren
pOutput->write( pInput, iOutputLength, iInputStartOffset );
}
return false;
}
bool CITAVDLCubicSplineInterpolation::Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset /*=0*/ ) const
{
// Eingabe validieren
assert( iInputStartOffset >= 2 );
assert( iOutputLength > 0 );
assert( pOutput->length() >= iOutputLength + iOutputOffset );
// Samples, die aufgeholt (>0) bzw. abgebremst (<0) werden mssen
int iDeltaDelay = ( iInputLength - iInputStartOffset ) - iOutputLength;
// Resamplingfaktor
float r = 1 + iDeltaDelay / ( float ) iOutputLength;
// --= Resampling =--
const float* pInputData = pInput->data();
float* pOutputData = pOutput->data();
// Return value
float* pInputSecondDerivatives = new float[ iInputLength ];
spline_cubic_set_equidistant( iInputLength, pInputData, pInputSecondDerivatives ); // Natural spline interpolation
for( int i = 0; i < iOutputLength; i++ )
{
// Position in der Eingabe berechnen (relativ zum Anfang mit Offset)
float fXInput = ( i + 1 )*r - 1;
if( i == iOutputLength - 1 )
assert( iOutputLength + iDeltaDelay == iInputLength - iInputStartOffset );
// Spline-Segment auswerten
float fYOutput = spline_cubic_val_equidistant( iInputLength, fXInput + iInputStartOffset, pInputData, pInputSecondDerivatives );
pOutputData[ i + iOutputOffset ] = fYOutput;
}
delete[] pInputSecondDerivatives;
return true;
}
bool CITAVDLWindowedSincInterpolation::Interpolate( const ITASampleBuffer* pInput, int iInputLength, int iInputStartOffset, ITASampleBuffer* pOutput, int iOutputLength, int iOutputOffset /*=0*/ ) const
{
// Eingabe validieren
assert( iInputStartOffset >= m_iWindowSize / 2 );
assert( iOutputLength > 0 );
assert( pOutput->length() >= iOutputLength + iOutputOffset );
// Samples, die aufgeholt (>0) bzw. abgebremst (<0) werden mssen
int iDeltaDelay = ( iInputLength - iInputStartOffset ) - iOutputLength;
// Resamplingfaktor
float r = 1 + iDeltaDelay / ( float ) iOutputLength;
// --= Resampling =--
const float* pInputData = pInput->data();
float* pOutputData = pOutput->data();
for( int i = 0; i < iOutputLength; i++ )
{
// Position in der Eingabe berechnen (relativ zum Anfang mit Offset)
float fXInput = ( i + 1 )*r - 1;
if( i == iOutputLength - 1 )
assert( iOutputLength + iDeltaDelay == iInputLength - iInputStartOffset );
// Fenster anwenden & Si-Rekonstruktion
pOutputData[ i + iOutputOffset ] = 0.0f;
for( int j = ( int ) ceil( fXInput - m_iWindowSize / 2 ); j < ( int ) floor( fXInput + m_iWindowSize / 2 ); j++ )
{
int iCenterIndex = j + iInputStartOffset;
assert( iCenterIndex >= 0 && iCenterIndex < iInputLength + m_iWindowSize / 2 ); // berlauf verhindern
float fAmplitude = pInputData[ iCenterIndex ];
// Relativer Abstand zum aktuellen Lesecursor in der Eingabe
float fSincArg = fXInput - iCenterIndex + iInputStartOffset;
assert( std::fabs( fSincArg ) <= ( float ) m_iWindowSize / 2 ); // Fensterberlauf verhindern
float fSincValue = ( fSincArg != 0.0f ) ? ( sin( ITAConstants::PI_F * fSincArg ) / ( ITAConstants::PI_F * fSincArg ) ) : 1.0f;
// Rect
//float fWindowValue = 1.0f;
// Hanning
//float fWindowValue = 0.5f * (1-cos(2*PI_F*fSincArg/(float) (m_iWindowSize/2)));
// Hamming
//float fWindowValue = 0.54f - 0.46f * cos(2*PI_F*fSincArg/(float) (m_iWindowSize/2));
// Lanczos
float a = m_iWindowSize / 2.0f;
float fWindowValue = ( fSincArg != 0.0f ) ? ( sin( ITAConstants::PI_F * fSincArg / a ) / ( ITAConstants::PI_F * fSincArg / a ) ) : 1.0f;
float fContribution = fWindowValue * fAmplitude * fSincValue;
pOutputData[ i + iOutputOffset ] += fContribution;
}
}
return true;
}
void spline_cubic_set_equidistant( const int n, const float* y, float* ypp )
{
assert( n > 2 );
// Temporrer Puffer
float* a = new float[ 3 * n - 1 ];
// Setup equations, save intermediate values to b
ypp[ 0 ] = 0.0f;
a[ 1 ] = 1.0f; // 1+0*3
a[ 3 ] = -1.0f; // 0+1*3
for( int i = 1; i < n - 1; i++ )
{
ypp[ i ] = y[ i + 1 ] - 2 * y[ i ] + y[ i - 1 ];
a[ 2 + ( i - 1 ) * 3 ] = 1.0f / 6.0f;
a[ 1 + i * 3 ] = 2.0f / 3.0f;
a[ 0 + ( i + 1 ) * 3 ] = 1.0f / 6.0f;
}
ypp[ n - 1 ] = 0.0f;
a[ 2 + ( n - 2 ) * 3 ] = -1.0f;
a[ 1 + ( n - 1 ) * 3 ] = 1.0f;
// Lineares Gleichungssystem lsen
a[ 1 + 1 * 3 ] = 5.0f / 6.0f;
ypp[ 1 ] = ypp[ 1 ] - ypp[ 0 ] / 6.0f;
float xmult;
for( int i = 2; i < n; i++ )
{
xmult = a[ 2 + ( i - 1 ) * 3 ] / a[ 1 + ( i - 1 ) * 3 ];
a[ 1 + i * 3 ] = a[ 1 + i * 3 ] - xmult * a[ 0 + i * 3 ];
ypp[ i ] = ypp[ i ] - xmult * ypp[ i - 1 ];
}
// Rckwrts einsetzen
for( int i = n - 2; 0 <= i; i-- )
ypp[ i ] = ( ypp[ i ] - a[ 0 + ( i + 1 ) * 3 ] * ypp[ i + 1 ] ) / a[ 1 + i * 3 ];
delete[] a;
}
float spline_cubic_val_equidistant( const int n, const float fX, const float* y, const float* ypp )
{
assert( n > fX );
int iXInt = ( int ) std::floor( fX );
float fFrac = fX - iXInt;
float y_interp = y[ iXInt ]
+ fFrac*( ( y[ iXInt + 1 ] - y[ iXInt ] ) - ( ypp[ iXInt + 1 ] / 6.0f + ypp[ iXInt ] / 3.0f )
+ fFrac * ( 0.5f * ypp[ iXInt ] + fFrac * ( ypp[ iXInt + 1 ] - ypp[ iXInt ] ) / 6.0f ) );
return y_interp;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment