Skip to content

GitLab

I
ITAPropagationModels
Commit cd2695d0 authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files
Options

Appyling new methods from ITAGeo (coordinates, directivities)

parent fbd1cc2a
Hide whitespace changes
Inline Side-by-side
Showing with 102 additions and 132 deletions
include/ITAPropagationModels/FilterEngine.h
View file @ cd2695d0
......@@ -73,7 +73,7 @@ namespace ITAPropagationModels
//! Standard constructor for a filter engine
CFilterEngine();
CFilterEngine( const ITAGeo::ECoordinateSystemConvention iConvention = ITAGeo::ECoordinateSystemConvention::ISO );
//! Default destructor of a filter engine
......@@ -145,6 +145,8 @@ namespace ITAPropagationModels
std::shared_ptr< const ITAGeo::Material::IManager > m_pMaterialManager;
std::shared_ptr< const ITAGeo::Directivity::IManager > m_pDirectivityManager;
ITAGeo::ECoordinateSystemConvention m_iCoordinateSystemConvention;
};
}
......
src/ITAPropagationModels/FilterEngine.cpp
View file @ cd2695d0
......@@ -28,8 +28,9 @@
using namespace ITAPropagationModels;
CFilterEngine::CFilterEngine()
:m_pMaterialManager( nullptr )
CFilterEngine::CFilterEngine(const ITAGeo::ECoordinateSystemConvention iCSC /* = ISO */)
: m_pMaterialManager(nullptr)
, m_iCoordinateSystemConvention(iCSC)
{
}
......@@ -37,60 +38,60 @@ CFilterEngine::~CFilterEngine()
{
}
int CFilterEngine::GetNumSensorChannels( const ITAGeo::CPropagationPathList& oPathList )
int CFilterEngine::GetNumSensorChannels(const ITAGeo::CPropagationPathList& oPathList)
{
//Check for correct structur of oPathList
if( !ValidateSameSensor( oPathList ) )
ITA_EXCEPT1( INVALID_PARAMETER, "The propagation path list has multiple sensor anchors or last anchor of paths is not a sensor." );
if (!ValidateSameSensor(oPathList))
ITA_EXCEPT1(INVALID_PARAMETER, "The propagation path list has multiple sensor anchors or last anchor of paths is not a sensor.");
//Cast the sensor that is always the last propagation anchor of a path
auto pSensor = std::dynamic_pointer_cast< ITAGeo::CSensor >( oPathList[ 0 ][ oPathList[ 0 ].size() - 1 ] );
auto pSensor = std::dynamic_pointer_cast<ITAGeo::CSensor>(oPathList[0][oPathList[0].size() - 1]);
return pSensor->GetNumChannels();
}
bool CFilterEngine::ValidateSameSensor( const ITAGeo::CPropagationPathList& oPathList )
bool CFilterEngine::ValidateSameSensor(const ITAGeo::CPropagationPathList& oPathList)
{
shared_ptr<ITAGeo::CSensor> pPeviousSensor = nullptr;
// Last anchor of the all paths must share the same sensor (pointer)
for( auto& oPath : oPathList )
for (auto& oPath : oPathList)
{
if( oPath.empty() )
ITA_EXCEPT_INVALID_PARAMETER( "Detected an empty path in path list during sensor validation, aborting" );
if (oPath.empty())
ITA_EXCEPT_INVALID_PARAMETER("Detected an empty path in path list during sensor validation, aborting");
auto pLastAnchor = oPath[ oPath.size() - 1 ];
auto pLastAnchor = oPath[oPath.size() - 1];
if( pLastAnchor->iAnchorType != ITAGeo::CPropagationAnchor::ACOUSTIC_SENSOR )
if (pLastAnchor->iAnchorType != ITAGeo::CPropagationAnchor::ACOUSTIC_SENSOR)
return false;
if( pPeviousSensor && pLastAnchor != pPeviousSensor )
if (pPeviousSensor && pLastAnchor != pPeviousSensor)
return false;
pPeviousSensor = std::dynamic_pointer_cast< ITAGeo::CSensor >( pLastAnchor );
pPeviousSensor = std::dynamic_pointer_cast<ITAGeo::CSensor>(pLastAnchor);
}
return true;
}
bool CFilterEngine::ValidateSameEmitter( const ITAGeo::CPropagationPathList& oPathList )
bool CFilterEngine::ValidateSameEmitter(const ITAGeo::CPropagationPathList& oPathList)
{
shared_ptr<ITAGeo::CEmitter> pPeviousEmitter = nullptr;
// Last anchor of the all paths must share the same sensor (pointer)
for( auto& oPath : oPathList )
for (auto& oPath : oPathList)
{
if( oPath.empty() )
ITA_EXCEPT_INVALID_PARAMETER( "Detected an empty path in path list during emitter validation, aborting" );
if (oPath.empty())
ITA_EXCEPT_INVALID_PARAMETER("Detected an empty path in path list during emitter validation, aborting");
auto pFirstAnchor = oPath[ 0 ];
auto pFirstAnchor = oPath[0];
if( pFirstAnchor->iAnchorType != ITAGeo::CPropagationAnchor::ACOUSTIC_EMITTER )
if (pFirstAnchor->iAnchorType != ITAGeo::CPropagationAnchor::ACOUSTIC_EMITTER)
return false;
if( pPeviousEmitter && pFirstAnchor != pPeviousEmitter )
if (pPeviousEmitter && pFirstAnchor != pPeviousEmitter)
return false;
pPeviousEmitter = std::dynamic_pointer_cast< ITAGeo::CEmitter >( pFirstAnchor );
pPeviousEmitter = std::dynamic_pointer_cast<ITAGeo::CEmitter>(pFirstAnchor);
}
return true;
......@@ -152,7 +153,7 @@ ITA_EXCEPT_NOT_IMPLEMENTED;
}
*/
void ITAPropagationModels::CFilterEngine::SetMaterialManager( std::shared_ptr< const ITAGeo::Material::IManager > pMaterialManager )
void ITAPropagationModels::CFilterEngine::SetMaterialManager(std::shared_ptr< const ITAGeo::Material::IManager > pMaterialManager)
{
m_pMaterialManager = pMaterialManager;
}
......@@ -162,7 +163,7 @@ std::shared_ptr< const ITAGeo::Material::IManager > ITAPropagationModels::CFilte
return m_pMaterialManager;
}
void ITAPropagationModels::CFilterEngine::SetDirectivityManager( std::shared_ptr< const ITAGeo::Directivity::IManager > pDirectivityManager )
void ITAPropagationModels::CFilterEngine::SetDirectivityManager(std::shared_ptr< const ITAGeo::Directivity::IManager > pDirectivityManager)
{
m_pDirectivityManager = pDirectivityManager;
}
......@@ -172,71 +173,72 @@ std::shared_ptr< const ITAGeo::Directivity::IManager > ITAPropagationModels::CFi
return m_pDirectivityManager;
}
void CFilterEngine::Generate( const ITAGeo::CPropagationPathList & oPathList, ITABase::CHDFTSpectra & oHDFTSpectra, bool* pbDFTDegreeTooSmall )
void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITABase::CHDFTSpectra & oHDFTSpectra, bool* pbDFTDegreeTooSmall)
{
if( pbDFTDegreeTooSmall )
if (pbDFTDegreeTooSmall)
*pbDFTDegreeTooSmall = false;
//Check for correct structur of oPathList
if( !ValidateSameSensor( oPathList ) )
ITA_EXCEPT1( INVALID_PARAMETER, "The propagation path list does not share the same sensor." );
if (!ValidateSameSensor(oPathList))
ITA_EXCEPT1(INVALID_PARAMETER, "The propagation path list does not share the same sensor.");
if( !ValidateSameEmitter( oPathList ) )
ITA_EXCEPT1( INVALID_PARAMETER, "The propagation path list has more than multiple sensor anchors or last anchor of paths is not a sensor." );
if (!ValidateSameEmitter(oPathList))
ITA_EXCEPT1(INVALID_PARAMETER, "The propagation path list has more than multiple sensor anchors or last anchor of paths is not a sensor.");
//Initialize the path spectra
m_pTempPropPathSpectra = make_unique<ITABase::CHDFTSpectra>( ( float ) oHDFTSpectra.GetSampleRate(), oHDFTSpectra.GetNumChannels(), oHDFTSpectra.GetDFTSize(), true );
m_pAccumulatedSpectra = make_unique<ITABase::CHDFTSpectra>( ( float ) oHDFTSpectra.GetSampleRate(), oHDFTSpectra.GetNumChannels(), oHDFTSpectra.GetDFTSize(), true );
m_pTempPropPathSpectra = make_unique<ITABase::CHDFTSpectra>((float)oHDFTSpectra.GetSampleRate(), oHDFTSpectra.GetNumChannels(), oHDFTSpectra.GetDFTSize(), true);
m_pAccumulatedSpectra = make_unique<ITABase::CHDFTSpectra>((float)oHDFTSpectra.GetSampleRate(), oHDFTSpectra.GetNumChannels(), oHDFTSpectra.GetDFTSize(), true);
double dScalarMagnitude; //For the combination of all scalar filter components
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctFactors; // For the combination of all third octave filter components
for( auto& oPath : oPathList )
for (auto& oPath : oPathList)
{
//Path length
double dPathLength = oPath.GetLength();
if( dPathLength == 0.0f )
ITA_EXCEPT_INVALID_PARAMETER( "Encountered a path with no geometrical distance, could not generate filter" );
if (dPathLength == 0.0f)
ITA_EXCEPT_INVALID_PARAMETER("Encountered a path with no geometrical distance, could not generate filter");
//Include sound energy decrease due to the length of the path @todo this might get more complex with diffraction
dScalarMagnitude /= dPathLength;
dScalarMagnitude = 1/dPathLength;
// Air attenuation according to ISO9613-1 (overrides third octaves, so no resetting required, do not change order!)
ITABase::ISO9613::AtmosphericAbsorption( oThirdOctFactors, ( float ) dPathLength, oEnvironmentState.dTemperature, oEnvironmentState.dHumidity );
ITABase::ISO9613::AtmosphericAbsorption(oThirdOctFactors, (float)dPathLength, oEnvironmentState.dTemperature, oEnvironmentState.dHumidity);
// Multiply the transmissionfactor of the air attenuation alpha: sqrt(1-|alpha|)
for( int iFrequencyIndex = 0; iFrequencyIndex < oThirdOctFactors.GetNumBands(); iFrequencyIndex++ )
oThirdOctFactors[ iFrequencyIndex ] *= sqrt( 1 - abs( oThirdOctFactors[ iFrequencyIndex ] ) );
for (int iFrequencyIndex = 0; iFrequencyIndex < oThirdOctFactors.GetNumBands(); iFrequencyIndex++)
oThirdOctFactors[iFrequencyIndex] *= sqrt(1 - abs(oThirdOctFactors[iFrequencyIndex]));
// Process path and assemble filter components domain by domain
for( size_t i = 0; i < oPath.size(); i++ )
for (size_t i = 0; i < oPath.size(); i++)
{
auto pAnchor( oPath[ i ] );
switch( pAnchor->iAnchorType )
auto pAnchor(oPath[i]);
switch (pAnchor->iAnchorType)
{
case( ITAGeo::CPropagationAnchor::ACOUSTIC_EMITTER ) :
case(ITAGeo::CPropagationAnchor::ACOUSTIC_EMITTER):
{
auto pEmitter = std::dynamic_pointer_cast< const ITAGeo::CEmitter >( pAnchor );
auto pEmitter = std::dynamic_pointer_cast<const ITAGeo::CEmitter>(pAnchor);
if( i == oPath.size() )
ITA_EXCEPT_INVALID_PARAMETER( "Propagation path invalid: detected an emitting anchor point at the end of the path." );
auto pNextAnchor( oPath[ i + 1 ] );
if (i == oPath.size())
ITA_EXCEPT_INVALID_PARAMETER("Propagation path invalid: detected an emitting anchor point at the end of the path.");
auto pNextAnchor(oPath[i + 1]);
if( pEmitter->pDirectivity )
if (pEmitter->pDirectivity)
{
ITAGeo::Coordinates::CSpherical_RightHanded_Degree oDirection;
switch( pEmitter->pDirectivity->GetFormat() )
const VistaVector3D v3Dir( pNextAnchor->v3InteractionPoint - pEmitter->v3InteractionPoint);
ITAGeo::Coordinates::CSpherical oDirection(v3Dir, m_iCoordinateSystemConvention);
switch (pEmitter->pDirectivity->GetFormat())
{
case( ITAGeo::Directivity::IDirectivity::OPENDAFF ) :
case(ITAGeo::Directivity::IDirectivity::OPENDAFF):
{
auto pDAFFDirectivity = std::dynamic_pointer_cast< const ITAGeo::Directivity::CDAFF_Format >( pEmitter->pDirectivity );
if( pDAFFDirectivity->GetDomain() == ITADomain::ITA_FREQUENCY_DOMAIN )
auto pDAFFDirectivity = std::dynamic_pointer_cast<const ITAGeo::Directivity::CDAFF_Format>(pEmitter->pDirectivity);
if (pDAFFDirectivity->GetDomain() == ITADomain::ITA_FREQUENCY_DOMAIN)
{
auto pDAFFDirectivityMS = std::dynamic_pointer_cast< const ITAGeo::Directivity::CDAFF_MagnitudeSpectrum >( pAnchor );
pDAFFDirectivityMS->MultiplyNearestNeighbourMagnitudeSpectrum( oDirection, oThirdOctFactors );
auto pDAFFDirectivityMS = std::dynamic_pointer_cast<const ITAGeo::Directivity::CDAFF_MagnitudeSpectrum>(pEmitter->pDirectivity);
pDAFFDirectivityMS->MultiplyNearestNeighbourMagnitudeSpectrum(oDirection, oThirdOctFactors);
}
}
}
......@@ -275,38 +277,38 @@ void CFilterEngine::Generate( const ITAGeo::CPropagationPathList & oPathList, IT
//DFT coordinates for the respective frequencies
vector<float> vfDFTCoordinatesIn;
int iHDFTSize = ( *m_pTempPropPathSpectra )[ 0 ]->GetSize();
int iHDFTSize = (*m_pTempPropPathSpectra)[0]->GetSize();
//Get the corresponding DFT coordinates from the third octave coordinates
for( auto& fFreq : ITABase::CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies() )
for (auto& fFreq : ITABase::CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies())
{
vfDFTCoordinatesIn.push_back( fFreq*( iHDFTSize - 1 ) / ( ( float ) oHDFTSpectra.GetSampleRate() / 2 ) );
vfDFTCoordinatesIn.push_back(fFreq*(iHDFTSize - 1) / ((float)oHDFTSpectra.GetSampleRate() / 2));
}
//Combine all magnitudes without the eventually existent directivity of the emitter/sensor
for( int iChannel = 0; iChannel < m_pTempPropPathSpectra->GetNumChannels(); iChannel++ )
for (int iChannel = 0; iChannel < m_pTempPropPathSpectra->GetNumChannels(); iChannel++)
{
if( iChannel == 0 )
if (iChannel == 0)
{
//Interpolate third octaves and multiply calculated values with third octave and scalar factors and add delay due to sound propagation
for( int i = 0; i < iHDFTSize; i++ )
for (int i = 0; i < iHDFTSize; i++)
{
// @todo JST/AER: I think this is not correctly using the spline interp, especially tval should be a base value of given frequency, not integer/index
const float fDFTSpectrum = spline_b_val( ( int ) vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[ 0 ], &oThirdOctFactors[ 0 ], ( float ) i );
const float fDFTSpectrum = spline_b_val((int)vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &oThirdOctFactors[0], (float)i);
// @todo: m_pTempPropPathSpectra is untouched -> unity until here ...
const complex<float> fSpectrumData( ( *m_pTempPropPathSpectra )[ 0 ]->GetData()[ 2 * i ], ( *m_pTempPropPathSpectra )[ 0 ]->GetData()[ 2 * i + 1 ] );
const complex<float> fSpectrumData((*m_pTempPropPathSpectra)[0]->GetData()[2 * i], (*m_pTempPropPathSpectra)[0]->GetData()[2 * i + 1]);
complex<float> cfMultipliedSpectra = fSpectrumData * fDFTSpectrum * ( float ) dScalarMagnitude;
complex<float> cfMultipliedSpectra = fSpectrumData * fDFTSpectrum * (float)dScalarMagnitude;
( *m_pTempPropPathSpectra )[ 0 ]->SetCoeff( i, cfMultipliedSpectra );
(*m_pTempPropPathSpectra)[0]->SetCoeff(i, cfMultipliedSpectra);
//Get the delay [samples] and calculate the phase shift according to the shift theorem
// DFT(x(k-Delta)) = exp(-j * omega_k * delta) * X(k) , with omega_k = 2 * pi * k / N
float fDelaySamples = float( dPathLength / oEnvironmentState.fSpeedOfSound * oHDFTSpectra.GetSampleRate() );
float fDelayPhase = -fDelaySamples * ITAConstants::TWO_PI_F * i / ( float ) oHDFTSpectra.GetDFTSize();
float fDelaySamples = float(dPathLength / oEnvironmentState.fSpeedOfSound * oHDFTSpectra.GetSampleRate());
float fDelayPhase = -fDelaySamples * ITAConstants::TWO_PI_F * i / (float)oHDFTSpectra.GetDFTSize();
( *m_pTempPropPathSpectra )[ 0 ]->SetPhasePreserveMagnitude( i, fDelayPhase );
(*m_pTempPropPathSpectra)[0]->SetPhasePreserveMagnitude(i, fDelayPhase);
}
......@@ -314,15 +316,15 @@ void CFilterEngine::Generate( const ITAGeo::CPropagationPathList & oPathList, IT
else
{
//Copy first channel to other channels
( *m_pTempPropPathSpectra )[ iChannel ]->Copy( ( *m_pTempPropPathSpectra )[ 0 ] );
(*m_pTempPropPathSpectra)[iChannel]->Copy((*m_pTempPropPathSpectra)[0]);
}
}
//Add current spectra to accumulated spectra
m_pAccumulatedSpectra->add( m_pTempPropPathSpectra.get() );
m_pAccumulatedSpectra->add(m_pTempPropPathSpectra.get());
}
oHDFTSpectra.CopyFrom( m_pAccumulatedSpectra.get() );
oHDFTSpectra.CopyFrom(m_pAccumulatedSpectra.get());
}
tests/ReceiverDirectivityTest/ReceiverDirectivityTest.cpp
View file @ cd2695d0
......@@ -53,7 +53,7 @@ int main( int, char** )
oPathFromRight.push_back(pSenderRight);
oPathFromRight.push_back(pReceiver);
CFilterEngine oFilterEngine;
CFilterEngine oFilterEngine( OPENGL );
int iFilterLengthSamples = (int) ceil(Utils::EstimateFilterLengthSamples(oPathFromLeft, fSampleRate));
ITABase::CHDFTSpectra oTransmissionFilterLeft(fSampleRate, pReceiver->GetNumChannels(), iFilterLengthSamples );
......
tests/SourceDirectivityTest/SourceDirectivityTest.cpp
View file @ cd2695d0
......@@ -18,82 +18,48 @@
#include <ITAPropagationModels/FilterEngine.h>
#include <ITAPropagationModels/Utils.h>
#include <ITAGeo/Base.h>
#include <ITAGeo/Material/Manager.h>
#include <ITAGeo/Directivity/Base.h>
#include <ITAGeo/Directivity/DAFF_ImpulseResponse.h>
#include <ITAGeo/Directivity/DAFF_MagnitudeSpectrum.h>
#include <ITAFFTUtils.h>
#include <ITAISO9613.h>
using namespace std;
using namespace ITABase;
using namespace ITAConstants;
using namespace ITAGeo;
using namespace ITAFFTUtils;
using namespace ITAPropagationModels;
int main( int, char** )
{
/*
|___x___|
/\
s \ r
\ /
___x___
|______|
*/
auto pSender = make_shared< CEmitter >( VistaVector3D( -2.0f, 0.0f, 0.0f ) );
auto pReceiver = make_shared< CSensor >( VistaVector3D( 2.0f, 0.0f, 0.0f ) );
auto pMaterialDirectory = std::make_shared< Material::CDirectory >( "./" );
auto pReflection = make_shared< CSpecularReflection >( VistaVector3D( 0.0f, 2.0f, 0.0f ) );
pReflection->pAcousticMaterial = pMaterialDirectory->GetMaterial( "stonewall" );
const float fSampleRate = 44.1e3f;
auto pW = std::make_shared< CITADiffractionOuterWedgeAperture >();
pW->v3AperturePoint.SetValues( 0.0f, 1.0f, 0.0f );
pW->v3MainWedgeFaceNormal.SetValues( -1.0f, 1.0f, .0f );
pW->v3MainWedgeFaceNormal.Normalize();
pW->v3OppositeWedgeFaceNormal.SetValues( 1.0f, 1.0f, .0f );
pW->v3OppositeWedgeFaceNormal.Normalize();
pW->v3VertextStart.SetValues( .0f, 1.0f, -1.0f );
pW->v3VertextEnd.SetValues( .0f, 1.0f, 1.0f );
CPropagationPath oPathDirect;
oPathDirect.push_back( pSender );
oPathDirect.push_back( pReceiver );
CPropagationPath oPathReflection;
oPathReflection.push_back( pSender );
oPathReflection.push_back( pReflection );
oPathReflection.push_back( pReceiver );
CPropagationPath oPathDiffraction;
oPathDiffraction.push_back( pSender );
oPathDiffraction.push_back( pW );
oPathDiffraction.push_back( pReceiver );
int main(int, char**)
{
shared_ptr< Directivity::CDAFF_MagnitudeSpectrum > pTrumpetDirectivity = make_shared< Directivity::CDAFF_MagnitudeSpectrum >();
pTrumpetDirectivity->LoadFromFile("Trumpet1.v17.ms.daff");
CPropagationPathList oPathList;
oPathList.push_back(oPathDirect);
oPathList.push_back(oPathReflection);
oPathList.push_back(oPathDiffraction);
auto pSender = make_shared< CEmitter >(VistaVector3D(2.0f, 1.3f, 2.0f));
pSender->pDirectivity = pTrumpetDirectivity;
CFilterEngine oFilterEngine;
oFilterEngine.SetMaterialManager( pMaterialDirectory );
auto pReceiver = make_shared< CSensor >(VistaVector3D(0.0f, 0.0f, 0.0f));
CPropagationPath oPath;
oPath.push_back(pSender);
oPath.push_back(pReceiver);
// Set filter length according to the maximum path length
const float fSpeedOfSound = DEFAULT_SPEED_OF_SOUND_F; //Approximation of speed of sound at ~20C
const float fSampleRate = 44.1e3f;
int iFilterLengthSamples = ( int ) ( oPathList.GetMaxLength() / fSpeedOfSound * fSampleRate );
iFilterLengthSamples = int( iFilterLengthSamples + 4096 );
CFilterEngine oFilterEngine(OPENGL);
ITABase::CHDFTSpectra oTransmissionFilter( fSampleRate, pReceiver->GetNumChannels(), iFilterLengthSamples );
int iFilterLengthSamples = (int)ceil(Utils::EstimateFilterLengthSamples(oPath, fSampleRate));
CHDFTSpectra oTransmissionFilter(fSampleRate, pReceiver->GetNumChannels(), iFilterLengthSamples);
bool bDFTDegreeTooSmallFlag;
oFilterEngine.Generate( oPathList, oTransmissionFilter, &bDFTDegreeTooSmallFlag );
oFilterEngine.Generate(oPath, oTransmissionFilter, &bDFTDegreeTooSmallFlag);
if( bDFTDegreeTooSmallFlag )
if (bDFTDegreeTooSmallFlag)
cerr << "DFT lengh too small, could not include all path completely into target filter" << endl;
ITAFFTUtils::Export( &oTransmissionFilter, "FilterEngineTest.wav" ); // Exports in time domain as impulse response (IR)
Export(&oTransmissionFilter, "SourceDirectivityTest.wav"); // Exports in time domain as impulse response (IR)
}
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