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

Trying to find bug in Burg algorithm

parent 3342241b
......@@ -44,10 +44,12 @@ namespace ITADSP
inline CFilterCoefficients( unsigned int uiOrder = 0, bool bARMAInit = false )
: iDesignAlgorithm( ITADSP::UNSPECIFIED )
, bIsARMA( bARMAInit )
{
Initialise( uiOrder, bARMAInit );
};
//! Initialize (sets coeffs to zero)
inline void Initialise( unsigned int uiOrder_, bool bARMAInit_ )
{
bIsARMA = bARMAInit_;
......@@ -58,8 +60,11 @@ namespace ITADSP
vfNumerator.resize( 1 );
vfDenominator.resize( uiOrder + 1 );
SetZero();
};
//! Initialize with given order (sets coeffs to zero)
inline void InitialiseOrder( unsigned int uiOrder_ )
{
uiOrder = uiOrder_;
......@@ -69,11 +74,13 @@ namespace ITADSP
vfNumerator.resize( 1 );
vfDenominator.resize( uiOrder + 1 );
SetZero();
};
inline void SetZero()
{
for( size_t i = 1; i < vfNumerator.size(); i++ )
for( size_t i = 0; i < vfNumerator.size(); i++ )
vfNumerator[ i ] = 0.0f;
if( vfDenominator.size() > 0 )
......
......@@ -10,26 +10,26 @@ void ITADSP::IIRFilterGenerator::Yulewalk( const ITABase::CFiniteImpulseResponse
int na = oCoeffs.vfDenominator.size();
int n = oIR.GetLength();
int n2 = floor((n + 1) / 2);
int n2 = floor( ( n + 1 ) / 2 );
int nb = na;
int nr = 4 * na;
//nt = 0:nr-1
//pick nr coefficients
ITASampleBuffer R(nr, false);
for (int i = 0; i < nr; i++)
R.GetData()[i] = oIR.GetData()[i] * (0.54 + 0.46*cos(ITAConstants::PI_F*i/(nr-1)));
ITASampleBuffer R( nr, false );
for( int i = 0; i < nr; i++ )
R.GetData()[ i ] = oIR.GetData()[ i ] * ( 0.54 + 0.46*cos( ITAConstants::PI_F*i / ( nr - 1 ) ) );
//window used t extract right wing of two sided covaianc sequence
ITASampleBuffer Rwindow(n, true);
Rwindow[0] = 0.5f;
for (int i = 1; i <= n2; i++ )
Rwindow[i] = 1.0f;
ITASampleBuffer Rwindow( n, true );
Rwindow[ 0 ] = 0.5f;
for( int i = 1; i <= n2; i++ )
Rwindow[ i ] = 1.0f;
//find roots of
float *A = new float(oCoeffs.vfDenominator.size());
float *A = new float( oCoeffs.vfDenominator.size() );
denf( A, R, na );
......@@ -37,20 +37,20 @@ void ITADSP::IIRFilterGenerator::Yulewalk( const ITABase::CFiniteImpulseResponse
}
void ITADSP::IIRFilterGenerator::denf(float *A, const ITASampleBuffer& R, int na) {
void ITADSP::IIRFilterGenerator::denf( float *A, const ITASampleBuffer& R, int na ) {
//output coefficients A of length na + 1, input vector R, and order na
float **out;
out = new float*[R.GetLength() - na - 1];
for (int i = 0; i < R.GetLength() - na - 1; i++) {
out[i] = new float[na];
out = new float*[ R.GetLength() - na - 1 ];
for( int i = 0; i < R.GetLength() - na - 1; i++ ) {
out[ i ] = new float[ na ];
}
toeplitz( out, R, na);
toeplitz( out, R, na );
for (int i = 0; i < R.GetLength() - na - 1; i++) {
for (int j = 0; j < na; j++) {
std::cout << out[i][j] << '\t';
for( int i = 0; i < R.GetLength() - na - 1; i++ ) {
for( int j = 0; j < na; j++ ) {
std::cout << out[ i ][ j ] << '\t';
}
std::cout << '\n';
}
......@@ -59,52 +59,52 @@ void ITADSP::IIRFilterGenerator::denf(float *A, const ITASampleBuffer& R, int na
void ITADSP::IIRFilterGenerator::toeplitz(float **out, const ITASampleBuffer &row, const ITASampleBuffer &column) {
void ITADSP::IIRFilterGenerator::toeplitz( float **out, const ITASampleBuffer &row, const ITASampleBuffer &column ) {
//return a toeplitz matrix out constructed from input vectors row and column
//initialise 2d array out to size [column][row]
out = new float*[column.GetLength()];
for (int i = 0; i < column.GetLength(); i++) {
out[i] = new float[row.GetLength()];
out = new float*[ column.GetLength() ];
for( int i = 0; i < column.GetLength(); i++ ) {
out[ i ] = new float[ row.GetLength() ];
}
int test = std::min(row.GetLength(), column.GetLength());
int test = std::min( row.GetLength(), column.GetLength() );
//set the upper diagonal (excluding leading edge)
for (int i = 1; i < row.GetLength(); i++) {
for( int j=0; j<test-i; j++ )
out[j][j+i] = row[i];
for( int i = 1; i < row.GetLength(); i++ ) {
for( int j = 0; j < test - i; j++ )
out[ j ][ j + i ] = row[ i ];
}
//set lower diagonal (including leading edge)
for (int i = 0; i < column.GetLength(); i++) {
for (int j = 0; j < test-i; j++)
out[j+i][j] = column[i];
for( int i = 0; i < column.GetLength(); i++ ) {
for( int j = 0; j < test - i; j++ )
out[ j + i ][ j ] = column[ i ];
}
return;
}
void ITADSP::IIRFilterGenerator::toeplitz(float **out, const ITASampleBuffer &in, int na) {
void ITADSP::IIRFilterGenerator::toeplitz( float **out, const ITASampleBuffer &in, int na ) {
//return a toeplitz matrix out constructed from input vectors row and column
/*
//initialise 2d array out to size [in.GetLength()-na-1][na]
out = new float*[in.GetLength() - na - 1];
for (int i = 0; i < in.GetLength() - na - 1; i++) {
out[i] = new float[na];
out[i] = new float[na];
}
*/
//set the upper diagonal (excluding leading edge)
for (int i = 1; i < na; i++) {
for (int j = 0; j < na - i; j++)
out[j][j + i] = in[na-i];
for( int i = 1; i < na; i++ ) {
for( int j = 0; j < na - i; j++ )
out[ j ][ j + i ] = in[ na - i ];
}
//set lower diagonal (including leading edge)
for (int i = 0; i < in.GetLength()-na-1; i++) {
for (int j = 0; j < na - i; j++)
out[j + i][j] = in[na+i-1];
for( int i = 0; i < in.GetLength() - na - 1; i++ ) {
for( int j = 0; j < na - i; j++ )
out[ j + i ][ j ] = in[ na + i - 1 ];
}
return;
......@@ -113,58 +113,63 @@ void ITADSP::IIRFilterGenerator::toeplitz(float **out, const ITASampleBuffer &in
void ITADSP::IIRFilterGenerator::Burg(const ITABase::CFiniteImpulseResponse& oIR, CFilterCoefficients& oCoeffs) {
void ITADSP::IIRFilterGenerator::Burg( const ITABase::CFiniteImpulseResponse& oIR, CFilterCoefficients& oCoeffs ) {
//oIR -> targe impulse response -> x, filter order = oCoeffs.uiOrder -> p
float k;
float k, k_numerator, k_denominator;
int buffer_length = oIR.GetLength() - 1;
ITASampleBuffer efp(buffer_length);
ITASampleBuffer ebp(buffer_length);
ITASampleBuffer efp_temp(buffer_length);
ITASampleBuffer efp( buffer_length );
ITASampleBuffer ebp( buffer_length );
ITASampleBuffer efp_temp( buffer_length );
oIR.read(efp.GetData(), buffer_length, 1); //read oIR 1:end
oIR.read(ebp.GetData(), buffer_length, 0); //read oIR 0:end - 1
oIR.read( efp.GetData(), buffer_length, 1 ); //read oIR 1:end
oIR.read( ebp.GetData(), buffer_length, 0 ); //read oIR 0:end - 1
std::vector<float> a_temp;
a_temp.resize(oCoeffs.uiOrder + 1);
for (auto it = oCoeffs.vfDenominator.begin(); it != oCoeffs.vfDenominator.end(); it++)
a_temp.resize( oCoeffs.uiOrder + 1 );
for( auto it = oCoeffs.vfDenominator.begin(); it != oCoeffs.vfDenominator.end(); it++ )
*it = 0.0f;
oCoeffs.vfNumerator[0] = InnerProduct( oIR.GetData(), oIR.GetData(), oIR.GetLength()) / oIR.GetLength();
oCoeffs.vfDenominator[0] = 1.0f;
oCoeffs.vfNumerator[ 0 ] = InnerProduct( oIR.GetData(), oIR.GetData(), oIR.GetLength() ) / oIR.GetLength();
oCoeffs.vfDenominator[ 0 ] = 1.0f;
oCoeffs.bIsARMA = false;
oCoeffs.iDesignAlgorithm = ITADSP::BURG; //record that the Burg algorithm was used to generate the coefficients
for (int m = 0; m < oCoeffs.uiOrder; m++) {
k = (-2 * InnerProduct(ebp.GetData(), efp.GetData()+m, buffer_length)) /
(InnerProduct(efp.GetData()+m,efp.GetData()+m, buffer_length) + InnerProduct(ebp.GetData(),ebp.GetData(), buffer_length));
for( int m = 0; m < oCoeffs.uiOrder; m++ )
{
k_numerator = ( -2 * InnerProduct( ebp.GetData(), efp.GetData() + m, buffer_length ) );
k_denominator = ( InnerProduct( efp.GetData() + m, efp.GetData() + m, buffer_length ) + InnerProduct( ebp.GetData(), ebp.GetData(), buffer_length ) );
if( k_denominator != 0.0f )
k = k_numerator / k_denominator;
else
k = 0.0f;
buffer_length--;
efp_temp = efp;
efp_temp.MulAdd(ebp, k, 0, m, buffer_length); //makes sure efp and abp are added in correct alignment
ebp.MulAdd(efp, k, m, 0, buffer_length);
efp_temp.MulAdd( ebp, k, 0, m, buffer_length ); //makes sure efp and abp are added in correct alignment
ebp.MulAdd( efp, k, m, 0, buffer_length );
efp = efp_temp;
for (int j = 0; j <= m; j++)
a_temp[j+1] = oCoeffs.vfDenominator[j + 1] + (k*oCoeffs.vfDenominator[m - j]);
for (int i = 1; i <= m+1; i++)
oCoeffs.vfDenominator[i] = a_temp[i];
oCoeffs.vfNumerator[0] = (1.0f - pow(k,2)) * oCoeffs.vfNumerator[0];
for( int j = 0; j <= m; j++ )
a_temp[ j + 1 ] = oCoeffs.vfDenominator[ j + 1 ] + ( k*oCoeffs.vfDenominator[ m - j ] );
for( int i = 1; i <= m + 1; i++ )
oCoeffs.vfDenominator[ i ] = a_temp[ i ];
oCoeffs.vfNumerator[ 0 ] = ( 1.0f - pow( k, 2 ) ) * oCoeffs.vfNumerator[ 0 ];
}
oCoeffs.vfNumerator[0] = sqrt( oCoeffs.vfNumerator[0] * oIR.GetLength() );
oCoeffs.vfNumerator[ 0 ] = sqrt( oCoeffs.vfNumerator[ 0 ] * oIR.GetLength() );
}
float ITADSP::IIRFilterGenerator::InnerProduct( const float *x, const float *y, const int length ) {
float ITADSP::IIRFilterGenerator::InnerProduct( const float *x, const float *y, const int length )
{
//calculate the inner product of two sample buffer inputs. Lenghts of the inputs must be equal
float out = 0.0f;
for (int i = 0; i < length; i++) {
out += x[i] * y[i];
}
for( int i = 0; i < length; i++ )
out += x[ i ] * y[ i ];
return out;
}
\ No newline at end of file
......@@ -42,7 +42,7 @@ void TestThirdOctaveFilterbankIIRBiquad()
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime identity magnitudes:" << sw.ToString() << endl;
cout << "Runtime identity magnitudes: " << sw.ToString() << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_Identity.wav", &y, g_dSampleRate );
......@@ -56,7 +56,7 @@ void TestThirdOctaveFilterbankIIRBiquad()
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime zero magnitude:" << sw.ToString() << endl;
cout << "Runtime zero magnitude: " << sw.ToString() << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_Zeros.wav", &y, g_dSampleRate );
......@@ -76,7 +76,7 @@ void TestThirdOctaveFilterbankIIRBiquad()
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime real magnitudes:" << sw.ToString() << endl;
cout << "Runtime real magnitudes: " << sw.ToString() << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBiquad_SomeBands.wav", &y, g_dSampleRate );
......@@ -93,8 +93,8 @@ void TestThirdOctaveFilterbankIIRBurg()
x[ 0 ] = 1.0f;
ITABase::CThirdOctaveGainMagnitudeSpectrum oMags;
oMags.SetIdentity();
oMags.SetIdentity();
pIIRFilterbank->SetMagnitudes( oMags, false );
ITASampleBuffer y( iSampleLength );
......@@ -103,7 +103,8 @@ void TestThirdOctaveFilterbankIIRBurg()
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime identity magnitudes:" << sw.ToString() << endl;
cout << "Runtime identity magnitudes: " << sw.ToString() << endl;
cout << "RMS (energy): " << y.RootMeanSquare( false ) << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_Identity.wav", &y, g_dSampleRate );
......@@ -111,13 +112,32 @@ void TestThirdOctaveFilterbankIIRBurg()
pIIRFilterbank->Clear();
oMags.SetIdentity();
oMags.Multiply( 0.5f );
pIIRFilterbank->SetMagnitudes( oMags, false );
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime identity half energy magnitudes: " << sw.ToString() << endl;
cout << "RMS (energy): " << y.RootMeanSquare( false ) << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_IdentityHalfEnergy.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;
cout << "Runtime zero magnitude: " << sw.ToString() << endl;
cout << "RMS (energy): " << y.RootMeanSquare( false ) << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_Zeros.wav", &y, g_dSampleRate );
......@@ -137,11 +157,28 @@ void TestThirdOctaveFilterbankIIRBurg()
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime real magnitudes:" << sw.ToString() << endl;
cout << "Runtime real magnitudes: " << sw.ToString() << endl;
cout << "RMS (energy): " << y.RootMeanSquare( false ) << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_SomeBands.wav", &y, g_dSampleRate );
// High pass
oMags.SetIdentity();
for( int i = 15; i >= 0; i-- )
oMags[ i ] = powf( 0.5f, 15 - i );
pIIRFilterbank->SetMagnitudes( oMags );
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime high pass magnitudes: " << sw.ToString() << endl;
cout << "RMS (energy): " << y.RootMeanSquare( false ) << endl;
y.Normalize();
writeAudiofile( "ITADSPThirdOctaveFilterbankTest_IIRBurg_HighPass.wav", &y, g_dSampleRate );
delete pIIRFilterbank;
}
......
......@@ -70,14 +70,16 @@ 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_unity_half = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_IdentityHalfEnergy.wav' );
iir_burg_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_SomeBands.wav' );
iir_burg_hp = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIRBurg_HighPass.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 = ita_merge( iir_burg_zeros, iir_burg_unity, iir_burg_unity_half, iir_burg_sb, iir_burg_hp );
iir_burg_all.channelNames = { 'zero', 'ident', 'ident half energy', 'some bands', 'high pass' };
iir_burg_all.pf
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