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

Fixing bug in minimum phase impl

parent 01486d0f
......@@ -21,9 +21,10 @@
#include <ITADSPDefinitions.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <ITAFFT.h>
#include <ITASampleBuffer.h>
#include <ITAStopWatch.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <string.h>
#include <vector>
......@@ -36,13 +37,6 @@
class ITA_DSP_API CITAThirdOctaveFIRFilterGenerator
{
public:
// Design methods
enum
{
LINEAR_PHASE = 0,
MINIMUM_PHASE
};
// Konstruktor
CITAThirdOctaveFIRFilterGenerator( const double dSampleRate, const int iFilterLength );
......@@ -61,27 +55,31 @@ public:
// Mittlere Laufzeit der Erzeuger-Methode zurckgeben
double GetAverageRuntime() const;
// Dateiname fr Debug-Ausgaben setzen
void SetDumpFilename( const std::string& sFilename );
// Filter erzeugen
void GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase = false);
void GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase = false) const;
//! Generate FIR filter from third octave (linear) magnitude spectrum
/**
* @param[in] oTOGainMagnitudes Target energetic spectrum
* @param[out] sbImpulseResponse Generated impulse response (time domain filter coefficients) (must be equal or longer than initial filter length of generator class)
* @param[out] bMinimumPhase Request minimum phase
*/
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, ITASampleBuffer& sbImpulseResponse, bool bMinimumPhase /*=false*/ ) const;
private:
double m_dSamplerate; // Abtastrate der Ausgabefilter [Hz]
int m_iFilterLength; // Lnge der Ausgabefilter
int m_iDFTCoeffs; // Anzahl symetrischer DFT-Koeffizienten
float m_fDeltaF; // DFT-Frequenzauflsung [Hz]
int m_iInputFreqs; // Anzahl Frequenzsttzstellen fr Interpolation (31 Terzen + 2 Rnder)
float* m_pfInputFreqs; // Frequenzsttzstellen der Eingabedaten (Terzen) @todo jst: std::vector
float* m_pfInputData; // Zwischenspeicher fr die Interpolation @todo jst: std::vector
float* m_ypp; // Interpolationsdaten
float* m_pfBuf1;
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;
ITAFFT m_ifft;
std::string m_sDumpFilename;
mutable ITAFFT m_ifft;
mutable ITAStopWatch m_sw;
};
......
......@@ -14,20 +14,20 @@ CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const doub
: m_dSamplerate( dSampleRate )
, m_iFilterLength( iFilterLength )
, m_ypp( nullptr )
, m_pfInputFreqs( nullptr )
, m_pfInputData( nullptr )
, m_pfBuf1( nullptr )
, m_pfInputFreqVector( nullptr )
, m_pfInputFreqData( nullptr )
, m_pfFreqDataInterpolated( 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
m_iNumInputFreqs = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
m_pfInputFreqVector = fm_falloc( m_iNumInputFreqs, true );
m_pfInputFreqVector[ 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_pfInputFreqVector[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
m_pfInputFreqVector[ m_iNumInputFreqs - 1 ] = ( float ) dSampleRate / 2; // Right margin: Nyquist frequency
m_pfInputData = fm_falloc( m_iInputFreqs, true );
m_pfInputFreqData = fm_falloc( m_iNumInputFreqs, true );
// DFT frequency bandwidth
m_fDeltaF = ( float ) dSampleRate / ( float ) iFilterLength;
......@@ -35,7 +35,7 @@ CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const doub
// Number of symetric DFT coefficients;
m_iDFTCoeffs = iFilterLength / 2 + 1;
m_pfBuf1 = fm_falloc( 2 * m_iDFTCoeffs, false );
m_pfFreqDataInterpolated = fm_falloc( 2 * m_iDFTCoeffs, false );
m_pfBuf2 = fm_falloc( iFilterLength, false );
m_pfWindow = fm_falloc( iFilterLength, false );
......@@ -46,16 +46,14 @@ CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const doub
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";
m_ifft.plan( ITAFFT::IFFT_C2R, iFilterLength, m_pfFreqDataInterpolated, m_pfBuf2 );
}
CITAThirdOctaveFIRFilterGenerator::~CITAThirdOctaveFIRFilterGenerator()
{
fm_free( m_pfInputFreqs );
fm_free( m_pfInputData );
fm_free( m_pfBuf1 );
fm_free( m_pfInputFreqVector );
fm_free( m_pfInputFreqData );
fm_free( m_pfFreqDataInterpolated );
fm_free( m_pfBuf2 );
fm_free( m_pfWindow );
}
......@@ -81,12 +79,14 @@ double CITAThirdOctaveFIRFilterGenerator::GetAverageRuntime() const
return m_sw.mean();
}
void CITAThirdOctaveFIRFilterGenerator::SetDumpFilename( const std::string& sFilename )
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, ITASampleBuffer& sbImpulseResponse, bool bMinimumPhase /*=false*/ ) const
{
m_sDumpFilename = sFilename;
if( sbImpulseResponse.GetLength() < m_iFilterLength )
ITA_EXCEPT_INVALID_PARAMETER( "Given impulse response is shorter than initial filter size of generator, cannot proceed" );
GenerateFilter( oTOGainMagnitudes, sbImpulseResponse.GetData(), bMinimumPhase );
}
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ )
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ ) const
{
m_sw.start();
......@@ -101,23 +101,26 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0.0f;
pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
if( bMinimumPhase )
pfFilterCoeffs[ 0 ] = 1.0f;
else
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 ] );
m_pfInputFreqData[ 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;
m_pfInputFreqData[ 0 ] = m_pfInputFreqData[ 1 ];
m_pfInputFreqData[ m_iNumInputFreqs - 1 ] = 0.0f;
// Initialize cubic spline interpolation
m_ypp = spline_cubic_set( m_iInputFreqs,
m_pfInputFreqs,
m_pfInputData,
m_ypp = spline_cubic_set( m_iNumInputFreqs,
m_pfInputFreqVector,
m_pfInputFreqData,
1, // Left boundary condition => 1st derivative m=0
0,
1, // Right boundary condition => 1st derivative m=0
......@@ -126,40 +129,41 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
const float fScale = 1 / ( float ) m_iFilterLength;
// No DC offset, ever!
m_pfBuf1[ 0 ] = 0;
m_pfBuf1[ 1 ] = 0;
m_pfFreqDataInterpolated[ 0 ] = 0;
m_pfFreqDataInterpolated[ 1 ] = 0;
if( bMinimumPhase ) {
if( bMinimumPhase )
{
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iInputFreqs,
m_pfInputFreqs,
float x = spline_cubic_val( m_iNumInputFreqs,
m_pfInputFreqVector,
i*m_fDeltaF,
m_pfInputData,
m_pfInputFreqData,
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;
m_pfFreqDataInterpolated[ 2 * i ] = x * fScale; //minimum phase
m_pfFreqDataInterpolated[ 2 * i + 1 ] = 0;
}
}
else {
else
{
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iInputFreqs,
m_pfInputFreqs,
float x = spline_cubic_val( m_iNumInputFreqs,
m_pfInputFreqVector,
i*m_fDeltaF,
m_pfInputData,
m_pfInputFreqData,
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;
m_pfFreqDataInterpolated[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
m_pfFreqDataInterpolated[ 2 * i + 1 ] = 0;
}
}
......
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