......@@ -3,7 +3,7 @@
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
* RWTH Aachen University, Germany, 2015-2020
*
* ----------------------------------------------------------------
* ____ __________ _______
......
#include <ITADSP/PD/JetEngine.h>
ITADSP::PD::CJetEngine::CJetEngine( double dSampleRate, float fRPMInit, bool bColdStart )
: m_dSampleRate( dSampleRate )
{
m_vfTurbineBaseModeFrequencies = { 3097.0f, 4495.0f, 5588.0f, 7414.0f, 11000.0f };
m_vfTurbineBaseModeAmplitudes = { 0.25f, 0.25f, 1.0f, 0.4f, 0.4f };
m_vfRPMRange = { 1000.0f, 5000.0f };
const float fCenterFrequency1 = 8000.0f;
const float fQ1 = 0.5f;
m_oForcedFlameBP1.setup( 1, float( dSampleRate ), fCenterFrequency1, fCenterFrequency1 / fQ1 );
const float fCutoffFrequency1 = 120.0f;
m_oForcedFlameHP1.setup( 1, float( dSampleRate ), fCutoffFrequency1 );
const float fLPCutoffFrequency1 = 11000.0f;
m_oJetEngineLP1.setup( 1, float( dSampleRate ), fLPCutoffFrequency1 );
const float fLPCutoffFrequency2 = 1.0f;
m_oJetEngineLP2.setup( 1, float( dSampleRate ), fLPCutoffFrequency2 );
oCtxAudio.vfTurbineModePhaseShift = { 0, 0, 0, 0, 0 };
oCtxAudio.vfTurbineModeFrequencies = { 0, 0, 0, 0, 0 };
SetRPM( fRPMInit );
if( !bColdStart )
{
// Run input low-pass for 3 secs (until 0.2 Hz LP died away on unit step function) to akkumulate
// of the default input, otherwise there is a run-up from zero to target RPM value
const float fNormalizedRPMInit = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMUpsample;
float* pfNormalizedRPMUpsampleAlias = &fNormalizedRPMUpsample;
for( int n = 0; n < ( int ) dSampleRate * 3; n++ )
{
fNormalizedRPMUpsample = fNormalizedRPMInit;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMUpsampleAlias );
}
}
};
ITADSP::PD::CJetEngine::~CJetEngine()
{
}
void ITADSP::PD::CJetEngine::SetRPM( float fRPM )
{
oCtxAudio.vfRPM = fRPM; // Direct write-through, no lock & swap implemented here
}
float ITADSP::PD::CJetEngine::GetRPMNormalized( float fRPM ) const
{
// Validate input range
assert( m_vfRPMRange.size() == 2 );
const float fValidRPM = std::min( std::max( m_vfRPMRange[ 0 ], fRPM ), m_vfRPMRange[ 1 ] );
const float fNormalizedRPM = std::max( 0.1f, fValidRPM / m_vfRPMRange[ 1 ] );
return fNormalizedRPM;
}
void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
{
/* Validate, normalize and up-sampled RPM and apply low-pass filter */
// Processed DSP coefficient with normalized input
const float fNormalizedRPMControl = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMControlUpscaled;
// We need this ref pointer for Dsp functions, but its ugly we dont use it apart from that
float* pfTempSampleAlias = &m_fTempSample;
float* pfNormalizedRPMControlAlias = &fNormalizedRPMControlUpscaled;
for( int n = 0; n < iNumSamples; n++ )
{
// Low-pass input
fNormalizedRPMControlUpscaled = fNormalizedRPMControl;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMControlAlias );
fNormalizedRPMControlUpscaled = *pfNormalizedRPMControlAlias;
/* Update DSP network coefficients */
float& fCurrentSample( pfOutputBuffer[ n ] );
// Update turbine frequencies
assert( m_vfTurbineBaseModeFrequencies.size() == oCtxAudio.vfTurbineModeFrequencies.size() );
for( int i = 0; i < m_vfTurbineBaseModeFrequencies.size(); i++ )
oCtxAudio.vfTurbineModeFrequencies[ i ] = m_vfTurbineBaseModeFrequencies[ i ] * fNormalizedRPMControlUpscaled;
// Update forced flame filters
const float fCenterFrequency2 = fNormalizedRPMControlUpscaled * fNormalizedRPMControlUpscaled * 150.0f;
m_oForcedFlameBP2.setup( 1, float( m_dSampleRate ), fCenterFrequency2, fCenterFrequency2 / 1.0f );
const float fCenterFrequency3 = fNormalizedRPMControlUpscaled * 12000.0f;
const float fQ3 = 0.6f;
m_oForcedFlameBP3.setup( 1, float( m_dSampleRate ), fCenterFrequency3, fCenterFrequency3 / fQ3 );
/* Process DSP */
// Forced flame
m_fTempSample = oRNG.GenerateFloat( -1.0f, 1.0f ); // noise~
m_oForcedFlameBP1.process( 1, &pfTempSampleAlias ); // bd~ 8000 0.5
m_oForcedFlameBP2.process( 1, &pfTempSampleAlias ); // vcf~ 0 1 (real part only = band pass)
m_oForcedFlameHP1.process( 1, &pfTempSampleAlias ); // hip~ 120
m_fTempSample *= 120.0f; // *~ 120
m_fTempSample = ( m_fTempSample < -1.0f ) ? -1 : ( ( m_fTempSample > 1 ) ? 1 : m_fTempSample ); // clip~ -1 1
m_oForcedFlameBP3.process( 1, &pfTempSampleAlias ); // vcf~ 0 0.6 (real part only = band pass)
fCurrentSample = m_fTempSample; // override buffer
/* Turbine (adds to output buffer) */
m_fTempSample = 0.0f;
assert( oCtxAudio.vfTurbineModeFrequencies.size() == oCtxAudio.vfTurbineModePhaseShift.size() );
assert( oCtxAudio.vfTurbineModeFrequencies.size() == m_vfTurbineBaseModeAmplitudes.size() );
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
const float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
float fAmplitude( m_vfTurbineBaseModeAmplitudes[ i ] );
const float t = fmodf( ITAConstants::TWO_PI_F_L * fFrequency * float( n ) / float( m_dSampleRate ), ITAConstants::TWO_PI_F );
// [MANUAL MODIFICATION] sine-sqared roll-in of amplitude for the frequencies, otherwise very synthetic
if( fNormalizedRPMControlUpscaled < 0.1 )
fAmplitude *= sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f ) * sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f );
m_fTempSample += fAmplitude * sin( t + fPhaseShift ); // all osc~ and *~
}
m_fTempSample = ( m_fTempSample < -0.9f ) ? -0.9f : ( ( m_fTempSample > 0.9f ) ? 0.9f : m_fTempSample ); // clip~ -0.9 0.9
/* Jet engine */
m_fTempSample *= 0.5f; // *~ 0.5
const float fManualBalance = 0.1f / 0.5f; // [MANUAL MODIFICATION] not included in patch
m_fTempSample = fCurrentSample + fManualBalance * m_fTempSample; // combine turbine and flame from jet engine patch (with a manual balance that sound better)
m_oJetEngineLP1.process( 1, &pfTempSampleAlias ); // ~lop 11000
fCurrentSample = 0.2f * m_fTempSample; // *~0.2
}
/* Update phases for turbine */
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
fPhaseShift = fmodf( ITAConstants::TWO_PI_F_L * fFrequency / float( m_dSampleRate ) * float( iNumSamples ) + fPhaseShift, ITAConstants::TWO_PI_F );
}
}
#include <ITADSP/ThirdOctaveFilterbank/IIRBurg.h>
#include <ITAIIRFilterEngine.h>
#include <ITAIIRFilterGenerator.h>
#include <ITAThirdOctaveFIRFilterGenerator.h>
#include <algorithm>
//#include "ITAIIRUtils.h" // rm me
using namespace ITADSP::ThirdOctaveFilterbank;
CIIRBurg::CIIRBurg( const double dSampleRate, const int iBlockLength, const int iFilterOrder /* = 4 */ )
: m_iBlockLength( iBlockLength )
{
m_pFilterEngine = new CITAIIRFilterEngine( iFilterOrder );
int iFIRGenLength = int( 4 * ceil( dSampleRate / ITABase::CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
m_pFilterGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, iFIRGenLength );
}
CIIRBurg::~CIIRBurg()
{
delete m_pFilterGenerator;
delete m_pFilterEngine;
}
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(), 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 );
}
void CIIRBurg::Clear()
{
m_pFilterEngine->ClearAccumulators();
}
void ITADSP::ThirdOctaveFilterbank::CIIRBurg::Process( const float* pfInputSamples, float* pfOutputSamples )
{
m_pFilterEngine->Process( pfInputSamples, pfOutputSamples, m_iBlockLength );
}
......@@ -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
#include <ITAIIRUtils.h>
#include <ITAException.h>
#include <fstream>
......@@ -41,6 +42,8 @@ void ITADSP::ExportIIRCoefficientsToJSON( const std::string& sJSONFilePath, cons
fsOut.close();
#else
ITA_EXCEPT_NOT_IMPLEMENTED( "Export function not implemented or no JSON support activated for ITA DSP library" );
ITA_EXCEPT1( NOT_IMPLEMENTED, "Export function not implemented or no JSON support activated for ITA DSP library" );
#endif
}
......@@ -13,21 +13,15 @@ using namespace ITABase;
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 = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
m_pfInputFreqs = fm_falloc( m_iInputFreqs, true );
m_pfInputFreqs[ 0 ] = 0; // Left margin
m_iNumInputFreqs = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
m_vfInputFreqVector.resize( m_iNumInputFreqs );
m_vfInputFreqVector[ 0 ] = 0; // Left margin
for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputFreqs[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
m_pfInputFreqs[ m_iInputFreqs - 1 ] = ( float ) dSampleRate / 2; // Right margin: Nyquist frequency
m_vfInputFreqVector[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
m_vfInputFreqVector[ m_iNumInputFreqs - 1 ] = ( float ) dSampleRate / 2; // Right margin: Nyquist frequency
m_pfInputData = fm_falloc( m_iInputFreqs, true );
m_vfInputFreqData.resize( m_iNumInputFreqs );
// DFT frequency bandwidth
m_fDeltaF = ( float ) dSampleRate / ( float ) iFilterLength;
......@@ -35,55 +29,42 @@ CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const doub
// 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";
m_vfFreqDataInterpolated.resize( 2 * m_iDFTCoeffs );
m_vfImpulseResponse.resize( iFilterLength );
m_ifft.plan( ITAFFT::IFFT_C2R, iFilterLength, &m_vfFreqDataInterpolated[0], &m_vfImpulseResponse[0] );
}
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;
}
double CITAThirdOctaveFIRFilterGenerator::GetSampleRate() const
{
return m_dSamplerate;
}
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 )
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, ITASampleBuffer& sbImpulseResponse, bool bMinimumPhase /*=false*/ ) const
{
m_sDumpFilename = sFilename;
if( sbImpulseResponse.GetLength() < m_iFilterLength )
ITA_EXCEPT_INVALID_PARAMETER( "Given impulse response is shorter than initial filter size of generator, cannot proceed" );
GenerateFilter( oTOGainMagnitudes, sbImpulseResponse.GetData(), bMinimumPhase );
}
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ )
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ ) const
{
m_sw.start();
if( oTOGainMagnitudes.IsZero() )
{
for( int i = 0; i < m_iFilterLength; i++ )
......@@ -95,23 +76,26 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
{
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = 0.0f;
pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
if( bMinimumPhase )
pfFilterCoeffs[ 0 ] = 1.0f;
else
pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
return;
}
// 1st step: Interpolate the magnitudes
for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
m_pfInputData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );
m_vfInputFreqData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );
// Bounndaries must be defined in a meaningful way.
m_pfInputData[ 0 ] = m_pfInputData[ 1 ];
m_pfInputData[ m_iInputFreqs - 1 ] = 0.0f;
m_vfInputFreqData[ 0 ] = m_vfInputFreqData[ 1 ];
m_vfInputFreqData[ m_iNumInputFreqs - 1 ] = 0.0f;
// Initialize cubic spline interpolation
m_ypp = spline_cubic_set( m_iInputFreqs,
m_pfInputFreqs,
m_pfInputData,
m_ypp = spline_cubic_set( m_iNumInputFreqs,
&m_vfInputFreqVector[0],
&m_vfInputFreqData[0],
1, // Left boundary condition => 1st derivative m=0
0,
1, // Right boundary condition => 1st derivative m=0
......@@ -120,59 +104,48 @@ void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOct
const float fScale = 1 / ( float ) m_iFilterLength;
// No DC offset, ever!
m_pfBuf1[ 0 ] = 0;
m_pfBuf1[ 1 ] = 0;
m_vfFreqDataInterpolated[ 0 ] = 0;
m_vfFreqDataInterpolated[ 1 ] = 0;
if( bMinimumPhase ) {
if( bMinimumPhase )
{
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iInputFreqs,
m_pfInputFreqs,
float x = spline_cubic_val( m_iNumInputFreqs,
&m_vfInputFreqVector[0],
i*m_fDeltaF,
m_pfInputData,
&m_vfInputFreqData[0],
m_ypp,
&fDummy,
&fDummy );
// Phase-shift by half the FFT-period
m_pfBuf1[ 2 * i ] = pow( x * fScale, 2 ) * m_iFilterLength; //minimum phase
m_pfBuf1[ 2 * i + 1 ] = 0;
m_vfFreqDataInterpolated[ 2 * i ] = x * fScale; //minimum phase
m_vfFreqDataInterpolated[ 2 * i + 1 ] = 0;
}
}
else {
else
{
for( int i = 1; i < m_iDFTCoeffs; i++ )
{
float x = spline_cubic_val( m_iInputFreqs,
m_pfInputFreqs,
float x = spline_cubic_val( m_iNumInputFreqs,
&m_vfInputFreqVector[0],
i*m_fDeltaF,
m_pfInputData,
&m_vfInputFreqData[0],
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;
m_vfFreqDataInterpolated[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
m_vfFreqDataInterpolated[ 2 * i + 1 ] = 0;
}
}
// 2nd step: Convert into time-domain (out-of-place C2R-IFFT)
m_ifft.execute();
m_ifft.execute(); // populates m_vfImpulseResponse
// 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();
// Copy to output
for( int i = 0; i < m_iFilterLength; i++ )
pfFilterCoeffs[ i ] = m_vfImpulseResponse[ i ];
}
......@@ -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;
......
......@@ -3,7 +3,7 @@
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
* RWTH Aachen University, Germany, 2015-2020
*
* ----------------------------------------------------------------
* ____ __________ _______
......@@ -16,4 +16,83 @@
*
*/
#include "VAThirdOctaveFilterbank.h"
#include "ITAThirdOctaveFilterbankFIR.h"
#include <ITAUPConvolution.h>
#include <algorithm>
CITAThirdOctaveFilterbankFIR::CITAThirdOctaveFilterbankFIR( const double dSampleRate, const int iBlockLength )
: m_pGenerator( nullptr )
, m_pConvolver( nullptr )
, m_iBlocklength( iBlockLength )
{
// [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_vfGeneratedFIR.resize( m_iGeneratorFilterLength );
m_pGenerator = new CITAThirdOctaveFIRFilterGenerator( dSampleRate, m_iGeneratorFilterLength );
m_iConvolutionTaps = ( std::max )( 128, iBlockLength );
m_pConvolver = new ITAUPConvolution( iBlockLength, m_iConvolutionTaps );
m_pConvolver->SetFilterExchangeFadingFunction( ITABase::FadingFunction::COSINE_SQUARE );
m_pConvolver->SetFilterCrossfadeLength( 32 );
SetIdentity( false );
}
CITAThirdOctaveFilterbankFIR::~CITAThirdOctaveFilterbankFIR()
{
delete m_pGenerator;
delete m_pConvolver;
}
void CITAThirdOctaveFilterbankFIR::SetIdentity( const bool bSmoothChangeover /*=true*/ )
{
int iLatency = m_pGenerator->GetLatency();
if( iLatency > m_iGeneratorFilterLength )
ITA_EXCEPT_INVALID_PARAMETER( "Latency exceeds filter length in FIR filterbank identity setter" );
// Set a dirac with given latency
fm_zero( &m_vfGeneratedFIR[ 0 ], m_iGeneratorFilterLength );
m_vfGeneratedFIR[ iLatency ] = 1;
// Window (with rect), if convolution filter is shorter than generated filter
int iShift = std::max( 0, iLatency - m_iConvolutionTaps / 2 );
int iLength = std::min( m_iConvolutionTaps, m_iGeneratorFilterLength - iShift );
ITAUPFilter* pFilter = m_pConvolver->RequestFilter();
pFilter->Load( &m_vfGeneratedFIR[ iShift ], iLength );
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_vfGeneratedFIR[ 0 ], false );
int iRange = ( std::min )( m_iGeneratorFilterLength, m_iConvolutionTaps );
int iShift = std::max( 0, m_pGenerator->GetLatency() - iRange / 2 );
int iLength = std::min( m_iConvolutionTaps, m_iGeneratorFilterLength - iShift );
ITAUPFilter* pFilter = m_pConvolver->RequestFilter();
pFilter->Load( &m_vfGeneratedFIR[ iShift ], iLength );
m_pConvolver->ExchangeFilter( pFilter, ( bSmoothChangeover ? ITABase::FadingFunction::COSINE_SQUARE : ITABase::FadingFunction::SWITCH ) );
pFilter->Release(); // Auto-release
}
int CITAThirdOctaveFilterbankFIR::GetLatency() const
{
return m_pGenerator->GetLatency();
}
void CITAThirdOctaveFilterbankFIR::Clear()
{
m_pConvolver->clear();
SetIdentity( true );
}
void CITAThirdOctaveFilterbankFIR::Process( const float* pfInputSamples, float* pfOutputSamples )
{
m_pConvolver->Process( pfInputSamples, m_iBlocklength, pfOutputSamples, m_iBlocklength );
}
\ No newline at end of file
......@@ -18,95 +18,9 @@ if( ITA_VISTA_BUILD_STATIC )
add_definitions( -DVISTATOOLS_STATIC -DVISTABASE_STATIC -DVISTAMATH_STATIC -DVISTAASPECTS_STATIC -DVISTAINTERPROCCOMM_STATIC )
endif( )
add_executable( ITADSPIIRFilterTest ITADSPIIRFilterTest.cpp )
target_link_libraries( ITADSPIIRFilterTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPIIRFilterTest )
vista_install( ITADSPIIRFilterTest )
vista_create_default_info_file( ITADSPIIRFilterTest )
set_property( TARGET ITADSPIIRFilterTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPIIRFilterGeneratorTest ITADSPIIRFilterGeneratorTest.cpp )
target_link_libraries( ITADSPIIRFilterGeneratorTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPIIRFilterGeneratorTest )
vista_install( ITADSPIIRFilterGeneratorTest )
vista_create_default_info_file( ITADSPIIRFilterGeneratorTest )
set_property( TARGET ITADSPIIRFilterGeneratorTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPBiquadTest ITADSPBiquadTest.cpp )
target_link_libraries( ITADSPBiquadTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPBiquadTest )
vista_install( ITADSPBiquadTest )
vista_create_default_info_file( ITADSPBiquadTest )
set_property( TARGET ITADSPBiquadTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPVariableDelayLineTest ITADSPVariableDelayLineTest.cpp )
target_link_libraries( ITADSPVariableDelayLineTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPVariableDelayLineTest )
vista_install( ITADSPVariableDelayLineTest )
vista_create_default_info_file( ITADSPVariableDelayLineTest )
set_property( TARGET ITADSPVariableDelayLineTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPSIMOVDLTest ITADSPSIMOVDLTest.cpp )
target_link_libraries( ITADSPSIMOVDLTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPSIMOVDLTest )
vista_install( ITADSPSIMOVDLTest )
vista_create_default_info_file( ITADSPSIMOVDLTest )
set_property( TARGET ITADSPSIMOVDLTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPThirdOctaveFilterGeneratorTest ITADSPThirdOctaveFilterGeneratorTest.cpp )
target_link_libraries( ITADSPThirdOctaveFilterGeneratorTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPThirdOctaveFilterGeneratorTest )
vista_install( ITADSPThirdOctaveFilterGeneratorTest )
vista_create_default_info_file( ITADSPThirdOctaveFilterGeneratorTest )
set_property( TARGET ITADSPThirdOctaveFilterGeneratorTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPSIMOVDLSourceInShoebox ITADSPSIMOVDLSourceInShoebox.cpp )
target_link_libraries( ITADSPSIMOVDLSourceInShoebox ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPSIMOVDLSourceInShoebox )
vista_install( ITADSPSIMOVDLSourceInShoebox )
vista_create_default_info_file( ITADSPSIMOVDLSourceInShoebox )
set_property( TARGET ITADSPSIMOVDLSourceInShoebox PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITADSPThirdOctaveFilterbankTest ITADSPThirdOctaveFilterbankTest.cpp )
target_link_libraries( ITADSPThirdOctaveFilterbankTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPThirdOctaveFilterbankTest )
vista_install( ITADSPThirdOctaveFilterbankTest )
vista_create_default_info_file( ITADSPThirdOctaveFilterbankTest )
set_property( TARGET ITADSPThirdOctaveFilterbankTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_executable( ITASIMOVariableDelayLineBaseTest ITASIMOVariableDelayLineBaseTest.cpp )
target_link_libraries( ITASIMOVariableDelayLineBaseTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITASIMOVariableDelayLineBaseTest )
vista_install( ITASIMOVariableDelayLineBaseTest )
vista_create_default_info_file( ITASIMOVariableDelayLineBaseTest )
set_property( TARGET ITASIMOVariableDelayLineBaseTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP" )
add_subdirectory( "vdl" )
add_subdirectory( "filtering" )
add_subdirectory( "generators" )
if( ITA_DSP_WITH_INTEL_MKL_SUPPORT )
add_subdirectory( "mkl" )
......@@ -115,3 +29,7 @@ endif( )
if( DSPFILTERS_FOUND OR TRUE )
add_subdirectory( "dspfilters" )
endif( )
if( ITA_DSP_WITH_PURE_DATA OR TRUE )
add_subdirectory( "pd" )
endif( )
%Import the output of the c++ file "ITADSPFilterTest", where an impulse
%response of a filter with coefficients b = [1 2 1], a = [1 0.56 0.78] is generated.
%Plot the generated impulse response
ITADSPIIRFilterTest_out = ita_read( 'ITADSPIIRFilterTest_out.wav' );
%ITADSPIIRFilterTest_out = ITADSPIIRFilterTest_out .* 0.004603998475022463843231435021152719855;
%ITADSPIIRFilterTest_out.pt
ITADSPIIRFilterTest_coeffs = jsondecode( fileread( 'ITADSPIIRFilterTest_coeffs.json' ) )
ITADSPIIRFilterTest_out_uninitialized = ita_read( 'ITADSPIIRFilterTest_out_uninitialized.wav' );
ITADSPIIRFilterTest_out_identity = ita_read( 'ITADSPIIRFilterTest_out_identity.wav' );
ITADSPIIRFilterTest_out_zero = ita_read( 'ITADSPIIRFilterTest_out_zero.wav' );
ITADSPIIRFilterTest_out_all = ita_merge( ITADSPIIRFilterTest_out_uninitialized, ITADSPIIRFilterTest_out_identity, ITADSPIIRFilterTest_out_zero );
ITADSPIIRFilterTest_out_all.channelNames = { 'Uninitialized filter engine out', 'Identity filter engine out', 'Zero''d filter entine out' };
ITADSPIIRFilterTest_out_all.signalType = 'energy';
ita_plot_freq( ITADSPIIRFilterTest_out_all, 'axis', [ 100 20e3 -7e3 0 ] )
\ No newline at end of file
{
"order" : 2,
"numerator" : {
"b0" : 1,
"b1" : 2,
"b2" : 1
},
"denominator" : {
"a1" : 1,
"a2" : 0.56,
"a3" : 0.78
},
"design_algorithm" : 0,
"is_ARMA" : true
}
\ No newline at end of file
%% 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();
int main( int, char** )
{
TestThirdOctaveFilterbankIIR();
TestThirdOctaveFilterbankFIR();
return 255;
}
void TestThirdOctaveFilterbankIIR()
{
const int iSampleLength = ( 1 << 17 );
CITAThirdOctaveFilterbank* pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::IIR_BIQUADS_ORDER10 );
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_IIR_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_IIR_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_IIR_SomeBands.wav", &y, g_dSampleRate );
delete pIIRFilterbank;
}
void TestThirdOctaveFilterbankFIR()
{
const int iSampleLength = ( 1 << 17 );
CITAThirdOctaveFilterbank* pIIRFilterbank = CITAThirdOctaveFilterbank::Create( g_dSampleRate, iSampleLength, CITAThirdOctaveFilterbank::FIR_SPLINE_LINEAR_PHASE );
ITASampleBuffer x( iSampleLength );
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_FIR_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_FIR_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_FIR_SomeBands.wav", &y, g_dSampleRate );
delete pIIRFilterbank;
}
%% ITADSPThirdOctaveFilterbankTest evaluation Matlab script
% 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.signalType = 'energy';
iir_ident.channelNames = { 'Identity' };
iir_ident.pf
test_iir_ident = ita_normalize_dat( ita_merge( ita_demosound, ita_convolve( ita_demosound, iir_ident ) ), 'allchannels', true );
test_iir_ident2 = ita_normalize_spk( ita_merge( ita_demosound, ita_convolve( ita_demosound, iir_ident ) ), 'allchannels', true );
test_iir_ident2.play
iir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_Zeros.wav' );
iir_zeros.channelNames = { 'Zeros' };
iir_zeros.pf
%
% fdatool
%
% SOS =
%
% 1.0000 0 -1.0000 1.0000 0.2606 0.8247
% 1.0000 0 -1.0000 1.0000 -0.8148 0.8398
% 1.0000 0 -1.0000 1.0000 -0.5621 0.6067
% 1.0000 0 -1.0000 1.0000 0.0527 0.5859
% 1.0000 0 -1.0000 1.0000 -0.2483 0.5095
%
% G
%
% G =
%
% 0.2843
% 0.2843
% 0.2544
% 0.2544
% 0.2452
% 1.0000
%
ITADSPThirdOctaveFilterbankTest_IIR_Zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_Zeros.wav' );
ITADSPThirdOctaveFilterbankTest_IIR_SomeBands = ita_read( 'ITADSPThirdOctaveFilterbankTest_IIR_SomeBands.wav' );
ITADSPThirdOctaveFilterbankTest_IIR_Zeros.channelNames = { 'Zero' };
ITADSPThirdOctaveFilterbankTest_IIR_Zeros.pf
ITADSPThirdOctaveFilterbankTest_IIR_SomeBands.pf
%% FIR
fir_zeros = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIR_Zeros.wav' );
fir_zeros.signalType = 'power';
fir_unity = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIR_Identity.wav' );
fir_unity.signalType = 'power';
fir_sb = ita_read( 'ITADSPThirdOctaveFilterbankTest_FIR_SomeBands.wav' );
fir_sb.signalType = 'power';
fir_sb_unity = fir_sb.ch( 1 );
for n = 2:fir_sb.nChannels
fir_sb_unity = fir_sb_unity + fir_sb.ch( n );
end
fir_all = ita_merge( fir_zeros, fir_unity, fir_sb );
fir_all.pf
......@@ -3,7 +3,7 @@
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
* RWTH Aachen University, Germany, 2015-2020
*
* ----------------------------------------------------------------
* ____ __________ _______
......@@ -20,13 +20,16 @@
#include <DspFilters/Butterworth.h>
#include <ITAAudioSample.h>
#include <ITAConstants.h>
#include <ITAFiniteImpulseResponse.h>
#include <ITAStopWatch.h>
#include <cassert>
#include <iostream>
#include <cmath>
#include <vector>
#include "VistaTools/VistaRandomNumberGenerator.h"
using namespace std;
using namespace Dsp;
......@@ -44,6 +47,8 @@ void manual_band_pass_1kHz();
void tri_band();
void downcast_test();
void octave_bands();
void live_coeff_change_test();
void very_low_low_pass();
int main( int, char** )
{
......@@ -57,6 +62,10 @@ int main( int, char** )
octave_bands();
live_coeff_change_test();
very_low_low_pass();
return 255;
}
......@@ -191,3 +200,84 @@ void octave_bands()
oHPF.process( iFilterLength, &pf );
oIR.Store( "NativeDspFiltersTest_IR_OctaveBands_HP.wav" );
}
void live_coeff_change_test()
{
float fMidCenterFrequency = 1200.0f;
float fQ = 2;
SimpleFilter< Butterworth::BandPass< iOrder >, iChannels > oBPFDynamic, oBPFStatic;
assert( fMidCenterFrequency > 0 );
oBPFDynamic.setup( iOrder, fSampleRate, fMidCenterFrequency, fMidCenterFrequency / fQ );
oBPFStatic.setup( iOrder, fSampleRate, fMidCenterFrequency, fMidCenterFrequency / fQ );
CITAAudioSample sbCoeffChangeTest( 3, int( fSampleRate * 10 ), fSampleRate );
ITAStopWatch sw;
for( int n = 0; n < sbCoeffChangeTest.GetLength(); n++ )
{
float fRandomSample = VistaRandomNumberGenerator::GetStandardRNG()->GenerateFloat( -1.0f, 1.0f );
sbCoeffChangeTest[ 0 ][ n ] = fRandomSample;
sbCoeffChangeTest[ 1 ][ n ] = fRandomSample;
float fModulatedCenterFrequency = fMidCenterFrequency + 600.0f * sin( ITAConstants::TWO_PI_F * n / float( sbCoeffChangeTest.GetLength() ) * 10 );
sbCoeffChangeTest[ 2 ][ n ] = ( fModulatedCenterFrequency - fMidCenterFrequency ) / 2.0f / fMidCenterFrequency;
sw.start();
oBPFDynamic.setup( iOrder, fSampleRate, fModulatedCenterFrequency, fModulatedCenterFrequency / fQ );
sw.stop(); // takes about 30us on Intel i7 3GHZ from 2015 in debug mode .. not entirely irrelevant :/
float* pfSampleAliasDynamic = &sbCoeffChangeTest[ 0 ][ n ];
oBPFDynamic.process( 1, &pfSampleAliasDynamic );
float* pfSampleAliasStatic = &sbCoeffChangeTest[ 1 ][ n ];
oBPFStatic.process( 1, &pfSampleAliasStatic );
}
cout << "SimpleFilter setup routine stats: " << sw << endl;
sbCoeffChangeTest.Store( "NativeDspFiltersTest_LiveCoeffChangeTest.wav" );
}
void very_low_low_pass()
{
float fVeryLowLowPassFrequencyCuttoff = 0.2f; // 0.2 Hz
SimpleFilter< Butterworth::LowPass< 1 >, 1 > oVLLPF;
oVLLPF.setup( 1, fSampleRate, fVeryLowLowPassFrequencyCuttoff );
CFiniteImpulseResponse oSF( iFilterLength * 30, fSampleRate );
oSF.SetUnitStepFunction( 0.1f );
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 2 * iFilterLength + n ] = 0.2f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 4 * iFilterLength + n ] = 0.3f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 6 * iFilterLength + n ] = 0.4f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 8 * iFilterLength + n ] = 0.5f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 10 * iFilterLength + n ] = 1.0f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 12 * iFilterLength + n ] = 0.0f;
oSF.Store( "NativeDspFiltersTest_SF_VeryLowLowPass_in.wav" );
for( size_t n = 0; n < oSF.GetLength(); n++ )
{
auto pf = oSF.GetData() + n;
oVLLPF.process( 1, &pf );
}
oSF.Store( "NativeDspFiltersTest_SF_VeryLowLowPass_out.wav" );
oVLLPF.reset();
CFiniteImpulseResponse oIR( iFilterLength * 30, fSampleRate );
oIR.SetDirac();
for( size_t n = 0; n < oIR.GetLength(); n++ )
{
auto pf = oIR.GetData() + n;
oVLLPF.process( 1, &pf);
}
oIR.Store( "NativeDspFiltersTest_IR_VeryLowLowPass.wav" );
}
......@@ -42,3 +42,23 @@ IR_Octaves_HP = ita_read( 'NativeDspFiltersTest_IR_OctaveBands_HP.wav' );
IR_Octaves = ita_merge( IR_Octaves_LP, IR_Octaves_BP1, IR_Octaves_BP2,IR_Octaves_BP3, IR_Octaves_BP4, IR_Octaves_BP5, IR_Octaves_BP6, IR_Octaves_BP7, IR_Octaves_BP8, IR_Octaves_HP );
IR_Octaves = ita_merge( IR_Octaves, IR_Octaves.sum );
IR_Octaves.pf
%% Time-variant filtering
LiveCoeffChangeTest = ita_read( 'NativeDspFiltersTest_LiveCoeffChangeTest.wav' );
LiveCoeffChangeTest.pt
LiveCoeffChangeTest.ch(1).play
%% Very low low pass
vllp_ir = ita_read( 'NativeDspFiltersTest_IR_VeryLowLowPass.wav' );
vllp_sf_in = ita_read( 'NativeDspFiltersTest_SF_VeryLowLowPass_in.wav' );
vllp_sf_out = ita_read( 'NativeDspFiltersTest_SF_VeryLowLowPass_out.wav' );
vllp = ita_merge( vllp_ir, vllp_sf_in, vllp_sf_out );
%vllp.channelUnits = 'energy';
vllp.channelNames = { 'Impulse response', 'Step function [in]', 'Step function [out]' };
title 'Very low pass filter (2 Hz) first order zero pole'
vllp.pt
cmake_minimum_required( VERSION 2.8 )
project( ITADSPFilteringTest )
list( APPEND CMAKE_MODULE_PATH "$ENV{VISTA_CMAKE_COMMON}" )
include( VistaCommon )
vista_use_package( ITABase REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITADSP REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITADataSources REQUIRED FIND_DEPENDENCIES )
if( ITA_CORE_LIBS_BUILD_STATIC )
add_definitions( -DITA_BASE_STATIC -DITA_DSP_STATIC -DITA_DATA_SOURCES_STATIC -DITA_CONVOLUTION_STATIC )
endif( )
if( ITA_VISTA_BUILD_STATIC )
add_definitions( -DVISTATOOLS_STATIC -DVISTABASE_STATIC -DVISTAMATH_STATIC -DVISTAASPECTS_STATIC -DVISTAINTERPROCCOMM_STATIC )
endif( )
add_executable( ITADSPThirdOctaveFilterbankTest ITADSPThirdOctaveFilterbankTest.cpp )
target_link_libraries( ITADSPThirdOctaveFilterbankTest ${VISTA_USE_PACKAGE_LIBRARIES} )
vista_configure_app( ITADSPThirdOctaveFilterbankTest )
vista_install( ITADSPThirdOctaveFilterbankTest )
vista_create_default_info_file( ITADSPThirdOctaveFilterbankTest )
set_property( TARGET ITADSPThirdOctaveFilterbankTest PROPERTY FOLDER "ITACoreLibs/Tests/ITADSP/filtering" )
add_subdirectory( "iir" )
#include <ITAThirdOctaveFilterbank.h>
#include <ITAAudiofileWriter.h>
#include <ITANumericUtils.h>
#include <ITAStopWatch.h>
#include <ITAStringUtils.h>
#include <ITASampleBuffer.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>
#ifdef WITH_JSON_SUPPORT
#include <ITABase/UtilsJSON.h>
#endif
#include <iostream>
using namespace std;
using namespace ITABase;
const int g_iSampleLength = ( 1 << 17 );
const int g_iBlockLength = 128;
const double g_dSampleRate = 44100;
const int g_iFIRGenLength = int( 4 * ceil( g_dSampleRate / CThirdOctaveMagnitudeSpectrum::GetCenterFrequencies()[ 0 ] ) );
ITAStopWatch sw;
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 oReflection;
ITABase::CThirdOctaveGainMagnitudeSpectrum oDiffraction;
void Run( CITAThirdOctaveFilterbank* pFilterBank, std::string sExportName );
int main( int, char** )
{
// Initialize
sbDirac[ 0 ] = 1.0f;
oZero.SetZero();
#ifdef WITH_JSON_SUPPORT
ITABase::Utils::JSON::Export( &oZero, "ZeroSpectrum.json" );
#endif
oIdentity.SetIdentity();
oIdentity.Multiply( 1.000001f );
#ifdef WITH_JSON_SUPPORT
ITABase::Utils::JSON::Export( &oIdentity, "IdentitySpectrum.json" );
#endif
oHighPass.SetIdentity();