Adding third-octave IIR filterbank using Burg coefficient design

CMakeLists.txt

@@ -28,6 +28,7 @@ include_directories( "include" )
 
 # sources
 set( ITADSPHeader
+	"include/ITADSP/ThirdOctaveFilterbank/IIRBurg.h"
 	"include/ITADSPDefinitions.h"
 	"include/ITASIMOVariableDelayLine.h"
 	"include/ITASIMOVariableDelayLineBase.h"
@@ -42,6 +43,7 @@ set( ITADSPHeader
 	"include/ITAThirdOctaveFIRFilterGenerator.h"
 	)
 set( ITADSPSources
+	"src/ITADSP/ThirdOctaveFilterbank/IIRBurg.cpp"
 	"src/ITASIMOVariableDelayLine.cpp"
 	"src/ITASIMOVariableDelayLineBase.cpp"
 	"src/ITABiquad.cpp"

include/ITAThirdOctaveFIRFilterGenerator.h

@@ -52,6 +52,9 @@ public:
 	// Filterlnge zurckgeben [Samples]
 	int GetFilterLength() const;
 
+	//! Returns sampling rate [Samples/s]
+	double GetSampleRate() const;
+
 	// Latenz in den Filtern zurckgeben [Samples]
 	int GetLatency() const;
 

include/ITAThirdOctaveFilterbank.h

@@ -40,7 +40,9 @@ public:
 	enum
 	{
 		FIR_SPLINE_LINEAR_PHASE = 0, //!< Linearphasiges FIR Filter mit Spline-Interpolation
-		IIR_BIQUADS_ORDER10			 //!< Rekursives IIR Filter umgesetzt durch kaskadierte Biquads
+		IIR_BIQUADS_ORDER10,			 //!< Rekursives IIR Filter umgesetzt durch kaskadierte Biquads
+		IIR_BURG_ORDER4, //!< Recursive 4th-order IIR filtering with a Burg coefficient design
+		IIR_BURG_ORDER10, //!< Recursive 10th-order IIR filtering with a Burg coefficient design
 	} FilterbankRealisationMethods;
 
 	//! Factory method zum erzeugen einer Filterbank

src/ITAThirdOctaveFIRFilterGenerator.cpp

@@ -64,6 +64,12 @@ int CITAThirdOctaveFIRFilterGenerator::GetFilterLength() const
 {
 	return m_iFilterLength;
 }
+
+double CITAThirdOctaveFIRFilterGenerator::GetSampleRate() const
+{
+	return m_dSamplerate;
+}
+
 int CITAThirdOctaveFIRFilterGenerator::GetLatency() const
 {
 	// Latency = Half DFT period (ceil)

src/ITAThirdOctaveFilterbank.cpp

@@ -2,6 +2,7 @@
 
 #include <ITAThirdOctaveFilterbankFIR.h>
 #include <ITAThirdOctaveFilterbankIIR.h>
+#include <ITADSP/ThirdOctaveFilterbank/IIRBurg.h>
 
 #include <cassert>
 
@@ -13,6 +14,10 @@ CITAThirdOctaveFilterbank* CITAThirdOctaveFilterbank::Create( const double dSamp
 		return new CITAThirdOctaveFilterbankFIR( dSampleRate, iBlockLength );
 	case IIR_BIQUADS_ORDER10:
 		return new CITAThirdOctaveFilterbankIIR( dSampleRate, iBlockLength );
+	case IIR_BURG_ORDER4:
+		return new ITADSP::ThirdOctaveFilterbank::CIIRBurg( dSampleRate, iBlockLength, 4 );
+	case IIR_BURG_ORDER10:
+		return new ITADSP::ThirdOctaveFilterbank::CIIRBurg( dSampleRate, iBlockLength, 10 );
 	default:
 		assert( false );
 		return nullptr;

tests/ITADSPThirdOctaveFilterbankTest.cpp

@@ -11,17 +11,19 @@ using namespace std;
 const double g_dSampleRate = 44100;
 const int g_iFilterLength = 512;
 
-void TestThirdOctaveFilterbankIIR();
+void TestThirdOctaveFilterbankIIRBiquad();
+void TestThirdOctaveFilterbankIIRBurg();
 void TestThirdOctaveFilterbankFIR();
 
 int main( int, char** )
 {
-	TestThirdOctaveFilterbankIIR();
+	TestThirdOctaveFilterbankIIRBiquad();
+	TestThirdOctaveFilterbankIIRBurg();
 	TestThirdOctaveFilterbankFIR();
 	return 255;
 }
 
-void TestThirdOctaveFilterbankIIR()
+void TestThirdOctaveFilterbankIIRBiquad()
 {
 	const int iSampleLength = ( 1 << 17 );
 	CITAThirdOctaveFilterbank* pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
@@ -43,7 +45,7 @@ void TestThirdOctaveFilterbankIIR()
 	cout << "Runtime identity magnitudes:" << sw.ToString() << endl;
 
 	y.Normalize();
-	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIR_Identity.wav", &y, g_dSampleRate );
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_Identity.wav", &y, g_dSampleRate );
 	y.Zero();
 
 	pIIRFilterbank->Clear();
@@ -57,7 +59,7 @@ void TestThirdOctaveFilterbankIIR()
 	cout << "Runtime zero magnitude:" << sw.ToString() << endl;
 
 	y.Normalize();
-	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIR_Zeros.wav", &y, g_dSampleRate );
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.wav", &y, g_dSampleRate );
 
 	pIIRFilterbank->Clear();
 
@@ -77,7 +79,68 @@ void TestThirdOctaveFilterbankIIR()
 	cout << "Runtime real magnitudes:" << sw.ToString() << endl;
 
 	y.Normalize();
-	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIR_SomeBands.wav", &y, g_dSampleRate );
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_SomeBands.wav", &y, g_dSampleRate );
+
+	delete pIIRFilterbank;
+}
+
+void TestThirdOctaveFilterbankIIRBurg()
+{
+	const int iSampleLength = ( 1 << 17 );
+	CITAThirdOctaveFilterbank* pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BURG_ORDER4 );
+
+	ITASampleBuffer x( iSampleLength, true );
+	x[ 0 ] = 1.0f;
+
+	ITABase::CThirdOctaveGainMagnitudeSpectrum oMags;
+	oMags.SetIdentity();
+
+	pIIRFilterbank->SetMagnitudes( oMags, false );
+
+	ITASampleBuffer y( iSampleLength );
+
+	ITAStopWatch sw;
+	sw.start();
+	pIIRFilterbank->Process( x.GetData(), y.GetData() );
+	sw.stop();
+	cout << "Runtime identity magnitudes:" << sw.ToString() << endl;
+
+	y.Normalize();
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_Identity.wav", &y, g_dSampleRate );
+	y.Zero();
+
+	pIIRFilterbank->Clear();
+
+	oMags.SetZero();
+	pIIRFilterbank->SetMagnitudes( oMags, false );
+
+	sw.start();
+	pIIRFilterbank->Process( x.GetData(), y.GetData() );
+	sw.stop();
+	cout << "Runtime zero magnitude:" << sw.ToString() << endl;
+
+	y.Normalize();
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_Zeros.wav", &y, g_dSampleRate );
+
+	pIIRFilterbank->Clear();
+
+	// dB values
+	oMags[ 3 ] = 3.0f;
+	oMags[ 10 ] = -3.0f;
+	oMags[ 11 ] = -6.0f;
+	oMags[ 12 ] = -9.0f;
+	oMags[ 13 ] = -9.0f;
+	oMags[ 14 ] = -7.0f;
+
+	pIIRFilterbank->SetMagnitudes( oMags );
+
+	sw.start();
+	pIIRFilterbank->Process( x.GetData(), y.GetData() );
+	sw.stop();
+	cout << "Runtime real magnitudes:" << sw.ToString() << endl;
+
+	y.Normalize();
+	writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_SomeBands.wav", &y, g_dSampleRate );
 
 	delete pIIRFilterbank;
 }

tests/ITADSPThirdOctaveFilterbankTest.m

@@ -3,7 +3,7 @@
 % Loads the results of filterbanks from the test ITADSPThirdOctaveFilterbankTest
 % Requires ITA Toolbox, see http://ita-toolbox.org
 
-iir_ident = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_Identity.wav' );
+iir_ident = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquad_Identity.wav' );
 %iir_ident.signalType = 'energy';
 iir_ident.channelNames = { 'Identity' };
 iir_ident.pf
@@ -13,7 +13,7 @@ test_iir_ident2 = ita_normalize_spk( ita_merge( ita_demosound, ita_convolve( ita
 test_iir_ident2.play
 
 
-iir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_Zeros.wav' );
+iir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.wav' );
 iir_zeros.channelNames = { 'Zeros' };
 iir_zeros.pf
 % 
@@ -40,14 +40,14 @@ iir_zeros.pf
 % 
 
 
-ITADSPThirdOctaveFilterbankTest_IIR_Zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_Zeros.wav' );
-ITADSPThirdOctaveFilterbankTest_IIR_SomeBands = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_SomeBands.wav' );
+ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.wav' );
+ITADSPThirdOctaveFilterbankTest_IIRBiquad_SomeBands = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquad_SomeBands.wav' );
 
 
-ITADSPThirdOctaveFilterbankTest_IIR_Zeros.channelNames = { 'Zero' };
-ITADSPThirdOctaveFilterbankTest_IIR_Zeros.pf
+ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.channelNames = { 'Zero' };
+ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.pf
 
-ITADSPThirdOctaveFilterbankTest_IIR_SomeBands.pf
+ITADSPThirdOctaveFilterbankTest_IIRBiquad_SomeBands.pf
 
 %% FIR
 fir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIR_Zeros.wav' );
@@ -66,3 +66,18 @@ fir_all = ita_merge( fir_zeros, fir_unity, fir_sb );
 
 fir_all.pf
 
+
+%% IIR Burg
+iir_burg_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_Zeros.wav' );
+iir_burg_unity = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_Identity.wav' );
+iir_burg_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_SomeBands.wav' );
+
+iir_burg_sb_unity = iir_burg_sb.ch( 1 );
+for n = 2:iir_burg_sb.nChannels
+    iir_burg_sb_unity = iir_burg_sb_unity + iir_burg_sb.ch( n );
+end
+
+iir_burg_all = ita_merge( iir_burg_zeros, iir_burg_unity, iir_burg_sb );
+
+iir_burg_all.pf
+