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

Adding third octave filter banks (FIR and IIR implementation)

parent 909df859
......@@ -8,6 +8,9 @@ include( VistaCommon )
# dependencies
vista_use_package( ITABase REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAFFT REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAConvolution REQUIRED FIND_DEPENDENCIES ) # required for FIR filtering (uses block convolver)
vista_use_package( TBB REQUIRED FIND_DEPENDENCIES )
vista_use_package( SPLINE REQUIRED FIND_DEPENDENCIES )
# includes
include_directories( "include" )
......@@ -18,10 +21,18 @@ set( ITADSPHeader
"include/ITADSPDefinitions.h"
"include/ITASIMOVariableDelayLine.h"
"include/ITABiquad.h"
"include/ITAThirdOctaveFilterbank.h"
"include/ITAThirdOctaveFilterbankFIR.h"
"include/ITAThirdOctaveFilterbankIIR.h"
"include/ITAThirdOctaveFIRFilterGenerator.h"
)
set( ITADSPSources
"src/ITASIMOVariableDelayLine.cpp"
"src/ITABiquad.cpp"
"src/ITAThirdOctaveFilterbank.cpp"
"src/ITAThirdOctaveFilterbankIIR.cpp"
"src/ITAThirdOctaveFIRFilterGenerator.cpp"
"src/ITAThirdOctaveFilterbankCoefficients.h"
)
......
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2017
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_THIRD_OCTAVE_FILTER_GENERATOR
#define IW_ITA_THIRD_OCTAVE_FILTER_GENERATOR
#include <ITADSPDefinitions.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <ITAFFT.h>
#include <ITAStopWatch.h>
#include <string.h>
#include <vector>
/**
* Diese Klasse erzeugt FIR-Filter (Impulsantworten) aus Terzband-Betragsspektren.
* Hierzu wird ein linear-phasiges FIR-Filter mittels Spline-Interpolation erzeugt.
* Die Filter erzeugen eine Latenz der halben Filterlnge.
*/
class ITA_DSP_API CITAThirdOctaveFIRFilterGenerator
{
public:
// Design methods
enum
{
LINEAR_PHASE = 0,
MINIMUM_PHASE
};
// Konstruktor
CITAThirdOctaveFIRFilterGenerator( const double dSampleRate, const int iFilterLength );
// Destruktor
~CITAThirdOctaveFIRFilterGenerator();
// Filterlnge zurckgeben [Samples]
int GetFilterLength() const;
// Latenz in den Filtern zurckgeben [Samples]
int GetLatency() const;
// 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 CITAThirdOctaveMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs );
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;
float* m_pfBuf2;
bool m_bWindow; // Apply window?
float* m_pfWindow;
ITAFFT m_ifft;
std::string m_sDumpFilename;
mutable ITAStopWatch m_sw;
};
#endif // IW_ITA_THIRD_OCTAVE_FILTER_GENERATOR
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2017
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_THIRD_OCTAVE_FILTERBANK
#define IW_ITA_THIRD_OCTAVE_FILTERBANK
#include <ITAThirdOctaveMagnitudeSpectrum.h>
//! Terzbandspektrum-Filterbank zur digitalen Filterung
/**
* Diese Oberklasse für Filterbänke, die zur Realisierung der Filterung von Terzbandspektren dient,
* definiert die Schnittstellen zur Nutzung einer solchen Filterbank.
*
* Sie wird durch die Realisierungsmethoden (#FilterbankRealisationMethods) mittels Factory Method
* erstellt und kann dann Eingabesamples im Process()-Schritt filtern.
*
* Die Datensätze der Terzbandspektren werden durch die Klasse \CVAThirdOctaveMagnitues verwaltet.
*
*/
class CITAThirdOctaveFilterbank
{
public:
//! Realisierungsmethoden
enum
{
FIR_SPLINE_LINEAR_PHASE = 0, //!< Linearphasiges FIR Filter mit Spline-Interpolation
IIR_BIQUADS_ORDER10 //!< Rekursives IIR Filter umgesetzt durch kaskadierte Biquads
} FilterbankRealisationMethods;
//! Factory method zum erzeugen einer Filterbank
/**
* \param dSamplerate Samplingrate
* \param iBlocklength Blockgröße
* \param iMethod Realisierungsmethode, eines von #FilterbankRealisationMethods (default: linearphasiges FIR Filter)
*
* \return Zeiger auf die erzeugte Terzfilterbank
*/
static CITAThirdOctaveFilterbank* Create( const double dSampleRate, const int iBlockLength, const int iMethod = FIR_SPLINE_LINEAR_PHASE );
//! Destruktor
virtual inline ~CITAThirdOctaveFilterbank() {};
//! Idealer Übertrager setzen
/**
* \param bSmoothChangeover Überbenden (default, true) oder direktes Umschalten (false)
*/
virtual inline void SetIdentity( const bool bSmoothChangeover = true )
{
CITAThirdOctaveMagnitudeSpectrum oIdentity;
SetGains( oIdentity, bSmoothChangeover );
};
//! Verstärkungsfaktoren setzen
/**
* \param oGains Neue Verstärkungsfaktoren
* \param bSmoothChangeover Überbenden (default, true) oder direktes Umschalten (false)
*/
virtual void SetGains( const CITAThirdOctaveMagnitudeSpectrum& oGains, const bool bSmoothChangeover = true ) = 0;
//! Latenz (Verzögerung) der Filterbank zurückgeben
/**
* \return Latenz (Verzögerung) in ganzzahligen Sample
*/
int GetLatency() const;
//! Löscht alle internen Puffer
virtual void Clear() = 0;
//! Filtert einen Block Samples (muss die angegebene Blocklänge haben, s.o.)
/**
* \pfInputSamples Eingabesamples (Anzahl = Blocklänge)
* \pfOutputSamples Ausgabesamples (Anzahl = Blocklänge)
*
* @todo jst: check if this couln't be decoupled from block-bound processing (like ITABiquad)
*/
virtual void Process( const float* pfInputSamples, float* pfOutputSamples ) = 0;
protected:
//! Standardkonstruktor (Deaktiviert, FactoryMethod Create() benutzen)
inline CITAThirdOctaveFilterbank() {};
};
#endif // IW_ITA_THIRD_OCTAVE_FILTERBANK
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2017
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_THIRD_OCTAVE_FILTERBANK_FIR
#define IW_ITA_THIRD_OCTAVE_FILTERBANK_FIR
#include <ITAThirdOctaveFilterbank.h>
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <ITAUPConvolution.h>
#include <ITAUPFilter.h>
#include <ITAFastMath.h>
#include <cassert>
//! Terzfilterbank mittels Dynamic Single-Channel Blockfalter (FIR Filter)
class ITA_DSP_API CITAThirdOctaveFilterbankFIR : public CITAThirdOctaveFilterbank
{
public:
inline CITAThirdOctaveFilterbankFIR( const double dSampleRate, const int iBlockLength )
: m_pfFilter( nullptr )
, m_pGenerator( nullptr )
, m_pConvolver( nullptr )
{
m_iBlocklength = iBlockLength;
// [fwe] Hier wird die Filterlänge für Directivities festgelegt
m_iFilterLength = 128;
m_pfFilter = fm_falloc( m_iFilterLength, true );
m_pGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, m_iFilterLength );
m_pConvolver = new ITAUPConvolution( iBlockLength, m_iFilterLength );
m_pConvolver->SetFilterExchangeMode( ITAUPConvolution::CROSSFADE_COSINE_SQUARE );
m_pConvolver->SetFilterCrossfadeLength( 32 );
SetIdentity( false );
}
inline virtual ~CITAThirdOctaveFilterbankFIR()
{
fm_free( m_pfFilter );
delete m_pGenerator;
delete m_pConvolver;
}
inline virtual void SetIdentity( const bool bSmoothChangeover )
{
ITAUPFilter* pFilter = m_pConvolver->RequestFilter();
int iLatency = m_pGenerator->GetLatency();
assert( iLatency < m_iFilterLength );
fm_zero( m_pfFilter, m_iFilterLength );
m_pfFilter[ iLatency ] = 1;
pFilter->Load( m_pfFilter, m_iFilterLength );
m_pConvolver->ExchangeFilter( pFilter, ( bSmoothChangeover ? ITAUPConvolution::AUTO : ITAUPConvolution::SWITCH ) );
pFilter->Release(); // Auto-release
}
inline virtual void SetGains( const CITAThirdOctaveMagnitudeSpectrum& oGains, const bool bSmoothChangeover = true )
{
m_pGenerator->GenerateFilter( oGains, m_pfFilter );
ITAUPFilter* pFilter = m_pConvolver->RequestFilter();
pFilter->Load( m_pfFilter, m_iFilterLength );
m_pConvolver->ExchangeFilter( pFilter, ( bSmoothChangeover ? ITAUPConvolution::AUTO : ITAUPConvolution::SWITCH ) );
pFilter->Release(); // Auto-release
}
inline int GetLatency() const
{
return m_pGenerator->GetLatency();
}
inline virtual void Clear()
{
m_pConvolver->clear();
SetIdentity( true );
}
inline virtual void Process( const float* pfInputSamples, float* pfOutputSamples )
{
m_pConvolver->Process( pfInputSamples, m_iBlocklength, pfOutputSamples, m_iBlocklength );
}
private:
int m_iBlocklength;
int m_iFilterLength;
float* m_pfFilter;
CITAThirdOctaveFIRFilterGenerator* m_pGenerator;
ITAUPConvolution* m_pConvolver;
};
#endif // IW_ITA_THIRD_OCTAVE_FILTERBANK_FIR
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2017
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_THIRD_OCTAVE_FILTERBANK_IIR
#define IW_ITA_THIRD_OCTAVE_FILTERBANK_IIR
#include <ITADSPDefinitions.h>
#include <ITAThirdOctaveFilterbank.h>
#include <ITABiquad.h>
#include <vector>
#include <tbb/concurrent_queue.h>
//! Terzfilterbank mittels Biquads (IIR Filter)
/**
* Diese Klasse realisiert eine Terzfilterbank (#CVAThirdOctaveFilterbank) mit der Methode der kaskadierten
* Biquads (#CVABiquad) für Verstärkungsfaktoren eines Terzbank-Betragsspektrums (#CVAThirdOctaveMagnitudes).
*
*/
class ITA_DSP_API CITAThirdOctaveFilterbankIIR : public CITAThirdOctaveFilterbank
{
public:
//! Konstruktor mit Samplerate und Blocklänge
/**
* \param dSamplerate Samplingrate
* \param iBlocklength Blocklänge
*/
CITAThirdOctaveFilterbankIIR( const double dSampleRate, const int iBlockLength );
//! Destruktor
virtual ~CITAThirdOctaveFilterbankIIR();
//! Filterlatenz in Samples zurückgeben
int GetLatency() const;
//! Verstärkungen (Gains) setzen
/**
* \param oGains Verstärkungsfaktoren
* \param bSmoothChangeover Wenn true, dann überblenden (default), sonst sofort internen Gain umschalten
*/
void SetGains( const CITAThirdOctaveMagnitudeSpectrum& oGains, const bool bSmoothChangeover = true );
//! Alle internen Zustände zurücksetzen (Akkumulatoren der Biquads)
void Clear();
//! Verarbeite Samples (Filtern)
/**
* \param pfInputSamples Eingabesamples (Anzahl = Blocklänge)
* \param pfOutputSamples Ausgabesamples (Anzahl = Blocklänge)
*/
void Process( const float* pfInputSamples, float* pfOutputSamples );
private:
//! Interne Datenklasse für das Verarbeiten eines neuen Gain Datensatzes
class GainUpdate
{
public:
CITAThirdOctaveMagnitudeSpectrum oGains; //!< Gain-Daten
int iBlendSamples; //!< Anzahl Samples zum Überblenden
};
double m_dSampleRate; //!< Samplingrate
int m_iBlockLength; //!< Blocklänge
int m_nBandsInternal; //!< Anzahl der internen Bänder
int m_nBiquadsPerBand; //!< Anzahl von Biqads pro Band
std::vector< CITABiquad > m_vBiquads; //!< Biquads, Zugriff: [Band][BiquadNummer]
tbb::strict_ppl::concurrent_queue< GainUpdate > m_qGains; //!< Liste von neuen Verstärkungsfaktoren
CITAThirdOctaveMagnitudeSpectrum m_oGainsInternal; //!< Interne Verstärkungsfaktoren
float* m_pfTempFilterBuf; //!< Zwischenpuffer für Filter
float* m_pfTempOutputBuf; //!< Zwischenpuffer für Ausgabe
};
#endif // IW_ITA_THIRD_OCTAVE_FILTERBANK_IIR
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <ITAConstants.h>
#include <ITAFastMath.h>
#include <ITANumericUtils.h>
#include <ITAStringUtils.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <spline.h>
// Debug:
#include <iomanip>
#include <iostream>
#include <fstream>
// Debugging flags
#define DUMP_INTERPOLATED_MAGNITUDES 0
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 = CITAThirdOctaveMagnitudeSpectrum::GetNumBands() + 2;
m_pfInputFreqs = fm_falloc( m_iInputFreqs, true );
m_pfInputFreqs[ 0 ] = 0; // Left margin
for( int i = 0; i < CITAThirdOctaveMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputFreqs[ i + 1 ] = CITAThirdOctaveMagnitudeSpectrum::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 CITAThirdOctaveMagnitudeSpectrum& oTOMagnitudes, float* pfFilterCoeffs )
{
m_sw.start();
// [fwe] Intermediate bugfix. Wrong filter for identity
const float EPSILON = 0.001F;
bool bZero = true;
bool bIdent = true;
for( int i = 0; i<CITAThirdOctaveMagnitudeSpectrum::GetNumBands(); i++ )
{
bZero &= ( oTOMagnitudes[ i ] > -EPSILON ) && ( oTOMagnitudes[ i ] < EPSILON );
bIdent &= ( oTOMagnitudes[ i ] > 1 - EPSILON ) && ( oTOMagnitudes[ i ] < 1 + EPSILON );
}
if( bZero )
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0;
return;
}
if( bIdent )
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0;
pfFilterCoeffs[ m_iFilterLength / 2 ] = 1;
return;
}
// 1st step: Interpolate the magnitudes
m_pfInputData[ 0 ] = 1;
memcpy( &m_pfInputData[ 1 ], &oTOMagnitudes.GetMagnitudes(), CITAThirdOctaveMagnitudeSpectrum::GetNumBands() * sizeof( float ) );
m_pfInputData[ m_iInputFreqs - 1 ] = 0;
// 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;
float fScale = 1 / ( float ) m_iFilterLength;
// No DC offset, ever!
m_pfBuf1[ 0 ] = 0;
m_pfBuf1[ 1 ] = 0;
#if (DUMP_INTERPOLATED_MAGNITUDES==1)
std::ofstream coeffs_file;
coeffs_file.open(m_sDumpFilename);
coeffs_file << "# 1/3 magnitudes: " << FloatArrayToString(pfThirdOctaveMagnitudes, CVAThirdOctaveMagnitudes::nBands, 3) << std::endl;
coeffs_file << "# Num DFT coeffs: " << m_iDFTCoeffs << std::endl;
coeffs_file << "# DFT bin bandwidth: " << m_fDeltaF << " Hz" << std::endl << std::endl;
#endif
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 );
#if (DUMP_INTERPOLATED_MAGNITUDES==1)
//std::setprecision(12) << std::ios::fixed <<
coeffs_file << (i*m_fDeltaF) << "\t"
<< x << std::endl;
#endif
// 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;
}
#if (DUMP_INTERPOLATED_MAGNITUDES==1)
coeffs_file.close();
#endif
// DEBUG:
//std::cout << FloatArrayToString(m_pfBuf1, m_iDFTCoeffs*2, 3) << std::endl;
// 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 ];
}
//writeAudiofile()
// TODO: Minimum-phase?
m_sw.stop();
}
#include <ITAThirdOctaveFilterbank.h>
#include <ITAThirdOctaveFilterbankFIR.h>
#include <ITAThirdOctaveFilterbankIIR.h>
#include <cassert>
CITAThirdOctaveFilterbank* CITAThirdOctaveFilterbank::Create( const double dSampleRate, const int iBlockLength, const int iMethod )
{
switch( iMethod )