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

Cleanup and moving from float array to float vector implementation

parent d2821402
......@@ -23,7 +23,6 @@
#include <ITAFFT.h>
#include <ITASampleBuffer.h>
#include <ITAStopWatch.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <string.h>
......@@ -56,7 +55,7 @@ public:
double GetAverageRuntime() const;
// Filter erzeugen
void GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase = false) const;
void GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase = false ) const;
//! Generate FIR filter from third octave (linear) magnitude spectrum
/**
......@@ -72,16 +71,12 @@ private:
int m_iDFTCoeffs; // Anzahl symetrischer DFT-Koeffizienten
float m_fDeltaF; // DFT-Frequenzauflsung [Hz]
int m_iNumInputFreqs; // Anzahl Frequenzsttzstellen fr Interpolation (31 Terzen + 2 Rnder)
float* m_pfInputFreqVector; // Frequenzsttzstellen der Eingabedaten (Terzen) @todo jst: std::vector
float* m_pfInputFreqData; // Zwischenspeicher fr die Interpolation @todo jst: std::vector
mutable float* m_ypp; // Interpolationsdaten
float* m_pfFreqDataInterpolated;
float* m_pfBuf2;
bool m_bWindow; // Apply window?
float* m_pfWindow;
std::vector< float > m_vfInputFreqVector; // Frequenzsttzstellen der Eingabedaten (Terzen) @todo jst: std::vector
mutable std::vector< float > m_vfInputFreqData; // Zwischenspeicher fr die Interpolation @todo jst: std::vector
mutable float* m_ypp; //! Cubic spline interplation data
mutable std::vector< float > m_vfFreqDataInterpolated;
std::vector< float > m_vfImpulseResponse;
mutable ITAFFT m_ifft;
mutable ITAStopWatch m_sw;
};
#endif // IW_ITA_THIRD_OCTAVE_FILTER_GENERATOR
......@@ -13,21 +13,15 @@ using namespace ITABase;
CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const double dSampleRate, const int iFilterLength )
: m_dSamplerate( dSampleRate )
, m_iFilterLength( iFilterLength )
, m_ypp( nullptr )
, m_pfInputFreqVector( nullptr )
, m_pfInputFreqData( nullptr )
, m_pfFreqDataInterpolated( nullptr )
, m_pfBuf2( nullptr )
, m_bWindow( false )
{
m_iNumInputFreqs = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
m_pfInputFreqVector = fm_falloc( m_iNumInputFreqs, true );
m_pfInputFreqVector[ 0 ] = 0; // Left margin
m_vfInputFreqVector.resize( m_iNumInputFreqs );
m_vfInputFreqVector[ 0 ] = 0; // Left margin
for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputFreqVector[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
m_pfInputFreqVector[ m_iNumInputFreqs - 1 ] = ( float ) dSampleRate / 2; // Right margin: Nyquist frequency
m_vfInputFreqVector[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
m_vfInputFreqVector[ m_iNumInputFreqs - 1 ] = ( float ) dSampleRate / 2; // Right margin: Nyquist frequency
m_pfInputFreqData = fm_falloc( m_iNumInputFreqs, true );
m_vfInputFreqData.resize( m_iNumInputFreqs );
// DFT frequency bandwidth
m_fDeltaF = ( float ) dSampleRate / ( float ) iFilterLength;
......@@ -35,27 +29,14 @@ CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const doub
// Number of symetric DFT coefficients;
m_iDFTCoeffs = iFilterLength / 2 + 1;
m_pfFreqDataInterpolated = 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_pfFreqDataInterpolated, m_pfBuf2 );
m_vfFreqDataInterpolated.resize( 2 * m_iDFTCoeffs );
m_vfImpulseResponse.resize( iFilterLength );
m_ifft.plan( ITAFFT::IFFT_C2R, iFilterLength, &m_vfFreqDataInterpolated[0], &m_vfImpulseResponse[0] );
}
CITAThirdOctaveFIRFilterGenerator::~CITAThirdOctaveFIRFilterGenerator()
{
fm_free( m_pfInputFreqVector );
fm_free( m_pfInputFreqData );
fm_free( m_pfFreqDataInterpolated );
fm_free( m_pfBuf2 );
fm_free( m_pfWindow );
}
int CITAThirdOctaveFIRFilterGenerator::GetFilterLength() const
......@@ -74,10 +55,6 @@ int CITAThirdOctaveFIRFilterGenerator::GetLatency() const
return uprdiv( m_iFilterLength, 2 );
}
double CITAThirdOctaveFIRFilterGenerator::GetAverageRuntime() const
{
return m_sw.mean();
}
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, ITASampleBuffer& sbImpulseResponse, bool bMinimumPhase /*=false*/ ) const
{
......@@ -88,8 +65,6 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ ) const
{
m_sw.start();
if( oTOGainMagnitudes.IsZero() )
{
for( int i = 0; i < m_iFilterLength; i++ )
......@@ -111,16 +86,16 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
// 1st step: Interpolate the magnitudes
for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputFreqData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );
m_vfInputFreqData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );
// Bounndaries must be defined in a meaningful way.
m_pfInputFreqData[ 0 ] = m_pfInputFreqData[ 1 ];
m_pfInputFreqData[ m_iNumInputFreqs - 1 ] = 0.0f;
m_vfInputFreqData[ 0 ] = m_vfInputFreqData[ 1 ];
m_vfInputFreqData[ m_iNumInputFreqs - 1 ] = 0.0f;
// Initialize cubic spline interpolation
m_ypp = spline_cubic_set( m_iNumInputFreqs,
m_pfInputFreqVector,
m_pfInputFreqData,
&m_vfInputFreqVector[0],
&m_vfInputFreqData[0],
1, // Left boundary condition => 1st derivative m=0
0,
1, // Right boundary condition => 1st derivative m=0
......@@ -129,24 +104,24 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
const float fScale = 1 / ( float ) m_iFilterLength;
// No DC offset, ever!
m_pfFreqDataInterpolated[ 0 ] = 0;
m_pfFreqDataInterpolated[ 1 ] = 0;
m_vfFreqDataInterpolated[ 0 ] = 0;
m_vfFreqDataInterpolated[ 1 ] = 0;
if( bMinimumPhase )
{
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iNumInputFreqs,
m_pfInputFreqVector,
&m_vfInputFreqVector[0],
i*m_fDeltaF,
m_pfInputFreqData,
&m_vfInputFreqData[0],
m_ypp,
&fDummy,
&fDummy );
// Phase-shift by half the FFT-period
m_pfFreqDataInterpolated[ 2 * i ] = x * fScale; //minimum phase
m_pfFreqDataInterpolated[ 2 * i + 1 ] = 0;
m_vfFreqDataInterpolated[ 2 * i ] = x * fScale; //minimum phase
m_vfFreqDataInterpolated[ 2 * i + 1 ] = 0;
}
}
else
......@@ -154,35 +129,23 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iNumInputFreqs,
m_pfInputFreqVector,
&m_vfInputFreqVector[0],
i*m_fDeltaF,
m_pfInputFreqData,
&m_vfInputFreqData[0],
m_ypp,
&fDummy,
&fDummy );
// Phase-shift by half the FFT-period: Negate all odd DFT coefficients
m_pfFreqDataInterpolated[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
m_pfFreqDataInterpolated[ 2 * i + 1 ] = 0;
m_vfFreqDataInterpolated[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
m_vfFreqDataInterpolated[ 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 ];
}
m_ifft.execute(); // populates m_vfImpulseResponse
// @todo: Minimum-phase?
m_sw.stop();
// Copy to output
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = m_vfImpulseResponse[ i ];
}
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