#include #include #include #include #include #include #include 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 m_pfInputData[0] = 1.0f; for (int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++) m_pfInputData[1 + i] = float(oTOGainMagnitudes[i]); m_pfInputData[m_iInputFreqs - 1] = 0.0f; // @todo jst: check if this is good // 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(); }