diff --git a/include/ITAPropagationModels/Maekawa.h b/include/ITAPropagationModels/Maekawa.h index b114ab17c91697b5b29be2dbd5a5da044aedb18f..2ce9dcc6ab1534b511d2203012c92402ce71f69f 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 8023b4df40ab4ef60174bcc752670d760e1860dc..146bccc858d9e94a66d173431bf8caf504852dba 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 873e8342ad3a6acc5e7527fa81263f421e89e411..0394e18b5bd7f3b7709014a648d4337817491269 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 72d08a1bd90570c93af62224674ebf0f6affb0ca..cf37fa0ffb006c9706a4f6b5cad89985b782499a 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