Commit cee1b125 authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files

Adding itaGeoPropagation class and moving functions

parent 9fd619b7
......@@ -4,8 +4,9 @@ function point_is_outside = point_outside_wedge( obj, point )
assert( ita_diffraction_point_is_of_dim3( point ) );
distance_main_plane = dot( point - obj.location, obj.main_face_normal );
distance_opposite_plane = dot( point - obj.location, obj.opposite_face_normal );
v = ( point - obj.location ) / norm( point - obj.location );
distance_main_plane = dot( v, obj.main_face_normal );
distance_opposite_plane = dot( v, obj.opposite_face_normal );
spatial_threshold = -obj.set_get_geo_eps * 10;
point_is_outside = ( distance_main_plane >= spatial_threshold ) || ...
......
......@@ -39,7 +39,7 @@ if in_shadow_zone
% From Handbook of Acoustics page 117 eq. 4.13 + inverted phase (pressure release first)
H_diffr = ( +1 ) * ( ( 10^(5/20) * sqrt( 2 * pi * N ) ./ tanh( sqrt( 2 * pi * N ) ) ).^(-1) .* H_dir );
else
H_diffr = zeros( size(frequencies) ); % diffraction field outside the shadow zone is not considered with Maekawa's method
H_diffr = zeros( size( frequencies ) ); % diffraction field outside the shadow zone is not considered with Maekawa's method
end
end
classdef itaGeoPropagation < handle
%ITAGEOPROPAGATION Summary of this class goes here
% Detailed explanation goes here
properties
pps;
fs = 44100;
c = 341.0;
diffraction_model = 'utd';
end
properties (Access = protected)
n = 2^15 + 1;
directivity_db = struct();
end
properties (Dependent)
freq_vec;
num_bins;
end
methods
function obj = itaGeoPropagation( fs, num_bins )
if nargin >= 1
obj.fs = fs;
end
if nargin >= 2
obj.n = num_bins;
end
end
function fs = get.fs( obj )
fs = obj.fs;
end
function num_bins = get.num_bins( obj )
num_bins = obj.n;
end
function f = get.freq_vec( obj )
% Returns frequency base vector
% taken from itaAudio (ITA-Toolbox)
if rem( obj.n, 2 ) == 0
f = linspace( 0, obj.fs / 2, obj.n )';
else
f = linspace( 0, obj.fs / 2 * ( 1 - 1 / ( 2 * obj.n - 1 ) ), obj.n )';
end
end
end
end
function [ directivity_id ] = load_directivity( obj, directivity_path, directivity_id_user )
%load_directivity loads a directivity and returns the id
[ ~, directivity_id, directivity_file_ext ] = fileparts( directivity_path );
if nargin >= 3
directivity_id = directivity_id_user;
end
if strcmpi( directivity_file_ext, '.daff' )
obj.directivity_db.( directivity_id ) = DAFF( directivity_path );
else
error( 'Could not load directivity, unrecognized file extension "%s"', directivity_file_ext )
end
end
function load_paths( obj, json_file_path )
json_txt = fileread( json_file_path );
json_struct = jsondecode( json_txt );
if ~isfield( json_struct, 'class' )
error( 'Could not import propagation path list from file %f, structure is missing field "class"', json_file_path )
end
assert( strcmpi( json_struct.class, 'propagation_path_list' ) || strcmpi( json_struct.class, 'PropagationPathList' ) )
if ~isfield( json_struct, 'propagation_paths' )
error( 'Could not import propagation path list from file %f, structure is missing field "propagation_paths"', json_file_path )
end
obj.pps = json_struct.propagation_paths;
end
function phase_by_delay = ita_propagation_delay( distance, c, fs, fft_degree )
%ITA_RPOPAGATION_SPREADING_LOSS Calculates spreading loss for different
%wave types for a given distance (straight line between emitter center
%point and sensor point)
if nargin ~= 4
error 'Expecting 4 arguments'
end
function phase_by_delay = phase_delay( obj, distance )
%PHASE_DELAY Calculates phase delay for the distance for all frequency bins
if distance <= 0
error 'Distance cannot be zero or negative'
end
delay_tf = itaAudio();
delay_tf.samplingRate = fs;
delay_tf.fftDegree = fft_degree;
f = delay_tf.freqVector;
lambda = c ./ f( 2:end ); % Wavelength
lambda = obj.c ./ obj.freq_vec( 2:end ); % Wavelength
k = 2 * pi ./ lambda; % Wavenumber
phase_by_delay = [ 0; exp( -1i .* k .* distance ) ]; % Note: DC value set to ZERO
......
function [ freq_data_linear ] = run( obj )
%RUN Calculates the transfer function (tf) of the superimposed (geometrical) propagation path in frequency domain
freq_data_linear = ones( obj.num_bins, 1 );
% Iterate over propagation paths, calculate transfer function and sum up
for n = 1:numel( obj.pps )
pp = obj.pps( n );
pp_tf = obj.tf( pp );
freq_data_linear = freq_data_linear + pp_tf;
end
end
function [ freq_data_linear ] = ita_propagation_tf( propagation_path, fs, fft_degree, c )
%ITA_PROPAGATION_PATH_TF Calculates the transfer function (tf)
%of a (geometrical) propagation path in frequency domain for a
%given sampling rate and fft degree. Provide speed of sound, see also
%ita_propagation_load_defaults
%
function [ freq_data_linear ] = tf( obj, pp )
%TFS Calculates the transfer functions (tfs) of the (geometrical) propagation paths in frequency domain
if ~isfield( propagation_path, 'propagation_anchors' )
if ~isfield( obj.pps, 'propagation_anchors' )
error( 'The propagation_path argument does not contain a field "propagation_anchors"' )
end
N = numel( propagation_path.propagation_anchors );
N = numel( pp.propagation_anchors );
if N < 2
error( 'Propagation path has less than two anchor points, cannot calculate a transfer function' )
end
freq_data_linear = ones( obj.num_bins, 1 );
prop_tfs = itaAudio;
prop_tfs.samplingRate = fs;
prop_tfs.fftDegree = fft_degree;
prop_tfs.freqData = ones( prop_tfs.nBins, N );
%lambda = obj.c ./ obj.freq_vec( 2:end )'; % Wavelength
%k = ( 2 * pi ) ./ lambda; % Wavenumber
lambda = c ./ prop_tfs.freqVector( 2:end )'; % Wavelength
k = (2 * pi) ./ lambda; % Wavenumber
distance = ita_propagation_path_length( propagation_path );
if distance / c > prop_tfs.trackLength
error 'Propagation path length too long, increase fft degree to generate transfer function for this propagation path'
distance = ita_propagation_path_length( pp );
if distance / obj.c > 2 * obj.num_bins / obj.fs
error 'Propagation path length too long, increase number of bins to generate transfer function for this propagation path'
end
phase_by_delay = ita_propagation_delay( distance, c, fs, fft_degree );
spreading_loss = ita_propagation_spreading_loss( distance, 'spherical' );
freq_data_linear = phase_by_delay .* spreading_loss;
incident_spreading_loss_applied = false;
for m = 1 : N
if N > 2
anchor = propagation_path.propagation_anchors{ m };
anchor = pp.propagation_anchors{ m };
else
anchor = propagation_path.propagation_anchors( m );
anchor = pp.propagation_anchors( m );
end
assert( strcmpi( anchor.class, 'propagation_anchor' ) )
......@@ -49,22 +38,33 @@ for m = 1 : N
if m == N
if N > 2
target_pos = propagation_path.propagation_anchors{ N - 1 }.interaction_point;
target_pos = pp.propagation_anchors{ m - 1 }.interaction_point;
else
target_pos = propagation_path.propagation_anchors( N - 1 ).interaction_point;
target_pos = pp.propagation_anchors( m - 1 ).interaction_point;
end
% If not already applied by a corresponding anchor type
% (i.e. diffraction), include incident field now
if ~incident_spreading_loss_applied
effective_source_distance = ita_propagation_effective_source_distance( pp, m );
phase_by_delay = obj.phase_delay( effective_source_distance );
spreading_loss = ita_propagation_spreading_loss( distance, 'spherical' );
freq_data_linear = freq_data_linear .* phase_by_delay .* spreading_loss;
incident_spreading_loss_applied = true;
end
else
if N > 2
target_pos = propagation_path.propagation_anchors{ m + 1 }.interaction_point;
target_pos = pp.propagation_anchors{ m + 1 }.interaction_point;
else
target_pos = propagation_path.propagation_anchors( m + 1 ).interaction_point;
target_pos = pp.propagation_anchors( m + 1 ).interaction_point;
end
end
target_position = target_pos - anchor.interaction_point; % Outgoing direction vector
prop_tfs.freqData( :, m ) = ita_propagation_directivity( anchor, target_position / norm( target_position ), fs, fft_degree );
target_position_relative = target_pos - anchor.interaction_point; % Outgoing direction vector
freq_data_linear = freq_data_linear .* obj.tf_directivity( anchor, target_position_relative / norm( target_position_relative ) );
case 'specular_reflection'
......@@ -72,46 +72,51 @@ for m = 1 : N
error( 'Detected a specular reflection at beginning or end of propagation path.' )
end
source_pos = propagation_path.propagation_anchors{ m - 1 }.interaction_point;
target_pos = propagation_path.propagation_anchors{ m + 1 }.interaction_point;
source_pos = pp.propagation_anchors{ m - 1 }.interaction_point;
target_pos = pp.propagation_anchors{ m + 1 }.interaction_point;
effective_source_position = anchor.interaction_point - source_pos;
target_position = target_pos - anchor.interaction_point;
target_position_relative = target_pos - anchor.interaction_point;
incident_direction_vec = effective_source_position / norm( effective_source_position );
emitting_direction_vec = target_position / norm( target_position );
emitting_direction_vec = target_position_relative / norm( target_position_relative );
freq_data_linear = freq_data_linear .* obj.tf_reflection( anchor, incident_direction_vec, emitting_direction_vec );
prop_tfs.freqData( :, m ) = ita_propagation_specular_reflection( anchor, incident_direction_vec, emitting_direction_vec, fs, fft_degree );
case { 'outer_edge_diffraction', 'inner_edge_diffraction' }
if m == 1 || m == N
error( 'Detected a diffraction at beginning or end of propagation path.' )
end
source_pos = propagation_path.propagation_anchors{ m - 1 }.interaction_point;
target_pos = propagation_path.propagation_anchors{ m + 1 }.interaction_point;
source_pos = pp.propagation_anchors{ m - 1 }.interaction_point;
target_pos = pp.propagation_anchors{ m + 1 }.interaction_point;
source_direction = ( source_pos - anchor.interaction_point ) / norm( source_pos - anchor.interaction_point );
target_direction = ( target_pos - anchor.interaction_point ) / norm( target_pos - anchor.interaction_point );
effective_source_distance = ita_propagation_effective_source_distance( propagation_path, m );
effective_target_distance = ita_propagation_effective_target_distance( propagation_path, m );
effective_source_distance = ita_propagation_effective_source_distance( pp, m );
effective_target_distance = ita_propagation_effective_target_distance( pp, m );
effective_source_position = anchor.interaction_point + source_direction * effective_source_distance;
effective_target_position = anchor.interaction_point + target_direction * effective_target_distance;
prop_tfs.freqData( :, m ) = ita_propagation_diffraction( anchor, effective_source_position, effective_target_position, 'utd', fs, fft_degree );
freq_data_linear = freq_data_linear .* obj.tf_diffraction( anchor, effective_source_position, effective_target_position );
% If not already applied by another diffraction edge, include incident field
if ~incident_spreading_loss_applied
phase_by_delay = obj.phase_delay( effective_source_distance );
spreading_loss = ita_propagation_spreading_loss( distance, 'spherical' );
freq_data_linear = phase_by_delay .* spreading_loss;
incident_spreading_loss_applied = true;
end
otherwise
sprintf( 'Detected unrecognized anchor type "%s", attempting to continue', anchor.anchor_type )
end
freq_data_linear = freq_data_linear .* prop_tfs.freqData( :, m );
end
end
function [ linear_freq_data ] = tf_diffraction( obj, anchor, effective_source_position, effective_receiver_position )
%TF_REFLECTION Returns the specular reflection transfer function for a diffraction point to an effective receiver point with
% an incident and emitting direction
if nargin < 5
diffraction_model = obj.diffraction_model;
end
if ~isfield( anchor, 'anchor_type' )
error( 'The anchor argument does not contain a field "anchor_type"' )
end
linear_freq_data = ones( obj.num_bins, 1 );
% Assemble wedge
n1 = anchor.main_wedge_face_normal( 1:3 )';
n2 = anchor.opposite_wedge_face_normal( 1:3 )';
loc = anchor.vertex_start( 1:3 )';
endPt = anchor.vertex_end( 1:3 )';
len = norm( endPt - loc );
aperture_dir = ( endPt - loc ) / len;
% check if wedge is a screen
if ~strcmpi( diffraction_model, 'btms' )
if abs( cross(n1, n2) ) < itaInfiniteWedge.set_get_geo_eps
w = itaSemiInfinitePlane( n1, loc, aperture_dir );
else
w = itaInfiniteWedge( n1, n2, loc );
end
else
w = itaFiniteWedge( n1, n2, loc, len );
w.aperture_direction = aperture_dir;
w.aperture_end_point = endPt;
end
% Legacy
if size( effective_source_position, 1 ) == 4
effective_source_position = effective_source_position( 1:3 )';
end
if size( effective_receiver_position, 1 ) == 4
effective_receiver_position = effective_receiver_position( 1:3 )';
end
% Return only the diffraction component and the subsequent
% sphere-cylinder-shaped wave spreading loss factor to next effective
% receiver (ignore effective source pos)
switch( diffraction_model )
case 'utd'
linear_freq_data = obj.tf_diffraction_utd( w, effective_source_position, effective_receiver_position );
case 'maekawa'
linear_freq_data = obj.tf_diffraction_maekawa( w, effective_source_position, effective_receiver_position );
case 'btms'
btms_ir = obj.tf_diffraction_btms( w, effective_source_position( 1:3 )', effective_receiver_position( 1:3 )' );
diffraction_dft = fft( btms_ir, obj.num_bins * 2 - 1 ); % odd DFT spectrum
diffraction_hdft = diffraction_dft( 1:ceil( obj.num_bins ) );
apex_point = w.get_aperture_point( effective_source_position( 1:3 )', effective_receiver_position( 1:3 )' );
eff_source_distance = norm( apex_point - effective_source_position );
spreading_loss = ita_propagation_spreading_loss( eff_source_distance );
phase_delay = obj.phase_delay( eff_source_distance );
normilization_tf = spreading_loss * phase_delay;
linear_freq_data = diffraction_hdft ./ normilization_tf;
otherwise
warnning 'Unknown diffraction model, returning neutral transfer function'
end
end
function [ linear_freq_data ] = tf_diffraction_maekawa( wedge, eff_source_pos, eff_receiver_pos )
%TF_DIFFRACTION_UTD Calculates the diffraction filter based on uniform
%theory of diffraction (with Kawai approximation).
[ H_diffr, ~ ] = ita_diffraction_maekawa( wedge, eff_source_pos( 1:3 ), eff_receiver_pos( 1:3 ), obj.freq_vec( 2:end ), obj.c );
eff_direct = norm( eff_source_pos( 1:3 ) - eff_source_pos( 1:3 ) );
H_eff = 1 ./ eff_direct;
linear_freq_data = [ 0; H_diffr ./ H_eff ]; % Normalizing incident field, it's taken care of by previous segment.
end
function [ linear_freq_data ] = tf_diffraction_utd( obj, wedge, source_pos, receiver_pos )
%TF_DIFFRACTION_UTD Calculates the diffraction filter based on uniform
%theory of diffraction (with Kawai approximation).
[ ~, D, A ] = ita_diffraction_utd( wedge, source_pos, receiver_pos, obj.freq_vec( 2:end ), obj.c );
linear_freq_data = [ 0; A .* D ];
end
function [ linear_freq_data ] = tf_directivity( obj, anchor, wave_front_direction )
%TF_DIRECTIVITY Returns the directivity transfer function for an anchor and
%a target (incoming or outgoing wave front direction relative (not in world coordinates!) to anchor point)
linear_freq_data = ones( obj.num_bins, 1 );
if ~isfield( anchor, 'directivity_id' )
return
end
if ~isfield( obj.directivity_db, anchor.directivity_id )
warning( 'Directivity id "%s" not found in database, skipping directivity tf calculation', anchor.directivity_id )
return
end
directivity_data = obj.directivity_db.( anchor.directivity_id );
if isa( directivity_data, 'DAFF' )
q_object = quaternion( anchor.orientation );
v = wave_front_direction( 1:3 ) / norm( wave_front_direction( 1:3 ) );
q_target = quaternion.rotateutov( [ 1 0 0 ], v );
q_combined = q_target * conj( q_object );
euler_angles = q_combined.EulerAngles( 'zxy' );
azi_deg = rad2deg( euler_angles( 1 ) );
ele_deg = rad2deg( euler_angles( 2 ) );
idx = directivity_data.nearest_neighbour_index( azi_deg, ele_deg );
if strcmpi( directivity_data.properties.contentType, 'ir' )
directivity_ir = directivity_data.record_by_index( idx )';
directivity_dft = fft( directivity_ir, obj.num_bins * 2 - 1 ); % odd DFT length
directivity_hdft = directivity_dft( 1:( ceil( obj.num_bins ) ) );
linear_freq_data = directivity_hdft;
else
warning( 'Unrecognized DAFF content type "%s" of directivity with id "%s"', directivity_data.properties.contentType, anchor.directivity_id )
end
else
warning( 'Unrecognized directivity format "%s" of directivity with id "%s"', class( directivity_data ), anchor.directivity_id )
end
end
function [ linear_freq_data ] = tf_reflection( obj, anchor, incident_direction_vec, emitting_direction_vec )
%TF_REFLECTION Returns the specular reflection transfer function for an reflection point anchor with
% an incident and emitting direction
linear_freq_data = ones( obj.num_bins, 1 );
if ~isfield( anchor, 'material_id' )
return
end
if ~isfield( obj.material_db, anchor.material_id )
warning( 'Material id "%s" not found in database, skipping reflection tf calculation', anchor.material_id )
return
end
material_data = obj.material_db.( anchor.material_id );
if isa( material_data, 'struct' )
else
warning( 'Unrecognized material format "%s" of material with id "%s"', class( material_data ), anchor.material_id )
end
end
......@@ -3,7 +3,7 @@ function [ freq_data_linear ] = ita_propagation_diffraction( anchor, effective_s
%of the edge diffraction for certain model and given effective in/out positions
%sampling rate and fft degree (defaults to fs = 44100 and fft_degree = 15)
%
if nargin < 3
error 'You are missing several input argument'
end
......@@ -53,6 +53,7 @@ else
w = itaInfiniteWedge( n1, n2, loc );
end
% Legacy
if size( effective_source_position, 1 ) == 4
effective_source_position = effective_source_position( 1:3 )';
......
function [ freq_data_linear ] = ita_propagation_directivity( anchor, target_direction, fs, fft_degree )
%ITA_PROPAGATION_DIRECTIVITY Calculates the directivity function
% in frequency domain for a given target direction, sampling rate and fft degree (defaults to fs = 44100 and fft_degree = 15)
% Can be used for emitter and sensor objects.
if nargin < 2
error 'You are missing second argument "target_direction"'
end
if nargin < 3
fs = 44100;
end
if nargin < 4
fft_degree = 15;
end
if ~isfield( anchor, 'anchor_type' )
error( 'The anchor argument does not contain a field "anchor_type"' )
end
warning( 'ita_propagation_directivity not implemented yet, returning neutral transfer function values' )
% @todo select target_direction data set from directivity database and interpolate in frequency
% domain, if necessary
directivity_tf = itaAudio();
directivity_tf.samplingRate = fs;
directivity_tf.fftDegree = fft_degree;
directivity_tf.freqData = ones( directivity_tf.nBins, 1 );
freq_data_linear = directivity_tf.freqData;
end
function [ pps_with_directivity ] = ita_propagation_paths_set_directivity( pps, anchor_type, directivity_id )
%
% ita_propagation_paths_set_directivity Assigns all given anchor types a
% directivity id
%%
% Example: pps_with_directivity = ita_propagation_paths_set_directivity( pps, 'emitter', 'Genelec8020' )
%
if nargin < 3
error( 'Too few arguments, need propagation path list and directivity id' )
end
if isempty( pps )