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

Tests progress

parent 6c9da095
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <ITAConstants.h>
#include <ITAFastMath.h>
#include <ITANumericUtils.h>
#include <ITAStringUtils.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#include <spline.h>
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();
if( oTOMagnitudes.IsZero() )
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0.0f;
return;
}
if( oTOMagnitudes.IsIdentity() )
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0.0f;
pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
return;
}
// 1st step: Interpolate the magnitudes
m_pfInputData[ 0 ] = 1.0f;
for( int i = 0; i < CITAThirdOctaveMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputData[ 1 + i ] = float( db10_to_ratio( oTOMagnitudes[ i ] ) );
m_pfInputData[ m_iInputFreqs - 1 ] = 0.0f; // @todo jst: check if this is good
// 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;
const float fScale = 1 / ( float ) m_iFilterLength;
// No DC offset, ever!
m_pfBuf1[ 0 ] = 0;
m_pfBuf1[ 1 ] = 0;
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 );
// 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;
}
// 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 ];
}
// @todo: Minimum-phase?
m_sw.stop();
}
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <ITAThirdOctaveFilterBank.h>
#include <ITAAudiofileWriter.h>
#include <ITAStringUtils.h>
#include <ITAStopWatch.h>
#include <ITASampleBuffer.h>
#include <ITASampleFrame.h>
#include <iostream>
using namespace std;
const double g_dSampleRate = 44100;
const int g_iFilterLength = 512;
void TestThirdOctaveFilterGeneratorFIRIdentity();
void TestThirdOctaveFilterGeneratorFIRZero();
void TestThirdOctaveFilterGeneratorFIRSingleBands();
void main( int, char** )
{
TestThirdOctaveFilterGeneratorFIRIdentity();
TestThirdOctaveFilterGeneratorFIRZero();
TestThirdOctaveFilterGeneratorFIRSingleBands();
}
void TestThirdOctaveFilterGeneratorFIRIdentity()
{
CITAThirdOctaveMagnitudeSpectrum oMags;
oMags.SetIdentity();
ITASampleBuffer oFilter( g_iFilterLength );
CITAThirdOctaveFIRFilterGenerator oFilterGenerator( g_dSampleRate, g_iFilterLength );
oFilterGenerator.GenerateFilter( oMags, oFilter.GetData() );
string sFilePath = "ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR.wav";
writeAudiofile( sFilePath, &oFilter, g_dSampleRate );
cout << "Exported result to " << sFilePath << endl;
}
void TestThirdOctaveFilterGeneratorFIRZero()
{
CITAThirdOctaveMagnitudeSpectrum oMags;
oMags.SetZero();
ITASampleBuffer oFilter( g_iFilterLength );
CITAThirdOctaveFIRFilterGenerator oFilterGenerator( g_dSampleRate, g_iFilterLength );
oFilterGenerator.GenerateFilter( oMags, oFilter.GetData() );
string sFilePath = "ITADSPThirdOctaveFilterGeneratorTest_Zero_FIR.wav";
writeAudiofile( sFilePath, &oFilter, g_dSampleRate );
cout << "Exported result to " << sFilePath << endl;
}
void TestThirdOctaveFilterGeneratorFIRSingleBands()
{
ITASampleFrame oFilter( CITAThirdOctaveMagnitudeSpectrum::GetNumBands(), g_iFilterLength, true );
CITAThirdOctaveFIRFilterGenerator oFilterGenerator( g_dSampleRate, g_iFilterLength );
CITAThirdOctaveMagnitudeSpectrum oMags;
for( int i = 0; i < CITAThirdOctaveMagnitudeSpectrum::GetNumBands(); i++ )
{
oMags.SetZero();
oMags[ i ] = 0.0f; // dB
oFilterGenerator.GenerateFilter( oMags, oFilter[ i ].GetData() );
}
string sFilePath = "ITADSPThirdOctaveFilterGeneratorTest_SingleBands_FIR.wav";
writeAudiofile( sFilePath, &oFilter, g_dSampleRate );
cout << "Exported result to " << sFilePath << endl;
}
%% ITADSPThirdOctaveFilterGeneratorTest evaluation Matlab script
% Loads the FIR filters resulting from the test ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR
% Requires ITA Toolbox, see http://ita-toolbox.org
ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR = ita_read( 'ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR.wav' );
ITADSPThirdOctaveFilterGeneratorTest_Zero_FIR = ita_read( 'ITADSPThirdOctaveFilterGeneratorTest_Zero_FIR.wav' );
ITADSPThirdOctaveFilterGeneratorTest_SingleBands_FIR = ita_read( 'ITADSPThirdOctaveFilterGeneratorTest_SingleBands_FIR.wav' );
ITADSPThirdOctaveFilterGeneratorTest_Zero_FIR.channelNames = { 'Zero' };
ITADSPThirdOctaveFilterGeneratorTest_Zero_FIR.pf
ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR.channelNames = { 'Identity' };
ITADSPThirdOctaveFilterGeneratorTest_Identity_FIR.pf
ITADSPThirdOctaveFilterGeneratorTest_SingleBands_FIR.pf
#include <ITAThirdOctaveFilterBank.h>
#include <ITAAudiofileWriter.h>
#include <ITAStopWatch.h>
#include <ITASampleBuffer.h>
#include <iostream>
using namespace std;
const double g_dSampleRate = 44100;
const int g_iFilterLength = 512;
void TestThirdOctaveFilterbankIIR();
void TestThirdOctaveFilterbankFIR() {};
void main( int, char** )
{
TestThirdOctaveFilterbankIIR();
TestThirdOctaveFilterbankFIR();
}
void TestThirdOctaveFilterbankIIR()
{
const int iSampleLength = ( 1 << 17 );
CITAThirdOctaveFilterbank* pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
ITASampleBuffer x( iSampleLength );
x[ 0 ] = 1.0f;
CITAThirdOctaveMagnitudeSpectrum oMags;
oMags.SetIdentity();
pIIRFilterbank->SetGains( oMags );
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( "bankout_zeros_iir.wav", &y, g_dSampleRate );
y.Zero();
delete pIIRFilterbank;
pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
oMags.SetZero();
pIIRFilterbank->SetGains( oMags );
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime zero magnitude:" << sw.ToString() << endl;
y.Normalize();
writeAudiofile( "bankout_identity_iir.wav", &y, g_dSampleRate );
delete pIIRFilterbank;
pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
oMags[ 10 ] = 0.01f;
oMags[ 11 ] = 0.001f;
oMags[ 12 ] = 0.050f;
oMags[ 13 ] = 0.30f;
oMags[ 14 ] = 0.80f;
oMags[ 3 ] = 2.0f;
pIIRFilterbank->SetGains( oMags );
sw.start();
pIIRFilterbank->Process( x.GetData(), y.GetData() );
sw.stop();
cout << "Runtime real magnitudes:" << sw.ToString() << endl;
y.Normalize();
writeAudiofile( "bankout_realmags_iir.wav", &y, g_dSampleRate );
delete pIIRFilterbank;
}
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