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

WIP

parent d23dd2bc
function ita_diffraction_visualize_scene( w, source_pos, target_pos )
function ita_diffraction_visualize_scene( w, source_pos, target_pos, is_opengl_cs )
%% ita_diffraction_visualize_scene Visualises a source-receiver-wedge scene
%
% Example: ita_diffraction_visualize_scene( wedge, source_pos, target_pos )
%
if nargin < 4
is_opengl_cs = false;
end
apx = w.get_aperture_point2( source_pos, target_pos );
d1 = norm( source_pos - apx );
d2 = norm( target_pos - apx );
......@@ -110,7 +114,25 @@ Z3 = [ source_pos( 3 ); apx( 3 ); target_pos( 3 ) ];
%% Polygon 3D plot
fill3( X, Y, Z, 'red', X2, Y2, Z2, 'yellow', X3, Y3, Z3, 'green', X4, Y4, Z4, 'blue' )
M1 = [ X, Y, Z ];
M2 = [ X2, Y2, Z2 ];
M3 = [ X3, Y3, Z3 ];
M41 = [ X4(:,1), Y4(:,1), Z4(:,1) ];
M42 = [ X4(:,2), Y4(:,2), Z4(:,2) ];
if is_opengl_cs
M1 = ita_openGL2Matlab( M1 );
M2 = ita_openGL2Matlab( M2 );
M3 = ita_openGL2Matlab( M3 );
M41 = ita_openGL2Matlab( M41 );
M42 = ita_openGL2Matlab( M42 );
end
figure
fill3( M1(:,1), M1(:,2), M1(:,3), 'red', ...
M2(:,1), M2(:,2), M2(:,3), 'yellow', ...
M3(:,1), M3(:,2), M3(:,3), 'green', ...
M41(:,1), M41(:,2), M41(:,3), 'blue', ...
M42(:,1), M42(:,2), M42(:,3), 'blue' )
end
%% Config
n1 = [1 1 0];
n2 = [-1 1 0];
loc = [0 0 2];
src = 5/sqrt(1) * [-1 0 0];
rcv_start_pos = 3/sqrt(2) * [ -1 -1 0];
inf_wdg = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc);
apex_point = inf_wdg.get_aperture_point(src, rcv_start_pos);
apex_dir = inf_wdg.aperture_direction;
delta = 0.05;
n1 = [ 1 1 0 ] / sqrt( 2 );
n2 = [ -1 1 0 ] / sqrt( 2 );
loc = [ 0 0 0 ];
source_pos = 5 * [ -1 0 0] / sqrt( 2 );
receiver_start_pos = 5 * [ -1 -1 0 ] / sqrt( 2 );
w = itaInfiniteWedge( n1, n2, loc );
apex_point = w.get_aperture_point( source_pos, receiver_start_pos );
apex_dir = w.aperture_direction;
%delta = 0.05;
c = 344; % Speed of sound
freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]';
alpha_d = linspace( 0, inf_wdg.opening_angle, 100 );
freq = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
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
rcv_positions = ita_diffraction_align_points_around_aperture( inf_wdg, rcv_start_pos, alpha_d, apex_point, ref_face );
rcv_positions = rcv_positions(2 : end-1, :);
recevier_positions = norm( receiver_start_pos ) * [ cos( alpha_d - pi/4 ); sin( alpha_d - pi/4 ); zeros( 1, numel( alpha_d ) ) ]';
% wavenumber
k = 2 * pi * freq ./ c;
......@@ -27,40 +31,45 @@ res_data_template.freqVector = freq;
%% Calculations
N = size( rcv_positions, 1 );
N = size( recevier_positions, 1 );
for n = 1 : N
current_rcv_pos = rcv_positions( n, 1:3 );
receiver_pos = recevier_positions( n, 1:3 );
r_dir = norm( current_rcv_pos - src );
dir_field_comp = 1 ./ r_dir .* exp( -1i .* k .* r_dir );
r_dir = norm( receiver_pos - source_pos );
H_direct_field = 1 ./ r_dir .* exp( -1i .* k .* r_dir );
rcv_in_shadow_zone = ita_diffraction_shadow_zone( inf_wdg, src, current_rcv_pos );
rcv_in_shadow_zone = ita_diffraction_shadow_zone( w, source_pos, receiver_pos );
% UTD total wave field
H_diffracted_field = ita_diffraction_utd( w, source_pos, receiver_pos, freq, c );
if rcv_in_shadow_zone
H_total_field = H_diffracted_field;
else
H_total_field = H_diffracted_field + H_direct_field;
end
att_sum = res_data_template;
att_sum.freqData = ita_diffraction_utd( inf_wdg, src, current_rcv_pos, freq, c );
if ~rcv_in_shadow_zone
att_sum.freqData = att_sum.freqData + dir_field_comp;
att_sum.freqData = H_total_field ./ H_direct_field; % Insertion loss
% UTD total wave field with approximation by Tsingos et. al.
[ H_diffracted_field_approx ] = ita_diffraction_utd_approx( w, source_pos, receiver_pos, freq, c );
if rcv_in_shadow_zone
H_total_field_approx = H_diffracted_field_approx;
else
H_total_field_approx = H_diffracted_field_approx + H_direct_field;
end
att_sum.freqData = att_sum.freqData ./ dir_field_comp;
% UTD combined with Approximation by Tsingos et. al.
att_sum_approx = res_data_template;
[ H ] = ita_diffraction_utd_approx( inf_wdg, src, current_rcv_pos, freq, c );
att_sum_approx.freqData = H;
if ~rcv_in_shadow_zone
att_sum_approx.freqData = att_sum_approx.freqData + dir_field_comp;
end
att_sum_approx.freqData = att_sum_approx.freqData ./ dir_field_comp;
att_sum_approx.freqData = H_total_field_approx ./ H_direct_field;
% Error caused by approximation
att_sum_error = res_data_template;
att_sum_error.freqData = att_sum_approx.freqData ./ att_sum.freqData;
att_sum_error = att_sum_approx / att_sum;
%att_sum_error.freqData = att_sum_approx.freqData ./ att_sum.freqData;
att_sum_diffr = res_data_template;
att_sum_diffr.freqData = ita_diffraction_utd( inf_wdg, src, current_rcv_pos, freq, c );
att_sum_diffr.freqData = att_sum_diffr.freqData ./ dir_field_comp;
att_sum_diffr.freqData = ita_diffraction_utd( w, source_pos, receiver_pos, freq, c );
att_sum_diffr.freqData = att_sum_diffr.freqData ./ H_direct_field;
if n == 1
att_sum_all = att_sum;
......@@ -77,7 +86,7 @@ for n = 1 : N
end
%% Tsingos paper plot
angles = rad2deg( alpha_d(2 : end-1) );
angles = rad2deg( linspace( alpha_d( end ), alpha_d( 1 ), numel( alpha_d ) ) ); % Reverse angle on axis
figure( 'units', 'normalized', 'outerposition', [0 0 1 1] );
subplot( 2, 2, 1 );
plot( angles, att_sum_all.freqData_dB' )
......@@ -85,7 +94,6 @@ title( 'Total wave field with UTD' )
legend( num2str( freq ), 'Location', 'southwest' )
xlabel( 'alpha_d in degree (shadow boundary at 225deg)' )
ylabel( 'dB SPL' )
xlim( [alpha_d(1), inf_wdg.opening_angle_deg] );
ylim( [-35, 10] );
grid on;
......@@ -95,7 +103,6 @@ title( 'Tsingos et al.: Total wave field plot with approximated UTD (Figure 6b)'
legend( num2str( freq ), 'Location', 'southwest' )
xlabel( 'alpha_d in degree (shadow boundary at 225deg)' )
ylabel( 'dB SPL' )
xlim( [angles(1), inf_wdg.opening_angle_deg] );
ylim( [-35, 10] );
grid on;
......@@ -105,7 +112,6 @@ title( 'Tsingos et al.: Error by approximation (Figure 6c)' )
legend( num2str( freq ), 'Location', 'southwest' )
xlabel( 'alpha_d in degree (shadow boundary at 225deg)' )
ylabel( 'dB SPL' )
xlim( [angles(1), inf_wdg.opening_angle_deg] );
ylim( [-35, 10] );
grid on;
......@@ -115,5 +121,4 @@ title( 'Diffracted field by UTD' )
legend( [num2str( freq ), repmat(' Hz', numel(freq),1)], 'Location', 'southeast' )
xlabel( 'theta_R []' )
ylabel( 'dB SPL' )
xlim( [angles(1), inf_wdg.opening_angle_deg] );
grid on;
%% Test transition at shadown zone of UTD diffraction
w = itaInfiniteWedge( [ 1 0 0 ], [ 0 0 -1 ], [ 0 0 0 ] ); % OpenGL coordinates
delta = 0.001;
s_shadow = [ ( -3 - delta ) 0 -3 ];
s_illuminated = [ ( -3 + delta ) 0 -3 ];
r = [ 3 0 3 ];
assert( ita_diffraction_shadow_zone( w, s_shadow, r ) )
assert( ~ita_diffraction_shadow_zone( w, s_illuminated, r ) )
utd_tf = itaAudio( 1 );
utd_tf.fftDegree = 11;
f = utd_tf.freqVector( 2:end );
c = 341;
[ H1, D1, A1 ] = ita_diffraction_utd( w, s_shadow, r, f, c );
utd_tf.freqData( :, 1 ) = [ 0 H1 ];
[ H2, D2, A2 ] = ita_diffraction_utd( w, s_illuminated, r, f, c );
utd_tf.freqData( :, 2 ) = [ 0 H2 ];
d = norm( r - s_illuminated );
c = 343;
k = 2 * pi * f / c;
H_direct = 1 ./ d .* exp( -1i .* k .* d );
utd_tf.freqData( :, 3 ) = [ 0 H1 ./ H_direct ];
utd_tf.freqData( :, 4 ) = [ 0 ( H2 + H_direct ) ./ H_direct ];
utd_tf.channelNames = { 'Diffracted field (shadow)', 'Diffracted field (illuminated)', 'Insertion loss (shadowed)', 'Insertion loss (illuminated)' };
utd_tf.pf
Supports Markdown
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