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

WIP, iir filtering looking good now

parent b23d8e09
......@@ -39,17 +39,17 @@ namespace ITADSP
{
ITA_DSP_API void Yulewalk( const ITABase::CFiniteImpulseResponse& oIR, CFilterCoefficients& oCoeffs );
ITA_DSP_API void Burg(const ITABase::CFiniteImpulseResponse& oIR, CFilterCoefficients& oCoeffs);
ITA_DSP_API void Burg( const ITABase::CFiniteImpulseResponse& oIR, CFilterCoefficients& oCoeffs );
ITA_DSP_API float InnerProduct(const float *x, const float *y, const int length);
ITA_DSP_API float InnerProduct( const float *x, const float *y, const int length );
ITA_DSP_API void numf();
ITA_DSP_API void denf( float *A, const ITASampleBuffer& R, int na );
ITA_DSP_API void toeplitz(float **out, const ITASampleBuffer &row, const ITASampleBuffer &column);
ITA_DSP_API void toeplitz( float **out, const ITASampleBuffer &row, const ITASampleBuffer &column );
ITA_DSP_API void toeplitz(float **out, const ITASampleBuffer &in, int na);
ITA_DSP_API void toeplitz( float **out, const ITASampleBuffer &in, int na );
}
}
......
......@@ -45,8 +45,9 @@ public:
private:
int m_iBlocklength;
int m_iFilterLength;
float* m_pfFilter;
int m_iConvolutionFilterLength;
int m_iGeneratorFilterLength;
std::vector< float > m_vfFilter;
CITAThirdOctaveFIRFilterGenerator* m_pGenerator;
ITAUPConvolution* m_pConvolver;
};
......
......@@ -4,6 +4,7 @@
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <algorithm>
//#include "ITAIIRUtils.h" // rm me
using namespace ITADSP::ThirdOctaveFilterbank;
......@@ -11,7 +12,8 @@ CIIRBurg::CIIRBurg( const double dSampleRate, const int iBlockLength, const int
: m_iBlockLength( iBlockLength )
{
m_pFilterEngine = new CITAIIRFilterEngine( iFilterOrder );
m_pFilterGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, (std::max)( 2 * iFilterOrder, 32 ) );
int iFIRGenLength = int( 4 * ceil( dSampleRate / ITABase::CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
m_pFilterGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, iFIRGenLength );
}
CIIRBurg::~CIIRBurg()
......@@ -22,12 +24,16 @@ CIIRBurg::~CIIRBurg()
void CIIRBurg::SetMagnitudes( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oGains, const bool )
{
ITABase::CFiniteImpulseResponse oIR( m_pFilterGenerator->GetFilterLength(), (float) m_pFilterGenerator->GetSampleRate(), true );
m_pFilterGenerator->GenerateFilter( oGains, oIR.GetData(), false );
ITABase::CFiniteImpulseResponse oIR( m_pFilterGenerator->GetFilterLength(), ( float ) m_pFilterGenerator->GetSampleRate(), true );
m_pFilterGenerator->GenerateFilter( oGains, oIR.GetData(), true );
//oIR.Store( "iirburg_fb_ir_input.wav" );
CFilterCoefficients oCoeffs( m_pFilterEngine->GetOrder() );
ITADSP::IIRFilterGenerator::Burg( oIR, oCoeffs );
//ITADSP::ExportIIRCoefficientsToJSON( "iirburg_fb_coeffs.json", oCoeffs );
m_pFilterEngine->SetCoefficients( oCoeffs );
}
......
......@@ -19,19 +19,21 @@
#include "ITAThirdOctaveFilterbankFIR.h"
#include <ITAUPConvolution.h>
#include <algorithm>
CITAThirdOctaveFilterbankFIR::CITAThirdOctaveFilterbankFIR( const double dSampleRate, const int iBlockLength )
: m_pfFilter( nullptr )
, m_pGenerator( nullptr )
: m_pGenerator( nullptr )
, m_pConvolver( nullptr )
{
m_iBlocklength = iBlockLength;
// [fwe] Hier wird die Filterlnge fr Directivities festgelegt
m_iFilterLength = 128;
m_pfFilter = fm_falloc( m_iFilterLength, true );
// [fwe] Hier wird die Filterlnge fr Directivities festgelegt, jst: dont take too short lengths as the spline interp requires zero DC offset and will cut low freq amplitudes
m_iGeneratorFilterLength = int( 4 * ceil( dSampleRate / ITABase::CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
m_vfFilter.resize( m_iGeneratorFilterLength );
m_pGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, m_iFilterLength );
m_pGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, m_iGeneratorFilterLength );
m_pConvolver = new ITAUPConvolution( iBlockLength, m_iFilterLength );
m_iConvolutionFilterLength = ( std::max )( 128, iBlockLength );
m_pConvolver = new ITAUPConvolution( iBlockLength, m_iConvolutionFilterLength );
m_pConvolver->SetFilterExchangeFadingFunction( ITABase::FadingFunction::COSINE_SQUARE );
m_pConvolver->SetFilterCrossfadeLength( 32 );
......@@ -40,7 +42,6 @@ CITAThirdOctaveFilterbankFIR::CITAThirdOctaveFilterbankFIR( const double dSample
CITAThirdOctaveFilterbankFIR::~CITAThirdOctaveFilterbankFIR()
{
fm_free( m_pfFilter );
delete m_pGenerator;
delete m_pConvolver;
}
......@@ -49,19 +50,22 @@ void CITAThirdOctaveFilterbankFIR::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 );
assert( iLatency < m_iGeneratorFilterLength );
fm_zero( &m_vfFilter[ 0 ], m_vfFilter.size() );
m_vfFilter[ iLatency ] = 1;
pFilter->Load( &m_vfFilter[ 0 ], ( std::min )( m_iGeneratorFilterLength, m_iConvolutionFilterLength ) );
m_pConvolver->ExchangeFilter( pFilter, ( bSmoothChangeover ? ITABase::FadingFunction::COSINE_SQUARE : ITABase::FadingFunction::SWITCH ) );
pFilter->Release(); // Auto-release
}
void CITAThirdOctaveFilterbankFIR::SetMagnitudes( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oMags, const bool bSmoothChangeover /*= true*/ )
{
m_pGenerator->GenerateFilter( oMags, m_pfFilter );
m_pGenerator->GenerateFilter( oMags, &m_vfFilter[0], false );
ITAUPFilter* pFilter = m_pConvolver->RequestFilter();
pFilter->Load( m_pfFilter, m_iFilterLength );
int iRange = ( std::min )( m_iGeneratorFilterLength, m_iConvolutionFilterLength );
int iShift = m_pGenerator->GetLatency() - iRange / 2;
assert( iShift > 0 && iShift + iRange < m_iGeneratorFilterLength );
pFilter->Load( &m_vfFilter[ iShift ], iRange );
m_pConvolver->ExchangeFilter( pFilter, ( bSmoothChangeover ? ITABase::FadingFunction::COSINE_SQUARE : ITABase::FadingFunction::SWITCH ) );
pFilter->Release(); // Auto-release
}
......
......@@ -5,6 +5,7 @@
#include <ITAStopWatch.h>
#include <ITAStringUtils.h>
#include <ITASampleBuffer.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#ifdef WITH_JSON_SUPPORT
#include <ITABase/UtilsJSON.h>
......@@ -13,20 +14,22 @@
#include <iostream>
using namespace std;
using namespace ITABase;
const int iSampleLength = ( 1 << 17 );
const int g_iSampleLength = ( 1 << 17 );
const int g_iBlockLength = 128;
const double g_dSampleRate = 44100;
const int g_iFilterLength = 512;
const int g_iFIRGenLength = int( 4 * ceil( g_dSampleRate / CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
ITAStopWatch sw;
ITASampleBuffer sbImpulseResponse( iSampleLength );
ITASampleBuffer sbDirac( iSampleLength, true );
ITASampleBuffer sbImpulseResponse( g_iSampleLength );
ITASampleBuffer sbDirac( g_iSampleLength, true );
ITABase::CThirdOctaveGainMagnitudeSpectrum oZero;
ITABase::CThirdOctaveGainMagnitudeSpectrum oIdentity;
ITABase::CThirdOctaveGainMagnitudeSpectrum oSomeBands;
ITABase::CThirdOctaveGainMagnitudeSpectrum oHighPass;
ITABase::CThirdOctaveGainMagnitudeSpectrum oDirectivity;
ITABase::CThirdOctaveGainMagnitudeSpectrum oAbsorption;
ITABase::CThirdOctaveGainMagnitudeSpectrum oReflection;
ITABase::CThirdOctaveGainMagnitudeSpectrum oDiffraction;
void Run( CITAThirdOctaveFilterbank* pFilterBank, std::string sExportName );
......@@ -67,10 +70,10 @@ int main( int, char** )
#endif
vector< float > vfAbsorptionSpectrumTO = { 0.02f, 0.02f, 0.02f, 0.0233333f, 0.0266667f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.0333333f, 0.0366667f, 0.04f, 0.0433333f, 0.0466667f, 0.05f, 0.0566667f, 0.0633333f, 0.07f, 0.07f, 0.07f, 0.07f, 0.0733333f, 0.0766667f, 0.08f, 0.08f };
for( int i = 0; i < oAbsorption.GetNumBands(); i++ )
oAbsorption.SetMagnitude( i, 1 - vfAbsorptionSpectrumTO[ i ] * vfAbsorptionSpectrumTO[ i ] );
for( int i = 0; i < oReflection.GetNumBands(); i++ )
oReflection.SetMagnitude( i, 1 - vfAbsorptionSpectrumTO[ i ] * vfAbsorptionSpectrumTO[ i ] );
#ifdef WITH_JSON_SUPPORT
ITABase::Utils::JSON::Export( &oAbsorption, "ReflectedSpectrum.json" );
ITABase::Utils::JSON::Export( &oReflection, "ReflectedSpectrum.json" );
#endif
vector< float > vfTrumpetDirectivitySpectrumTO = { 0.860280990600586f, 0.860280632972717f, 0.860280632972717f, 0.860280632972717f, 0.860280632972717f, 0.860280632972717f, 0.860280632972717f, 0.860280692577362f, 0.860280692577362f, 0.863133668899536f, 0.752475023269653f, 0.674751520156860f, 0.499466955661774f, 0.738280415534973f, 0.653443872928619f, 0.507191777229309f, 0.533296763896942f, 0.503476321697235f, 0.376767426729202f, 0.353374809026718f, 0.269741356372833f, 0.207140043377876f, 0.153062343597412f, 0.112099960446358f, 0.127615734934807f, 0.0946486070752144f, 0.0785422623157501f, 0.0600289255380631f, 0.0488252453505993f, 0.0387985333800316f, 0.0315645076334477f };
......@@ -78,27 +81,29 @@ int main( int, char** )
#ifdef WITH_JSON_SUPPORT
ITABase::Utils::JSON::Export( &oDirectivity, "DirectivitySpectrum.json" );
#endif
vector< float > vfDiffractionSpectrumTO = { 0.111934553579764f, 0.106460935614618f, 0.100699704548526f, 0.0947001721104573f, 0.0890866491204974f, 0.0835569704783756f, 0.0775507540050590f, 0.0722382340567499f, 0.0670519961369095f, 0.0621838133597414f, 0.0565935658858452f, 0.0519303873398500f, 0.0473436952392266f, 0.0428843783633814f, 0.0389753239378729f, 0.0351914143864374f, 0.0315689677042525f, 0.0284535398945892f, 0.0255957249048243f, 0.0227266136854582f, 0.0203857703689526f, 0.0182702796747406f, 0.0162996829305633f, 0.0144787461518097f, 0.0129580533457382f, 0.0115033291355433f, 0.0102513560003608f, 0.00917063160436257f, 0.00820336236165088f, 0.00725137721071123f, 0.00648611579946720f };
oDiffraction.SetMagnitudes( vfDiffractionSpectrumTO );
#ifdef WITH_JSON_SUPPORT
ITABase::Utils::JSON::Export( &oDiffraction, "DiffractionSpectrum.json" );
#endif
auto pIIRFilterbankBiquadOrder10 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
auto pIIRFilterBurgOrder10 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, g_iBlockLength, CITAThirdOctaveFilterbank::IIR_BURG_ORDER10 );
Run( pIIRFilterBurgOrder10, "IIRBurgOrder10" );
delete pIIRFilterBurgOrder10;
auto pIIRBurgOrder4 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, g_iBlockLength, CITAThirdOctaveFilterbank::IIR_BURG_ORDER4 );
Run( pIIRBurgOrder4, "IIRBurgOrder4" );
delete pIIRBurgOrder4;
auto pIIRFilterbankBiquadOrder10 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, g_iBlockLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
Run( pIIRFilterbankBiquadOrder10, "IIRBiquadOrder10" );
delete pIIRFilterbankBiquadOrder10;
auto pFIRConv = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::FIR_SPLINE_LINEAR_PHASE );
auto pFIRConv = CITAThirdOctaveFilterbank::Create( g_dSampleRate, g_iBlockLength, CITAThirdOctaveFilterbank::FIR_SPLINE_LINEAR_PHASE );
Run( pFIRConv, "FIRSplineLinearPhase" );
delete pFIRConv;
auto pIIRFilterBurgOrder10 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BURG_ORDER10 );
Run( pIIRFilterBurgOrder10, "IIRBurgOrder10" );
delete pIIRFilterBurgOrder10;
//auto pIIRBurgOrder4 = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BURG_ORDER10 );
//Run( pIIRBurgOrder4, "IIRBurgOrder4" );
return 255;
}
......@@ -106,69 +111,76 @@ int main( int, char** )
void Run( CITAThirdOctaveFilterbank* pFilterBank, std::string sExportName )
{
double dT;
cout << " ### " << sExportName << " ### " << endl;
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oDirectivity );
sw.start();
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime directivity magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( sExportName + "_Dir.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oIdentity, false );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i * g_iBlockLength );
dT = sw.stop();
cout << "Runtime identity magnitudes: \t" << timeToString( dT ) << " (" << dT / sbImpulseResponse.GetLength() * g_dSampleRate * 100.0f << "% of budget)" << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_Identity.wav", &sbImpulseResponse, g_dSampleRate );
writeAudiofile( sExportName + "_Identity.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oZero, false );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime zero magnitude:\t\t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_Zeros.wav", &sbImpulseResponse, g_dSampleRate );
writeAudiofile( sExportName + "_Zeros.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oSomeBands );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime real magnitudes:\t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_SomeBands.wav", &sbImpulseResponse, g_dSampleRate );
writeAudiofile( sExportName + "_SomeBands.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oHighPass );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime high-pass magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_HighPass.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oAbsorption );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
dT = sw.stop();
cout << "Runtime absorption magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_Absorb.wav", &sbImpulseResponse, g_dSampleRate );
writeAudiofile( sExportName + "_HighPass.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oDirectivity );
pFilterBank->SetMagnitudes( oReflection );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime directivity magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_Dir.wav", &sbImpulseResponse, g_dSampleRate );
cout << "Runtime reflection magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( sExportName + "_Reflect.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
pFilterBank->Clear();
pFilterBank->SetMagnitudes( oDiffraction );
sw.start();
pFilterBank->Process( sbDirac.GetData(), sbImpulseResponse.GetData() );
for( int i = 0; i < g_iSampleLength / g_iBlockLength; i++ )
pFilterBank->Process( sbDirac.GetData() + i*g_iBlockLength, sbImpulseResponse.GetData() + i*g_iBlockLength );
dT = sw.stop();
cout << "Runtime diffraction magnitudes: \t" << timeToString( dT ) << endl;
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_" + sExportName + "_Diffr.wav", &sbImpulseResponse, g_dSampleRate );
writeAudiofile( sExportName + "_Diffr.wav", &sbImpulseResponse, g_dSampleRate );
sbImpulseResponse.Zero();
cout << endl;
......
......@@ -6,48 +6,59 @@
% We use energy/power signal types as the impulse responses will give
% appropriate results of flat spectrum after filtering (a la insertion loss)
addpath( '../generators' )
%% FIR spline interpolation linear phase
fir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_Zeros.wav' );
fir_unity = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_Identity.wav' );
fir_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_SomeBands.wav' );
fir_hp = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_HighPass.wav' );
fir_abs = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_Absorb.wav' );
fir_dir = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_Dir.wav' );
fir_dif = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIRSplineLinearPhase_Diffr.wav' );
%% IIR Burg
iir_burg_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_Zeros.wav' );
iir_burg_unity = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_Identity.wav' );
iir_burg_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_SomeBands.wav' );
iir_burg_hp = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_HighPass.wav' );
iir_burg_abs = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_Absorb.wav' );
iir_burg_dir = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_Dir.wav' );
iir_burg_dif = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurgOrder10_Diffr.wav' );
fir_zeros = ita_read( 'FIRSplineLinearPhase_Zeros.wav' );
fir_unity = ita_read( 'FIRSplineLinearPhase_Identity.wav' );
fir_sb = ita_read( 'FIRSplineLinearPhase_SomeBands.wav' );
fir_hp = ita_read( 'FIRSplineLinearPhase_HighPass.wav' );
fir_refl = ita_read( 'FIRSplineLinearPhase_Reflect.wav' );
fir_dir = ita_read( 'FIRSplineLinearPhase_Dir.wav' );
fir_dif = ita_read( 'FIRSplineLinearPhase_Diffr.wav' );
%% IIR Burg order 4
iir_burg4_zeros = ita_read( 'IIRBurgOrder4_Zeros.wav' );
iir_burg4_unity = ita_read( 'IIRBurgOrder4_Identity.wav' );
iir_burg4_sb = ita_read( 'IIRBurgOrder4_SomeBands.wav' );
iir_burg4_hp = ita_read( 'IIRBurgOrder4_HighPass.wav' );
iir_burg4_refl = ita_read( 'IIRBurgOrder4_Reflect.wav' );
iir_burg4_dir = ita_read( 'IIRBurgOrder4_Dir.wav' );
iir_burg4_dif = ita_read( 'IIRBurgOrder4_Diffr.wav' );
%% IIR Burg order 10
iir_burg10_zeros = ita_read( 'IIRBurgOrder10_Zeros.wav' );
iir_burg10_unity = ita_read( 'IIRBurgOrder10_Identity.wav' );
iir_burg10_sb = ita_read( 'IIRBurgOrder10_SomeBands.wav' );
iir_burg10_hp = ita_read( 'IIRBurgOrder10_HighPass.wav' );
iir_burg10_refl = ita_read( 'IIRBurgOrder10_Reflect.wav' );
iir_burg10_dir = ita_read( 'IIRBurgOrder10_Dir.wav' );
iir_burg10_dif = ita_read( 'IIRBurgOrder10_Diffr.wav' );
%% IIR BiQuads
iir_bq_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_Zeros.wav' );
iir_bq_unity = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_Identity.wav' );
iir_bq_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_SomeBands.wav' );
iir_bq_hp = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_HighPass.wav' );
iir_bq_abs = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_Absorb.wav' );
iir_bq_dir = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_Dir.wav' );
iir_bq_dif = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBiquadOrder10_Diffr.wav' );
iir_bq_zeros = ita_read( 'IIRBiquadOrder10_Zeros.wav' );
iir_bq_unity = ita_read( 'IIRBiquadOrder10_Identity.wav' );
iir_bq_sb = ita_read( 'IIRBiquadOrder10_SomeBands.wav' );
iir_bq_hp = ita_read( 'IIRBiquadOrder10_HighPass.wav' );
iir_bq10_refl = ita_read( 'IIRBiquadOrder10_Reflect.wav' );
iir_bq10_dir = ita_read( 'IIRBiquadOrder10_Dir.wav' );
iir_bq_dif = ita_read( 'IIRBiquadOrder10_Diffr.wav' );
%% JSON import
json_zeros = load_result_from_json_file( 'ZeroSpectrum.json' );
json_ident = load_result_from_json_file( 'IdentitySpectrum.json' );
json_hp = load_result_from_json_file( 'HighPassSpectrum.json' );
json_sb = load_result_from_json_file( 'SomeBandsSpectrum.json' );
json_refl = load_result_from_json_file( 'ReflectedSpectrum.json' );
json_dir = load_result_from_json_file( 'DirectivitySpectrum.json' );
json_diffr = load_result_from_json_file( 'DiffractionSpectrum.json' );
json_zeros = ita_result_from_json( 'ZeroSpectrum.json' );
json_ident = ita_result_from_json( 'IdentitySpectrum.json' );
json_hp = ita_result_from_json( 'HighPassSpectrum.json' );
json_sb = ita_result_from_json( 'SomeBandsSpectrum.json' );
json_refl = ita_result_from_json( 'ReflectedSpectrum.json' );
json_dir = ita_result_from_json( 'DirectivitySpectrum.json' );
json_diffr = ita_result_from_json( 'DiffractionSpectrum.json' );
%% Plots
cn = { 'Original input spectrum', 'FIR Spline-Interpolation Linear-Phase', 'IIR Burg Design Algorithm', 'IIR BiQuad filterbank structure' };
cn = { 'Target energetic spectrum', 'FIR spline-interpolation linear-phase (128 taps)', 'IIR Burg Design Algorithm (order 4)', 'IIR Burg Design Algorithm (order 10)', 'IIR Biquad filterbank (order 10)' };
if 0
cp_zeros = ita_merge( fir_zeros, iir_burg_zeros, iir_bq_zeros );
cp_zeros = ita_merge( fir_zeros, iir_burg4_zeros, iir_burg10_zeros, iir_bq_zeros );
cp_zeros.signalType = 'energy';
cp_zeros_r = ita_spk2frequencybands( cp_zeros );
json_zero_r = cp_zeros_r.ch(1);
......@@ -59,7 +70,7 @@ if 0
end
if 0
cp_ident = ita_merge( fir_unity, iir_burg_unity, iir_bq_unity );
cp_ident = ita_merge( fir_unity, iir_burg4_unity, iir_burg10_unity, iir_bq_unity );
cp_ident.signalType = 'energy';
cp_ident_r = ita_spk2frequencybands( cp_ident );
json_zero_r = cp_ident_r.ch(1);
......@@ -71,7 +82,7 @@ if 0
end
if 0
cp_sb = ita_merge( fir_sb, iir_burg_sb, iir_bq_sb );
cp_sb = ita_merge( fir_sb, iir_burg4_sb, iir_burg10_sb, iir_bq_sb );
cp_sb.signalType = 'energy';
cp_sb_r = ita_spk2frequencybands( cp_sb );
json_sb_r = cp_sb_r.ch(1);
......@@ -83,7 +94,7 @@ if 0
end
if 0
cp_hp = ita_merge( fir_hp, iir_burg_hp, iir_bq_hp );
cp_hp = ita_merge( fir_hp, iir_burg4_hp, iir_burg10_hp, iir_bq_hp );
cp_hp.signalType = 'energy';
cp_hp_r = ita_spk2frequencybands( cp_hp );
json_hp_r = cp_hp_r.ch(1);
......@@ -96,7 +107,7 @@ end
if 1
cp_dir = ita_merge( fir_dir, iir_burg_dir, iir_bq_dir );
cp_dir = ita_merge( fir_dir, iir_burg4_dir, iir_burg10_dir, iir_bq10_dir );
cp_dir.signalType = 'energy';
cp_dir_r = ita_spk2frequencybands( cp_dir );
json_dir_r = cp_dir_r.ch(1);
......@@ -109,22 +120,16 @@ end
%% Helper functions
function ita_result_t = load_result_from_json_file( file_path )
json_t = jsondecode( fileread( file_path ) );
ita_result_t = itaResult();
ita_result_t.domain = 'freq';
fn = fields( json_t.spectrum );
ita_result_t.freqVector = zeros( numel( fn ), 1 );
ita_result_t.freqData = zeros( numel( fn ), 1 );
for n = 1:numel( fn )
ita_result_t.freqVector( n ) = json_t.spectrum.( fn{ n } ).center_frequency;
ita_result_t.freqData( n ) = json_t.spectrum.( fn{ n } ).value;
if 0
cp_refl = ita_merge( fir_refl, iir_burg4_refl, iir_burg10_refl, iir_bq10_refl );
cp_refl.signalType = 'energy';
cp_refl_r = ita_spk2frequencybands( cp_refl );
json_refl_r = cp_refl_r.ch(1);
json_refl_r.freqData = json_refl.freqData;
cp_refl_r = ita_merge( json_refl_r, cp_refl_r );
cp_refl_r.channelNames = cn;
cp_refl_r.pf
title( 'Reflection spectrum' )
end
end
......@@ -48,7 +48,7 @@ using namespace std;
using namespace ITABase;
const float g_fSampleRate = 44100;
const int g_iFIRGenLength = int( 4 * ceil( g_fSampleRate / CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );;
const int g_iFIRGenLength = int( 4 * ceil( g_fSampleRate / CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
// Absorption filter
vector< float > vfAbsorptionSpectrumTO = { 0.02f, 0.02f, 0.02f, 0.0233333f, 0.0266667f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.03f, 0.0333333f, 0.0366667f, 0.04f, 0.0433333f, 0.0466667f, 0.05f, 0.0566667f, 0.0633333f, 0.07f, 0.07f, 0.07f, 0.07f, 0.0733333f, 0.0766667f, 0.08f, 0.08f };
......@@ -77,7 +77,8 @@ int main( int, char** )
void DesignFromSpectrum()
{
ITADSP::CFilterCoefficients oBurgCoeffs( 4, false );
ITADSP::CFilterCoefficients oBurgCoeffsOrder4( 4, false );
ITADSP::CFilterCoefficients oBurgCoeffsOrder10( 10, false );
// Generate an impulse response from the given spectra
ITABase::CThirdOctaveGainMagnitudeSpectrum oMags;
......@@ -93,9 +94,11 @@ void DesignFromSpectrum()
#endif
oIRGenerator.GenerateFilter( oMags, oIRTarget.GetData(), true );
oIRTarget.Store( "ReflectionIR_min_phase.wav" );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffs );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder4 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder10 );
#ifdef WITH_JSON_SUPPORT
ITADSP::ExportIIRCoefficientsToJSON( "ReflectionIIRCoeffs.json", oBurgCoeffs );
ITADSP::ExportIIRCoefficientsToJSON( "ReflectionIIRCoeffsOrder4.json", oBurgCoeffsOrder4 );
ITADSP::ExportIIRCoefficientsToJSON( "ReflectionIIRCoeffsOrder4.json", oBurgCoeffsOrder10 );
#endif
// Directivity spectrum
......@@ -105,10 +108,13 @@ void DesignFromSpectrum()
#endif
oIRGenerator.GenerateFilter( oMags, oIRTarget.GetData(), true );
oIRTarget.Store( "DirectivityIR_min_phase.wav" );
oBurgCoeffs.InitialiseOrder( 4 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffs );
oBurgCoeffsOrder4.InitialiseOrder( 4 );
oBurgCoeffsOrder10.InitialiseOrder( 10 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder4 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder10 );
#ifdef WITH_JSON_SUPPORT
ITADSP::ExportIIRCoefficientsToJSON( "DirectivityIIRCoeffs.json", oBurgCoeffs );
ITADSP::ExportIIRCoefficientsToJSON( "DirectivityIIRCoeffsOrder4.json", oBurgCoeffsOrder4 );
ITADSP::ExportIIRCoefficientsToJSON( "DirectivityIIRCoeffsOrder10.json", oBurgCoeffsOrder10 );
#endif
// Diffraction spectrum
......@@ -118,10 +124,13 @@ void DesignFromSpectrum()
#endif
oIRGenerator.GenerateFilter( oMags, oIRTarget.GetData(), true );
oIRTarget.Store( "DiffractionIR_min_phase.wav" );
oBurgCoeffs.InitialiseOrder( 4 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffs );
oBurgCoeffsOrder4.InitialiseOrder( 4 );
oBurgCoeffsOrder10.InitialiseOrder( 10 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder4 );
ITADSP::IIRFilterGenerator::Burg( oIRTarget, oBurgCoeffsOrder10 );
#ifdef WITH_JSON_SUPPORT
ITADSP::ExportIIRCoefficientsToJSON( "DiffractionIIRCoeffs.json", oBurgCoeffs );
ITADSP::ExportIIRCoefficientsToJSON( "DiffractionIIRCoeffsOrder4.json", oBurgCoeffsOrder4 );
ITADSP::ExportIIRCoefficientsToJSON( "DiffractionIIRCoeffsOrder10.json", oBurgCoeffsOrder10 );
#endif
}
......
......@@ -5,17 +5,17 @@ reflection_ir.channelNames = { 'C++ IR from spectrum generator' };
excitation_ir = ita_generate_impulse( 'fftDegree', reflection_ir.fftDegree );
reflection_coeffs = jsondecode( fileread( 'ReflectionIIRCoeffs.json' ) );
reflection_iir = dsp.IIRFilter;
reflection_iir.Numerator = [ reflection_coeffs.numerator.b0 0 0 0 0 ];
reflection_iir.Denominator = [ reflection_coeffs.denominator.a1 ...
reflection_coeffs.denominator.a2 ...
reflection_coeffs.denominator.a3 ...
reflection_coeffs.denominator.a4 ...
reflection_coeffs.denominator.a5 ];
reflection_coeffs4 = jsondecode( fileread( 'ReflectionIIRCoeffsOrder4.json' ) );
reflection_iir_o4 = dsp.IIRFilter;
reflection_iir_o4.Numerator = [ reflection_coeffs4.numerator.b0 0 0 0 0 ];
reflection_iir_o4.Denominator = [ reflection_coeffs4.denominator.a1 ...
reflection_coeffs4.denominator.a2 ...
reflection_coeffs4.denominator.a3 ...
reflection_coeffs4.denominator.a4 ...
reflection_coeffs4.denominator.a5 ];
reflection_ir_filtering = reflection_ir;