Skip to content

GitLab

I
ITAPropagationModels
Commit 31464c3e authored by Armin Erraji's avatar Armin Erraji
Browse files
Options

Merge branch 'develop' of https://git.rwth-aachen.de/ita/ITAPropagationModels into develop

parents b21c8c92 e3a3e49f
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1891 additions and 386 deletions
CMakeLists.txt
View file @ 31464c3e
......@@ -20,15 +20,13 @@ set( ITAPropagationModelsHeader
"include/ITAPropagationModels/Definitions.h"
"include/ITAPropagationModels/DiffractionFilter.h"
"include/ITAPropagationModels/FilterEngine.h"
"include/ITAPropagationModels/FilterGenerator.h"
"include/ITAPropagationModels/Kirchhoff.h"
"include/ITAPropagationModels/Maekawa.h"
"include/ITAPropagationModels/Svensson.h"
"include/ITAPropagationModels/UTD.h"
)
set( ITAPropagationModelsSources
"src/ITAPropagationModels/FilterEngine.cpp"
"src/ITAPropagationModels/FilterGenerator.cpp"
"src/ITAPropagationModels/Kirchhoff.cpp"
"src/ITAPropagationModels/Maekawa.cpp"
"src/ITAPropagationModels/Svensson.cpp"
"src/ITAPropagationModels/UTD.cpp"
)
......
benchmarks/DiffractionModelsBenchmark.cpp
View file @ 31464c3e
......@@ -21,7 +21,7 @@
#include <ITAPropagationModels/Base.h>
#include <ITAPropagationModels/Kirchhoff.h>
#include <ITAPropagationModels/Maekawa.h>
#include <ITAPropagationModels/UTD.h>
#include <ITAPropagationModels/Svensson.h>
......@@ -47,7 +47,7 @@
using namespace std;
using namespace ITAPropagationModels;
float g_fSpeedOfSound = 344.0f;
float g_fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F;
float g_fSampleRate = 44.1e3;
int g_iFilterLength = 64;
......@@ -82,44 +82,47 @@ void benchmark_single_wedge()
oReceiver->qOrient = VistaQuaternion( 0, 1, 0, 0 );
oReceiver->vPos.SetValues( -.4f, .6f, .1f );
auto oWedge = std::make_shared< ITAGeo::CITADiffractionWedgeAperture >();
oWedge->v3MainWedgeFaceNormal = VistaVector3D( 1.0f, 1.0f, .0f ).GetNormalized();
oWedge->v3OppositeWedgeFaceNormal = VistaVector3D( -1.0f, 1.0f, .0f ).GetNormalized();
oWedge->v3VertextStart.SetValues( .0f, 1.0f, -1.0f );
oWedge->v3VertextEnd.SetValues( .0f, 1.0f, 1.0f );
auto pWedge = std::make_shared< ITAGeo::CITADiffractionOuterWedgeAperture >();
pWedge->v3MainWedgeFaceNormal = VistaVector3D( 1.0f, 1.0f, .0f ).GetNormalized();
pWedge->v3OppositeWedgeFaceNormal = VistaVector3D( -1.0f, 1.0f, .0f ).GetNormalized();
pWedge->v3VertextStart.SetValues( .0f, 1.0f, -1.0f );
pWedge->v3VertextEnd.SetValues( .0f, 1.0f, 1.0f );
ITAGeoUtils::CalculateDiffractionAperturePoint( oWedge->GetEdgeRay(), oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, oWedge->v3AperturePoint );
ITAGeoUtils::CalculateDiffractionAperturePoint( pWedge->GetEdgeRay(), oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, pWedge->v3AperturePoint );
if( !ITAGeoUtils::IsDiffractionAperturePointInRange( oWedge->v3VertextStart, oWedge->v3VertextEnd, oWedge->v3AperturePoint ) )
if( !ITAGeoUtils::IsDiffractionAperturePointInRange( pWedge->v3VertextStart, pWedge->v3VertextEnd, pWedge->v3AperturePoint ) )
ITA_EXCEPT_INVALID_PARAMETER( "Aperture point over wedge is outside edge vertices" );
if( !oWedge->IsOutsideWedge( oSource->v3InteractionPoint ) )
if( !pWedge->IsOutsideWedge( oSource->v3InteractionPoint ) )
ITA_EXCEPT_INVALID_PARAMETER( "Source point is inside wedge (solid part)" );
if( !oWedge->IsOutsideWedge( oReceiver->v3InteractionPoint ) )
if( !pWedge->IsOutsideWedge( oReceiver->v3InteractionPoint ) )
ITA_EXCEPT_INVALID_PARAMETER( "Receiver point is inside wedge (solid part)" );
ITAGeo::CPropagationPath oPropPath;
oPropPath.push_back( oSource );
oPropPath.push_back( oWedge );
oPropPath.push_back( pWedge );
oPropPath.push_back( oReceiver );
cout << "Path: " << oPropPath << endl;
// Kirchhoff
if( Kirchhoff::IsApplicable( oPropPath ) )
if( Maekawa::IsApplicable( oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, pWedge ) )
ITA_EXCEPT_INVALID_PARAMETER( "Kirchhoff model can not be applied for this propagation path." );
ITAStopWatch sw;
cout << "Starting Kirchoff model benchmark:" << endl;
double dStartTime = ITAClock::getDefaultClock()->getTime();
ITABase::CHDFTSpectrum oTF( g_fSampleRate, g_iFilterLength, true );
for( int i = 0; i < 1e8; i++ )
{
sw.start();
Kirchhoff::CalculateDiffractionFilter( oPropPath, oTF, g_fSpeedOfSound );
double dDir, dDet;
Maekawa::GetDirectLengthAndDetourLength( oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, pWedge, dDir, dDet );
Maekawa::CalculateDiffractionFilter( dDir, dDet, oTF, g_fSpeedOfSound );
sw.stop();
}
double dStopTime = ITAClock::getDefaultClock()->getTime();
......@@ -136,7 +139,7 @@ void benchmark_single_wedge()
for( int i = 0; i < 1e6; i++ )
{
sw.start();
UTD::CalculateDiffractionFilter(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, oWedge, oTF, iMethod, g_fSpeedOfSound);
UTD::CalculateDiffractionFilter(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, pWedge, oTF, iMethod, g_fSpeedOfSound);
sw.stop();
}
dStopTime = ITAClock::getDefaultClock()->getTime();
......@@ -149,14 +152,14 @@ void benchmark_single_wedge()
sw.reset();
cout << "Starting Biot-Tolstoy-Medwin-Svensson model benchmark:" << endl;
dStartTime = ITAClock::getDefaultClock()->getTime();
const float fMinDelay = oWedge->GetMinimumWavefrontDelayTime(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, g_fSpeedOfSound);
const float fMaxDelay = oWedge->GetMaximumWavefrontDelayTime(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, g_fSpeedOfSound);
const float fMinDelay = pWedge->GetMinimumWavefrontDelayTime(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, g_fSpeedOfSound);
const float fMaxDelay = pWedge->GetMaximumWavefrontDelayTime(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, g_fSpeedOfSound);
const int iIRLength = ( int ) ceil( float( fMaxDelay - fMinDelay ) * g_fSampleRate );
ITABase::CFiniteImpulseResponse oIR( iIRLength, g_fSampleRate, true );
for( int i = 0; i < 1e4; i++ )
{
sw.start();
Svensson::CalculateDiffractionIR(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, oWedge, oIR, g_fSpeedOfSound);
Svensson::CalculateDiffractionIR(oSource->v3InteractionPoint, oReceiver->v3InteractionPoint, pWedge, oIR, g_fSpeedOfSound);
sw.stop();
}
dStopTime = ITAClock::getDefaultClock()->getTime();
......
include/ITAPropagationModels/FilterEngine.h
View file @ 31464c3e
......@@ -25,14 +25,18 @@
// ITA includes
#include <ITAGeo/Base.h>
#include <ITAGeo/Material//Manager.h>
#include <ITAHDFTSpectra.h>
#include <ITAHDFTSpectrum.h>
// STL includes
#include <vector>
namespace ITAPropagationModels
{
using namespace std;
//! Transfer function filter generator for propagation paths
/**
* Generates transfer functions in the frequency-domain that
......@@ -65,12 +69,37 @@ namespace ITAPropagationModels
//! Applies all acoustic models like reflection diffraction models
void ApplyAcousticModels( ITAGeo::CPropagationPathList& oPathList );
void ApplyDiffractionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = 1 );
void ApplyReflectionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = 1 );
//!Applies diffraction models
/**
* @param[out] oPathList Propagation path list
* @param[in] iModel Acoustic material model type. If value is set to -1, the default material model is used.
*/
void ApplyDiffractionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = -1);
//!Applies reflection models
/**
* @param[out] oPathList Propagation path list
* @param[in] iModel Acoustic material model type. If value is set to -1, the default material model is used.
*/
void ApplyReflectionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = -1 );
void ApplyTransmissionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = 1 ); // @todo
//!Applies emitter models
/**
* @param[out] oPathList Propagation path list
* @param[in] iModel Acoustic material model type. If value is set to -1, the default material model is used.
*/
void ApplyEmitterModel(ITAGeo::CPropagationPathList& oPathList, int iModel = -1);
//! Generate multi-channel propagation path (for multi-channel receiver directivitiy)
//!Applies sensor models
/**
* @param[out] oPathList Propagation path list
* @param[in] iModel Acoustic material model type. If value is set to -1, the default material model is used.
*/
void ApplySensorModel(ITAGeo::CPropagationPathList& oPathList, int iModel = -1);
void ApplyTransmissionModel( ITAGeo::CPropagationPathList& oPathList, int iModel = -1 ); // @todo
//! Generate multi-channel propagation path (for multi-channel receiver directivity)
/**
* @todo AER
*/
......@@ -79,15 +108,44 @@ namespace ITAPropagationModels
//! Generate single-channel propagation path (for single-channel receiver directivity)
inline void Generate( const ITAGeo::CPropagationPathList& oPathList, ITABase::CHDFTSpectrum& oFilter )
{
const std::vector< ITABase::CHDFTSpectrum* > vpSpectra = { &oFilter };
const vector< ITABase::CHDFTSpectrum* > vpSpectra = { &oFilter };
ITABase::CHDFTSpectra oTF( vpSpectra );
Generate( oPathList, oTF );
};
//! Sets a connection to the material manager
void SetMaterialManager( const ITAGeo::CMaterialManager* pMaterialManager );
// Returns pointer to material manager or null
const ITAGeo::CMaterialManager* GetMaterialManager() const;
private:
ITABase::CHDFTSpectra m_oAccumulatedSpectra; //!< Gathered propagation paths from list
ITABase::CHDFTSpectra m_oTempPropPathSpectra; //!< Single prop-path spectra
unique_ptr<ITABase::CHDFTSpectra> m_pAccumulatedSpectra; //!< Gathered propagation paths from list
unique_ptr<ITABase::CHDFTSpectra> m_pTempPropPathSpectra; //!< Single prop-path spectra
static struct m_DefaultValues //!< Default values
{
static const int iReflectionModel = ITAGeo::IAcousticMaterial::SCALAR;
static const int iDiffractionModel = ITAGeo::IAcousticMaterial::THIRD_OCTAVE;
static const int iEmitterModel = ITAGeo::IAcousticMaterial::SCALAR;
static const int iSensorModel = ITAGeo::IAcousticMaterial::SCALAR;
};
const ITAGeo::CMaterialManager* m_pMaterialManager;
const double m_dHumidity = 80.0;
const double m_dTemperature = 20.0;
const float m_fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F;
static struct m_DefaultDiffractionModel //!< Default values for diffractions
{
int iModel = 1;
};
};
}
......
include/ITAPropagationModels/FilterGenerator.h deleted 100644 → 0
View file @ b21c8c92
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2018
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef INCLUDE_WATCHER_ITA_PROPAGATION_MODELS_FILTER_GENERATOR
#define INCLUDE_WATCHER_ITA_PROPAGATION_MODELS_FILTER_GENERATOR
#include "Base.h"
#include "Definitions.h"
#include "DiffractionFilter.h"
// ITA includes
#include <ITAGeo/Base.h>
#include <ITAHDFTSpectra.h>
#include <ITAHDFTSpectrum.h>
// STL includes
#include <vector>
namespace ITAPropagationModels
{
//! Transfer function filter generator for propagation paths
/**
* Generates transfer functions in the frequency-domain that
* can be used as filters, e.g. for auralization.
*/
class ITA_PROPAGATION_MODELS_API CTFGenerator
{
public:
//! Construct a generator with predetermined output channels, e.g. for binaural or SH filters
/**
* @param[in] iNumOutputChannels Number of output channels
* @param[in] iDFTSize Length of the DFT used to construct and combine filter components
* @param[in] fSamplingRate Filer sampling rate
*/
CTFGenerator( const int iNumOutputChannels, const int iDFTSize = 128, const float fSamplingRate = 44100.0f );
//! Cunstruct a generator with matching properties of a target spectrum
inline CTFGenerator( const ITABase::CHDFTSpectrum& pTransferFunction )
: CTFGenerator( 1, pTransferFunction.GetDFTSize(), pTransferFunction.GetSampleRate() )
{
};
//! Cunstruct a generator with matching properties of multi-channel target spectra
CTFGenerator( const ITABase::CHDFTSpectra& pTransferFunctions )
: CTFGenerator( pTransferFunctions.GetNumChannels(), pTransferFunctions.GetDFTSize(), pTransferFunctions.GetSampleRate() )
{
};
//! Generate a single-channel transfer function
inline void GenerateTF( ITABase::CHDFTSpectrum& pTransferFunction )
{
const std::vector< ITABase::CHDFTSpectrum* > vpSpectra = { &( pTransferFunction ) };
ITABase::CHDFTSpectra oTF( vpSpectra );
GenerateTF( oTF );
};
//! Generate a multi-channel transfer function
inline void GenerateTF( ITABase::CHDFTSpectra& pTransferFunctions );
private:
};
}
#endif // INCLUDE_WATCHER_ITA_PROPAGATION_MODELS_FILTER_GENERATOR
include/ITAPropagationModels/Svensson.h
View file @ 31464c3e
......@@ -49,7 +49,7 @@ namespace ITAPropagationModels
*
* @sa GetMinimumWavefrontDelayTime, GetMaximumWavefrontDelayTime
*/
ITA_PROPAGATION_MODELS_API bool CalculateDiffractionIR( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeAperture > pApex, ITABase::CFiniteImpulseResponse& oEffectiveDiffractionIR, const float fSpeedOfSound = 344.0f );
ITA_PROPAGATION_MODELS_API bool CalculateDiffractionIR( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex, ITABase::CFiniteImpulseResponse& oEffectiveDiffractionIR, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F );
}
}
......
include/ITAPropagationModels/UTD.h
View file @ 31464c3e
......@@ -40,6 +40,9 @@ namespace ITAPropagationModels
{
namespace UTD
{
using namespace std;
//! UTD diffraction coefficient calculation methods (only approximation implemented yet!)
enum UTD_DIFFRACTION_ALGORITHM
{
......@@ -60,10 +63,24 @@ namespace ITAPropagationModels
* @param[in] v3TargetPos target position (outside wedge)
* @param[in] pApex Aperture point including wedge parameters (angles, etc)
* @param[in] fWaveNumber Wave number ( k )
* @param[out] cfCoeff Complex-valued diffraction coefficient ( D( n, k,\rho, r, \phi_i, \alpha_i, \alpha_d ) )
* @param[out] cfFactor Complex-valued diffraction coefficient ( D( n, k,\rho, r, \phi_i, \alpha_i, \alpha_d ) )
* @param[in] iMethod Approximation algorithm
*/
ITA_PROPAGATION_MODELS_API bool CalculateDiffractionCoefficient( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeAperture > pApex, const float fWaveNumber, std::complex< float >& cfCoeff, const int iMethod = UTD_APPROX_KAWAI_KOUYOUMJIAN );
bool ITA_PROPAGATION_MODELS_API CalculateDiffractionCoefficient( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex, const float fFrequency, complex< float >& cfCoeff, const int iMethod = UTD_APPROX_KAWAI_KOUYOUMJIAN, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F );
//! Calculate diffraction factor based on the Maekawa detour method
/**
* This function calculates the Maekawa diffraction factor based on detour length around a wedge obstacle
* and the direct path length with and without any obstacles.
* The diffraction factor contains the influence of the wedge without the influence of the propagation path length.
*
* @param[in] dPropagationLengthFreeField Direct path length [m]
* @param[in] dPropagationLengthOverApex Detour path length [m]
* @param[in] fFrequency Frequency [Hz]
* @param[in] fSpeedOfSound Speed of sound [m/s]
*
*/
bool ITA_PROPAGATION_MODELS_API CalculateDiffractionFactor(const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex, const float fFrequency, complex< float >& cfFactor, const int iMethod = UTD_APPROX_KAWAI_KOUYOUMJIAN, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F );
//! Calculates the UTD diffraction filter in frequency domain for a given constellation
/**
......@@ -78,7 +95,7 @@ namespace ITAPropagationModels
* @param[in] iMethod Approximation algorithm
* @param[in] fSpeedOfSound Sound speed
*/
ITA_PROPAGATION_MODELS_API void CalculateDiffractionFilter( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeAperture > pApex, ITABase::CHDFTSpectrum& oTF, const int iMethod = UTD_APPROX_KAWAI_KOUYOUMJIAN, const float fSpeedOfSound = 344.0f );
void ITA_PROPAGATION_MODELS_API CalculateDiffractionFilter( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex, ITABase::CHDFTSpectrum& oTF, const int iMethod = UTD_APPROX_KAWAI_KOUYOUMJIAN, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F );
}
}
......
src/ITAPropagationModels/FilterEngine.cpp
View file @ 31464c3e
This diff is collapsed.
src/ITAPropagationModels/FilterGenerator.cpp deleted 100644 → 0
View file @ b21c8c92
#include <ITAPropagationModels/FilterGenerator.h>
// ITA includes
#include <ITAException.h>
#include <ITAGeo/Halfedge/MeshModel.h>
// Vista includes
#include <VistaMath/VistaGeometries.h>
// STL includes
#include <assert.h>
#include <stdio.h>
using namespace ITAGeo;
using namespace ITAPropagationModels;
ITAPropagationModels::CTFGenerator::CTFGenerator( const int iNumOutputChannels, const int iDFTSize /*= 128 */, const float fSamplingRate )
{
}
void ITAPropagationModels::CTFGenerator::GenerateTF( ITABase::CHDFTSpectra& pTransferFunctions )
{
}
src/ITAPropagationModels/Kirchhoff.cpp deleted 100644 → 0
View file @ b21c8c92
#include <ITAPropagationModels/Kirchhoff.h>
#include <ITAException.h>
using namespace ITAGeo;
using namespace ITAPropagationModels;
bool Kirchhoff::IsApplicable( const CPropagationPath& oPropPath )
{
if( oPropPath.size() < 3 )
ITA_EXCEPT1( INVALID_PARAMETER, "Propagation path needs at least 3 anchors for Kirchhoff method" );
auto pAnchorDirectLazy( oPropPath[ 0 ] );
if( pAnchorDirectLazy->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "First geo propagation anchor is invalid, please purge first." );
auto pAnchorTail( oPropPath[ oPropPath.size() - 1 ] );
if( pAnchorTail->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Last geo propagation anchor is invalid, please purge first." );
// We don't care for first and last anchor's type (except invalid)
for( size_t i = 1; i < oPropPath.size() - 1; i++ )
{
auto pAnchorPrev( oPropPath[ i - 1 ] );
auto pAnchorCur( oPropPath[ i ] );
auto pAnchorNext( oPropPath[ i + 1 ] );
if( pAnchorCur->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Geo propagation path has invalid anchors, please purge first." );
if( pAnchorCur->iAnchorType == CPropagationAnchor::DIFFRACTION_APEX )
{
auto pWedge = std::dynamic_pointer_cast< CITADiffractionWedgeAperture>( pAnchorCur );
// Check if prev and next anchors are occluded by wedge
if( !pWedge->IsOccluding( pAnchorPrev, pAnchorNext ) )
return false;
}
}
return true;
}
void Kirchhoff::GetDirectLengthAndDetourLength( const CPropagationPath& oPropPath, double& dPropagationLengthDirect, double& dPropagationLengthDetour )
{
if( oPropPath.size() < 3 )
ITA_EXCEPT1( INVALID_PARAMETER, "Propagation path needs at least 3 anchors for Kirchhoff method" );
auto pAnchorDirectLazy( oPropPath[ 0 ] );
if( pAnchorDirectLazy->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "First geo propagation anchor is invalid, please purge." );
auto pAnchorTail( oPropPath[ oPropPath.size() - 1 ] );
if( pAnchorTail->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Last geo propagation anchor is invalid, please purge." );
// We don't care for first and last anchor's type (except invalid)
dPropagationLengthDirect = 0.0f;
dPropagationLengthDetour = 0.0f;
for( size_t i = 1; i < oPropPath.size(); i++ )
{
auto pAnchorPrev( oPropPath[ i - 1 ] );
auto pAnchorCur( oPropPath[ i ] );
if( pAnchorCur->iAnchorType == CPropagationAnchor::INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Geo propagation path has invalid anchors, please purge first." );
// Detour is everything ...
const VistaVector3D vIntermediateLine = pAnchorCur->v3InteractionPoint - pAnchorPrev->v3InteractionPoint;
const double dIntermediateLength = vIntermediateLine.GetLength();
dPropagationLengthDetour += dIntermediateLength;
// ... but direct path skips diffraction anchors
if( pAnchorCur->iAnchorType != CPropagationAnchor::DIFFRACTION_APEX || pAnchorCur == pAnchorTail )
{
const VistaVector3D vIntermediateDirectLine = pAnchorCur->v3InteractionPoint - pAnchorDirectLazy->v3InteractionPoint;
const double dIntermediateDirectLine = vIntermediateDirectLine.GetLength();
dPropagationLengthDirect += dIntermediateDirectLine;
pAnchorDirectLazy = pAnchorCur;
}
}
// Zero-length path would be crude, make assertion here to give user a hint that something is wrong with his input
assert( dPropagationLengthDetour > 0.0f );
assert( dPropagationLengthDirect > 0.0f );
return;
}
double Kirchhoff::KirchhoffDelta( const double dPropagationLengthDirect, const double dPropagationLengthDetour )
{
if( dPropagationLengthDirect <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction direct path length can not be zero or negative" );
if( dPropagationLengthDetour <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction detour path length can not be zero or negative" )
if( dPropagationLengthDetour < dPropagationLengthDirect )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction detour path length can not be smaller than direct path" );
const double dKirchhoffDelta = dPropagationLengthDetour - dPropagationLengthDirect; // Difference in path length, no magic here
return dKirchhoffDelta;
}
void Kirchhoff::CalculateDiffractionFilter( const double dDirect, const double dDetour, ITABase::CHDFTSpectrum& oTransferFunction, const float fSpeedOfSound /*= 344.0f*/ )
{
assert( fSpeedOfSound > 0.0f );
if( fSpeedOfSound <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Speed of sound can't be zero or negative" );
// Per definition, DC value is real and 1.0
oTransferFunction.SetCoeffRI( 0, 1.0f, 0.0f );
const int iDFTSize = oTransferFunction.GetSize();
for( int i = 1; i < iDFTSize - 1; i++ )
{
const double dFrequency = oTransferFunction.GetFrequencyResolution() * i;
const double dFresnelNumber = 2.0f * Kirchhoff::KirchhoffDelta( dDirect, dDetour ) * dFrequency / fSpeedOfSound; // Sometimes referred to as 'N'
oTransferFunction.SetCoeffRI( i, float( dFresnelNumber ) );
}
}
src/ITAPropagationModels/Maekawa.cpp
View file @ 31464c3e
#include <ITADiffractionMaekawa.h>
#include <ITAPropagationModels/Maekawa.h>
#include <ITAException.h>
bool IsMaekawaApplicable( const CITAGeoPropagationPath& oPropPath )
{
if( oPropPath.size() < 3 )
ITA_EXCEPT1( INVALID_PARAMETER, "Propagation path needs at least 3 anchors for Maekawa method" );
#include <cassert>
const CITAGeoPropagationAnchor* pAnchorDirectLazy( oPropPath[ 0 ] );
if( pAnchorDirectLazy->iAnchorType == CITAGeoPropagationAnchor::ITA_ANCHOR_INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "First geo propagation anchor is invalid, please purge first." );
using namespace ITAGeo;
using namespace ITAPropagationModels;
const CITAGeoPropagationAnchor* pAnchorTail( oPropPath[ oPropPath.size() - 1 ] );
if( pAnchorTail->iAnchorType == CITAGeoPropagationAnchor::ITA_ANCHOR_INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Last geo propagation anchor is invalid, please purge first." );
bool Maekawa::IsApplicable( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex )
{
return pApex->IsOccluding( v3SourcePos, v3TargetPos );
}
// We don't care for first and last anchor's type (except invalid)
void Maekawa::GetDirectLengthAndDetourLength( const VistaVector3D& v3SourcePos, const VistaVector3D& v3TargetPos, std::shared_ptr< const ITAGeo::CITADiffractionWedgeApertureBase > pApex, double& dPropagationLengthFreeField, double& dPropagationLengthOverApex )
{
dPropagationLengthOverApex = ( pApex->v3InteractionPoint - v3SourcePos ).GetLength() + ( v3TargetPos - pApex->v3InteractionPoint ).GetLength();
dPropagationLengthFreeField = ( v3TargetPos - v3SourcePos ).GetLength();
double dIncrementalAngle = 0.0f;
assert( dPropagationLengthOverApex > 0.0f );
assert( dPropagationLengthFreeField > 0.0f );
}
for( size_t i = 1; i < oPropPath.size() - 1; i++ )
{
const CITAGeoPropagationAnchor* pAnchorPrev( oPropPath[ i - 1 ] );
const CITAGeoPropagationAnchor* pAnchorCur( oPropPath[ i ] );
const CITAGeoPropagationAnchor* pAnchorNext( oPropPath[ i + 1 ] );
double Maekawa::KirchhoffDelta( const double dPropagationLengthFreeField, const double dPropagationLengthOverApex )
{
if( dPropagationLengthFreeField <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction direct path length can not be zero or negative" );
if( dPropagationLengthOverApex <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction detour path length can not be zero or negative" )
if( pAnchorCur->iAnchorType == CITAGeoPropagationAnchor::ITA_ANCHOR_INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Geo propagation path has invalid anchors, please purge first." );
if( dPropagationLengthOverApex < dPropagationLengthFreeField )
ITA_EXCEPT1( INVALID_PARAMETER, "Diffraction detour path length can not be smaller than direct path" );
if( pAnchorCur->iAnchorType == CITAGeoPropagationAnchor::ITA_EDGE_DIFFRACTION )
{
const CITADiffractionWedge* pWedge( ( const CITADiffractionWedge* ) pAnchorCur );
dIncrementalAngle += pWedge->GetOpeningAngleRad();
const double dKirchhoffDelta = dPropagationLengthOverApex - dPropagationLengthFreeField; // Difference in path length, no magic here
return dKirchhoffDelta;
}
if( dIncrementalAngle > ITAConstants::TWO_PI_D )
return false;
void Maekawa::CalculateDiffractionFilter( const double dDirect, const double dDetour, ITABase::CHDFTSpectrum& oTransferFunction, const float fSpeedOfSound /*= ITAConstants::SPEED_OF_SOUND_F*/ )
{
assert( fSpeedOfSound > 0.0f );
if( fSpeedOfSound <= 0.0f )
ITA_EXCEPT1( INVALID_PARAMETER, "Speed of sound can't be zero or negative" );
// Per definition, DC value (N = 0) is set to or 10^(5/20) / r_dir (5dB loss as well as loss due to detour length)
float fCoefficient = pow(10.0f, 5.0f / 20.0f)/dDirect;
oTransferFunction.SetCoeffRI(0, fCoefficient, 0.0f);
}
const int iDFTSize = oTransferFunction.GetSize();
for( int i = 1; i < iDFTSize; i++ )
{
const float fFrequency = oTransferFunction.GetFrequencyResolution() * i;
fCoefficient = (float) CalculateDiffractionFilterCoefficient( dDirect, dDetour, fFrequency, fSpeedOfSound );
oTransferFunction.SetCoeffRI( i, fCoefficient );
//Set wave number and phase shift
const float k = ITAConstants::TWO_PI_F_L * fFrequency / fSpeedOfSound;
oTransferFunction.SetPhasePreserveMagnitude(i, - k * (float) dDetour);
}
return true;
}
void MaekawaDetourAndLength( const CITAGeoPropagationPath& oPropPath, double& dPropagationLengthDetour, double& dPropagationLengthDirect )
double Maekawa::CalculateDiffractionFilterCoefficient( const double dPropagationLengthFreeField, const double dPropagationLengthOverApex, const float fFrequency, const float fSpeedOfSound /*= ITAConstants::SPEED_OF_SOUND_F */ )
{
if( oPropPath.size() < 3 )
ITA_EXCEPT1( INVALID_PARAMETER, "Propagation path needs at least 3 anchors for Maekawa method" );
const CITAGeoPropagationAnchor* pAnchorDirectLazy( oPropPath[ 0 ] );
if( pAnchorDirectLazy->iAnchorType == CITAGeoPropagationAnchor::ITA_ANCHOR_INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "First geo propagation anchor is invalid, please purge first." );
//Difference between detour length and direct length
const double dDifference = dPropagationLengthOverApex - dPropagationLengthFreeField;
const CITAGeoPropagationAnchor* pAnchorTail( oPropPath[ oPropPath.size() - 1 ] );
if( pAnchorTail->iAnchorType == CITAGeoPropagationAnchor::ITA_ANCHOR_INVALID )
ITA_EXCEPT1( INVALID_PARAMETER, "Last geo propagation anchor is invalid, please purge first." );
//Fresnel number N according to Maekawa //TODO: Equation number
const double dN = 2.0 * dDifference / fSpeedOfSound * fFrequency;
// We don't care for first and last anchor's type (except invalid)
//Calculation of the diffraction filter coefficient according to Kurze, Anderson
const double dTemp = sqrt(ITAConstants::TWO_PI_D * dN);
const double dFilterCoefficient = pow(pow(10.0, 5.0 / 20.0)* dTemp / tanh(dTemp), -1) / dPropagationLengthFreeField;
dPropagationLengthDetour = 0.0f;
dPropagationLengthDirect = 0.0f;