Implemented and tested artificial reverb with 1, 3 or 10 reverberation time values

parent af276202
......@@ -65,6 +65,10 @@
#include <tbb/atomic.h>
#include <tbb/concurrent_queue.h>
#include <DspFilters/Dsp.h>
using namespace ITAConstants;
const double g_dMinReverberationTime = 0.25f; // Minimum allowed reverberation time [s]
const double g_dMinRoomVolume = 3 * 2 * 1; // Minimum allowed room volume [m^3]
const double g_dMinRoomSurfaceArea = 2 * ( 3 * 2 + 3 + 2 ); // Minimum allowed room surface area [m^2]
......@@ -210,6 +214,67 @@ private:
std::atomic< int > m_iState;
};
const int iOrder = 10;
const int iChannels = 1;
int iFilterLength = ( int ) pow( 2, 12 );
class CLowPassFilter : public CVABinauralArtificialReverbAudioRenderer::CFilterBase
{
public:
inline CLowPassFilter( float fLowerFrequencyCuttoff, float fSampleRate )
{
oFilterEngine.setup( iOrder, fSampleRate, fLowerFrequencyCuttoff );
};
inline void Process( int iNumSamples, float* pfInOut, bool bAutoReset = true )
{
if( bAutoReset )
oFilterEngine.reset();
oFilterEngine.process( iNumSamples, &pfInOut );
};
Dsp::SimpleFilter< Dsp::Butterworth::LowPass< iOrder >, iChannels > oFilterEngine;
};
class CHighPassFilter : public CVABinauralArtificialReverbAudioRenderer::CFilterBase
{
public:
inline CHighPassFilter( float fLowerFrequencyCuttoff, float fSampleRate )
{
oFilterEngine.setup( iOrder, fSampleRate, fLowerFrequencyCuttoff );
};
inline void Process( int iNumSamples, float* pfInOut, bool bAutoReset = true )
{
if( bAutoReset )
oFilterEngine.reset();
oFilterEngine.process( iNumSamples, &pfInOut );
};
Dsp::SimpleFilter< Dsp::Butterworth::HighPass< iOrder >, iChannels > oFilterEngine;
};
class CBandPassFilter : public CVABinauralArtificialReverbAudioRenderer::CFilterBase
{
public:
inline CBandPassFilter( float fCenterFrequency, float fBandwidth, float fSampleRate )
{
oFilterEngine.setup( iOrder, fSampleRate, fCenterFrequency, fBandwidth );
};
inline void Process( int iNumSamples, float* pfInOut, bool bAutoReset = true )
{
if( bAutoReset )
oFilterEngine.reset();
oFilterEngine.process( iNumSamples, &pfInOut );
};
Dsp::SimpleFilter< Dsp::Butterworth::BandPass< iOrder >, iChannels > oFilterEngine;
};
// --= Renderer =--
......@@ -251,7 +316,24 @@ void CVABinauralArtificialReverbAudioRenderer::Init( const CVAStruct& oArgs )
{
CVAConfigInterpreter conf( oArgs );
conf.OptNumber( "ReverberationTime", m_dReverberationTime, 0.71f );
std::string sReverberationTimes;
conf.OptString( "ReverberationTimes", sReverberationTimes );
m_vdReverberationTimes = StringToDoubleVec( sReverberationTimes );
std::string sReverberationTime;
conf.OptString( "ReverberationTime", sReverberationTime );
if( sReverberationTimes.empty() && !sReverberationTime.empty() )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Deprecated ReverberationTime used in configuration, please switch to ReverberationTimes and provide a list of e.g. 3 values (low, mid, high). Approximating reverb." );
double dReverberationTime = StringToDouble( sReverberationTime );
m_vdReverberationTimes = { dReverberationTime * sqrt( 2 ), dReverberationTime, dReverberationTime / sqrt( 2 ) };
}
DesignEQFilters();
conf.OptNumber( "RoomVolume", m_dRoomVolume, 10.0f *10.0f * 2.0f );
conf.OptNumber( "RoomSurfaceArea", m_dRoomSurfaceArea, 4 * ( 10.0f * 2.0f ) + 2 * ( 2.0f * 2.0f ) );
conf.OptNumber( "SoundPowerCorrectionFactor", m_dSoundPowerCorrectionFactor, 0.05 );
......@@ -264,7 +346,7 @@ void CVABinauralArtificialReverbAudioRenderer::Init( const CVAStruct& oArgs )
conf.OptNumber( "TimeSlotResolution", m_dTimeSlotResolution, 0.005f );
conf.OptNumber( "MaxReflectionDensity", m_dMaxReflectionDensity, 12000.0f );
conf.OptNumber( "ScatteringCoefficient", m_dScatteringCoefficient, 0.1f );
std::string sVLDInterpolationAlgorithm;
conf.OptString( "VDLSwitchingAlgorithm", sVLDInterpolationAlgorithm, "linear" );
sVLDInterpolationAlgorithm = toLowercase( sVLDInterpolationAlgorithm );
......@@ -282,10 +364,14 @@ void CVABinauralArtificialReverbAudioRenderer::Init( const CVAStruct& oArgs )
else
ITA_EXCEPT1( INVALID_PARAMETER, "Unrecognized interpolation algorithm '" + sVLDInterpolationAlgorithm + "' in BinauralFreefieldAudioRendererConfig" );
if( m_dReverberationTime < g_dMinReverberationTime )
for( size_t i = 0; i < m_vdReverberationTimes.size(); i++ )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << m_dReverberationTime << " s too small, minimum is " << g_dMinReverberationTime << " s." );
m_dReverberationTime = g_dMinReverberationTime;
auto dReverberationTime = m_vdReverberationTimes[ i ];
if( dReverberationTime < g_dMinReverberationTime )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << dReverberationTime << " s too small, minimum is " << g_dMinReverberationTime << " s." );
m_vdReverberationTimes[ i ] = g_dMinReverberationTime;
}
}
if( m_dRoomSurfaceArea < g_dMinRoomSurfaceArea )
......@@ -303,6 +389,37 @@ void CVABinauralArtificialReverbAudioRenderer::Init( const CVAStruct& oArgs )
return;
}
void CVABinauralArtificialReverbAudioRenderer::DesignEQFilters()
{
if( m_vdReverberationTimes.size() == 3 )
{
m_vpFilter.clear();
m_vpFilter.push_back( std::make_shared< CLowPassFilter >( 330.0f, GetSampleRate() ) );
m_vpFilter.push_back( std::make_shared< CBandPassFilter >( 2300.0f, 4000.0f, GetSampleRate() ) );
m_vpFilter.push_back( std::make_shared< CHighPassFilter >( 3900.0f, GetSampleRate() ) );
}
else if( m_vdReverberationTimes.size() == 10 )
{
m_vpFilter.clear();
float fLowPassFrequencyCuttoff = OCTAVE_CENTER_FREQUENCIES_ISO_F[ 0 ] * sqrt( 2.0f );
m_vpFilter.push_back( std::make_shared< CLowPassFilter >( fLowPassFrequencyCuttoff, GetSampleRate() ) );
for( size_t i = 1; i < OCTAVE_CENTER_FREQUENCIES_ISO_F.size() - 1; i++ )
{
float fMidCenterFrequency = OCTAVE_CENTER_FREQUENCIES_ISO_F[ i ];
float fMidBandWidth = ( OCTAVE_CENTER_FREQUENCIES_ISO_F[ i + 1 ] - OCTAVE_CENTER_FREQUENCIES_ISO_F[ i - 1 ] ) / 2.15f;
m_vpFilter.push_back( std::make_shared< CBandPassFilter >( fMidCenterFrequency, fMidBandWidth, GetSampleRate() ) );
}
float fHighPassFrequencyCuttoff = OCTAVE_CENTER_FREQUENCIES_ISO_F[ 0 ] / sqrt( 2.0f );
m_vpFilter.push_back( std::make_shared< CHighPassFilter >( fHighPassFrequencyCuttoff, GetSampleRate() ) );
}
else
{
VA_ERROR( "BinauralArtificialReverbAudioRenderer", "Current implementation only accepts 3 reverb times for low, mid & high frequencies or octave resolution (8 values)" );
}
}
CVAObjectInfo CVABinauralArtificialReverbAudioRenderer::GetObjectInfo() const
{
CVAObjectInfo oInfo;
......@@ -332,7 +449,7 @@ CVAStruct CVABinauralArtificialReverbAudioRenderer::CallObject( const CVAStruct&
if( sValue == "STATUS" )
{
VA_PRINT( "RoomReverberationTime: " << m_dReverberationTime << " seconds" );
VA_PRINT( "RoomReverberationTimes: " << DoubleVecToString( m_vdReverberationTimes ) << " seconds" );
VA_PRINT( "MaxReverbFilterLengthSamples: " << m_iMaxReverbFilterLengthSamples << " samples (" << m_iMaxReverbFilterLengthSamples / GetSampleRate() << " seconds)" );
VA_PRINT( "RoomVolume: " << m_dRoomVolume << " m^3" );
VA_PRINT( "RoomSurfaceArea: " << m_dRoomSurfaceArea << " m^2" );
......@@ -356,25 +473,28 @@ CVAStruct CVABinauralArtificialReverbAudioRenderer::CallObject( const CVAStruct&
( pStruct->GetDatatype() == CVAStructValue::BOOL ) )
{
double dRoomReverberationTime = *pStruct;
if( dRoomReverberationTime < g_dMinReverberationTime )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << m_dReverberationTime << "s too small, minimum is " << g_dMinReverberationTime << "s." );
m_dReverberationTime = g_dMinReverberationTime;
}
else
{
m_dReverberationTime = dRoomReverberationTime;
}
VA_INFO( "BinauralArtificialReverbAudioRenderer", "Set reverberation time to " << m_dReverberationTime );
for( size_t i = 0; i < m_vdReverberationTimes.size(); i++ )
if( dRoomReverberationTime < g_dMinReverberationTime )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << dRoomReverberationTime << "s too small, minimum is " << g_dMinReverberationTime << "s." );
m_vdReverberationTimes[ i ] = g_dMinReverberationTime;
}
else
{
m_vdReverberationTimes[ i ] = dRoomReverberationTime;
}
VA_INFO( "BinauralArtificialReverbAudioRenderer", "Set reverberation time to " << dRoomReverberationTime );
if( dRoomReverberationTime * GetSampleRate() > m_iMaxReverbFilterLengthSamples )
VA_WARN( "BinauralArtificialReverbAudioRenderer",
"Reverberation time greater than target filter length, reverb will be truncated at "
<< m_iMaxReverbFilterLengthSamples / GetSampleRate() << " seconds"
);
}
if( m_dReverberationTime * GetSampleRate() > m_iMaxReverbFilterLengthSamples )
VA_WARN( "BinauralArtificialReverbAudioRenderer",
"Reverberation time greater than target filter length, reverb will be truncated at "
<< m_iMaxReverbFilterLengthSamples / GetSampleRate() << " seconds"
);
m_bForceARUpdateOnce = true;
DesignEQFilters();
UpdateArtificialReverbPaths();
return oReturn;
......@@ -462,7 +582,9 @@ CVAStruct CVABinauralArtificialReverbAudioRenderer::CallObject( const CVAStruct&
std::string sValue = toUppercase( *pStruct );
if( sValue == "ROOMREVERBERATIONTIME" )
oReturn[ "Return" ] = m_dReverberationTime;
oReturn[ "Return" ] = m_vdReverberationTimes[ 0 ];
else if( sValue == "ROOMREVERBERATIONTIMES" )
oReturn[ "Return" ] = DoubleVecToString( m_vdReverberationTimes );
else if( sValue == "ROOMVOLUME" )
oReturn[ "Return" ] = m_dRoomVolume;
else if( sValue == "ROOMSURFACEAREA" )
......@@ -567,7 +689,7 @@ void CVABinauralArtificialReverbAudioRenderer::UpdateScene( CVASceneState* pNewS
UpdateArtificialReverbPaths( false ); // Unforced update validation
if( m_pCurSceneState )
if( m_pCurSceneState )
m_pCurSceneState->RemoveReference();
m_pCurSceneState = m_pNewSceneState;
......@@ -581,6 +703,62 @@ ITADatasource* CVABinauralArtificialReverbAudioRenderer::GetOutputDatasource()
return this;
}
CVAStruct CVABinauralArtificialReverbAudioRenderer::GetParameters( const CVAStruct& oInArgs ) const
{
CVAStruct oRet;
std::vector< float > vf;
for( auto d : m_vdReverberationTimes )
vf.push_back( ( float ) d );
oRet[ "room_reverberation_times" ] = CVAStructValue( ( void* ) &vf[ 0 ], vf.size() * sizeof( float ) );
return oRet;
}
void CVABinauralArtificialReverbAudioRenderer::SetParameters( const CVAStruct& oInArgs )
{
if( oInArgs.HasKey( "room_reverberation_times" ) )
{
const auto& oReverbTimes( oInArgs[ "room_reverberation_times" ] );
if( oReverbTimes.IsDouble() )
{
double dReverbTime = oReverbTimes;
if( dReverbTime < g_dMinReverberationTime )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << dReverbTime << " s too small, minimum is " << g_dMinReverberationTime << " s." );
dReverbTime = g_dMinReverberationTime;
}
m_vdReverberationTimes = { dReverbTime * sqrt( 2 ), dReverbTime, dReverbTime / sqrt( 2 ) };
DesignEQFilters();
UpdateArtificialReverbPaths( true );
}
else if( oReverbTimes.IsData() )
{
int iNumReverbTimes = int( oReverbTimes.GetDataSize() / sizeof( float ) );
if( iNumReverbTimes == 3 || iNumReverbTimes == 10 )
{
m_vdReverberationTimes.resize( iNumReverbTimes );
const auto pfReverbTimes = ( const float* ) oReverbTimes.GetData();
for( size_t i = 0; i < iNumReverbTimes; i++ )
if( pfReverbTimes[ i ] < g_dMinReverberationTime )
{
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Requested reverberation time of " << pfReverbTimes[ i ] << " s too small, minimum is " << g_dMinReverberationTime << " s." );
m_vdReverberationTimes[ i ] = g_dMinReverberationTime;
}
else
m_vdReverberationTimes[ i ] = pfReverbTimes[ i ];
}
else
VA_EXCEPT2( INVALID_PARAMETER, "Reverberation times must containt 1, 3 (low, mid, high) or 8 (octave) values" );
DesignEQFilters();
UpdateArtificialReverbPaths( true );
}
}
}
void CVABinauralArtificialReverbAudioRenderer::ManageArtificialReverbPaths( const CVASceneStateDiff* pDiff )
{
// ber aktuelle Pfade iterieren und gelschte markieren
......@@ -636,7 +814,7 @@ void CVABinauralArtificialReverbAudioRenderer::ManageArtificialReverbPaths( cons
}
// ber aktuelle Quellen und Hrer iterieren und gelschte entfernen
for( std::vector<int>::const_iterator cit = pDiff->viDelSoundSourceIDs.begin(); cit != pDiff->viDelSoundSourceIDs.end(); ++cit )
for( std::vector<int>::const_iterator cit = pDiff->viDelSoundSourceIDs.begin(); cit != pDiff->viDelSoundSourceIDs.end(); ++cit )
{
DeleteSource( *cit );
}
......@@ -1226,14 +1404,15 @@ void CVABinauralArtificialReverbAudioRenderer::UpdateArtificialReverbFilter( Lis
// local copies
const double dRoomSurfaceArea = m_dRoomSurfaceArea;
const double dRoomVolume = m_dRoomVolume;
const double dReverberationTime = m_dReverberationTime;
assert( dReverberationTime > 0.0f );
const double dSoundPowerCorrectionFactor = m_dSoundPowerCorrectionFactor;
if( ceil( dReverberationTime * GetSampleRate() ) > m_iMaxReverbFilterLengthSamples )
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Target filter length smaller than requested reverberation time, truncating." )
for( auto dReverberationTime : m_vdReverberationTimes )
{
assert( dReverberationTime > 0.0f );
ITASampleFrame sfBRIR( 2, m_iMaxReverbFilterLengthSamples, true ); // 0=L, 1=R, TODO member var?
if( ceil( dReverberationTime * GetSampleRate() ) > m_iMaxReverbFilterLengthSamples )
VA_WARN( "BinauralArtificialReverbAudioRenderer", "Target filter length smaller than requested reverberation time, truncating." )
}
if( pListener->pNewHRIR != nullptr )
{
......@@ -1243,169 +1422,193 @@ void CVABinauralArtificialReverbAudioRenderer::UpdateArtificialReverbFilter( Lis
}
CVADirectivityDAFFHRIR* pHRIR = ( CVADirectivityDAFFHRIR* ) pListener->pCurHRIR;
if( !pHRIR )
{
VA_ERROR( "BinauralArtificialReverbAudioRenderer", "Could not cast receiver HRIR to type CVADirectivityDAFFHRIR, please use a DAFF file with IR content." );
}
int iHRIRLength = pHRIR->GetProperties()->iFilterLength;
ITASampleFrame sfTempHRIR( 2, iHRIRLength, true ); // TODO member?
// Some prior assertions
assert( dRoomSurfaceArea > 0.0f );
assert( dRoomVolume > 0.0f );
assert( dReverberationTime > 0.0f );
assert( m_dTimeSlotResolution > 0.0f );
// START LAS -----------------------------------------------------------------------------------------
srand( 667 ); // fixed poisson sequence
const double dAverageLengthRoom = dRoomVolume / dRoomSurfaceArea;
const double dMeanfreePath = 4 * dAverageLengthRoom;
if( m_vdReverberationTimes.size() != m_vpFilter.size() )
{
VA_ERROR( "BinauralArtificialReverbAudioRenderer", "EQ filters not equalized, aborting reverb calculation" );
}
const double dScatterReflectionFactor = 0.75f + ( m_dScatteringCoefficient / 4.0f );
// START LAS -----------------------------------------------------------------------------------------
const double dSamplesPerTimeSlot = GetSampleRate() * m_dTimeSlotResolution;
ITASampleFrame sfBRIR_temp( 2, m_iMaxReverbFilterLengthSamples, true ); // 0=L, 1=R, TODO member var? Temporal buffer per frequency
ITASampleFrame sfBRIR( 2, m_iMaxReverbFilterLengthSamples, true ); // Accumulated reverbs
// get time of last image source
const double dTimeLastImageSource = dMeanfreePath / m_pCore->oHomogeneousMedium.dSoundSpeed;
for( size_t r = 0; r < m_vdReverberationTimes.size(); r++ )
{
const double& dReverberationTime( m_vdReverberationTimes[ r ] );
assert( dReverberationTime > 0.0f );
// calc critical distances and initial reverberation energies
const double dTimeConst = -13.816f / dReverberationTime;
assert( dTimeConst != 0.0f );
srand( 667 ); // fixed poisson sequence
const double dAreaNormIntegral = ( exp( dTimeConst * double( m_iMaxReverbFilterLengthSamples ) ) - exp( dTimeConst * dTimeLastImageSource ) ) / dTimeConst;
assert( dAreaNormIntegral > 0.0f );
sfBRIR_temp.zero();
// calc critical distance
double dAverageEquivalentAbsorptionArea = 0.163f * dRoomVolume / dReverberationTime;
double dCriticalDistance = sqrt( dAverageEquivalentAbsorptionArea / ( 16.0f * ITAConstants::PI_D ) );
assert( dCriticalDistance > 0.0f );
const double dAverageLengthRoom = dRoomVolume / dRoomSurfaceArea;
const double dMeanFreePath = 4 * dAverageLengthRoom;
// calc initial energy (formula see diploma thesis Lukas Aspck)
// correction factor as envelope curve has approx 3 times the energy the synth ir has
double dInitialEnergy = 3.0f / pow( dCriticalDistance, 2 ) / m_pCore->oHomogeneousMedium.dSoundSpeed / dAreaNormIntegral;
assert( dInitialEnergy > 0.0f );
const double dScatterReflectionFactor = 0.75f + ( m_dScatteringCoefficient / 4.0f );
double dCurrentScaleFactor; // PROBLEM
const double dSamplesPerTimeSlot = GetSampleRate() * m_dTimeSlotResolution;
//############ CREATE NOISE SEQUENCES WITH TIME DOMAIN WEIGHTING ############//
// get time of last image source
const double dTimeLastImageSource = dMeanFreePath / m_pCore->oHomogeneousMedium.dSoundSpeed;
// new loop over fixed time slots
int iBlockLength = int( dSamplesPerTimeSlot );
int iNumberOfTimeSlots = int( m_iMaxReverbFilterLengthSamples / double( iBlockLength ) ) + 1;
// calc critical distances and initial reverberation energies
const double dTimeConst = -13.816f / dReverberationTime;
assert( dTimeConst != 0.0f );
for( int i = 0; i < int( iNumberOfTimeSlots ); i++ )
{
int& iTimeSlot( i );
// calculate reflection density for current block
const double dAreaNormIntegral = ( exp( dTimeConst * double( m_iMaxReverbFilterLengthSamples ) ) - exp( dTimeConst * dTimeLastImageSource ) ) / dTimeConst;
assert( dAreaNormIntegral > 0.0f );
// calculates reflection density according to kuttruffs equation
double dTime = ( iTimeSlot * iBlockLength + 1 ) / GetSampleRate();
double dCurrentReflectionDensity = 4.0f * ITAConstants::PI_D * m_pCore->oHomogeneousMedium.dSoundSpeed* m_pCore->oHomogeneousMedium.dSoundSpeed* m_pCore->oHomogeneousMedium.dSoundSpeed * dTime * dTime / dRoomVolume;
// calc critical distance
double dAverageEquivalentAbsorptionArea = 0.163f * dRoomVolume / dReverberationTime;
double dCriticalDistance = sqrt( dAverageEquivalentAbsorptionArea / ( 16.0f * ITAConstants::PI_D ) );
assert( dCriticalDistance > 0.0f );
// if the density exceeds maximumReflectionDensity it will be kept constant at this value
if( dCurrentReflectionDensity > m_dMaxReflectionDensity )
dCurrentReflectionDensity = m_dMaxReflectionDensity;
// calc initial energy (formula see diploma thesis Lukas Aspck)
// correction factor as envelope curve has approx 3 times the energy the synth ir has
double dInitialEnergy = 3.0f / pow( dCriticalDistance, 2 ) / m_pCore->oHomogeneousMedium.dSoundSpeed / dAreaNormIntegral;
assert( dInitialEnergy > 0.0f );
// calculate reflection density per sample (average scattering also taken into account)
double dCurrentReflectionDensityPerSample = dCurrentReflectionDensity * dScatterReflectionFactor / GetSampleRate();
double dCurrentScaleFactor; // PROBLEM
// calculate number of diracs in current timeslot
int iNumberDiracsTimeSlot = int( dCurrentReflectionDensityPerSample * iBlockLength );
//############ CREATE NOISE SEQUENCES WITH TIME DOMAIN WEIGHTING ############//
// calculate current Poisson Scale Factor in TimeSlot
double fCurrentScaleFactorConstantPoissonEnergy = ( iNumberDiracsTimeSlot > 1 ) ? sqrt( 1.0f / double( iNumberDiracsTimeSlot ) ) : 1.0f;
// new loop over fixed time slots
int iBlockLength = int( dSamplesPerTimeSlot );
int iNumberOfTimeSlots = int( m_iMaxReverbFilterLengthSamples / double( iBlockLength ) ) + 1;
// set all diracs in current time slot
for( int j = 0; j < iNumberDiracsTimeSlot; j++ )
for( int i = 0; i < int( iNumberOfTimeSlots ); i++ )
{
// choose random position within timeSlot
int iCurrentDiracPosition = iTimeSlot * iBlockLength + int( iBlockLength * rand() / double( RAND_MAX ) );
// calculate absolute position (needed to choose if audible or not)
double dAbsoluteTimePosition = iCurrentDiracPosition / GetSampleRate(); // + fDelayTimeDirectSound;
// get random angles
double dAziAngle = rand() / RAND_MAX * 360.0f;
//double dAziAngleVistaImpl = VistaRandomNumberGenerator::GetStandardRNG()->GenerateDouble( 0.0f, 360.0f );
double dEleAngle = ( 1.0f - 2.0f * ( rand() / RAND_MAX ) ) * 90.0f;
int& iTimeSlot( i );
// calculate reflection density for current block
// get respective HRIRs
assert( iHRIRLength == sfTempHRIR.length() );
sfTempHRIR.zero(); // nicht ntig, wenn assert gilt
int iIdx = -1;
pHRIR->GetNearestNeighbour( float( dAziAngle ), float( dEleAngle ), &iIdx ); // view direction
pHRIR->GetHRIRByIndex( &sfTempHRIR, iIdx, 1.0f );
// calculates reflection density according to kuttruffs equation
double dTime = ( iTimeSlot * iBlockLength + 1 ) / GetSampleRate();
double dCurrentReflectionDensity = 4.0f * ITAConstants::PI_D * m_pCore->oHomogeneousMedium.dSoundSpeed* m_pCore->oHomogeneousMedium.dSoundSpeed* m_pCore->oHomogeneousMedium.dSoundSpeed * dTime * dTime / dRoomVolume;
// get sign for dirac pulse
int iRandomSign = 2 * ( rand() % 2 ) - 1;
// if the density exceeds maximumReflectionDensity it will be kept constant at this value
if( dCurrentReflectionDensity > m_dMaxReflectionDensity )
dCurrentReflectionDensity = m_dMaxReflectionDensity;
// Beispiel mit Vista:
//double dRandomSignVistaImpl = ( int( VistaRandomNumberGenerator::GetStandardRNG()->GenerateInt31() ) > 0 ) ? 1.0f : -1.0f;
// calculate reflection density per sample (average scattering also taken into account)
double dCurrentReflectionDensityPerSample = dCurrentReflectionDensity * dScatterReflectionFactor / GetSampleRate();
bool bEarlyReflectionsPhase = ( dAbsoluteTimePosition < dTimeLastImageSource );
bool bLateReflectionPhase = ( dAbsoluteTimePosition > dTimeLastImageSource ) &&
( iCurrentDiracPosition + iHRIRLength < m_iMaxReverbFilterLengthSamples ); // PROBLEM, warum dieses Abbruchkrit. hier?
// calculate number of diracs in current timeslot
int iNumberDiracsTimeSlot = int( dCurrentReflectionDensityPerSample * iBlockLength );
// calculate current Poisson Scale Factor in TimeSlot
double fCurrentScaleFactorConstantPoissonEnergy = ( iNumberDiracsTimeSlot > 1 ) ? sqrt( 1.0f / double( iNumberDiracsTimeSlot ) ) : 1.0f;
if( bEarlyReflectionsPhase )
// set all diracs in current time slot
for( int j = 0; j < iNumberDiracsTimeSlot; j++ )
{
// phase during early reflections
// choose random position within timeSlot
int iCurrentDiracPosition = iTimeSlot * iBlockLength + int( iBlockLength * rand() / double( RAND_MAX ) );
// get Scatter coefficient and calc number of early scattered impulses
int iNumberOfScatteredPulses = int( m_dScatteringCoefficient * 10 );
// calculate absolute position (needed to choose if audible or not)
double dAbsoluteTimePosition = iCurrentDiracPosition / GetSampleRate(); // + fDelayTimeDirectSound;
assert( dTimeLastImageSource > 0.0f );
// calc energy of current Sample: apply weighting function (phase ends at dTimeLastImageSource)
double dSpectralEnergyThisSample =
dInitialEnergy // energy direct sound
* ( dAbsoluteTimePosition / dTimeLastImageSource ) // weighting function (IS Energy)
* exp( -13.816f * dAbsoluteTimePosition / dReverberationTime ); // reverb
// get random angles
double dAziAngle = rand() / RAND_MAX * 360.0f;
//double dAziAngleVistaImpl = VistaRandomNumberGenerator::GetStandardRNG()->GenerateDouble( 0.0f, 360.0f );
double dEleAngle = ( 1.0f - 2.0f * ( rand() / RAND_MAX ) ) * 90.0f;
// insert scaled dirac in each frequency band
dCurrentScaleFactor = iRandomSign * sqrt( dSpectralEnergyThisSample );
// get respective HRIRs
assert( iHRIRLength == sfTempHRIR.length() );
sfTempHRIR.zero(); // nicht ntig, wenn assert gilt
int iIdx = -1;
pHRIR->GetNearestNeighbour( float( dAziAngle ), float( dEleAngle ), &iIdx ); // view direction
pHRIR->GetHRIRByIndex( &sfTempHRIR, iIdx, 1.0f );
// insert scattered energy
for( int k = 0; k < iNumberOfScatteredPulses; k++ )
{
int iScatterDelayInSamples = int( ( ( k + 1 ) / double( iNumberOfScatteredPulses ) ) * ( 2.0f * dAverageLengthRoom / m_pCore->oHomogeneousMedium.dSoundSpeed * GetSampleRate() ) + rand() % 30 );
// get sign for dirac pulse
int iRandomSign = 2 * ( rand() % 2 ) - 1;
// Beispiel mit Vista:
//double dRandomSignVistaImpl = ( int( VistaRandomNumberGenerator::GetStandardRNG()->GenerateInt31() ) > 0 ) ? 1.0f : -1.0f;
assert( dSpectralEnergyThisSample > 0.0f );
double dCurrentScatterFactor = iRandomSign *
( m_dScatteringCoefficient / pow( 2.0f, k ) ) *
sqrt( dSpectralEnergyThisSample ) *
fCurrentScaleFactorConstantPoissonEnergy;
bool bEarlyReflectionsPhase = ( dAbsoluteTimePosition < dTimeLastImageSource );
bool bLateReflectionPhase = ( dAbsoluteTimePosition > dTimeLastImageSource ) &&
( iCurrentDiracPosition + iHRIRLength < m_iMaxReverbFilterLengthSamples ); // PROBLEM, warum dieses Abbruchkrit. hier?
int iDestOffset = iCurrentDiracPosition + iScatterDelayInSamples;
if( iDestOffset + iHRIRLength > sfBRIR.length() )
continue;
//float* pfDestCh1 = sfBRIR[ 0 ].data();
sfBRIR.muladd_frame( sfTempHRIR, float( dCurrentScatterFactor ), 0, iDestOffset, iHRIRLength );
if( bEarlyReflectionsPhase )
{
// phase during early reflections
// get Scatter coefficient and calc number of early scattered impulses
int iNumberOfScatteredPulses = int( m_dScatteringCoefficient * 10 );
assert( dTimeLastImageSource > 0.0f );
// calc energy of current Sample: apply weighting function (phase ends at dTimeLastImageSource)
double dSpectralEnergyThisSample =
dInitialEnergy // energy direct sound
* ( dAbsoluteTimePosition / dTimeLastImageSource ) // weighting function (IS Energy)
* exp( -13.816f * dAbsoluteTimePosition / dReverberationTime ); // reverb
// insert scaled dirac in each frequency band
dCurrentScaleFactor = iRandomSign * sqrt( dSpectralEnergyThisSample );
// insert scattered energy
for( int k = 0; k < iNumberOfScatteredPulses; k++ )
{
int iScatterDelayInSamples = int( ( ( k + 1 ) / double( iNumberOfScatteredPulses ) ) * ( 2.0f * dAverageLengthRoom / m_pCore->oHomogeneousMedium.dSoundSpeed * GetSampleRate() ) + rand() % 30 );
assert( dSpectralEnergyThisSample > 0.0f );
double dCurrentScatterFactor = iRandomSign *
( m_dScatteringCoefficient / pow( 2.0f, k ) ) *
sqrt( dSpectralEnergyThisSample ) *
fCurrentScaleFactorConstantPoissonEnergy;
int iDestOffset = iCurrentDiracPosition + iScatterDelayInSamples;
if( iDestOffset + iHRIRLength > sfBRIR_temp.length() )
continue;
//float* pfDestCh1 = sfBRIR[ 0 ].data();
sfBRIR_temp.muladd_frame( sfTempHRIR, float( dCurrentScatterFactor ), 0, iDestOffset, iHRIRLength );
}
}