From 4482fe3f28d30efe613957a0c210010b0e214fb1 Mon Sep 17 00:00:00 2001 From: Armin Erraji Date: Mon, 17 Dec 2018 15:25:42 +0100 Subject: [PATCH] Added functions CalculateDiffractionFilter() and CalculateDiffractionFilterCoefficient() for the calculation of the Diffraction transfer function and the diffraction coefficient, as well as the function CalculateDiffractionFactor() for the calculation of the (frequency dependent) diffraction factor without the influence of the path length. --- include/ITAPropagationModels/Maekawa.h | 22 ++++++++++++--- src/ITAPropagationModels/Maekawa.cpp | 39 ++++++++++++++++++++------ tests/Maekawa/MaekawaTest.cpp | 7 ++++- tests/Maekawa/MaekawaTest.m | 6 +++- 4 files changed, 60 insertions(+), 14 deletions(-) diff --git a/include/ITAPropagationModels/Maekawa.h b/include/ITAPropagationModels/Maekawa.h index b114ab1..2ce9dcc 100644 --- a/include/ITAPropagationModels/Maekawa.h +++ b/include/ITAPropagationModels/Maekawa.h @@ -108,8 +108,22 @@ namespace ITAPropagationModels //! Calculate diffraction filter based on the Maekawa detour method /** - * This function calculates the Maekawa diffraction filter based on detour length around a wedge obstacle - * and the direct path length with and withont any obstacles. + * This function calculates the Maekawa diffraction filter coefficient based on detour length around a wedge obstacle + * and the direct path length with and without any obstacles. + * + * @param[in] dPropagationLengthDirect Direct path length [m] + * @param[in] dPropagationLengthDetour Detour path length [m] + * @param[in] fFrequency Frequency [Hz] + * @param[in] fSpeedOfSound Speed of sound [m/s] + * + */ + double ITA_PROPAGATION_MODELS_API CalculateDiffractionFilterCoefficient( const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, 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] dPropagationLengthDirect Direct path length [m] * @param[in] dPropagationLengthDetour Detour path length [m] @@ -117,7 +131,7 @@ namespace ITAPropagationModels * @param[in] fSpeedOfSound Speed of sound [m/s] * */ - float ITA_PROPAGATION_MODELS_API CalculateDiffractionFilterCoefficient( const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, const float fSpeedOfSound = 344.0f ); + double ITA_PROPAGATION_MODELS_API CalculateDiffractionFactor(const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F); //! Calculate diffraction filter (frequency domain) based on the Maekawa detour method @@ -131,7 +145,7 @@ namespace ITAPropagationModels * @param[in] fSpeedOfSound Speed of sound [m/s] * */ - void ITA_PROPAGATION_MODELS_API CalculateDiffractionFilter( const double dPropagationLengthDirect, const double dPropagationLengthDetour, ITABase::CHDFTSpectrum& oTransferFunction, const float fSpeedOfSound = 344.0f ); + void ITA_PROPAGATION_MODELS_API CalculateDiffractionFilter( const double dPropagationLengthDirect, const double dPropagationLengthDetour, ITABase::CHDFTSpectrum& oTransferFunction, const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F); } } diff --git a/src/ITAPropagationModels/Maekawa.cpp b/src/ITAPropagationModels/Maekawa.cpp index 8023b4d..146bccc 100644 --- a/src/ITAPropagationModels/Maekawa.cpp +++ b/src/ITAPropagationModels/Maekawa.cpp @@ -36,26 +36,49 @@ double Maekawa::KirchhoffDelta( const double dPropagationLengthDirect, const dou return dKirchhoffDelta; } -void Maekawa::CalculateDiffractionFilter( const double dDirect, const double dDetour, ITABase::CHDFTSpectrum& oTransferFunction, const float fSpeedOfSound /*= 344.0f*/ ) +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 is real and 1.0 - oTransferFunction.SetCoeffRI( 0, 1.0f, 0.0f ); + // 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 - 1; i++ ) + for( int i = 1; i < iDFTSize; i++ ) { const float fFrequency = oTransferFunction.GetFrequencyResolution() * i; - const float fCoefficient = CalculateDiffractionFilterCoefficient( dDirect, dDetour, fFrequency, fSpeedOfSound ); + 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); } } -float ITA_PROPAGATION_MODELS_API ITAPropagationModels::Maekawa::CalculateDiffractionFilterCoefficient( const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, const float fSpeedOfSound /*= 344.0f */ ) +double Maekawa::CalculateDiffractionFilterCoefficient( const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, const float fSpeedOfSound /*= ITAConstants::SPEED_OF_SOUND_F */ ) { - // @todo daniel armin - ITA_EXCEPT_NOT_IMPLEMENTED; + //Difference between detour length and direct length + const double dDifference = dPropagationLengthDetour - dPropagationLengthDirect; + + //Fresnel number N according to Maekawa //TODO: Equation number + const double dN = 2.0 * dDifference / fSpeedOfSound * fFrequency; + + //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) / dPropagationLengthDirect; + + + return dFilterCoefficient; } + + +double Maekawa::CalculateDiffractionFactor(const double dPropagationLengthDirect, const double dPropagationLengthDetour, const float fFrequency, const float fSpeedOfSound /*= ITAConstants::SPEED_OF_SOUND_F */) +{ + //Return the diffraction factor without the influence of the propagation path length + return dPropagationLengthDetour * CalculateDiffractionFilterCoefficient(dPropagationLengthDirect, dPropagationLengthDetour, fFrequency, fSpeedOfSound); +} \ No newline at end of file diff --git a/tests/Maekawa/MaekawaTest.cpp b/tests/Maekawa/MaekawaTest.cpp index 873e834..0394e18 100644 --- a/tests/Maekawa/MaekawaTest.cpp +++ b/tests/Maekawa/MaekawaTest.cpp @@ -34,7 +34,7 @@ using namespace ITAGeo; using namespace ITAPropagationModels; static float g_fSampleRate = 44.1e3f; -static int g_iFilterLength = 128; +static int g_iFilterLength = 1280; //! Tests running the Maekawa model functions /** @@ -80,6 +80,11 @@ int main( int, char** ) Maekawa::CalculateDiffractionFilter( dPropLengthDirect, dPropLengthDetour, oDiffractionTF ); cout << "Maekawa filter generation calculation time: " << timeToString( sw.stop() ) << endl; + if (Maekawa::IsApplicable(pS->v3InteractionPoint, pR->v3InteractionPoint, pW)) + { + oDiffractionTF.SetZero(); + } + cout << oDiffractionTF.ToString() << endl; ITAFFTUtils::Export( &oDiffractionTF, "MaekawaTestIR.wav" ); diff --git a/tests/Maekawa/MaekawaTest.m b/tests/Maekawa/MaekawaTest.m index 72d08a1..cf37fa0 100644 --- a/tests/Maekawa/MaekawaTest.m +++ b/tests/Maekawa/MaekawaTest.m @@ -1,6 +1,6 @@ %% Load C++ implementation result maekawa_ir_cpp = ita_read( 'MaekawaTestIR.wav' ); - +maekawa_ir_cpp.signalType = 'energy'; %% Run ITA-Toolbox diffraction simulation n1 = [ -1 1 0 ]; @@ -20,3 +20,7 @@ maekawa_ir.freqData = ita_diffraction_maekawa( w, s, r, maekawa_ir_cpp.freqVecto %% Combined result maekawa_eval = ita_merge( maekawa_ir, maekawa_ir_cpp ); + +%% Plot result + +maekawa_eval.pt \ No newline at end of file -- 2.24.1