Commit 3fc57de9 authored by jonasstienen's avatar jonasstienen
Browse files

New diffraction functions and tests for UTD method

parent 63d9a58a
function [ f ] = ita_diffraction_utd_illumination_fadeout( wedge, source_pos, receiver_pos, reflection_at_main_face )
%ita_diffraction_utd_illumination_fadeout Returns a scalar factor to fade out
%last diffraction component between shadow and half way towards reflection boundary
assert( ~ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ) )
assert( ~ita_diffraction_reflection_zone( wedge, source_pos, receiver_pos, reflection_at_main_face ) )
assert( ita_diffraction_utd_illumination_sign( wedge, source_pos, receiver_pos, reflection_at_main_face ) == -1 )
if reflection_at_main_face
vw = wedge.main_face_normal;
else
vw = wedge.opposite_face_normal;
end
vs = source_pos - wedge.location;
vr = receiver_pos - wedge.location;
distance_plane2source = dot( vs / norm( vs ) , vw / norm( vw ) );
distance_plane2receiver = dot( vr / norm( vr ) , vw / norm( vw ) );
s = 1 - abs( distance_plane2source / distance_plane2receiver ); % 0 ... 1
f = cos( s * pi / 2 ) * cos( s * pi / 2 ); % cosine-square fading from 1 to 0
end
function [ diffr_comp_sign ] = ita_diffraction_utd_illumination_sign( wedge, source_pos, receiver_pos, reflection_at_main_face )
%ITA_DIFFRACTION_UTD_ILLUMINATION_SIGN Returns -1 if source is closer to
%shadow boundary and 1 if source is closer to reflection boundary
assert( ~ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ) )
assert( ~ita_diffraction_reflection_zone( wedge, source_pos, receiver_pos, reflection_at_main_face ) )
if reflection_at_main_face
vw = wedge.main_face_normal;
else
vw = wedge.opposite_face_normal;
end
vs = source_pos - wedge.location;
distance_plane2source = dot( vs / norm( vs ) , vw / norm( vw ) );
if distance_plane2source > 0
diffr_comp_sign = 1; % closer to reflection boundary
else
diffr_comp_sign = -1; % closer to shadow boundary
end
end
%% 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 = [ 20, 50, 100, 200, 400 ]';
%f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 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 = 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, false ) % opposite face reflection
H = H + H_r;
end
H_1 = H_d_1;
if ~ita_diffraction_shadow_zone( w, s, r3 )
H_1 = H_1 + 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 = zeros( numel( H_d_1 ), 1 );
if ita_diffraction_shadow_zone( w, s, r3 )
H_1_half = H_d_1; % Shadow zone yes diffraction only
else
H_1_half = 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 = H_1_half + sf .* H_d_1; % closer to shadow zone diffraction yes
end
else
humbug = 1;
end
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' );
h = plot( rad2deg( alpha_d ), 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 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
%% Export
set( fg, 'papersize', [ 16*1.5 9*1.5 ], 'units', 'centimeters' )
print( 'ita_diffraction_plot_utd_total_field_stienen_diss_with_ground', '-fillpage', '-dpdf' )
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