Commit d4053ba3 authored by jonasstienen's avatar jonasstienen
Browse files

Improving propagation diffraction handling and path management. modified path...

Improving propagation diffraction handling and path management. modified path identifier pattern! Now more clear.
parent c6459cf3
......@@ -8,40 +8,48 @@ assert( numel( receiver_pos ) == 3 )
%% Calculations
% Based on a line-plane intersection
source_receiver_vec = receiver_pos - source_pos;
source_receiver_dir = source_receiver_vec / norm( source_receiver_vec );
edge_dir = obj.edge_direction / norm( obj.edge_direction );
if norm( source_receiver_vec ) ~= 0
if norm( source_receiver_vec ) == 0
warning( '@todo auxiliar plane must be created differently if source and receiver positions are equal. Trying to continue.' )
end
% Based on a line-plane intersection
source_receiver_dir = source_receiver_vec / norm( source_receiver_vec );
% Auxilary plane spanned by source_receiver_dir and aux_plane_dir (closest
% line between aperture vector and source-receiver-vector
aux_plane_vec = cross( source_receiver_dir, edge_dir );
% Auxilary plane spanned by source_receiver_dir and aux_plane_dir (closest
% line between edge vector and source-receiver-vector
aux_plane_vec = cross( source_receiver_dir, edge_dir );
if norm( aux_plane_vec ) ~= 0
% Directions are not parallel, a closest point must exist
aux_plane_dir = aux_plane_vec / norm( aux_plane_vec );
aux_plane_normal = cross( source_receiver_dir, aux_plane_dir );
% Determine intersection of line (aperture) and auxiliary plane
lambda_divisor = dot( aux_plane_normal, edge_dir );
assert( lambda_divisor ~= 0 )
d = dot( aux_plane_normal, obj.location ); % Distance to origin
lambda = ( d - dot( aux_plane_normal, source_pos ) );% / lambda_divisor; % .. or receiver_pos
ap = obj.location + lambda * edge_dir;
else
% Directions are parallel, project middle point between source & target
% onto aperture
% Project middle point onto aperture
lambda = dot( ( source_pos + source_receiver_vec ./ 2 ) - obj.location, edge_dir );
ap = obj.location + lambda * edge_dir;
end
if norm( aux_plane_vec ) ~= 0
% Directions are no parallel, a closest point must exist
aux_plane_dir = aux_plane_vec / norm( aux_plane_vec );
aux_plane_normal = cross( source_receiver_dir, aux_plane_dir );
% Determine intersection of line (aperture) and auxiliary plane
lambda_divisor = dot( aux_plane_normal, edge_dir );
assert( lambda_divisor ~= 0 )
d = dot( aux_plane_normal, obj.location ); % Distance to origin
lambda = ( d - dot( aux_plane_normal, source_pos ) );% / lambda_divisor; % .. or receiver_pos
ap = obj.location + lambda * edge_dir;
else
% Directions are parallel, project middle point between souce & target
% onto aperture
% Special case where source & target coincide.
% Project middle point onto aperture
lambda = dot( ( source_pos + source_receiver_vec ./ 2 )' - obj.location, edge_dir );
% Project position on edge.
lambda = dot( source_pos - obj.location, edge_dir );
ap = obj.location + lambda * edge_dir;
end
......
function apx = apex_point_approx( obj, source_pos, receiver_pos, spatial_precision )
%apex_point_approx approximates the shortest path by a minimisation
%approach based on the Euclidean distance.
if nargin < 4
spatial_precision = 1e-3; % mm
end
......@@ -7,27 +8,25 @@ function apx = apex_point_approx( obj, source_pos, receiver_pos, spatial_precisi
ap_start = obj.location;
ap_dir = obj.edge_direction;
S_on_ap = orthogonal_projection( ap_start, ap_dir, source_pos ); %project the source to the aperture
S_t = (S_on_ap - ap_start) / ap_dir;
S_t = S_t(1,1); %find the parametric distance along the aperture
R_on_ap = orthogonal_projection( ap_start, ap_dir, receiver_pos ); %same as above but for receiver
R_t = (R_on_ap - ap_start) / ap_dir;
R_t = R_t(1,1);
S_on_ap = orthogonal_projection( ap_start, ap_dir, source_pos ); % project the source to the edge
S_t = norm( S_on_ap - ap_start ); % parametric distance along the aperture
R_on_ap = orthogonal_projection( ap_start, ap_dir, receiver_pos ); % same as above but for receiver
R_t = norm( R_on_ap - ap_start ); % same
start_t = min( S_t, R_t ); %start the optimisation at whichever of the projected source/ receiver comes first on the aperture
end_t = max( S_t, R_t ); %finish at the other
opts_t = optimset( 'TolX', spatial_precision, 'TolFun', spatial_precision, 'FunValCheck', 'on' );
t = fminbnd(@(t)total_path_distance(t, source_pos, receiver_pos, ap_start, ap_dir), start_t, end_t, opts_t );
t = fminbnd( @(t)total_path_distance( t, source_pos, receiver_pos, ap_start, ap_dir ), start_t, end_t, opts_t );
apx = ap_start + t .* ap_dir; %using the optimised parameter, find the aperture position
end
function dist = total_path_distance(t, source_pos, receiver_pos, start, dir)
P = start + (t*dir); %P = point on the aperture
dist = norm(P - source_pos) + norm(receiver_pos - P); %given point on aperture, source and receiver positions, calculate the distance traveled
function dist = total_path_distance( t, source_pos, receiver_pos, start, dir )
P = start + ( t * dir ); % P = point on the edge
dist = norm( P - source_pos ) + norm( receiver_pos - P ); % given point on aperture, source and receiver positions, calculate the distance traveled
end
function point_on_line = orthogonal_projection(line_point,line_dir,point)
point_on_line = line_point + dot(point-line_point,line_dir) / dot(line_dir,line_dir) * line_dir;
function point_on_line = orthogonal_projection( line_point, line_dir, point )
point_on_line = line_point + dot( point - line_point, line_dir ) / dot( line_dir, line_dir ) * line_dir;
end
\ No newline at end of file
%% init
%params
c = 344;
fs = 44100;
fftDegree = 12;
% wedges
n1 = [ 1, 1, 0];
n2 = [-1, 1, 0];
apexStart = [0, 1, -4];
apexEnd = [0, 1, 4];
apexLen = norm(apexEnd -apexStart);
apexDir = (apexEnd -apexStart) / apexLen;
infScreen = itaSemiInfinitePlane(n1, apexStart, apexDir);
infWedge = itaInfiniteWedge(n1, n2, apexStart);
finWedge = itaFiniteWedge(n1, n2, apexStart, apexLen);
% interaction points
src = [-3, 0, 0];
rcv = [ 3, 0, 0];
% result variables
diffr_tf = itaAudio();
diffr_tf.fftDegree = fftDegree;
diffr_tf.samplingRate = fs;
diffr_tf_maekawa = diffr_tf;
diffr_tf_maekawa.channelNames = {'maekawa'};
diffr_tf_utd = diffr_tf;
diffr_tf_utd.channelNames = {'utd'};
diffr_tf_btms = diffr_tf;
diffr_tf_btms.channelNames = {'btms'};
%% diffraction
diffr_tf_maekawa.freqData = ita_diffraction_maekawa(infScreen, src, rcv, diffr_tf_maekawa.freqVector, c);
diffr_tf_utd.freqData = ita_diffraction_utd(infWedge, src, rcv, diffr_tf_maekawa.freqVector, c);
diffr_tf_btms.timeData = ita_diffraction_btms(finWedge, src, rcv, fs, diffr_tf_maekawa.nSamples, c);
diffr_tf = ita_merge(diffr_tf_maekawa, diffr_tf_utd, diffr_tf_btms);
diffr_tf_norm = ita_normalize_spk(diffr_tf, 'allchannels');
diffr_tf_norm.channelNames = [{'normalized maekawa'}; {'normalized utd'}; {'normalized btms'}];
diffr_tf = ita_merge(diffr_tf, diffr_tf_norm);
%% plot
diffr_tf.pf;
title('diffraction filters at simple rectangular wedge');
ylim auto
%% Total sound field of a difraction problem at a rectagular wedge using UTD
% Author: Jonas Stienen
%% Scene
n1 = [ -1 1 0 ] / sqrt( 2 ); % main face
n2 = [ +1 1 0 ] / sqrt( 2 );
loc = [ 0 0 0 ];
w = itaInfiniteWedge( n1, n2, loc );
s = 5 * [ -1 0 0 ];
s_1 = 5 * [ 0 -1 0 ]; % image source n1
r1 = 5 * [ -1 +1 0 ] / sqrt( 2 );
a = w.apex_point( s, r1 );
ad = w.edge_direction;
c = 344; % speed of sound
f = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz', '800 Hz', ...
'1.6 kHz', '3.2 kHz', '6.4 kHz', '12.8 kHz', '24 kHz' };
num_angles = 1000;
alpha_d_start = w.opening_angle;
alpha_d_end = 0;
alpha_d = linspace( alpha_d_start, alpha_d_end, num_angles );
% Set different receiver positions rotated around the aperture
R = norm( r1 ) * [ cos( alpha_d - pi/4 ); sin( alpha_d - pi/4 ); zeros( 1, numel( alpha_d ) ) ]';
% wavenumber
k = 2 * pi * f ./ c;
% store results in ita audio format
res_data_template = itaResult();
res_data_template.resultType = 'signal';
res_data_template.freqVector = f;
%% Sound field simulation
N = size( R, 1 );
for n = 1 : N
r = R( n, 1:3 );
r_dir = norm( r - s );
H_i = 1 ./ r_dir .* exp( -1i .* k .* r_dir ); % direct sound component
r_1_dir = norm( r - s_1 );
H_r = 1 ./ r_1_dir .* exp( -1i .* k .* r_1_dir ); % reflected sound component
H_d = ita_diffraction_utd( w, s, r, f, c ); % diffracted sound component
H = H_d; % Total sound field ... diffracted component always present
if ~ita_diffraction_shadow_zone( w, s, r )
H = H + H_i;
end
if ita_diffraction_reflection_zone( w, s, r, true ) % main face reflection
H = H + H_r;
end
att_sum = res_data_template;
att_sum.freqData = H ./ H_i; % Insertion loss
if n == 1
att_sum_all = att_sum;
else
att_sum_all = ita_merge( att_sum_all, att_sum );
end
end
%% Plot
fg = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );
angles = rad2deg( linspace( alpha_d( end ), alpha_d( 1 ), numel( alpha_d ) ) ); % Reverse angle on axis
h = plot( angles, flip( att_sum_all.freqData_dB', 2 ) );
title( 'Transfer function nomalized to free-field transmission' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )
lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;
xlim( [ -5 275 ] );
ylim( [ -35 10 ] );
grid on;
lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
h( i ).LineWidth = lw( i );
end
%% Export
set( fg, 'papersize', [ 16*2 9*2 ], 'units', 'centimeters' )
print( 'ita_diffraction_plot_utd_total_field_stienen_diss', '-fillpage', '-dpdf' )
%% Total sound field of a difraction problem at a rectagular wedge using UTD
% Author: Jonas Stienen
%% Scene
n1 = [ 0 1 0 ]; % main face 1
n2 = [ 1 0 0 ]; % opposite face 2 / facade
n3 = [ 0 1 0 ]; % ground plane / face 3
loc = [ 0 0 0 ];
w = itaInfiniteWedge( n1, n2, loc );
r = [ 2 -2 0 ];
r2 = [ -2 -2 0 ]; % image receiver 2
r3 = [ 2 -4 0 ]; % image receiver 3
r23 = [ -2 -4 0 ]; % image receiver 23 / 32
s1 = 5 * [ 0 -1 0 ];
a = w.apex_point( s1, r );
ad = w.edge_direction;
c = 344; % speed of sound
%f = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
%f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz', '800 Hz', ...
%'1.6 kHz', '3.2 kHz', '6.4 kHz', '12.8 kHz', '24 kHz' };
f = [ 50, 500, 5000 ]';
f_legend = { '50 Hz', '500 Hz', '5000 Hz' };
num_angles = 1000;
alpha_d_start = 0;
alpha_d_end = pi;
alpha_d = linspace( alpha_d_start, alpha_d_end, num_angles );
% Set different receiver positions rotated around the aperture
S = norm( s1 ) * [ - cos( alpha_d ); sin( alpha_d ); zeros( 1, numel( alpha_d ) ) ]';
% wavenumber
k = 2 * pi * f ./ c;
% store results in ita audio format
res_data_template = itaResult();
res_data_template.resultType = 'signal';
res_data_template.freqVector = f;
%% Sound field simulation
N = size( S, 1 );
for n = 1 : N
s = S( n, 1:3 );
% Zero order
r_dir = norm( s - r );
H_i = 1 ./ r_dir .* exp( -1i .* k .* r_dir ); % direct sound component
r_1_dir = norm( s - r2 );
H_r = 1 ./ r_1_dir .* exp( -1i .* k .* r_1_dir ); % reflected sound component
H_d = ita_diffraction_utd( w, s, r, f, c ); % diffracted sound component
% Ground reflection image
r_dir_1 = norm( s - r3 );
H_i_1 = 1 ./ r_dir_1 .* exp( -1i .* k .* r_dir_1 ); % direct sound component
r_1_dir_1 = norm( s - r23 );
H_r_1 = 1 ./ r_1_dir_1 .* exp( -1i .* k .* r_1_dir_1 ); % reflected sound component
[ H_d_1 , D_half ] = ita_diffraction_utd( w, s, r3, f, c ); % diffracted sound component
% zero and first order components
H_0 = H_d; % Total sound field ... diffracted component always present
if ~ita_diffraction_shadow_zone( w, s, r )
H_0 = H_0 + H_i;
end
if ita_diffraction_reflection_zone( w, s, r, false ) % opposite face reflection
H_0 = H_0 + H_r;
end
H_1 = H_d_1;
H_1_no_2 = H_d_1; % Missing 2d order reflection component
if ~ita_diffraction_shadow_zone( w, s, r3 )
H_1 = H_1 + H_i_1;
H_1_no_2 = H_1_no_2 + H_i_1;
end
if ita_diffraction_reflection_zone( w, s, r3, false ) % opposite face reflection
H_1 = H_1 + H_r_1;
end
H_1_half_dfo = zeros( numel( H_d_1 ), 1 );
if ita_diffraction_shadow_zone( w, s, r3 )
H_1_half_dfo = H_d_1; % Shadow zone yes diffraction only
else
H_1_half_dfo = H_i_1; % Direct sound yes
if ~ita_diffraction_reflection_zone( w, s, r3, false )
if ita_diffraction_utd_illumination_sign( w, s, r3, false ) == -1
sf = ita_diffraction_utd_illumination_fadeout( w, s, r3, false );
H_1_half_dfo = H_1_half_dfo + sf .* H_d_1; % closer to shadow zone diffraction yes
end
else
humbug = 1;
end
end
att_sum_2o = res_data_template;
att_sum_2o.freqData = ( H_0 + H_1 ) ./ H_i;
if n == 1
att_sum_all_2o = att_sum_2o;
else
att_sum_all_2o = ita_merge( att_sum_all_2o, att_sum_2o );
end
att_sum_2o_cut = res_data_template;
att_sum_2o_cut.freqData = ( H_0 + H_1_no_2 ) ./ H_i; % TF w/o 2. order refletion
if n == 1
att_sum_all_2o_cut = att_sum_2o_cut;
else
att_sum_all_2o_cut = ita_merge( att_sum_all_2o_cut, att_sum_2o_cut );
end
att_sum_dfo = res_data_template;
att_sum_dfo.freqData = ( H_1_half_dfo ) ./ H_i; % Diffraction fade-out
if n == 1
att_sum_all_dfo = att_sum_dfo;
else
att_sum_all_dfo = ita_merge( att_sum_all_dfo, att_sum_dfo );
end
att_sum_1ho = res_data_template;
att_sum_1ho.freqData = ( H_0 + H_1_half_dfo ) ./ H_i; % TF w/o 2. order reflection and fade-out
if n == 1
att_sum_all_1ho = att_sum_1ho;
else
att_sum_all_1ho = ita_merge( att_sum_all_1ho, att_sum_1ho );
end
end
%% Plot 1
psize_cm = [ 16*1 9*0.6 ];
% First order
fg1o = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );
h = plot( rad2deg( alpha_d ), flip( att_sum_all_2o.freqData_dB', 2 ) );
title( 'TF / DS: ground reflection with 2nd order reflection' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )
lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;
xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;
lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
h( i ).LineWidth = lw( i );
end
set( fg1o, 'papersize', psize_cm, 'units', 'centimeters' )
print( 'jst_diss_utd_ground_w_o2_refl', '-fillpage', '-dpdf' )
% Second order
fg2o = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );
h = plot( rad2deg( alpha_d ), flip( att_sum_all_2o_cut.freqData_dB', 2 ) );
title( 'TF / DS: ground reflection w/o 2nd order reflection' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )
lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;
xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;
lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
h( i ).LineWidth = lw( i );
end
set( fg2o, 'papersize', psize_cm, 'units', 'centimeters' )
print( 'jst_diss_utd_ground_wo_o2_refl', '-fillpage', '-dpdf' )
% Diffraction contribution w/o reflection component of diffraction method
fgfdo = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );
h = plot( rad2deg( alpha_d ), flip( att_sum_all_dfo.freqData_dB', 2 ) );
title( 'TF / DS: diffraction component SB-RB fade-out' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )
lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;
xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;
lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
h( i ).LineWidth = lw( i );
end
set( fgfdo, 'papersize', psize_cm, 'units', 'centimeters' )
print( 'jst_diss_utd_ground_no_2o_refl_fade_out', '-fillpage', '-dpdf' )
% 1st order w/o reflection component of diffraction method
fg1ho = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );
h = plot( rad2deg( alpha_d ), flip( att_sum_all_1ho.freqData_dB', 2 ) );
title( 'TF / DS: ground reflection with SB-RB fade-out' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )
lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;
xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;
lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
h( i ).LineWidth = lw( i );
end
set( fg1ho, 'papersize', psize_cm, 'units', 'centimeters' )
print( 'jst_diss_utd_ground_w_o1h_refl_fade_out', '-fillpage', '-dpdf' )
function [ distance ] = ita_propagation_effective_source_distance( propagation_path, anchor_idx )
%ITA_PROPAGATION_EFFECTIVE_SOURCE_DISTANCE Returns the backwards distance from given anchor index to
% the previous anchor type that provides a field value, e.g. a sensor or
% diffraction item. Integrates distance when one or multiple specular
% reflections are ahead.
% the previous anchor type that provides a field value, e.g. an emitter. Integrates distance when reflections or diffractions are ahead.
%
if nargin < 2
......@@ -28,7 +26,9 @@ for m = anchor_idx : -1 : 2
else
anchor = propagation_path.propagation_anchors( m - 1 );
end
if strcmpi( anchor.anchor_type, 'specular_reflection' )
if strcmpi( anchor.anchor_type, 'specular_reflection' ) || ...
strcmpi( anchor.anchor_type, 'inner_edge_diffraction' ) || ...
strcmpi( anchor.anchor_type, 'outer_edge_diffraction' )
if isa( propagation_path.propagation_anchors, 'cell' )
current_segment_vec = propagation_path.propagation_anchors{ m - 2 }.interaction_point - anchor.interaction_point;
else
......