Commit 100dbb8d authored by Armin Erraji's avatar Armin Erraji

The Filter engine now calculates the diffraction factor with the Maekawa method.

parent 4482fe3f
#include <ITAPropagationModels/FilterGenerator.h>
#include <ITAPropagationModels/FilterGenerator.h>
// ITA includes
#include <ITAException.h>
#include <ITAGeo/Halfedge/MeshModel.h>
#include <ITAPropagationModels/UTD.h>
#include <ITAPropagationModels/Maekawa.h>
#include <ITAISO9613.h>
// Vista includes
......@@ -114,13 +115,15 @@ void CFilterEngine::ApplyDiffractionModel(ITAGeo::CPropagationPathList & oPathLi
for (auto fFrequency : vfFrequencies)
{
const float fWaveNumberRad = 2.0f * ITAConstants::PI_F * fFrequency / 344.0f; //TODO: Change speed of sound value
complex <float> cfDiffractionCoeff;
UTD::CalculateDiffractionCoefficient(pLastAnchor->v3InteractionPoint, pNextAnchor->v3InteractionPoint, pApex, fWaveNumberRad, cfDiffractionCoeff, UTD::UTD_APPROX_KAWAI_KOUYOUMJIAN);
//TODO:delete just for testing
cfDiffractionCoeff = 0.7f;
vfAbsorptionCoeffs.push_back(1.0f - abs(cfDiffractionCoeff));
double dDiffractionFactor;
//TODO:Change diffraction calculation, just for testing
double dDirectLength, dDetourLength;
Maekawa::GetDirectLengthAndDetourLength(pLastAnchor->v3InteractionPoint, pNextAnchor->v3InteractionPoint, pApex, dDirectLength, dDetourLength);
dDiffractionFactor= Maekawa::CalculateDiffractionFactor(dDirectLength, dDetourLength, fFrequency, ITAConstants::SPEED_OF_SOUND_F);
//Absorption coefficient alpha = 1 - |D|^2
vfAbsorptionCoeffs.push_back(1.0f - pow(dDiffractionFactor,2));
}
pAcousticMaterial->oAbsorptionCoefficients.SetMagnitudes(vfAbsorptionCoeffs);
......@@ -242,8 +245,6 @@ void ITAPropagationModels::CFilterEngine::ApplySensorModel(ITAGeo::CPropagationP
}
}
//TODO: delete pragma warning after implementation of function
#pragma warning( push )
#pragma warning( disable : 4100)
......@@ -274,16 +275,21 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
float fScalarMagnitude = 1.0f;
//For the combination of all third octave filter components
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctMagnitudes;
oThirdOctMagnitudes.SetIdentity();
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctFactors;
oThirdOctFactors.SetIdentity();
//Include sound energy decrease due to the length of the path
fScalarMagnitude /= pow(fPathLength,2);
//Subtract air attenuation according to ISO9613-1
fScalarMagnitude /= fPathLength;
//Air attenuation according to ISO9613-1
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctAbsorption;
ITABase::ISO9613::AtmosphericAbsorption(oThirdOctAbsorption, fPathLength, m_dTemperature, m_dHumidity);
oThirdOctMagnitudes.Sub(oThirdOctAbsorption);
//Multiply the transmissionfactor of the air attenuation alpha: sqrt(1-|alpha|)
for (int iFrequencyIndex = 0; iFrequencyIndex < oThirdOctFactors.GetNumBands(); iFrequencyIndex++)
{
oThirdOctFactors[iFrequencyIndex] *= sqrt(1 - abs(oThirdOctAbsorption[iFrequencyIndex]));
}
//Set scalar and third octave filter components of the propagation anchors
for (auto& pAnchor : oPath)
......@@ -307,12 +313,11 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
//Multiply the scale values of the respective bands
for (int i = 0; i < pThirdOctaveMaterial->GetNumBands(); i++)
{
oThirdOctMagnitudes[i] *= abs(1 - pThirdOctaveMaterial->oAbsorptionCoefficients[i]);
oThirdOctFactors[i] *= pThirdOctaveMaterial->GetTransmissionFactor(i);
}
}
}
//Reset the HDFT spectra of the current path
m_pTempPropPathSpectra->SetUnity();
......@@ -334,7 +339,7 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
//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++)
{
const float fDFTSpectrum = spline_b_val( (int)vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &oThirdOctMagnitudes[0], (float) i);
const float fDFTSpectrum = spline_b_val( (int)vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &oThirdOctFactors[0], (float) i);
const complex<float> fSpectrumData((*m_pTempPropPathSpectra)[0]->GetData()[2*i], (*m_pTempPropPathSpectra)[0]->GetData()[2*i+1]);
......
......@@ -76,10 +76,11 @@ int main(int, char**)
oPathDiffraction.push_back(pReceiver);
CPropagationPathList oPathList;
//oPathList.push_back(oPathDirect);
//oPathList.push_back(oPathReflection);
//oPathList.push_back(oPathDiffraction);
oPathList.Load("PropagationPathListExample.json");
oPathList.push_back(oPathDirect);
oPathList.push_back(oPathReflection);
oPathList.push_back(oPathDiffraction);
//oPathList.Load("PropagationPathListExample.json");
CFilterEngine oFilter;
......@@ -90,17 +91,17 @@ int main(int, char**)
//Set filter length according to the maximum path length
const float fSpeedOfSound = 344.0f; //Approximation of speed of sound at ~20C
const float fSpeedOfSound = ITAConstants::SPEED_OF_SOUND_F; //Approximation of speed of sound at ~20C
const float fSampleRate = 44.1e3f;
const int iFilterLength = (int)(oPathList.GetMaxLength()/ fSpeedOfSound * fSampleRate );
int iFilterLength = (int)(oPathList.GetMaxLength()/ fSpeedOfSound * fSampleRate );
iFilterLength =(int) iFilterLength * 1.5;
ITABase::CHDFTSpectra oSpectra(fSampleRate, pReceiver->iNumChannels, iFilterLength);
oFilter.Generate(oPathList, oSpectra);
ITAFFTUtils::Export(oSpectra[0], "ImageSource_Example.wav");
ITAFFTUtils::Export(oSpectra[0], "FilterEngine_Example.wav");
}
%% Example with one diffraction and one reflaction component
ThreePath_Example = ita_read('ThreePaths_Example.wav');
ita_plot(ThreePath_Example)
%% Image source example
ImageSource_Example = ita_read('ImageSource_Example.wav');
ImageSource_Example.pf
\ No newline at end of file
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