Skip to content

GitLab

Commit 7edae76d authored by Armin Erraji's avatar Armin Erraji
Browse files

In the function Generate(), added calculation of the sound energy decrease due... 

In the function Generate(), added calculation of the sound energy decrease due to the path length as well as the phase shift.
parent 8fa0729c
Hide whitespace changes
Inline Side-by-side
Showing with 126 additions and 60 deletions
benchmarks/DiffractionModelsBenchmark.cpp
View file @ 7edae76d
......@@ -115,6 +115,7 @@ void benchmark_single_wedge()
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++ )
{
......
include/ITAPropagationModels/FilterEngine.h
View file @ 7edae76d
......@@ -129,9 +129,10 @@ namespace ITAPropagationModels
};
const double m_dHumidity = 80.0;
const double m_dTemperature = 20.0;
const float m_fSpeedOfSound = 344.0f;
const float m_fSampleRate = 44.1e3;
const int m_iFilterLength = 64;
static struct m_DefaultDiffractionModel //!< Default values for diffractions
......
src/ITAPropagationModels/FilterEngine.cpp
View file @ 7edae76d
......@@ -4,6 +4,7 @@
#include <ITAException.h>
#include <ITAGeo/Halfedge/MeshModel.h>
#include <ITAPropagationModels/UTD.h>
#include <ITAISO9613.h>
// Vista includes
#include <VistaMath/VistaGeometries.h>
......@@ -83,21 +84,27 @@ void CFilterEngine::ApplyAcousticModels(ITAGeo::CPropagationPathList & oPathList
void CFilterEngine::ApplyDiffractionModel(ITAGeo::CPropagationPathList & oPathList, int iModel)
{
int iDiffractionModel;
if (iModel == -1)
iDiffractionModel = m_DefaultValues::iDiffractionModel;
else
iDiffractionModel = iModel;
for (auto& oPath : oPathList)
{
//Iterate over all anchors except the first and last anchor
for (size_t i = 1; i <= oPath.size() - 2; i++)
{
auto pLastAnchor = oPath[i - 1];
auto pCurrentAnchor = oPath[i];
auto pNextAnchor = oPath[i + 1];
auto& pLastAnchor = oPath[i - 1];
auto& pCurrentAnchor = oPath[i];
auto& pNextAnchor = oPath[i + 1];
//If the anchor is a reflection and the acoustic material is not set yet, set it to default values
if ((pCurrentAnchor->iAnchorType == ITAGeo::CPropagationAnchor::DIFFRACTION_APEX) && (pCurrentAnchor->pAcousticMaterial == NULL))
{
auto pApex = std::dynamic_pointer_cast<ITAGeo::CITADiffractionWedgeAperture>(pCurrentAnchor);
if (m_DefaultValues::iDiffractionModel == ITAGeo::IAcousticMaterial::THIRD_OCTAVE)
if (iDiffractionModel == ITAGeo::IAcousticMaterial::THIRD_OCTAVE)
{
auto pAcousticMaterial = make_shared<ITAGeo::CThirdOctaveMaterial>();
......@@ -107,12 +114,12 @@ void CFilterEngine::ApplyDiffractionModel(ITAGeo::CPropagationPathList & oPathLi
for (auto fFrequency : vfFrequencies)
{
const float fWaveNumberRad = 2.0f * ITAConstants::PI_F * fFrequency / 344.0; //TODO: Change speed of sound value
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.7;
cfDiffractionCoeff = 0.7f;
vfAbsorptionCoeffs.push_back(1.0f - abs(cfDiffractionCoeff));
}
......@@ -131,6 +138,12 @@ void CFilterEngine::ApplyDiffractionModel(ITAGeo::CPropagationPathList & oPathLi
void CFilterEngine::ApplyReflectionModel(ITAGeo::CPropagationPathList & oPathList, int iModel)
{
int iReflectionModel;
if (iModel == -1)
iReflectionModel = m_DefaultValues::iReflectionModel;
else
iReflectionModel = iModel;
for (auto& oPath : oPathList)
{
for (auto& pAnchor : oPath)
......@@ -138,7 +151,7 @@ void CFilterEngine::ApplyReflectionModel(ITAGeo::CPropagationPathList & oPathLis
//If the anchor is a reflection and the acoustic material is not set yet, set it to default values
if ((pAnchor->iAnchorType == ITAGeo::CPropagationAnchor::SPECULAR_REFLECTION) && (pAnchor->pAcousticMaterial == NULL))
{
if (m_DefaultValues::iReflectionModel == ITAGeo::IAcousticMaterial::SCALAR)
if (iReflectionModel == ITAGeo::IAcousticMaterial::SCALAR)
{
//Set material with default values(current reflection factor: 0.7. See constructor of CScalarMaterial)
pAnchor->pAcousticMaterial = make_shared<ITAGeo::CScalarMaterial>();
......@@ -154,6 +167,12 @@ void CFilterEngine::ApplyReflectionModel(ITAGeo::CPropagationPathList & oPathLis
void ITAPropagationModels::CFilterEngine::ApplyEmitterModel(ITAGeo::CPropagationPathList & oPathList, int iModel)
{
int iEmitterModel;
if (iModel == -1)
iEmitterModel = m_DefaultValues::iEmitterModel;
else
iEmitterModel = iModel;
//Apply the emitter acoustic material. The emitter is always the first anchor of the propagation paths.
for (auto& oPath : oPathList)
{
......@@ -162,7 +181,7 @@ void ITAPropagationModels::CFilterEngine::ApplyEmitterModel(ITAGeo::CPropagation
{
if (oPath[0]->pAcousticMaterial == NULL)
{
if (m_DefaultValues::iEmitterModel == ITAGeo::IAcousticMaterial::SCALAR)
if (iEmitterModel == ITAGeo::IAcousticMaterial::SCALAR)
{
//Set material with factor 1.0
auto pAcousticMaterial = make_shared<ITAGeo::CScalarMaterial>();
......@@ -193,6 +212,12 @@ void ITAPropagationModels::CFilterEngine::ApplySensorModel(ITAGeo::CPropagationP
if (!HasSameSensorAnchor(oPathList))
ITA_EXCEPT1(INVALID_PARAMETER, "The propagation path list has multiple sensor anchors or last anchor of paths is not a sensor.");
int iSensorModel;
if (iModel == -1)
iSensorModel = m_DefaultValues::iSensorModel;
else
iSensorModel = iModel;
//Because, the sensor is always the same anchor, only the first reference must be checked and eventually set
auto& oPath = oPathList[0];
......@@ -201,7 +226,7 @@ void ITAPropagationModels::CFilterEngine::ApplySensorModel(ITAGeo::CPropagationP
//Apply the sensor acoustic material. The sensor is always the last anchor of the propagation paths
if (oPath[oPath.size() - 1]->pAcousticMaterial == NULL)
{
if (m_DefaultValues::iSensorModel == ITAGeo::IAcousticMaterial::SCALAR)
if (iSensorModel == ITAGeo::IAcousticMaterial::SCALAR)
{
//Set material with factor 1.0
auto pAcousticMaterial = make_shared<ITAGeo::CScalarMaterial>();
......@@ -210,17 +235,25 @@ void ITAPropagationModels::CFilterEngine::ApplySensorModel(ITAGeo::CPropagationP
oPath[oPath.size() - 1]->pAcousticMaterial = pAcousticMaterial;
}
else //TODO: add cases like if (m_DefaultValues::iEmitterModel == ITAGeo::IAcousticMaterial::THIRD_OCTAVE)
else //TODO: add cases like if (m_DefaultValues::iSensorModel == ITAGeo::IAcousticMaterial::THIRD_OCTAVE)
{
ITA_EXCEPT_NOT_IMPLEMENTED;
}
}
}
//TODO: delete pragma warning after implementation of function
#pragma warning( push )
#pragma warning( disable : 4100)
void CFilterEngine::ApplyTransmissionModel(ITAGeo::CPropagationPathList & oPathList, int iModel)
{
ITA_EXCEPT_NOT_IMPLEMENTED;
}
//TODO: delete pragma warning pop after implementation of function
#pragma warning( pop )
void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITABase::CHDFTSpectra & oHDFTSpectra)
{
......@@ -228,85 +261,109 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
if (!HasSameSensorAnchor(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);
for (auto& oPath : oPathList)
{
//At first, combine all scalar filter components
//Path length
float fPathLength = (float)oPath.GetLength();
//For the combination of all scalar filter components
float fScalarMagnitude = 1.0f;
//For the combination of all third octave filter components
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctMagnitudes;
oThirdOctMagnitudes.SetIdentity();
//Include sound energy decrease due to the length of the path
fScalarMagnitude *= 1.0f / pow(fPathLength,2);
//Subtract air attenuation according to ISO9613-1
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctAbsorption;
ITABase::ISO9613::AtmosphericAbsorption(oThirdOctAbsorption, fPathLength, m_dTemperature, m_dHumidity);
oThirdOctMagnitudes.Sub(oThirdOctAbsorption);
//Set scalar and third octave filter components of the propagation anchors
for (auto& pAnchor : oPath)
{
if (pAnchor->pAcousticMaterial->GetType() == ITAGeo::IAcousticMaterial::SCALAR)
{
//Type cast of the acoustic material
auto pScalarMaterial = std::dynamic_pointer_cast<ITAGeo::CScalarMaterial>(pAnchor->pAcousticMaterial);
auto pScalarMaterial = dynamic_pointer_cast<ITAGeo::CScalarMaterial>(pAnchor->pAcousticMaterial);
if(pAnchor->iAnchorType == ITAGeo::CPropagationAnchor::SPECULAR_REFLECTION)
fScalarMagnitude *= abs(pScalarMaterial->cdReflectionFactor);
fScalarMagnitude *= (float) abs(pScalarMaterial->cdReflectionFactor);
else
fScalarMagnitude *= abs(pScalarMaterial->cdTransmissionFactor);
fScalarMagnitude *= (float) abs(pScalarMaterial->cdTransmissionFactor);
}
}
//Thereafter, combine all third octave filter components
vector<float> vfThirdOctMagnitudes;
for (auto& pAnchor : oPath)
{
if (pAnchor->pAcousticMaterial->GetType() == ITAGeo::IAcousticMaterial::THIRD_OCTAVE)
{
//Type cast of the acoustic material
auto pThirdOctaveMaterial = std::dynamic_pointer_cast<ITAGeo::CThirdOctaveMaterial>(pAnchor->pAcousticMaterial);
//At the first anchor with third octave material, set all start values to the prior calculated scalar magnitude
//The first and the last values are set to zero
if (vfThirdOctMagnitudes.size() == 0)
{
vfThirdOctMagnitudes.push_back(0.0f);
for (int i = 0; i < pThirdOctaveMaterial->GetNumBands();i++)
{
vfThirdOctMagnitudes.push_back(fScalarMagnitude);
}
vfThirdOctMagnitudes.push_back(0.0f);
}
auto pThirdOctaveMaterial = dynamic_pointer_cast<ITAGeo::CThirdOctaveMaterial>(pAnchor->pAcousticMaterial);
//Multiply the scale values of the respective bands
for (int i = 0; i < pThirdOctaveMaterial->GetNumBands(); i++)
{
vfThirdOctMagnitudes[i+1] *= abs(1 - pThirdOctaveMaterial->oAbsorptionCoefficients[i]);
oThirdOctMagnitudes[i] *= abs(1 - pThirdOctaveMaterial->oAbsorptionCoefficients[i]);
}
}
}
//Set the HDFT spectrum of the current path
ITABase::CHDFTSpectrum oHDFTSpectrum(oHDFTSpectra.GetSampleRate(), oHDFTSpectra.GetDFTSize(), true);
//Reset the HDFT spectra of the current path
m_pTempPropPathSpectra->SetUnity();
//DFT coordinates for the respective frequencies
vector<float> vfDFTCoordinatesIn;
vector<complex<float>> vcfDFTSpectrum;
const int iDFTSize = oHDFTSpectrum.GetSize();
int iHDFTSize = (*m_pTempPropPathSpectra)[0]->GetSize();
//Add zero, because the third octave spectrum doesn't contain 0 Hz
vfDFTCoordinatesIn.push_back(0);
//Get the corresponding DFT coordinates from the third octave coordinates
for (auto& fFreq : ITABase::CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies())
{
vfDFTCoordinatesIn.push_back(fFreq*(iDFTSize - 1) / (oHDFTSpectra.GetSampleRate() / 2));
vfDFTCoordinatesIn.push_back(fFreq*(iHDFTSize - 1) / ( (float) oHDFTSpectra.GetSampleRate() / 2));
}
//Add iDFTSize-1, because the third octave spectrum doesn't contain the threshold value at the half samplerate
vfDFTCoordinatesIn.push_back(iDFTSize - 1);
//Get the DFT spectrum by spline interpolation and set them to the HDFTSpectrum object
for (int i = 0; i < iDFTSize; i++)
//Combine all magnitudes without the eventually existent directivity of the emitter/sensor
for (int iChannel = 0; iChannel < m_pTempPropPathSpectra->GetNumChannels(); iChannel++)
{
vcfDFTSpectrum.push_back(spline_b_val(vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &vfThirdOctMagnitudes[0], i));
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++)
{
const float fDFTSpectrum = spline_b_val( (int)vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &oThirdOctMagnitudes[0], (float) i);
const complex<float> fSpectrumData((*m_pTempPropPathSpectra)[0]->GetData()[2*i], (*m_pTempPropPathSpectra)[0]->GetData()[2*i+1]);
complex<float> cfMultipliedSpectra = fSpectrumData * fDFTSpectrum * fScalarMagnitude;
(*m_pTempPropPathSpectra)[iChannel]->SetCoeff(i, cfMultipliedSpectra);
oHDFTSpectrum.SetCoeff(i, vcfDFTSpectrum[i]);
//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 = fPathLength / m_fSpeedOfSound * (float)oHDFTSpectra.GetSampleRate();
float fDelayPhase = -fDelaySamples* ITAConstants::TWO_PI_F * i / oHDFTSpectra.GetDFTSize();
(*m_pTempPropPathSpectra)[iChannel]->SetPhasePreserveMagnitude(i, fDelayPhase);
}
}
else
{
//Copy first channel to other channels
(*m_pTempPropPathSpectra)[iChannel]->Copy((*m_pTempPropPathSpectra)[0]);
}
}
m_pTempPropPathSpectra->
//Add current spectra to accumulated spectra
m_pAccumulatedSpectra->add(m_pTempPropPathSpectra.get());
}
ITA_EXCEPT_NOT_IMPLEMENTED;
oHDFTSpectra.CopyFrom(m_pAccumulatedSpectra.get());
}
tests/FilterEngine/FilterEngineTest.cpp
View file @ 7edae76d
......@@ -20,6 +20,8 @@
#include <ITAPropagationModels/FilterEngine.h>
#include <ITAGeo/Base.h>
#include <ITAISO9613.h>
using namespace std;
using namespace ITAConstants;
using namespace ITAGeo;
......@@ -27,9 +29,6 @@ using namespace ITAPropagationModels;
int main(int, char**)
{
const float fSpeedOfSound = 344.0f;
const float fSampleRate = 44.1e3;
const int iFilterLength = 64;
// _______
......@@ -81,8 +80,16 @@ int main(int, char**)
oFilter.ApplyAcousticModels(oPathList);
//Set filter length according to the maximum path length
const float fSpeedOfSound = 344.0f; //Approximation of speed of sound at ~20C
const float fSampleRate = 44.1e3;
const int iFilterLength = (int)(oPathList.GetMaxLength()/ fSpeedOfSound * fSampleRate ) + 1024;
ITABase::CHDFTSpectra oSpectra(fSampleRate, pReceiver->iNumChannels, iFilterLength);
oFilter.Generate(oPathList, oSpectra);
}
\ 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