WIP on filter generation, added some more configuration options for propagation filtering

parent a11c732b
......@@ -7,6 +7,7 @@ include( VistaCommon )
# dependencies
vista_use_package( ITABase REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAFFT REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAGeo REQUIRED FIND_DEPENDENCIES )
......
......@@ -50,7 +50,7 @@ namespace ITAPropagationModels
public:
//!< Material resolutions structure (used if not determined by material)
struct MaterialResolution
struct TMaterialResolution
{
ITAGeo::Material::IMaterial::Type iReflection = ITAGeo::Material::IMaterial::SCALAR; //!< Default reflection component resolution
ITAGeo::Material::IMaterial::Type iDiffraction = ITAGeo::Material::IMaterial::THIRD_OCTAVE; //!< Default diffraction component resolution
......@@ -58,7 +58,7 @@ namespace ITAPropagationModels
} oMaterialResolution;
//! Static environment state description
struct StaticEnvironmentState
struct TStaticEnvironmentState
{
double dHumidity = ITAConstants::DEFAULT_HUMIDITY_D;
double dTemperature = ITAConstants::DEFAULT_TEMPERATURE_D;
......@@ -71,11 +71,19 @@ namespace ITAPropagationModels
ISO9613,
};
enum EGeometricalSpreading
{
PLANE_WAVE = 1,
CYLINDRIC_WAVE,
SPHERICAL_WAVE,
};
//! Propagation model configuration description
struct PropagationModelConfiguration
struct TPropagationModelConfiguration
{
int iAtmosphericAbsorption = EAtmoshpericAbsorption::ISO9613;
int iDiffractionModel = DiffractionModels::MAEKAWA_DETOUR_LAW;
int iGeometricalSpreading = SPHERICAL_WAVE;
} oPropagationModelConfiguration; //!< Propagation model configuration
......
......@@ -14,6 +14,7 @@
#include <ITAException.h>
#include <ITAISO9613.h>
#include <ITAFiniteImpulseResponse.h>
#include <ITAFFTUtils.h>
// Vista includes
#include <VistaMath/VistaGeometries.h>
......@@ -21,10 +22,12 @@
// STL includes
#include <assert.h>
#include <stdio.h>
#include <complex>
//Spline
#include <spline.h>
#include "../../ITAFFT/include/ITAFFTUtils.h"
using namespace ITAPropagationModels;
......@@ -195,7 +198,8 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
double dScalarMagnitude; //For the combination of all scalar filter components
ITABase::CThirdOctaveFactorMagnitudeSpectrum oThirdOctFactors; // For the combination of all third octave filter components
oThirdOctFactors.SetIdentity();
ITABase::CFiniteImpulseResponse oIR(oHDFTSpectra.GetDFTSize(), (float)oHDFTSpectra.GetSampleRate(), false);
ITABase::CMultichannelFiniteImpulseResponse oMultiChannelIR( oHDFTSpectra.GetNumChannels(), oHDFTSpectra.GetDFTSize(), ( float ) oHDFTSpectra.GetSampleRate(), false );
//ITABase::CHDFTSpectra oMultichannelIR( oHDFTSpectra.GetDFTSize(), ( float ) oHDFTSpectra.GetSampleRate(), false );
for (auto& oPath : oPathList)
{
......@@ -278,7 +282,8 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
else if (pDAFFDirectivity->GetDomain() == ITADomain::ITA_TIME_DOMAIN)
{
auto pDAFFDirectivityIR = std::dynamic_pointer_cast<const ITAGeo::Directivity::CDAFF_ImpulseResponse>(pSensor->pDirectivity);
pDAFFDirectivityIR->GetNearestNeighbourImpulseResponse(oDirection, oIR);
pDAFFDirectivityIR->GetNearestNeighbourImpulseResponse(oDirection, oMultiChannelIR);
oMultiChannelIR.CyclicShift( pDAFFDirectivityIR->GetMeanTimeOfArrival() );
}
}
//case(ITAGeo::Directivity::IDirectivity::???):
......@@ -307,53 +312,39 @@ void CFilterEngine::Generate(const ITAGeo::CPropagationPathList & oPathList, ITA
}
}
//Reset the HDFT spectra of the current path
m_pTempPropPathSpectra->SetUnity();
// Combine domains
//DFT coordinates for the respective frequencies
vector<float> vfDFTCoordinatesIn;
int iHDFTSize = (*m_pTempPropPathSpectra)[0]->GetSize();
//Get the corresponding DFT coordinates from the third octave coordinates
for (auto& fFreq : ITABase::CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies())
{
vfDFTCoordinatesIn.push_back(fFreq*(iHDFTSize - 1) / ((float)oHDFTSpectra.GetSampleRate() / 2));
}
ITAFFTUtils::Convert( oMultiChannelIR, m_pTempPropPathSpectra.get() );
//Combine all magnitudes without the eventually existent directivity of the emitter/sensor
for (int iChannel = 0; iChannel < m_pTempPropPathSpectra->GetNumChannels(); iChannel++)
{
auto pTempPathSpectrum = (*m_pTempPropPathSpectra)[ iChannel ];
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++)
for( int i = 0; i < pTempPathSpectrum->GetSize(); i++ )
{
// @todo JST/AER: I think this is not correctly using the spline interp, especially tval should be a base value of given frequency, not integer/index
const float fDFTSpectrum = spline_b_val((int)vfDFTCoordinatesIn.size(), &vfDFTCoordinatesIn[0], &oThirdOctFactors[0], (float)i);
// @todo: m_pTempPropPathSpectra is untouched -> unity until here ...
const complex<float> fSpectrumData((*m_pTempPropPathSpectra)[0]->GetData()[2 * i], (*m_pTempPropPathSpectra)[0]->GetData()[2 * i + 1]);
const float fThirdOctaveInterp = spline_b_val( oThirdOctFactors.GetNumBands(), &( oThirdOctFactors.GetCenterFrequencies()[ 0 ] ), &( oThirdOctFactors[ 0 ] ), m_pTempPropPathSpectra->GetFrequencyOfBin( i ) );
complex<float> cfMultipliedSpectra = fSpectrumData * fDFTSpectrum * (float)dScalarMagnitude;
std::complex<float> cfMultipliedSpectra = pTempPathSpectrum->GetCoeff( i ) * fThirdOctaveInterp * ( float ) dScalarMagnitude;
(*m_pTempPropPathSpectra)[0]->SetCoeff(i, cfMultipliedSpectra);
pTempPathSpectrum->SetCoeff( i, cfMultipliedSpectra );
//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 = float(dPathLength / oEnvironmentState.fSpeedOfSound * oHDFTSpectra.GetSampleRate());
float fDelayPhase = -fDelaySamples * ITAConstants::TWO_PI_F * i / (float)oHDFTSpectra.GetDFTSize();
(*m_pTempPropPathSpectra)[0]->SetPhasePreserveMagnitude(i, fDelayPhase);
pTempPathSpectrum->SetPhasePreserveMagnitude( i, fDelayPhase );
}
}
else
{
//Copy first channel to other channels
(*m_pTempPropPathSpectra)[iChannel]->Copy((*m_pTempPropPathSpectra)[0]);
pTempPathSpectrum->Copy( ( *m_pTempPropPathSpectra )[ 0 ] );
}
}
//Add current spectra to accumulated spectra
......
......@@ -54,6 +54,8 @@ int main( int, char** )
oPathFromRight.push_back(pReceiver);
CFilterEngine oFilterEngine( OPENGL );
oFilterEngine.oPropagationModelConfiguration.iAtmosphericAbsorption = CFilterEngine::EAtmoshpericAbsorption::NONE;
oFilterEngine.oPropagationModelConfiguration.iGeometricalSpreading = CFilterEngine::EGeometricalSpreading::PLANE_WAVE;
int iFilterLengthSamples = (int) ceil(Utils::EstimateFilterLengthSamples(oPathFromLeft, fSampleRate));
ITABase::CHDFTSpectra oTransmissionFilterLeft(fSampleRate, pReceiver->GetNumChannels(), iFilterLengthSamples );
......
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