ITAThirdOctaveFIRFilterGenerator.cpp

#include <ITAThirdOctaveFIRFilterGenerator.h>

#include <ITAConstants.h>
#include <ITAFastMath.h>
#include <ITANumericUtils.h>
#include <ITAStringUtils.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>

#include <spline.h>

using namespace ITABase;

CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const double dSampleRate, const int iFilterLength )
	: m_dSamplerate( dSampleRate )
	, m_iFilterLength( iFilterLength )
	, m_ypp( nullptr )
	, m_pfInputFreqs( nullptr )
	, m_pfInputData( nullptr )
	, m_pfBuf1( nullptr )
	, m_pfBuf2( nullptr )
	, m_bWindow( false )
{
	m_iInputFreqs = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
	m_pfInputFreqs = fm_falloc( m_iInputFreqs, true );
	m_pfInputFreqs[ 0 ] = 0;	// Left margin
	for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
		m_pfInputFreqs[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
	m_pfInputFreqs[ m_iInputFreqs - 1 ] = ( float ) dSampleRate / 2;	// Right margin: Nyquist frequency

	m_pfInputData = fm_falloc( m_iInputFreqs, true );

	// DFT frequency bandwidth
	m_fDeltaF = ( float ) dSampleRate / ( float ) iFilterLength;

	// Number of symetric DFT coefficients;
	m_iDFTCoeffs = iFilterLength / 2 + 1;

	m_pfBuf1 = fm_falloc( 2 * m_iDFTCoeffs, false );
	m_pfBuf2 = fm_falloc( iFilterLength, false );
	m_pfWindow = fm_falloc( iFilterLength, false );

	// Windowing function (Hann window)
	float c = 2 * ITAConstants::PI_F / ( float ) ( m_iFilterLength - 1 );
	for( int i = 0; i < m_iFilterLength; i++ )
	{
		m_pfWindow[ i ] = 0.5F * ( 1 - cos( c * i ) );
	}

	m_ifft.plan( ITAFFT::IFFT_C2R, iFilterLength, m_pfBuf1, m_pfBuf2 );

	m_sDumpFilename = "interpolated_magnitudes.csv";
}

CITAThirdOctaveFIRFilterGenerator::~CITAThirdOctaveFIRFilterGenerator()
{
	fm_free( m_pfInputFreqs );
	fm_free( m_pfInputData );
	fm_free( m_pfBuf1 );
	fm_free( m_pfBuf2 );
	fm_free( m_pfWindow );
}

int CITAThirdOctaveFIRFilterGenerator::GetFilterLength() const
{
	return m_iFilterLength;
}
int CITAThirdOctaveFIRFilterGenerator::GetLatency() const
{
	// Latency = Half DFT period (ceil)
	return uprdiv( m_iFilterLength, 2 );
}

double CITAThirdOctaveFIRFilterGenerator::GetAverageRuntime() const
{
	return m_sw.mean();
}

void CITAThirdOctaveFIRFilterGenerator::SetDumpFilename( const std::string& sFilename )
{
	m_sDumpFilename = sFilename;
}

void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ )
{
	m_sw.start();

	if( oTOGainMagnitudes.IsZero() )
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = 0.0f;
		return;
	}

	if( oTOGainMagnitudes.IsIdentity() )
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = 0.0f;
		pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
		return;
	}

	// 1st step: Interpolate the magnitudes

	for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
		m_pfInputData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );

	// Bounndaries must be defined in a meaningful way.
	m_pfInputData[ 0 ] = m_pfInputData[ 1 ];
	m_pfInputData[ m_iInputFreqs - 1 ] = 0.0f;

	// Initialize cubic spline interpolation
	m_ypp = spline_cubic_set( m_iInputFreqs,
		m_pfInputFreqs,
		m_pfInputData,
		1, // Left boundary condition => 1st derivative m=0
		0,
		1, // Right boundary condition => 1st derivative m=0
		0 );
	float fDummy;
	const float fScale = 1 / ( float ) m_iFilterLength;

	// No DC offset, ever!
	m_pfBuf1[ 0 ] = 0;
	m_pfBuf1[ 1 ] = 0;


	if( bMinimumPhase ) {
		for( int i = 1; i < m_iDFTCoeffs; i++ )
		{
			float x = spline_cubic_val( m_iInputFreqs,
				m_pfInputFreqs,
				i*m_fDeltaF,
				m_pfInputData,
				m_ypp,
				&fDummy,
				&fDummy );

			// Phase-shift by half the FFT-period
			m_pfBuf1[ 2 * i ] = pow( x * fScale, 2 ) * m_iFilterLength; //minimum phase
			m_pfBuf1[ 2 * i + 1 ] = 0;
		}
	}
	else {
		for( int i = 1; i < m_iDFTCoeffs; i++ )
		{
			float x = spline_cubic_val( m_iInputFreqs,
				m_pfInputFreqs,
				i*m_fDeltaF,
				m_pfInputData,
				m_ypp,
				&fDummy,
				&fDummy );

			// Phase-shift by half the FFT-period: Negate all odd DFT coefficients
			m_pfBuf1[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
			m_pfBuf1[ 2 * i + 1 ] = 0;
		}
	}

	// 2nd step: Convert into time-domain (out-of-place C2R-IFFT)
	m_ifft.execute();

	// 3rd (optional) step: Hann window in the time-domain (optional)
	if( m_bWindow )
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = m_pfBuf2[ i ] * m_pfWindow[ i ];
	}
	else
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = m_pfBuf2[ i ];
	}

	// @todo: Minimum-phase?

	m_sw.stop();
}