Commit 7cc0c7e1 authored by Philipp Schäfer's avatar Philipp Schäfer
Browse files

- added class to run pidgeon simulations

- added classes for combined atmospheric / urban simulation
parent c7e92280
function [outTF_combined,outTF_separated,atmoPath,urbanPaths] = calcTransferFunct(obj,sourcePosition,receiverPosition,virtualSourceModus)
%RUN runs combined simulation
% Syntax: [outTF_combined,outTF_separated,atmoPath,urbanPaths] = ...
% calcTransferFunct(obj,sourcePosition,receiverPosition,virtualSourceModus)
if nargin < 4
virtualSourceModus = obj.virtualSourceModus;
end
%% run ray tracing
[atmoPath, receiverMissed] = obj.runRayTracing(sourcePosition,receiverPosition);
if receiverMissed
warning('Ray tracing result may be inaccurate.')
end
virtualSourcePosition = obj.findVirtualSource(atmoPath,receiverPosition,virtualSourceModus);
%% run urban propagation
urbanPaths = obj.urbanPropagation.run(virtualSourcePosition,receiverPosition);
%% calculate transfer function
[outTF_combined,outTF_separated] = obj.getTF(atmoPath,urbanPaths,receiverPosition,virtualSourcePosition);
end
function [fileListing] = getAuralizationData(obj,sourcePositions,receiverPosition,samplingRate,outputFolder,freqVector,forceCalculation)
%RUN Generates necessary data for auralization and saves individual frame
%updates
if nargin < 5
outFolder = fullfile(cd,'results');
end
if nargin < 6
freqVector = ita_ANSI_center_frequencies;
end
if nargin < 7
forceCalculation = true;
end
if size(sourcePositions,2) ~= 3
error('Source Positions must be defined as nx3-matrix')
end
nFrames = size(sourcePositions,1);
if ~exist(outputFolder,'dir')
mkdir(outputFolder)
end
%% run rayTracing
wbRT = itaWaitbar(nFrames,'Running Ray Tracing...');
for idFrame = 1:nFrames
atmoPathFile{idFrame} = ['AtmoPath_Frame_',num2str(idFrame),'.mat'];
if ~isfile(fullfile(outputFolder,atmoPathFile{idFrame})) || forceCalculation
[atmoPath, receiverMissed] = obj.runRayTracing(sourcePositions(idFrame,:),receiverPosition);
if receiverMissed
warning('Eigenrays could not be determined. Using previous path.')
load(fullfile(outputFolder,atmoPathFile{idFrame-1}));
end
save(fullfile(outputFolder,atmoPathFile{idFrame}),'atmoPath')
end
load(fullfile(outputFolder,atmoPathFile{idFrame}));
if obj.filter_delays
orig_prop_delay(idFrame) = atmoPath.t(end);
else
prop_delay(idFrame) = atmoPath.t(end);
end
wbRT.inc
end
wbRT.close
%% determine virtual sources
% low-pass filter propagation delays
if obj.filter_delays && size(sourcePositions,1) > 1
diff_prop_delay = lowPassFilter(diff(orig_prop_delay),samplingRate,'lowpassiir',samplingRate/6,samplingRate/3,1e-4);
for idFrame = 1:nFrames
prop_delay(idFrame) = orig_prop_delay(1) + sum(diff_prop_delay(1:idFrame-1));
end
elseif size(sourcePositions,1) == 1
prop_delay = orig_prop_delay;
end
for idFrame = 1:nFrames
load(fullfile(outputFolder,atmoPathFile{idFrame}));
% virtualSourcePosition(idFrame,:) = obj.findVirtualSource(atmoPath,receiverPosition);
distance = prop_delay(idFrame)*obj.atmosphere.c(receiverPosition(3));
spreading_loss(idFrame) = atmoPath.spreadingLoss;
virtualSourcePosition(idFrame,:) = ...
obj.findVirtualSource(atmoPath,receiverPosition,'delay',distance);
end
% low-pass filter virtual source positions and spreading loss
if size(virtualSourcePosition,1) > 1
spreading_loss = filterSpreadingLoss(spreading_loss,samplingRate);
virtualSourcePosition = lowPassFilter(virtualSourcePosition,samplingRate);
end
%% run urban propagation
obj.urbanPropagation.outFilePath = [outputFolder,'\'];
wbUP = itaWaitbar(nFrames,'Running urban propagation...');
for idFrame = 1:nFrames
urbanPathsFile{idFrame} = ['UrbanPaths_Frame_',num2str(idFrame)];
if ~isfile(fullfile(outputFolder,[urbanPathsFile{idFrame},'.json'])) || forceCalculation
obj.urbanPropagation.resultName = urbanPathsFile{idFrame};
obj.runUrbanPropagation(virtualSourcePosition(idFrame,:),receiverPosition);
end
wbUP.inc
end
wbUP.close
%% get and save parameters
TFGenerator = itaCauTFGenerator;
% set atmospheric properties
TFGenerator.constSpeed = obj.atmosphere.c(receiverPosition(3));
TFGenerator.constHumidity = obj.atmosphere.humidity;
TFGenerator.constTemp = obj.atmosphere.T(receiverPosition(3));
TFGenerator.constPressure = obj.atmosphere.p0(receiverPosition(3));
TFGenerator.virtualSourceModus = obj.virtualSourceModus;
wb = itaWaitbar(nFrames,'Calculating DSP parameters...');
for idFrame = 1:nFrames
file_name = sprintf( 'frame%05idata.mat', idFrame );
file_path = fullfile( outputFolder, file_name );
if ~isfile(file_path) || forceCalculation
if idFrame == 1
oldUrbanPaths = [];
end
load(fullfile(outputFolder,atmoPathFile{idFrame}));
urbanPaths = ita_propagation_load_paths([outputFolder,filesep,urbanPathsFile{idFrame},'.json']);
urbanPaths = obj.urbanPropagation.fixPaths(urbanPaths);
% get parameters of atmoPath
[~,~,airAbsorption,propDelay] = atmoTFParameters(atmoPath,obj.atmosphere,freqVector);
TFGenerator.atmoSpreadingLoss = spreading_loss(idFrame);
TFGenerator.atmoPropagationDelay = prop_delay(idFrame);%propDelay;
TFGenerator.atmoAirAbsorption = airAbsorption';
TFGenerator.urbanPaths = urbanPaths;
% determine DSP parameters
[frameData.freqMagnitudes,frameData.resGain,frameData.delay,~,~] = ...
TFGenerator.getAuralizationParameters(freqVector,obj.urbanPropagation.MaxReflectionOrder);
[frameData.srcPhi,frameData.srcTheta,frameData.src_orientation] = anglesAtSrc(urbanPaths,atmoPath);
% save parameters
[paths_update,all_paths_data,source_pos,receiver_pos,srcAngles] = ...
getPathUpdate(urbanPaths,frameData,oldUrbanPaths);
[oldUrbanPaths,hashTable] = ita_propagation_paths_add_identifiers(urbanPaths);
save( file_path, 'paths_update', 'source_pos', 'receiver_pos', 'all_paths_data', 'srcAngles','hashTable' );
end
wb.inc
end
wb.close
fileListing = dir([outputFolder, '/frame*data.mat']);
end
function [outTF_combined,outTF_separated] = getTF(this,atmoPath,urbanPaths,receiverPosition,virtSourcePos,receiver_height)
%CALCTF calculates the combined transfer function
TFGenerator = itaCauTFGenerator;
if nargin < 6
receiver_height = receiverPosition(3);
end
% set atmospheric properties
TFGenerator.constSpeed = this.atmosphere.c(receiver_height);
TFGenerator.constHumidity = this.atmosphere.humidity;
TFGenerator.constTemp = this.atmosphere.T(receiver_height);
TFGenerator.constPressure = this.atmosphere.p0(receiver_height);
% read out TF properties
freqVector = TFGenerator.getFreqVector;
[pathLength,spreadingLoss,airAbsorption,propDelay] = atmoTFParameters(atmoPath,this.atmosphere,freqVector);
TFGenerator.atmoSpreadingLoss = spreadingLoss;
TFGenerator.atmoPropagationDelay = propDelay;
TFGenerator.atmoAirAbsorption = airAbsorption';
TFGenerator.virtualSourceModus = this.virtualSourceModus;
TFGenerator.urbanPaths = urbanPaths;
[outTF_combined,outTF_separated] = TFGenerator.run(virtSourcePos,receiverPosition);
end
classdef itaCauSimulation
%ITA_CAU_AURALIZATION interface class for the combined atmospheric and
%urban auralization
properties
virtualSourceModus = 'delay';
windDirection = [-1 0 0];
atmosphere
rayTracer
urbanPropagation
cityModel
filter_delays = true;
end
%% initialization
methods
function obj = itaCauSimulation()
%ITA_CAU_AURALIZATION Construct an instance of this class
obj.rayTracer = initializeRayTracer;
obj.atmosphere = initializeAtmosphere([-1 0 0]);
obj.urbanPropagation = initializeUrbanPropagation;
end
end
%% public functions
methods(Access = public)
function [atmoPaths,receiverMissed] = runRayTracing(obj,sourcePosition,receiverPosition)
%RUNRAYTRACING runs atmospheric ray tracing
try
atmoPaths = obj.rayTracer.FindEigenrays(obj.atmosphere,sourcePosition,receiverPosition);
minDistancePoint = atmoPaths(1).r.cart(atmoPaths(1).numPoints,:);
distToReceiver = norm(minDistancePoint-receiverPosition);
catch
obj.rayTracer.maxAngleForGeomSpreading = 0.1;
try
atmoPaths = obj.rayTracer.FindEigenrays(obj.atmosphere,sourcePosition,receiverPosition);
minDistancePoint = atmoPaths(1).r.cart(atmoPaths(1).numPoints,:);
distToReceiver = norm(minDistancePoint-receiverPosition);
catch
atmoPaths = [];
distToReceiver = 1e9;
end
obj.rayTracer.maxAngleForGeomSpreading = 0.01;
end
% check if receiver is hit with sufficient accuracy
receiverMissed = 0;
distThreshold = 1;
if distToReceiver > distThreshold
warning(['Warning: Receiver not reached by eigenray. Minimum distance to receiver: ',num2str(distToReceiver),' m.'])
receiverMissed = 1;
end
end
function [urbanPaths] = runUrbanPropagation(obj,sourcePosition,receiverPosition)
urbanPaths = obj.urbanPropagation.run(sourcePosition,receiverPosition);
end
function [outTF, atmoPath, virtSource, urbanPaths] = run(obj,srcPos,recPos)
%RUN runs combined simulation and returns transfer function, atmospheric
%and urban paths and the virtual source position
% run ray tracing
[atmoPath, receiverMissed] = obj.runRayTracing(srcPos,recPos);
if receiverMissed
error('Eigenrays could not be determined.')
end
virtSource = obj.findVirtualSource(atmoPath,recPos);
% run urban propagation
urbanPaths = obj.urbanPropagation.run(virtSource,recPos);
% determine transfer function
outTF = obj.getTF(atmoPath,urbanPaths,recPos,virtSource);
end
[fileListing] = ...
getAuralizationData(obj,sourcePositions,receiverPosition,samplingRate,outputFolder,freqVector,forceCalculation)
[outTF_combined,outTF_separated,atmoPath,urbanPaths] = calcTransferFunct(obj,sourcePosition,receiverPosition)
[freqMagnitudes,resGain,delay,urbanPaths,srcPhi,srcTheta,atmoPath,srcOrientation] ...
= runUpdateStep(obj,sourcePosition,receiverPosition,freqVector,virtualSourceModus)
end
%% hidden function
methods(Access = public, Hidden = true)
[outTF_combined,outTF_separated] = getTF(obj,atmoPath,urbanPaths,receiverPosition,virtSourcePos,receiver_height)
end
%% private functions
methods(Access = private)
function [virtSource] = findVirtualSource(obj,atmoPath,receiverPosition,virtualSourceModus,distOrAtt)
if nargin < 4
virtualSourceModus = obj.virtualSourceModus;
end
if nargin < 5
distOrAtt = [];
end
[currentDirection, currentDelay, currentAttenuation, ~] = extractEndOfRay(atmoPath);
if ~isempty(distOrAtt)
% use input for distance/ attenuation
currentAttenuation = distOrAtt;
currentDistance = distOrAtt;
else
currentDistance = currentDelay*obj.atmosphere.c(receiverPosition(3));
end
% make sure virtual source doesn't end up below ground
if currentDirection(3) > 0
currentDirection(3) = 0;
currentDirection = currentDirection/norm(currentDirection);
end
if strcmp(virtualSourceModus,'spreading loss') || strcmp(virtualSourceModus,'attenuation')
virtSource = ita_cau_virtual_sound_source(receiverPosition, currentDirection, currentAttenuation,'spreading loss');
else
virtSource = ita_cau_virtual_sound_source(receiverPosition, currentDirection, currentDistance,'run time');
end
% set minimum height to 0
if virtSource(3) < 0
warning('Height of virtual source corrected to minimum of 0 m')
virtSource(3) = 0;
end
end
end
%% set functions
methods
function obj = set.windDirection(obj,input)
if ~isnumeric(input) || ~isequal(size(input),[1 3])
error('Wind direction must be defined as 1x3 numeric vector.')
end
obj.windDirection = input;
obj.atmosphere = initializeAtmosphere(input);
end
function obj = set.cityModel(obj,input)
obj.cityModel = input;
fileDir = dir(input);
obj.urbanPropagation.geometryFilePath = [fileDir.folder,'\'];
obj.urbanPropagation.modelName = erase(fileDir.name,'.skp');
end
end
end
function [outPhi,outTheta,srcOrientation] = anglesAtSrc(urbanPaths,atmoPath)
%ANGLESATSRC calculates output angles of an urban path at the source
%position
nPaths = numel(urbanPaths);
for idPath = 1:nPaths
nAnchors = numel(urbanPaths(idPath).propagation_anchors);
if nAnchors ~= 2
% add difference between urban and direct path to ray angle
directPath = urbanPaths(idPath).propagation_anchors{nAnchors}.interaction_point(1:3)' -...
urbanPaths(idPath).propagation_anchors{1}.interaction_point(1:3)';
directAngles = calcAngles(directPath);
firstPathSegment = urbanPaths(idPath).propagation_anchors{2}.interaction_point(1:3)' -...
urbanPaths(idPath).propagation_anchors{1}.interaction_point(1:3)';
segmentAngles = calcAngles(firstPathSegment);
deltaAngles = segmentAngles - directAngles;
else
deltaAngles = zeros(1,2);
end
outPhi(idPath) = atmoPath.phi + deltaAngles(1);
outTheta(idPath) = atmoPath.theta + deltaAngles(2);
end
%% determine src orientation
view = [0 0 1];
up = [0 0 1];
phi = atmoPath.phi/180*pi;
theta = atmoPath.theta/180*pi;
R_phi = [cos(phi) -sin(phi) 0;...
sin(phi) cos(phi) 0;...
0 0 1];
R_theta = [cos(theta) 0 sin(theta);...
0 1 0;...
-sin(theta) 0 cos(theta)];
view = R_phi*R_theta*view';
srcOrientation = [view';up];
end
function [pathLength,spreadingLoss,airAbsorption,propDelay ] = atmoTFParameters(eigenRay,atmosphere,freqVector)
%ATMOTF extracts parameters of the atmospheric transfer function from an
%eigenray
if nargin < 2
atmosphere = GetDefaultAtmosphere;
end
if nargin < 3
freqVector = ita_ANSI_center_frequencies([20 10000], 3);
end
%% read out spreading loss, delay and path length
spreadingLoss = eigenRay.spreadingLoss;
propDelay = eigenRay.t(eigenRay.numPoints);
pathLength = eigenRay.pathLength;
if spreadingLoss * pathLength >= 5
%singularity -> adjust spreading loss to 1/r-law
spreadingLoss = 1/pathLength;
end
%% calculate atmospheric absorption
alphaDB = zeros(size(freqVector));
for idxR = 1:(eigenRay.numPoints-1)
r1 = eigenRay.r.cart(idxR,:);
r2 = eigenRay.r.cart(idxR+1,:);
z = abs(r1(3));
l = norm(r2-r1);
alphaDB = alphaDB + atmosphere.attenuation(z,freqVector)*l;
end
airAbsorption = 10.^(-alphaDB/20);
end
function [angleVector] = calcAngles(directionVector)
%CALCANGLES This functions calculates the incident angles for a set
% of direction vectors as [azimuth,elevation]
numAngles = size(directionVector,1);
% refVector = [1, 0 ,0];
angleVector = zeros(numAngles,2);
for id=1:numAngles
directionVector(id,:) = directionVector(id,:)/norm(directionVector(id,:));
% angleVector(id) = atan2d(norm(cross(directionTmp,refVector)),dot(directionTmp,refVector));
end
[phiTmp,thetaTmp, ~] = cart2sph(directionVector(:,1),directionVector(:,2),directionVector(:,3));
thetaTmp = thetaTmp/pi*180;
phiTmp = phiTmp/pi*180;
% define theta in interval [0°, 90°]
thetaTmp = 90 - thetaTmp;
% define phi in interval [0°, 360°]
if phiTmp < 0
phiTmp = phiTmp + 360;
elseif phiTmp == 360
phiTmp = 0;
end
angleVector = [phiTmp,thetaTmp];
end
function [outDirection, outTime, outAttenuation, outDistance, minDistancePoint] = extractEndOfRay(atmoRay)
% outputs extract run time, propagation direction and spreading loss at
% last point of the ray
% Syntax: [outDirection, outTime, outAttenuation, outDistance, minDistancePoint] = extractEndOfRay(atmoRay)
try
outTime = atmoRay.tMin;
catch
outTime = atmoRay.t(atmoRay.numPoints);
end
minIdx = atmoRay.closestIdxToTime(outTime);
outDistance = atmoRay.pathLength;
outDirection = atmoRay.n.cart(minIdx,:);
minDistancePoint = atmoRay.r.cart(minIdx,:);
try
outAttenuation = atmoRay.spreadingLoss;
catch
disp('Warning: Ray Area at receiver could not be determined! Switched to propagation delay mode.')
outAttenuation = 1/(outTime*340);
end
%
end
function [outData] = filterSpreadingLoss(inData,samplingRate)
%FILTERSPREADINGLOSS removes single peaks from spreading loss data and
%applies low pass filter
%% remove singular peak
max_limit = prctile(inData,99);
outData = min(inData,max_limit);
%% aplly LPF
outData = lowPassFilter(outData,samplingRate,'lowpassiir',samplingRate/6,samplingRate/3,1e-12,80);
end
function [resPoint,relIndex] = findExactDistanceRayPoint(rayPoint1,rayPoint2,targetPoint,reqDistance)
%FINDEXACTDISTACERAYPOINT interpolates between two points of a ray to find
%the point with the exact distance to the targetPoint
%==========================================================================
% resPoint : result of interpolation
% relIndex : index of resPoint divided by number of interpolation
% points; can be used to e.g. estimate the travel time
% between rayPoint1 and resPoint
%--------------------------------------------------------------------------
% rayPoint1,
% rayPoint2 : two closest rayPoints to exact distance point
% targetPoint : target (e.g. receiver) to which the exact distance is
% to be maintained
% reqDistance : required distance to the target point
%==========================================================================
numInterpolationPoints = 100; % resolution for the interpolation
intStep = (rayPoint2-rayPoint1)/numInterpolationPoints;
for id=1:numInterpolationPoints
interpolatedPoints(id,:) = rayPoint1 + id*intStep;