Moving interpolation routines (from VDL implementation) to ITABase and making...

Moving interpolation routines (from VDL implementation) to ITABase and making them available for public use

CMakeLists.txt

@@ -65,6 +65,7 @@ set( ITABaseHeader
 	"include/ITAHDFTSpectra.h"
 	"include/ITAHDFTSpectrum.h"
 	"include/ITAFunctors.h"
+	"include/ITAInterpolation.h"
 	"include/ITALog.h"
 	"include/ITAMagnitudeSpectrum.h"
 	"include/ITANumericUtils.h"
@@ -98,6 +99,7 @@ set( ITABaseSources
 	"src/ITAFileSystemUtils.cpp"
 	"src/ITAHDFTSpectra.cpp"
 	"src/ITAHDFTSpectrum.cpp"
+	"src/ITAInterpolation.cpp"
 	"src/ITALog.cpp"
 	"src/ITAMagnitudeSpectrum.cpp"
 	"src/ITANumericUtils.cpp"

include/ITAInterpolation.h

+/*
+ * ----------------------------------------------------------------
+ *
+ *		ITA core libs
+ *		(c) Copyright Institute of Technical Acoustics (ITA)
+ *		RWTH Aachen University, Germany, 2015-2017
+ *
+ * ----------------------------------------------------------------
+ *				    ____  __________  _______
+ *				   //  / //__   ___/ //  _   |
+ *				  //  /    //  /    //  /_|  |
+ *				 //  /    //  /    //  ___   |
+ *				//__/    //__/    //__/   |__|
+ *
+ * ----------------------------------------------------------------
+ *
+ */
+
+#ifndef IW_ITA_INTERPOLATION
+#define IW_ITA_INTERPOLATION
+
+#include <ITABaseDefinitions.h>
+
+#include <string>
+
+class ITASampleBuffer;
+
+//! Abstract interface for interpolation algorithms
+/**
+* For equidistant sampling.
+*
+*/
+class IITAVDLInterpolationRoutine
+{
+public:
+	//! Human-readable name
+	virtual std::string GetName() const = 0;
+
+	//! Anzahl an Sttzwerten, die linksseitig (Anfang) bzw. rechtseitig (Ende) bentigt werden
+	/**
+	* Auf dem berlappungsbereich wird die Interpolation angewendet, auch wenn sie
+	* auf diesen keine gltigen Daten liefern kann. Erst ab dem Offset iLeft und bis
+	* zum Ende abzglich iRight werden in der Interpolationsroutine korrekte Daten
+	* abgelegt.
+	*
+	* @param[out] iLeft Anzahl Samples zur linken Seite (call-by-reference)
+	* @param[out] 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, const int iInputLength, const int iInputStartOffset, ITASampleBuffer* pOutput, const int iOutputLength, const int iOutputOffset = 0 ) const = 0;
+};
+
+//! Linear interpolation (fast)
+class ITA_BASE_API 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, const int iInputLength, const int iInputStartOffset, ITASampleBuffer* pOutput, const int iOutputLength, const int iOutputOffset = 0 ) const;
+};
+
+//! Cubic-Spline interpolation (efficient)
+class ITA_BASE_API 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, const int iInputLength, const int iInputStartOffset, ITASampleBuffer* pOutput, const int iOutputLength, const int iOutputOffset = 0 ) const;
+};
+
+//! Windowed Sinc-Interpolation (costly)
+class ITA_BASE_API CITAVDLWindowedSincInterpolation : public IITAVDLInterpolationRoutine
+{
+public:
+	inline CITAVDLWindowedSincInterpolation( const int iWindowSize = 13 ) : m_iWindowSize( iWindowSize )
+	{
+		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, const int iInputLength, const int iInputStartOffset, ITASampleBuffer* pOutput, const int iOutputLength, const int iOutputOffset = 0 ) const;
+
+private:
+	int m_iWindowSize;
+};
+
+#endif // IW_ITA_INTERPOLATION

src/ITAInterpolation.cpp

+#include <ITAInterpolation.h>
+
+#include <ITAConstants.h>
+#include <ITASampleBuffer.h>
+
+#include <cassert>
+// 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;
+}