From 7b2d10d8c7019fda92a39b35e4c5746757f8a1cc Mon Sep 17 00:00:00 2001 From: "Dipl.-Ing. Jonas Stienen" Date: Fri, 16 Nov 2018 10:36:07 +0100 Subject: [PATCH] Adding diffraction simulation scripts and wedge classes --- .../@itaFiniteWedge/itaFiniteWedge.m | 66 ++++ .../longer_branch_on_aperture.m | 16 + .../@itaFiniteWedge/point_on_aperture.m | 26 ++ .../@itaFiniteWedge/point_outside_wedge.m | 21 + .../get_angle_from_point_to_aperture.m | 45 +++ .../get_angle_from_point_to_wedge_face.m | 41 ++ .../@itaInfiniteWedge/get_aperture_point.m | 72 ++++ .../@itaInfiniteWedge/itaInfiniteWedge.m | 202 ++++++++++ .../point_facing_main_side.m | 44 +++ .../@itaInfiniteWedge/point_on_aperture.m | 25 ++ .../@itaInfiniteWedge/point_outside_wedge.m | 21 + .../@itaPlaneWedge/itaPlaneWedge.m | 11 + .../itaSemiInfinitePlane.m | 12 + .../diffraction/Svensson_BTM_test.m | 85 ++++ .../diffraction/UTD_approx_test.m | 239 ++++++++++++ .../diffraction/biot_tolstoy_test.m | 28 ++ .../diffraction/btm_total_field_plot.m | 55 +++ .../diffraction/diffraction_TF_IR.m | 227 +++++++++++ .../ita_align_points_around_aperture.m | 89 +++++ .../ita_diffraction_biot_tolstoy_jst.m | 19 + .../ita_diffraction_btm_infinite_wedge.m | 108 ++++++ .../ita_diffraction_btm_total_field_plots.m | 232 +++++++++++ .../diffraction/ita_diffraction_btms.m | 367 ++++++++++++++++++ .../diffraction/ita_diffraction_btms_approx.m | 149 +++++++ .../ita_diffraction_kawai_approx_fresnel.m | 22 ++ ...ta_diffraction_kawai_approx_fresnel_plot.m | 17 + .../diffraction/ita_diffraction_maekawa.m | 67 ++++ .../ita_diffraction_maekawa_approx.m | 76 ++++ ...ta_diffraction_maekawa_total_field_plots.m | 205 ++++++++++ .../diffraction/ita_diffraction_shadow_zone.m | 85 ++++ .../ita_diffraction_svensson_plots.m | 96 +++++ .../ita_diffraction_tsingos_utd_approx_plot.m | 64 +++ ...iffraction_tsingos_utd_total_field_plots.m | 91 +++++ .../diffraction/ita_diffraction_utd.m | 203 ++++++++++ .../ita_diffraction_utd_approximated.m | 57 +++ ...fraction_utd_cot_singularity_approx_test.m | 147 +++++++ .../diffraction/ita_diffraction_utd_test.m | 37 ++ .../diffraction/ita_reflection_zone.m | 153 ++++++++ .../diffraction/ita_rotation_rodrigues.m | 32 ++ .../VirtualAcoustics/diffraction/maekawa.m | 88 +++++ 40 files changed, 3640 insertions(+) create mode 100644 applications/VirtualAcoustics/diffraction/@itaFiniteWedge/itaFiniteWedge.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaFiniteWedge/longer_branch_on_aperture.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_on_aperture.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_outside_wedge.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_aperture.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_wedge_face.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_aperture_point.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/itaInfiniteWedge.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_facing_main_side.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_on_aperture.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_outside_wedge.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaPlaneWedge/itaPlaneWedge.m create mode 100644 applications/VirtualAcoustics/diffraction/@itaSemiInfinitePlane/itaSemiInfinitePlane.m create mode 100644 applications/VirtualAcoustics/diffraction/Svensson_BTM_test.m create mode 100644 applications/VirtualAcoustics/diffraction/UTD_approx_test.m create mode 100644 applications/VirtualAcoustics/diffraction/biot_tolstoy_test.m create mode 100644 applications/VirtualAcoustics/diffraction/btm_total_field_plot.m create mode 100644 applications/VirtualAcoustics/diffraction/diffraction_TF_IR.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_align_points_around_aperture.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_biot_tolstoy_jst.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_btm_infinite_wedge.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_btm_total_field_plots.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_btms.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_btms_approx.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel_plot.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_approx.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_total_field_plots.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_shadow_zone.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_svensson_plots.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_approx_plot.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_total_field_plots.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_utd.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_utd_approximated.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_utd_cot_singularity_approx_test.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_diffraction_utd_test.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_reflection_zone.m create mode 100644 applications/VirtualAcoustics/diffraction/ita_rotation_rodrigues.m create mode 100644 applications/VirtualAcoustics/diffraction/maekawa.m diff --git a/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/itaFiniteWedge.m b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/itaFiniteWedge.m new file mode 100644 index 0000000..7377276 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/itaFiniteWedge.m @@ -0,0 +1,66 @@ +classdef itaFiniteWedge < itaInfiniteWedge + + properties (Access = protected) + al % aperture length + sp % start point of aperture + ep % end point of aperture + end + + properties (Dependent) + length + aperture_start_point + aperture_end_point + end + + methods + function obj = itaFiniteWedge( main_face_normal, opposite_face_normal, location, length, wedge_type ) + % Creates an finite wedge with a location in 3-dim space + % Starting point of wedge aperture is the wedge location, end + % point is defined by length and direction of aperture + % Length must be a scalar value greater zero + if nargin < 5 + wedge_type = 'wedge'; + end + obj@itaInfiniteWedge( main_face_normal, opposite_face_normal, location, wedge_type ); + if numel( length ) > 1 || length <= 0 + error 'Length must be a scalar value greater zero' + end + obj.al = length; + obj.sp = obj.l; + + n_scaled = cross( obj.main_face_normal, obj.opposite_face_normal ); + if ~norm( n_scaled ) + warning 'aperture enp point could not be determined since aperture direction is not defined. Please set aperture direction manually.' + else + obj.ep = obj.sp + length * obj.aperture_direction; + end + end + + + function obj = set.aperture_end_point( obj, length_of_aperture ) + % Sets aperture direction manually (in case of linear + % dependent normals) + if norm( cross( obj.n1, obj.n2 ) ) + error 'Aperture is already fixed and cannot be modified' + end + if length_of_aperture <= 0 + error 'Aperture length must be a valid scalar > 0' + end + obj.ep = length_of_aperture; + end + + function l = get.length( obj ) + l = obj.al; + end + + function sp = get.aperture_start_point( obj ) + % 3-dim starting point of finite aperture + sp = obj.sp; + end + + function ep = get.aperture_end_point( obj ) + % 3-dim end point of finite aperture + ep = obj.ep; + end + end +end diff --git a/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/longer_branch_on_aperture.m b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/longer_branch_on_aperture.m new file mode 100644 index 0000000..b0acc55 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/longer_branch_on_aperture.m @@ -0,0 +1,16 @@ +function direction = longer_branch_on_aperture( obj, apex_point ) +%Returns the direction vector of the longer one of two sections on the aperture devided +%by the aperture point +% direction: vector(normalized) pointing in direction of longer section from given aperture +% point + +d1 = norm( obj.aperture_start_point - apex_point ); +d2 = norm( obj.aperture_end_point - apex_point ); +if d2 >= d1 + direction = obj.aperture_direction; +else + direction = -obj.aperture_direction; +end + +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_on_aperture.m b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_on_aperture.m new file mode 100644 index 0000000..37858bc --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_on_aperture.m @@ -0,0 +1,26 @@ +function b = point_on_aperture( obj, point ) +% Returns true if point is on aperture of the wedge and between start and +% end point of the aperture +dim = size( point ); +if dim(1) ~= 3 && dim(2) ~= 3 + error( 'Point must be of dimension 3' ); +end + +b = zeros( numel(point) / 3, 1 ); +norms = sqrt( sum( (point - obj.location).^2, 2 ) ); +condition1 = norms < obj.set_get_geo_eps; +if any( condition1 ) + b = b | condition1; +end +dir1 = ( point - obj.location ) ./ norms; +dir2 = ( point - obj.aperture_end_point ) ./ sqrt( sum( (point - obj.aperture_end_point).^2, 2 ) ); +dir1_norms = sqrt( sum( (dir1 - obj.aperture_direction).^2, 2 ) ); +dir2_norms = sqrt( sum( (dir2 + obj.aperture_direction).^2, 2 ) ); +condition2 = (dir1_norms < obj.set_get_geo_eps | dir2_norms < obj.set_get_geo_eps) & norms < obj.length; +if any( condition2 ) + b = b | condition2; +end +b = all(b); + +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_outside_wedge.m b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_outside_wedge.m new file mode 100644 index 0000000..7bf5fac --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaFiniteWedge/point_outside_wedge.m @@ -0,0 +1,21 @@ +function b = point_outside_wedge( obj, point ) + % Returns true if the point is outside the solid structure of + % the finite wedge + + dim = size( point ); + if dim(2) ~= 3 + if dim(1) ~= 3 + error( 'Point(s) must be of dimension 3') + end + point = point'; + dim = size( point ); + end + + dist_from_main_face = sum( (point - obj.location) .* obj.main_face_normal, 2 ); + dist_from_opposite_face = sum( (point - obj.location) .* obj.opposite_face_normal, 2 ); + + b = zeros( dim(1), 1 ); + mask = ( dist_from_main_face < -obj.set_get_geo_eps ) & ( dist_from_opposite_face < -obj.set_get_geo_eps ); + b(~mask) = true; + +end diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_aperture.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_aperture.m new file mode 100644 index 0000000..57994f1 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_aperture.m @@ -0,0 +1,45 @@ +function alpha = get_angle_from_point_to_aperture( obj, field_point, point_on_aperture ) +%Returns angle (radiant) between the ray from field point to aperture point and the +%aperture of the wedge. +% output angle alpha: 0 <= alpha <= pi/2 + +%% assertions +dim_fp = size(field_point); +dim_pa = size(point_on_aperture); +if dim_fp(2) ~= 3 + if dim_fp(1) ~= 3 + error( 'Field point must be a row vector of dimension 3' ); + end + field_point = field_point'; + dim_fp = size(field_point); +end +if dim_pa(2) ~= 3 + if dim_pa(1) ~= 3 + error( 'Point on Aperture must be a row vector of dimension 3.' ); + end + point_on_aperture = point_on_aperture'; + dim_pa = size(point_on_aperture); +end +if dim_fp(1) ~= 1 && dim_pa(1) ~= 1 + if dim_fp(1) ~= dim_pa(1) + error( 'Use same number of field points and points on aperture' ); + end +end +if any( ~obj.point_outside_wedge( field_point ) ) + error( 'Field point must be outside wedge' ); +end +if any( ~obj.point_on_aperture( point_on_aperture ) ) + error( 'No point on aperture found' ) +end + +%% begin +norms = sqrt( sum( (point_on_aperture - field_point).^2, 2 ) ); +dir_vec = ( point_on_aperture - field_point ) ./ norms; +alpha = acos( sum( dir_vec .* obj.aperture_direction, 2 ) ); +mask = alpha > pi/2; +if any(mask) + alpha(mask) = pi - alpha(mask); +end + +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_wedge_face.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_wedge_face.m new file mode 100644 index 0000000..fef79dd --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_angle_from_point_to_wedge_face.m @@ -0,0 +1,41 @@ +function theta = get_angle_from_point_to_wedge_face( obj, point_position, reference_face ) +%GET_ANGLE_FROM_SOURCE_TO_WEDGE_FACE Returns angle (radiant) between +% given point and a wedge face. +% input: point_pos arbitrary field point outside the wedge +% reference_face reference wedge face which the angle will +% be denoted from +% output: theta azimuth angle (radiant) in cylinder coordinates with +% the aperture as z axis + +%% Assertions +dim = size( point_position ); +if dim(2) ~= 3 + if dim(1) ~= 3 + error( 'Field points have to be of dimension 3' ); + end + point_position = point_position'; + dim = size( point_position ); +end +% if any( (reference_face ~= 'main face') & (reference_face ~= 'opposite face') ) +% error( 'Reference face must be either main face or opposite face. You may also use itaInfiniteWedge.get_source_facing_side().' ) +% end +if any( ~obj.point_outside_wedge( point_position ) ) + error('Point(s) must be outside the wedge'); +end + +%% Begin +e_z = repmat( obj.aperture_direction, dim(1), 1 ); +e_y = zeros( dim(1), 3 ); +e_y( reference_face, : ) = repmat( obj.main_face_normal, sum( reference_face ), 1 ); +e_y( ~reference_face, : ) = repmat( obj.opposite_face_normal, sum( ~reference_face ), 1 ); +e_x = cross( e_z, e_y, 2 ); + +% Calculate angle between incedent ray from source to aperture point and source facing wedge +% side +x_i = dot( point_position - obj.location, e_x, 2 ); % coordinates in new coordinate system +y_i = dot( point_position - obj.location, e_y, 2 ); +theta = atan2( y_i, x_i ); +theta( theta < 0 ) = theta( theta < 0 ) + 2*pi; + +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_aperture_point.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_aperture_point.m new file mode 100644 index 0000000..db30af5 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/get_aperture_point.m @@ -0,0 +1,72 @@ +function ap = get_aperture_point( obj, source_pos, receiver_pos ) + % Returns aperture point on wedge (closest point on wedge + % between source and receiver) + + %% Verification + dim_src = size( source_pos ); + dim_rcv = size( receiver_pos ); + if dim_src(2) ~= 3 + if dim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3') + end + source_pos = source_pos'; + dim_src = size( source_pos ); + end + if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Receiver point(s) must be of dimension 3') + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); + end + if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) + end + if dim_src(1) > dim_rcv(1) + dim_n = dim_src(1); + S = source_pos; + R = repmat( receiver_pos, dim_n, 1 ); + elseif dim_src(1) < dim_rcv(1) + dim_n = dim_rcv(1); + S = repmat( source_pos, dim_n, 1 ); + R = receiver_pos; + else + dim_n = dim_src(1); + S = source_pos; + R = receiver_pos; + end + + %% Variables + L = repmat( obj.location, dim_n, 1 ); + Apex_Dir = repmat( obj.aperture_direction, dim_n, 1 ); + + %% Calculations + SR = R - S; + norm_of_SR = Norm( SR ); + mask = norm_of_SR ~= 0; + SR_dir = SR(mask, :) ./ norm_of_SR(mask); + + % initialize result vector + ap = zeros( dim_n, 3 ); + + % Auxilary plane spanned by SR and aux_plane_dir + aux_plane_dir = cross( SR_dir, Apex_Dir(mask, :), 2 ) ./ Norm( cross( SR_dir, Apex_Dir(mask, :), 2 ) ); + aux_plane_normal = cross( SR_dir, aux_plane_dir, 2 ) ./ Norm( cross( SR_dir, aux_plane_dir, 2 ) ); + + % Distance of intersection of auxiliary plane and aperture direction + % from aperture location + % aux plane: dot( (x - source_point), aux_plane_normal) = 0 + % aperture line: x = location + dist * aperture_direction + dist = dot( S(mask, :) - L(mask, :), aux_plane_normal, 2 ) ./ dot( Apex_Dir(mask, :), aux_plane_normal, 2 ); + ap(mask, :) = L(mask, :) + dist .* Apex_Dir(mask, :); + + % In case receiver and source have same position + if any( norm_of_SR == 0 ) + dist = dot( R(~mask, :) - L(~mask, :), Apex_Dir(~mask, :), 2 ); + ap(~mask, :) = L(~mask, :) + dist * Apex_Dir(~mask, :); + end +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/itaInfiniteWedge.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/itaInfiniteWedge.m new file mode 100644 index 0000000..d6c5993 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/itaInfiniteWedge.m @@ -0,0 +1,202 @@ +classdef itaInfiniteWedge + + properties (Access = protected) + n1 % 3-dim normal vector of main face (internal) + n2 % 3-dim normal vector of opposite face (internal) + ad % 3-dim aperture direction vector (internal) + l % Internal location variable + wt % type of wedge (internal) + bc_hard % Internal boundary condition (hard = true) + end + + properties (Dependent) + main_face_normal % 3-dim normal vector of main face (normalized) + opposite_face_normal % 3-dim normal vector of opposite face (normalized) + aperture_direction % 3-dim normal vector of aperture direction (normalized) + location % Location of wedge (somewhere along aperture) + opening_angle % Angle from main to opposite face in propagation medium / air (radiants) + wedge_angle % Angle from main to opposite face in solid medium of wedge (radiants) + wedge_type % 'wedge' for opening angles > pi or 'corner' for opening angles < pi + boundary_condition % boundary condition of the wedge faces (hard or soft) + end + + methods + function obj = itaInfiniteWedge( main_face_normal, opposite_face_normal, location, wedge_type ) + % Create a wedge by a main face normal and an opposite face + % normal + % main_face_normal: Main face normal (3-dim) + % opposite_face_normal: Opposite face normal (3-dim) + % location: Point on aperture which defines + % location of the wedge in 3_dim sapce + % wedge_type: use 'wedge' for opening angles > pi (default) and + % 'corner' for opening angles < pi + % Note: 3-dim direction vectors will be normalized automatically + % + if nargin < 4 + wedge_type = 'wedge'; + end + if ~isequal( wedge_type, 'wedge' ) && ~isequal( wedge_type, 'corner' ) + error( 'Invalid wedge type. Use either wedge or corner' ) + end + if numel( main_face_normal ) ~= 3 + error 'Main face normal has to be a 3-dim vector' + end + if numel( opposite_face_normal ) ~= 3 + error 'Opposite face normal has to be a 3-dim vector' + end + if numel(location) ~= 3 + error( 'Location must be of dimension 3') + end + + obj.n1 = main_face_normal; + obj.n2 = opposite_face_normal; + obj.l = location; + obj.wt = wedge_type; + obj.bc_hard = true; + + if ~obj.validate_normals + warning 'Normalized face normals' + obj.n1 = main_face_normal ./ norm( main_face_normal ); + obj.n2 = opposite_face_normal ./ norm( opposite_face_normal ); + end + + n_scaled = cross( obj.main_face_normal, obj.opposite_face_normal ); + if ~norm( n_scaled ) + warning 'Normals are linear dependent and aperture direction could not be determined. Please set aperture direction manually.' + else + obj.ad = n_scaled ./ norm( n_scaled ); + end + end + + function n = get.main_face_normal( obj ) + n = obj.n1; + end + + function n = get.opposite_face_normal( obj ) + n = obj.n2; + end + + function n = get.aperture_direction( obj ) + % Returns normalized direction of aperture. Vectors main face normal, opposite face normal and aperture direction + % form a clockwise system. + if isempty( obj.ad ) + error 'Invalid wedge, aperture direction not set and face normals are linear dependent' + end + n = obj.ad; + end + + function obj = set.aperture_direction( obj, aperture_direction ) + % Sets aperture direction manually (in case of linear + % dependent normals) + if norm( cross( obj.n1, obj.n2 ) ) + error 'Aperture of linear independent normals is fixed can not be modified' + end + if ~norm( aperture_direction ) + error 'Aperture vector must be a valid direction' + end + if norm( aperture_direction ) ~= 1 + warning ' Normalizing aperture direction' + aperture_direction = aperture_direction / norm( aperture_direction ); + end + if ~( dot( aperture_direction, obj.n1 ) == 0 && dot( aperture_direction, obj.n2 ) == 0 ) + error 'Invalid aperture direction, vector must be perpendicular to face normals' + end + obj.ad = aperture_direction; + end + + function beta = get.wedge_angle( obj ) + % Returns angle from main to opposite face through solid medium + % of the wedge (radiant) + if isequal( obj.wt, 'wedge' ) + s = 1; + elseif isequal( obj.wt, 'corner' ) + s = -1; + end + beta = pi - s * acos(dot(obj.main_face_normal, obj.opposite_face_normal)); + end + + function beta_deg = wedge_angle_deg( obj ) + % Get the wedge angle angle in degree + beta_deg = rad2deg( obj.wedge_angle ); + end + + function theta = get.opening_angle( obj ) + % Returns angle from main face to opposite face through propagation medium / + % air (radiant) + theta = 2 * pi - obj.wedge_angle; + end + + function theta_deg = opening_angle_deg( obj ) + % Get the wedge opening angle in degree + theta_deg = rad2deg( obj.opening_angle ); + end + + function l = get.location( obj ) + l = obj.l; + end + + function b = validate_normals( obj ) + % Returns true, if the normals of the faces are both normalized + b = true; + if norm( obj.main_face_normal ) ~= 1 || norm( obj.opposite_face_normal ) ~= 1 + b = false; + end + end + + function wt = get.wedge_type( obj ) + wt = obj.wt; + end + + function bc = get.boundary_condition( obj ) + if obj.bc_hard + bc = 'hard'; + else + bc = 'soft'; + end + end + + function obj = set.boundary_condition( obj, bc ) + if ischar(bc) + if strcmpi('hard', bc) + obj.bc_hard = true; + elseif strcmpi('soft', bc) + obj.bc_hard = false; + else + error('boundary condition must be "hard" or "soft"!'); + end + else + error('boundary condtion must be of type character!'); + end + end + + function obj = set.bc_hard( obj, b ) + obj.bc_hard = b; + end + + function bc = is_boundary_condition_hard( obj ) + bc = obj.bc_hard; + end + + function bc = is_boundary_condition_soft( obj ) + bc = ~obj.bc_hard; + end + end + + + methods (Static) + function current_eps = set_get_geo_eps( new_eps ) + % Controls and returns the geometrical calculation precision value for + % deciding e.g. if a point is inside or outside a wedge + % (defaults to Matlab eps, but should be set for instance to + % millimeter (1e-3) or micrometer (1e-6). + persistent geo_eps; + if nargin > 0 + geo_eps = new_eps; + end + if isempty( geo_eps ) + geo_eps = eps; % Default eps from Matlab double precision + end + current_eps = geo_eps; + end + end +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_facing_main_side.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_facing_main_side.m new file mode 100644 index 0000000..6d32be7 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_facing_main_side.m @@ -0,0 +1,44 @@ +function res = point_facing_main_side( obj, point ) +% Returns string including the face normal of source +% facing wedge side +% output: 'n1' -> source is facing main face of wedge +% 'n2' -> source is facing opposite face of wedge + +%% Assertions +dim = size( point ); +if dim(2) ~= 3 + if dim(1) ~= 3 + error( 'Position vectors must be of dimension 3!' ); + end + point = point'; + dim = size( point ); +end +if any( ~obj.point_outside_wedge( point ) ) + error( 'Source position(s) must be outside the wedge!' ); +end + +%% begin +% Coordinate transformation with aperture location as origin +e_z = repmat( obj.aperture_direction, dim(1), 1 ); +e_y1 = repmat( obj.main_face_normal, dim(1), 1 ); +e_x1 = cross( e_y1, e_z, 2 ); % direction of main face of the wedge +e_y2 = repmat( obj.opposite_face_normal, dim(1), 1 ); +e_x2 = cross( e_z, e_y2, 2 ); + +% Calculate angle between incedent ray from source to aperture point and source facing wedge +% side +x_i1 = dot( point - obj.location, e_x1, 2 ); % coordinates in new coordinate system +y_i1 = dot( point - obj.location, e_y1, 2 ); +temp1 = atan2( y_i1, x_i1 ); +temp1( temp1 < 0 ) = temp1( temp1 < 0 ) + 2*pi; + +x_i2 = dot( point - obj.location, e_x2, 2 ); % coordinates in new coordinate system +y_i2 = dot( point - obj.location, e_y2, 2 ); +temp2 = atan2( y_i2, x_i2 ); +temp2( temp2 < 0 ) = temp2( temp2 < 0 ) + 2*pi; + +res = true( dim(1), 1 ); +res( temp2 < temp1 ) = false; + +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_on_aperture.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_on_aperture.m new file mode 100644 index 0000000..beb95ba --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_on_aperture.m @@ -0,0 +1,25 @@ +function b = point_on_aperture( obj, point ) +% Returns true if point is on aperture of the wedge + +dim = size( point ); +if dim(1) ~= 3 && dim(2) ~= 3 + error( 'Point must be of dimension 3' ); +end + +b = zeros( numel(point) / 3, 1 ); +norms = sqrt( sum( (point - obj.location).^2, 2 ) ); +condition1 = norms < obj.set_get_geo_eps; +if any( condition1 ) + b = b | condition1; +end +dir1 = ( point - obj.location ) ./ norms; +dir2 = ( point - (obj.location + 10 * obj.aperture_direction) ) ./ sqrt( sum( (point - (obj.location + 10 * obj.aperture_direction)).^2, 2 ) ); +dir1_norms = sqrt( sum( (abs(dir1) - abs(obj.aperture_direction)).^2, 2) ); +dir2_norms = sqrt( sum( (abs(dir2) - abs(obj.aperture_direction)).^2, 2) ); +condition2 = (dir1_norms < obj.set_get_geo_eps) | (dir2_norms < obj.set_get_geo_eps); +if any( condition2 ) + b = b | condition2; +end +b = all(b); +end + diff --git a/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_outside_wedge.m b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_outside_wedge.m new file mode 100644 index 0000000..f2d49a8 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaInfiniteWedge/point_outside_wedge.m @@ -0,0 +1,21 @@ +function b = point_outside_wedge( obj, point ) + % Returns true if point is outside the solid structure of + % the infinite wedge + + dim = size( point ); + if dim(2) ~= 3 + if dim(1) ~= 3 + error( 'Point(s) must be of dimension 3') + end + point = point'; + dim = size( point ); + end + + dist_from_main_face = sum( (point - obj.location) .* obj.main_face_normal, 2 ); + dist_from_opposite_face = sum( (point - obj.location) .* obj.opposite_face_normal, 2 ); + + b = false( dim(1), 1 ); + mask = ( dist_from_main_face < -obj.set_get_geo_eps ) & ( dist_from_opposite_face < -obj.set_get_geo_eps ); + b(~mask) = true; + +end diff --git a/applications/VirtualAcoustics/diffraction/@itaPlaneWedge/itaPlaneWedge.m b/applications/VirtualAcoustics/diffraction/@itaPlaneWedge/itaPlaneWedge.m new file mode 100644 index 0000000..7ae31aa --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaPlaneWedge/itaPlaneWedge.m @@ -0,0 +1,11 @@ +classdef itaPlaneWedge < itaInfiniteWedge + %ITAPLANEWEDGE special case of itaInfiniteWedge + % Case of wedge with opening angle of pi + + methods + function obj = itaPlaneWedge( plane_normal, location, aperture_direction ) + obj@itaInfiniteWedge( plane_normal, plane_normal, location ); + obj.aperture_direction = aperture_direction; + end + end +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/@itaSemiInfinitePlane/itaSemiInfinitePlane.m b/applications/VirtualAcoustics/diffraction/@itaSemiInfinitePlane/itaSemiInfinitePlane.m new file mode 100644 index 0000000..ca9c3d2 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/@itaSemiInfinitePlane/itaSemiInfinitePlane.m @@ -0,0 +1,12 @@ +classdef itaSemiInfinitePlane < itaInfiniteWedge + %ITASEMIINFINITEPLANE Special case of itaInfiniteWedge + % Case of a wedge with opening angle of 2pi. + + methods + function obj = itaSemiInfinitePlane( main_face_normal, location, aperture_direction ) + obj@itaInfiniteWedge( main_face_normal, -main_face_normal, location ); + obj.aperture_direction = aperture_direction; + end + end +end + diff --git a/applications/VirtualAcoustics/diffraction/Svensson_BTM_test.m b/applications/VirtualAcoustics/diffraction/Svensson_BTM_test.m new file mode 100644 index 0000000..3583d5e --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/Svensson_BTM_test.m @@ -0,0 +1,85 @@ +%% Config +n1 = [1, 1, 0]; +n2 = [-1, 1, 0]; +loc = [0, 0, -3]; +src = [-4, -2, 0]; +rcv = [5, -2, 0]; +p = [0, 0, 2]; +len = 18; +sampling_rate = 44100; +filter_length = 3534; + +infw = itaInfiniteWedge(n1, n2, loc); +finw = itaFiniteWedge(n1, n2, loc, len); + +ref = infw.point_facing_main_side( src ); % reference wedge + +%% EDB variables +fs = sampling_rate; +closwedang = finw.wedge_angle; +rs = sqrt( src(1)^2 + src(2)^2 ); +thetas = infw.get_angle_from_point_to_wedge_face( src, ref ); +zs = 0; +rr = sqrt( rcv(1)^2 + rcv(2)^2 ); +thetar = infw.get_angle_from_point_to_wedge_face( rcv, ref ); +zr = rcv(3) - src(3); +zw = [ loc(3) - src(3), len - src(3) ]; +Method = 'New'; + +%% Filter + +[ir,initdelay,singularterm] = EDB2wedge1st_int( fs, closwedang, rs, thetas, zs, rr, thetar, zr, zw, Method ); +res1 = ir(ir ~= 0); +res2 = ita_diffraction_btm_finite_wedge( finw, src, rcv, fs, filter_length ); + + +% figure(); +% plot(res1); +% title( 'EDB2 toolbox' ); +% xlim( [0, 1200] ); +% grid on; + +figure(); +plot( [res1, res2.timeData] ); +title( 'BTM_{finw}' ); +xlim( [0, 1200] ); +legend( 'EDB Toolbox', 'own Impl' ); +grid on; + +res3 = ita_diffraction_btm_infinite_wedge( infw, src, rcv, fs, filter_length ); +figure(); +plot(res3.timeData); +title( 'BTM_{infw}' ); +xlim( [0, 1200] ); +grid on; + +ratio1 = res1 ./ res2.timeData; +ratio2 = res2.timeData ./ res3.timeData; +figure(); +plot([ratio1, ratio2]); +title( 'Ratio between EDB2 toolbox and BTM_{finw}' ); +xlim( [0, 1200] ); +legend( 'svensson / btm_{finw}', 'btm_{finw} / btm_{infw}' ); +grid on; + +% Res = [res1, res2.timeData, ratio]'; +% +% figure(); +% yyaxis left +% plot( Res(1 : 2, :) ); +% ylabel( 'amplitude' ); +% legend( 'res1', 'res2' ); +% yyaxis right +% plot( Res(3, :) ); +% ylabel( 'ratio: res1/res2' ); + +% A = itaAudio; +% A.timeData = res1(1:1024); +% +% B = res2; +% +% A.pt; +% B.pt; +% A.pf; +% B.pf; + diff --git a/applications/VirtualAcoustics/diffraction/UTD_approx_test.m b/applications/VirtualAcoustics/diffraction/UTD_approx_test.m new file mode 100644 index 0000000..7f9295e --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/UTD_approx_test.m @@ -0,0 +1,239 @@ +% Config +n1 = [-1 1 0]; +n2 = [1 1 0]; +loc = [0 0 2]; +src = [-1 0 0]; +rcv = [1 1 0]; + +w = itaInfiniteWedge(n1, n2, loc); +apex_point = w.get_aperture_point(src, rcv); +apex_dir = w.aperture_direction; + +% freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]; +freq = 2000; + +delta = 0.01; +alpha_d = (5/4 * pi) - 2/8 * pi : delta : (5/4 * pi) + 2/8 * pi; +% Set different receiver positions rotated around the aperture +rcv_positions = zeros(numel(alpha_d), 3); +rcv_positions(1, :) = rcv; + +% Coordinate transformation +n3 = apex_dir; +n2 = w.main_face_normal; +n1 = cross( n2, n3 ); + +rho = norm( rcv - apex_point ); +z = rcv(3); + +rcv_pos_cylindrical = zeros(numel(alpha_d), 3); +for j = 1 : numel(alpha_d) + rcv_pos_cylindrical(j, :) = [rho, alpha_d(j), z]; +end + +for j = 2 : numel(rcv_positions(:, 1)) + rcv_positions(j, 1) = rcv_pos_cylindrical(j, 1) * cos( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 2) = rcv_pos_cylindrical(j, 1) * sin( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 3) = - rcv_pos_cylindrical(j, 3); +end + + +%% Variables +source_apex_direction = ( apex_point - src ) / norm( apex_point - src ); +rho = norm( apex_point - src ); % Distance of source to aperture point +r = norm( rcv - apex_point ); % Distance of receiver to aperture point + +c = 344; % Speed of Sound + +% Coordinate transformation with aperture point as origin +z = apex_dir; +y = w.main_face_normal; +x = cross( y, z ); % direction of main face of the wedge + +% Calculate angle between incedent ray from source to aperture point and source facing wedge +% side +u_i = dot( src - apex_point, x ); % coordinates in new coordinate system +v_i = dot( src - apex_point, y ); +alpha_i = atan2( v_i, u_i ); + +if alpha_i < 0 + alpha_i = alpha_i + pi; +end + + +% Angle between incedent ray from source to aperture point and aperture +% direction +theta_i = acos( dot( source_apex_direction, apex_dir ) ); +if theta_i > pi/2 + theta_i = pi - theta_i; +end + +n = w.opening_angle / pi; % Variable dependend on opening angle of the wedge + +lambda = c ./ freq; % Wavelength +k = (2 * pi) ./ lambda; % Wavenumber + +%% Diffraction Coefficient +A = sqrt( rho / ( r * (rho + r) ) ); +L = (rho * r) / (rho + r) * (sin( theta_i ))^2; + +epsilon = zeros( 1, numel(alpha_d) ); +data_1 = zeros( 2, numel( rcv_positions( :, 1 ) ) ); +data_2 = zeros( 2, numel( rcv_positions( :, 1 ) ) ); +data_3 = zeros( 2, numel( rcv_positions( :, 1 ) ) ); +small = 0.14; + +% D1_1 = cot( (pi + (alpha_d - alpha_i)) ./ (2 * n) ) * F( k * L * a_plus(alpha_d - alpha_i, n) ); +% D1_2 = n * exp( 1i * pi/4 ) * ... +% ( sqrt( 2 * pi .* k * L) * sgn( eps_plus(alpha_d - alpha_i, n) ) - ... +% 2 .* k * L * eps_plus(alpha_d - alpha_i, n) * exp( 1i * pi/4 ) ); +for j = 1 : numel(epsilon) + D1_1 = cot( (pi + (alpha_d(j) - alpha_i)) ./ (2 * n) ) * F( k * L * a_plus(alpha_d(j) - alpha_i, n) ); + magn_1 = zeros(1, numel(D1_1)); + for n = 1 : numel(magn_1) + magn_1(n) = 20 * log10( norm( D1_1(n) ) ); + end + D1_2 = n * exp( 1i * pi/4 ) * ... + ( sqrt( 2 * pi .* k * L) * sgn( eps_plus(alpha_d(j) - alpha_i, n) ) - ... + 2 .* k * L * eps_plus(alpha_d(j) - alpha_i, n) * exp( 1i * pi/4 ) ); + magn_2 = zeros(1, numel(D1_2)); + for m = 1 : numel(magn_2) + magn_2(m) = 20 * log10( norm( D1_2(m) ) ); + end + data_1(1, j) = magn_1; + data_1(2, j) = magn_2; + epsilon(j) = (alpha_d(j) - alpha_i) - 2 * pi * n * N_plus(alpha_d(j) - alpha_i, n) + pi; +end + +for i = 1 : numel(epsilon) + if norm(epsilon(i)) < small + data_3(1, i) = data_1(2, i); + else + data_3(1, i) = data_1(1, i); + end +end + +yyaxis left +plot( epsilon, data_1 ); +ylabel( 'magnitude in dB' ); +yyaxis right +plot( epsilon, alpha_d .* (360 / (2*pi)) ); +ylabel( 'alpha_d(°)' ); + +xlabel( 'eps' ); +legend( 'cot...', 'sqrt...', 'alpha_d' ); + +plot(epsilon, data_3(1, :)); + +for j = 1 : numel(epsilon) + D2_1 = cot( (pi - (alpha_d(j) - alpha_i)) ./ (2 * n) ) * F( k * L * a_minus(alpha_d(j) - alpha_i, n) ); + magn_1 = zeros(1, numel(D2_1)); + for n = 1 : numel(magn_1) + magn_1(n) = 20 * log10( norm( D2_1(n) ) ); + end + D2_2 = n * exp( 1i * pi/4 ) * ... + ( sqrt( 2 * pi .* k * L) * sgn( eps_minus(alpha_d(j) - alpha_i, n) ) - ... + 2 .* k * L * eps_minus(alpha_d(j) - alpha_i, n) * exp( 1i * pi/4 ) ); + magn_2 = zeros(1, numel(D2_2)); + for m = 1 : numel(magn_2) + magn_2(m) = 20 * log10( norm( D2_2(m) ) ); + end + data_2(1, j) = magn_1; + data_2(2, j) = magn_2; + epsilon(j) = 2 * pi * n * N_minus(alpha_d(j) - alpha_i, n) + pi - (alpha_d(j) - alpha_i); +end + +for i = 1 : numel(epsilon) + if norm(epsilon(i)) < small + data_3(2, i) = data_2(2, i); + else + data_3(2, i) = data_2(1, i); + end +end + +figure; +yyaxis left +plot( epsilon, cot( (pi - (2*pi*n*N_minus(alpha_d(j) - alpha_i, n) + pi - epsilon))/(2*n) ) ); + + +yyaxis left +plot( epsilon, data_2 ); +ylabel( 'magnitude in dB' ); +yyaxis right +plot( epsilon, alpha_d .* (360 / (2*pi)) ); +ylabel( 'alpha_d(°)' ); + +xlabel( 'eps' ); +legend( 'cot...', 'sqrt...', 'alpha_d' ); + +plot(epsilon, data_3(2, :)); + + + +%% Auxiliary Functions + +function res = a_plus( beta , n) + N = N_plus( beta, n ); + res = 2 * ( cos( ( 2 * pi * n * N - beta ) / 2 ) ) ^ 2; +end + +function res = a_minus( beta, n ) + N = N_minus(beta, n); + res = 2 * ( cos( ( 2 * pi * n * N - beta ) / 2 ) ) ^ 2; +end + +% Transition function with Kawai approximation +function out = F( X ) + out = zeros( 1, numel( X ) ); + for m = 1 : numel( out ) + if X(m) < 0.8 + out(m) = sqrt( pi * X(m) ) * ( 1 - ( sqrt(X(m)) / ( 0.7 * sqrt(X(m)) * 1.2 ) ) ) * exp( 1i * pi/4 * ( 1 - sqrt( X(m) / ( X(m) + 1.4 ) ) ) ); + else + out(m) = ( 1 - ( 0.8 / ( X(m) + 1.25 ) .^ 2 ) ) * exp( 1i * pi/4 * ( 1 - sqrt( X(m) / ( X(m) + 1.4 ) ) ) ); + end + end +end + +function res = eps_plus( beta, n ) + N = N_plus( beta, n ); + res = beta - 2 * pi * n * N + pi; +end + +function res = eps_minus( beta, n ) + N = N_minus( beta, n ); + res = 2 * pi * n * N + pi - beta; +end + +function res = sgn(x) + if x > 0 + res = 1; + else + res = -1; + end +end + +function N = N_plus( beta, n ) + x = (pi + beta) / (2 * pi * n); + a = floor(x); + b = floor(x) + 1; + d_a = norm(x - a); + d_b = norm(x - b); + if d_a < d_b + N = a; + else + N = b; + end +end + +function N = N_minus( beta, n ) + x = (beta - pi) / (2 * pi * n); + a = floor(x); + b = ceil(x); + d_a = norm(x - a); + d_b = norm(x - b); + if d_a < d_b + N = a; + else + N = b; + end +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/biot_tolstoy_test.m b/applications/VirtualAcoustics/diffraction/biot_tolstoy_test.m new file mode 100644 index 0000000..4e40400 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/biot_tolstoy_test.m @@ -0,0 +1,28 @@ +n1 = [1, 1, 0]; +n2 = [-1, 1, 0]; +loc = [0, 0, -3]; +w = itaInfiniteWedge(n1, n2, loc); +src = [-8, -1, -1]; +rcv = [8, -1, 1]; +ref = w.get_source_facing_side( src ); % reference wedge +p = [0, 0, 2]; +len = 18; +wt = 'corner'; + +finw = itaFiniteWedge(n1, n2, loc, len); +finw.point_on_aperture(p) + +res_finw = ita_diffraction_btm_finite_wedge(finw, src, rcv) +res_finw.pt; +res_finw.pf; + +p2 = [0, -1.0001, 5]; +finw.point_outside_wedge(p2) + +f_s = 44100; % sampling rate +l = 1024; % filter length + +res = ita_diffraction_btm_infinite_wedge(w, src, rcv, f_s, l) +res.pt; +% xlim([0, 0.006]); +res.pf; \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/btm_total_field_plot.m b/applications/VirtualAcoustics/diffraction/btm_total_field_plot.m new file mode 100644 index 0000000..4157e30 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/btm_total_field_plot.m @@ -0,0 +1,55 @@ +%% Config +n1 = [1 1 0]; +n2 = [-1 1 0]; +loc = [0 0 2]; +src = 5/sqrt(1) * [-1 0 0]; +rcv_start_pos = 5/sqrt(2) * [ 1 1 0]; +len = 15; +w = itaFiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc, len); +apex_point = w.get_aperture_point(src, rcv_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( pi, 3 * pi / 2, 500 ); + +% Set different receiver positions rotated around the aperture +rcv_positions = zeros(numel(alpha_d), 3); +rcv_positions(1, :) = rcv_start_pos; + +% Coordinate transformation +n3 = apex_dir; +n2 = w.main_face_normal; +n1 = cross( n2, n3 ); + +rho = norm( rcv_start_pos - apex_point ); +z = rcv_start_pos(3); + +rcv_pos_cylindrical = zeros(numel(alpha_d), 3); +for i = 1 : numel(alpha_d) + rcv_pos_cylindrical(i, :) = [rho, alpha_d(i), z]; +end + +for j = 2 : numel(rcv_positions(:, 1)) + rcv_positions(j, 1) = rcv_pos_cylindrical(j, 1) * cos( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 2) = rcv_pos_cylindrical(j, 1) * sin( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 3) = - rcv_pos_cylindrical(j, 3); +end + +%% Calculations +att_sum = itaAudio; +k = 2 * pi * freq ./ c; + +% UTD total wave field +for j = 1 : numel(rcv_positions(:, 1)) + r_dir = norm( rcv_positions(j, :) - src ); + E_dir = ( 1 / r_dir * exp( -1i .* k * r_dir ) )'; + att = ita_diffraction_btm_finite_wedge(w, src, rcv_positions(j, :), ); + % adding part of the incidence field if receiver not shadowed by wedge + if ~ita_diffraction_shadow_zone(w, src, rcv_positions(j, :)) + att.freqData = ( att.freqData + E_dir ); + end + att.freqData = att.freqData ./ E_dir; + att_sum = ita_merge( att_sum, att ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/diffraction_TF_IR.m b/applications/VirtualAcoustics/diffraction/diffraction_TF_IR.m new file mode 100644 index 0000000..79969c7 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/diffraction_TF_IR.m @@ -0,0 +1,227 @@ +%% Config +% wedge +n1 = [1, 1, 0]; +n2 = [-1, 1, 0]; +loc = [0, 0, -3]; +wedge = itaInfiniteWedge(n1, n2, loc); + +% fin wedge +len = 20; +finw = itaFiniteWedge(n1, n2, loc, len); + +% screen +n3 = [1, 0, 0]; +apex_dir = [0, 0, 1]; +screen = itaSemiInfinitePlane(n3, loc, apex_dir); + +% fin screen +fins = itaFiniteWedge(n3 + [0, 0.0001, 0], -n3 + [0, 0.0001, 0], loc, len); +% fins.aperture_direction = apex_dir; + +% source and receiver +src = 3 * [-1, 0, 0]; +angle1 = deg2rad(15); +angle2 = deg2rad(-15); +rcv1 = 3 * [cos(angle1), sin(angle1), 0]; +rcv2 = 3 * [cos(angle2), sin(angle2), 0]; + +% filter parameters +resolution = 5000; +% freq = linspace(20, 24000, resolution); +freq = ita_ANSI_center_frequencies; +c = 344; +f_s = 44100; +filter_length = 1024; +apex_point = wedge.get_aperture_point( src, rcv1 ); +R0 = norm( apex_point - src ) + norm( rcv1 - apex_point ); +tau0 = R0 / c; +tau = 0 : 1/f_s : tau0 + ( (filter_length - 1) * 1/f_s); + +ref_face_wedge = wedge.point_facing_main_side(src); +ref_face_screen = screen.point_facing_main_side(src); + +%% Filters +% Maekawa +% res_maekawa = itaResult; +res_maekawa1 = itaResult; +res_maekawa2 = itaResult; +res_maekawa3 = itaResult; +res_maekawa4 = itaResult; + +res_maekawa1.freqVector = freq; +res_maekawa2.freqVector = freq; +res_maekawa3.freqVector = freq; +res_maekawa4.freqVector = freq; + +res_maekawa1.freqData = ita_diffraction_maekawa(wedge, src, rcv1, freq, c); +res_maekawa2.freqData = ita_diffraction_maekawa(wedge, src, rcv2, freq, c); +res_maekawa3.freqData = ita_diffraction_maekawa(screen, src, rcv1, freq, c); +res_maekawa4.freqData = ita_diffraction_maekawa(screen, src, rcv2, freq, c); + +res_maekawa = ita_merge( res_maekawa1, res_maekawa2, res_maekawa3, res_maekawa4 ); + +% Maekawa approx +% res_maekawa_approx = itaResult; +res_maekawa_approx1 = itaResult; +res_maekawa_approx2 = itaResult; +res_maekawa_approx3 = itaResult; +res_maekawa_approx4 = itaResult; + +res_maekawa_approx1.freqVector = freq; +res_maekawa_approx2.freqVector = freq; +res_maekawa_approx3.freqVector = freq; +res_maekawa_approx4.freqVector = freq; + +res_maekawa_approx1.freqData = ita_diffraction_maekawa_approx(wedge, src, rcv1, freq, c); +res_maekawa_approx2.freqData = ita_diffraction_maekawa_approx(wedge, src, rcv2, freq, c); +res_maekawa_approx3.freqData = ita_diffraction_maekawa_approx(screen, src, rcv1, freq, c); +res_maekawa_approx4.freqData = ita_diffraction_maekawa_approx(screen, src, rcv2, freq, c); + +res_maekawa_approx = ita_merge( res_maekawa_approx1, res_maekawa_approx2, res_maekawa_approx3, res_maekawa_approx4 ); + +% UTD +% res_utd = itaResult; +res_utd1 = itaResult; +res_utd2 = itaResult; +res_utd3 = itaResult; +res_utd4 = itaResult; + +res_utd1.freqVector = freq; +res_utd2.freqVector = freq; +res_utd3.freqVector = freq; +res_utd4.freqVector = freq; + +res_utd1.freqData = ita_diffraction_utd(wedge, src, rcv1, freq, c); +res_utd2.freqData = ita_diffraction_utd(wedge, src, rcv2, freq, c); +res_utd3.freqData = ita_diffraction_utd(screen, src, rcv1, freq, c); +res_utd4.freqData = ita_diffraction_utd(screen, src, rcv2, freq, c); + +res_utd = ita_merge( res_utd1, res_utd2, res_utd3 ,res_utd4 ); + +% UTD approx +res_utd_approx1 = itaResult; +res_utd_approx2 = itaResult; +res_utd_approx3 = itaResult; +res_utd_approx4 = itaResult; + +res_utd_approx1.freqVector = freq; +res_utd_approx2.freqVector = freq; +res_utd_approx3.freqVector = freq; +res_utd_approx4.freqVector = freq; + +% res_utd_approx1.freqData = ita_diffraction_utd_approximated(wedge, src, rcv1, freq, c); +% res_utd_approx2.freqData = ita_diffraction_utd_approximated(wedge, src, rcv2, freq, c); +% res_utd_approx3.freqData = ita_diffraction_utd_approximated(screen, src, rcv1, freq, c); +% res_utd_approx4.freqData = ita_diffraction_utd_approximated(screen, src, rcv2, freq, c); + +res_utd_approx = ita_merge( res_utd_approx1, res_utd_approx2, res_utd_approx3, res_utd_approx4 ); + +% btms +res_btms1 = itaAudio; +res_btms2 = itaAudio; +res_btms3 = itaAudio; +res_btms4 = itaAudio; + +res_btms1.signalType = 'energy'; +res_btms2.signalType = 'energy'; +res_btms3.signalType = 'energy'; +res_btms4.signalType = 'energy'; + +res_btms1.samplingRate = f_s; +res_btms2.samplingRate = f_s; +res_btms3.samplingRate = f_s; +res_btms4.samplingRate = f_s; + +res_btms1.nSamples = filter_length; +res_btms2.nSamples = filter_length; +res_btms3.nSamples = filter_length; +res_btms4.nSamples = filter_length; + +res_btms1.timeData = ita_diffraction_btm_finite_wedge(finw, src, rcv1, tau, c, true ); +res_btms2.timeData = ita_diffraction_btm_finite_wedge(finw, src, rcv2, tau, c ); +res_btms3.timeData = ita_diffraction_btm_finite_wedge(fins, src, rcv1, tau, c ); +res_btms4.timeData = ita_diffraction_btm_finite_wedge(fins, src, rcv2, tau, c ); + +res_btms = ita_merge( res_btms1, res_btms2, res_btms3, res_btms4 ); + +%% Plots +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); + +str_lgnd = [ "wedge, illuminated"; "wedge, shadowed"; "screen, illuminated"; "screen, shadowed" ]; + +% maekawa +subplot( 2, 2, 1 ); +semilogx( freq', res_maekawa1.freqData_dB, '--', freq', res_maekawa2.freqData_dB, freq', res_maekawa3.freqData_dB, '--', freq', res_maekawa4.freqData_dB ); +title( 'Transfer function Maekawa method' ); +legend( str_lgnd ); +xlabel( 'f [Hz]' ); +ylabel( 'H(f) [dB]' ); +xlim( [freq(1), freq(end)] ); +ylim( [-45, -15] ); +grid on; + +% maekawa approx +subplot( 2, 2, 2 ); +semilogx( freq', res_maekawa_approx1.freqData_dB, '--', freq', res_maekawa_approx2.freqData_dB, freq', res_maekawa_approx3.freqData_dB, '--', freq', res_maekawa_approx4.freqData_dB ); +title( 'Transfer function Maekawa approx' ) +legend( str_lgnd ); +xlabel( 'f [Hz]' ); +ylabel( 'H(f) [dB]' ); +ylim( [-45, -15] ); +xlim( [freq(1), freq(end)] ); +grid on; + +% utd +subplot( 2, 2, 3 ); +semilogx( freq', res_utd1.freqData_dB, '--', freq', res_utd2.freqData_dB, freq', res_utd3.freqData_dB, '--', freq', res_utd4.freqData_dB ); +title( 'Transfer function UTD' ); +legend( str_lgnd ); +xlabel( 'f [Hz]' ); +ylabel( 'H(f) [dB]' ); +ylim( [-45, -15] ); +xlim( [freq(1), freq(end)] ); +grid on; + +% utd approx +subplot( 2, 2, 4 ); +semilogx( freq', res_utd_approx1.freqData_dB, '--', freq', res_utd_approx2.freqData_dB, freq', res_utd_approx3.freqData_dB, '--', freq', res_utd_approx4.freqData_dB ); +title( 'Transfer function UTD approx' ); +legend( str_lgnd ); +xlabel( 'f [Hz]' ); +ylabel( 'H(f) [dB]' ); +ylim( [-45, -15] ); +xlim( [freq(1), freq(end)] ); +grid on; + +% btms +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); + +subplot( 2, 2, 1); +semilogx( tau', res_btms1.timeData, '--', tau', res_btms2.timeData, tau', res_btms3.timeData, '--', tau', res_btms4.timeData ); +title( 'Impulse response BTMS' ); +legend( str_lgnd ); +xlabel( 't [s]' ); +ylabel( 'h(t)' ); +% ylim( [-45, -15] ); +xlim( [tau(1), tau(120)] ); +grid on; + +subplot(2, 2, 2); +semilogx( res_btms1.freqVector, res_btms1.freqData_dB, '--', res_btms1.freqVector, res_btms2.freqData_dB, res_btms1.freqVector, res_btms3.freqData_dB, '--', res_btms1.freqVector, res_btms4.freqData_dB ); +title( 'Transfer function BTMS' ); +legend( str_lgnd ); +xlabel( 'f [Hz]' ); +ylabel( 'H(f)' ); +ylim( [-45, -15] ); +xlim( [res_btms1.freqVector(1), res_btms1.freqVector(end)] ); +grid on; + + +%% more plots +% res_maekawa2.pf; +% res_maekawa_approx2.pf; +% +% res_utd.pf; +% res_utd_approx.pf; +% +% res_btms.pt; diff --git a/applications/VirtualAcoustics/diffraction/ita_align_points_around_aperture.m b/applications/VirtualAcoustics/diffraction/ita_align_points_around_aperture.m new file mode 100644 index 0000000..20452f1 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_align_points_around_aperture.m @@ -0,0 +1,89 @@ +function positions = ita_align_points_around_aperture( wedge, field_point, angles, point_of_rotation, ref_face ) +% ITA_ALIGN_POINTS_AROUND_APERTURE rotates field point around the aperture +% of the wedge around the point of rotation with the vector of angles respective to the reference face. +% Returns a set of 3-dim vectors, corresponding to the number of angles, in +% cartesian coordinates: dim( positions ) = #angles x 3, +% +% wedge: instance of class itaInfiniteWedge +% field_point: 3-dim point outside the wedge +% angles: field_point will be rotated by the angle(s) of this +% vector. angles cannot be higher than opening angle of +% the wedge. (radiant) +% point_of_rotation: 3-dim point on aperture the field point will be +% rotated around. +% ref_face: reference face angles are correspindingly denoted from +% (boolean: true if reference face is main face) + +%% Assertions +if ~( isa( wedge, 'itaInfiniteWedge' ) ) + error( 'wedge must an instance of class itaInfiniteWedge' ); +end +dim_fp = size( field_point ); +dim_pr = size( point_of_rotation ); +dim_an = size( angles ); +if dim_fp(1) ~= 1 + if dim_fp(2) ~= 1 || dim_fp(1) ~= 3 + error( 'field point must be of dimension 3' ); + end + field_point = field_point'; +end +if ~wedge.point_outside_wedge( field_point ) + error( 'field_point must be outside the wedge!' ); +end +if dim_pr(1) ~= 1 + if dim_pr(2) ~= 1 || dim_pr(1) ~= 3 + error( 'point of rotation must be of dimension 3' ); + end + point_of_rotation = point_of_rotation'; +end +if ~wedge.point_on_aperture( point_of_rotation ) + error( 'point of rotation must be on aperture of the wedge!' ); +end +if dim_an(1) ~= 1 + if dim_an(2) ~= 1 + error( 'angles have to be of dimension n x 1 or 1 x n!' ); + end + angles = angles'; +end +if any( angles > wedge.opening_angle ) + error( 'angles must be of smaller value than the opening angle of the wedge' ); +end + +%% Calculations +% initialize different positions rotated around the aperture +positions = zeros(numel( angles ), 3); +positions(1, :) = field_point; +apex_dir = wedge.aperture_direction; + +% Coordinate transformation +switch ref_face + case true + e3 = apex_dir; + e2 = wedge.main_face_normal; + e1 = cross( e2, e3 ); + phase_shift_x = acos( dot( e1, [1, 0, 0] ) ); + phase_shift_y = acos( dot( e2, [0, 1, 0] ) ); + case false + e3 = -apex_dir; + e2 = wedge.opposite_face_normal; + e1 = cross( e2, e3 ); + phase_shift_x = acos( dot( e1, [1, 0, 0] ) ); + phase_shift_y = - acos( dot( e2, [0, 1, 0] ) ); +end + +rho = repmat( norm( field_point - point_of_rotation ), numel( angles ), 1 ); +z = repmat( field_point(3), numel( angles ), 1 ); + +% Set field position in cylindrical coordinates +pos_cylindrical = zeros( numel( angles ), 3 ); +pos_cylindrical( :, 1 ) = rho; +pos_cylindrical( :, 2 ) = angles'; +pos_cylindrical( :, 3 ) = z; + +% Coordinate transformation to cartesian +positions( :, 1 ) = pos_cylindrical( :, 1 ) .* cos( pos_cylindrical( :, 2 ) + phase_shift_x ); +positions( :, 2 ) = pos_cylindrical( :, 1 ) .* sin( pos_cylindrical( :, 2 ) + phase_shift_y ); +positions( :, 3 ) = - pos_cylindrical( :, 3 ); + +end + diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_biot_tolstoy_jst.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_biot_tolstoy_jst.m new file mode 100644 index 0000000..e56b5d2 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_biot_tolstoy_jst.m @@ -0,0 +1,19 @@ +function att = ita_diffraction_biot_tolstoy_jst( wedge, source_pos, receiver_pos, sampling_rate, filter_length_samples ) + +if nargin < 4 + sampling_rate = 44100; +end +if nargin < 5 + filter_length_samples = 1024; +end + +T = 1 / sampling_rate; + +att = itaAudio(); +att.samplingRate = sampling_rate; +att.nSamples = filter_length_samples; + +Tau = 1:T:filter_length_samples; +att.timeData = - s * 1 / sinh( Tau ); + +end diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_infinite_wedge.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_infinite_wedge.m new file mode 100644 index 0000000..378d49f --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_infinite_wedge.m @@ -0,0 +1,108 @@ +function att = ita_diffraction_btm_infinite_wedge( wedge, source_pos, receiver_pos, sampling_rate, filter_length_samples, boundary_condition ) +%ITA_DIFFRATION_BIOT_TOLSTOY Summary of this function goes here +% Detailed explanation goes here + +%% Assertions +if nargin < 4 + sampling_rate = 44100; +end +if nargin < 5 + filter_length_samples = 1024; +end +if nargin < 6 + boundary_condition = 'hard'; +end +if ~( isequal(boundary_condition, 'hard') || isequal(boundary_condition, 'soft') ) + error('invalid boundary condition. Use hard or soft.'); +end +if ~isa(wedge, 'itaInfiniteWedge') + error('Invalid wedge input. use itaInfiniteWedge.') +end +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +if dim_src(2) ~= 3 + if sim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3' ) + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Reciever point(s) must be of dimension 3' ) + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +if dim_src(1) > dim_rcv(1) + dim_n = dim_src(1); + S = source_pos; + R = repmat( receiver_pos, dim_n, 1 ); +elseif dim_src(1) < dim_rcv(1) + dim_n = dim_rcv(1); + S = repmat( source_pos, dim_n, 1 ); + R = receiver_pos; +else + dim_n = dim_src(1); + S = source_pos; + R = receiver_pos; +end + + +%% Variables +Apex_Point = wedge.get_aperture_point( S, R ); +ref_face = wedge.point_facing_main_side( S ); +theta_S = wedge.get_angle_from_point_to_wedge_face( S, ref_face ); +theta_R = wedge.get_angle_from_point_to_wedge_face( R, ref_face ); +theta_i = wedge.get_angle_from_point_to_aperture( S, Apex_Point ); +theta_w = wedge.opening_angle; +SA = Norm( Apex_Point - S ); +AR = Norm( R - Apex_Point ); +r_S = SA .* sin(theta_i); +r_R = AR .* sin(theta_i); +z_R = (SA + AR) .* cos(theta_i); +ny = pi / theta_w; +L_0 = sqrt( (r_S + r_R).^2 + z_R.^2 ); +c = 343.21; % Speed of sound at 20°C + +att = itaAudio(); +att.signalType = 'energy'; +att.samplingRate = sampling_rate; +att.nSamples = filter_length_samples; +T = 1 / sampling_rate; +tau_0 = min( L_0 ) / c; +tau = ( tau_0 + T ) : T : ( tau_0 + T * (filter_length_samples) ); + +%% Calculation +Theta_S = repmat( theta_S, 1, numel(tau) ); +Theta_R = repmat( theta_R, 1, numel(tau) ); +Tau = repmat( tau, dim_n, 1 ); +eta = acosh( ( (c .* Tau).^2 - (r_S.^2 + r_R.^2 + z_R.^2) ) ./ ( 2 * r_S .* r_R ) ); + +phi1 = pi + Theta_S + Theta_R; +phi2 = pi + Theta_S - Theta_R; +phi3 = pi - Theta_S + Theta_R; +phi4 = pi - Theta_S - Theta_R; + +beta_pp = sin( ny * phi1 ) ./ ( cosh( ny .* eta ) - cos( ny * phi1 ) ); +beta_pm = sin( ny * phi2 ) ./ ( cosh( ny .* eta ) - cos( ny * phi2 ) ); +beta_mp = sin( ny * phi3 ) ./ ( cosh( ny .* eta ) - cos( ny * phi3 ) ); +beta_mm = sin( ny * phi4 ) ./ ( cosh( ny .* eta ) - cos( ny * phi4 ) ); + +switch boundary_condition + case 'hard' + beta = beta_pp + beta_pm + beta_mp + beta_mm; + case 'soft' + beta = -beta_pp + beta_pm + beta_mp - beta_mm; +end + +att.timeData = ( ( -(( c * ny ) / ( 2 * pi )) * beta ) ./ ( r_S .* r_R .* sinh( eta ) ) )' * T; + +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_total_field_plots.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_total_field_plots.m new file mode 100644 index 0000000..79d72fc --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_btm_total_field_plots.m @@ -0,0 +1,232 @@ +%% Config +n1 = [1 1 0]; +n2 = [-1 1 0]; +% apex_dir = [0 0 1]; +loc = [0 0 -3]; +len = 8; +src = 3/sqrt(1) * [-1 0 0]; +rcv_start_pos = 3/sqrt(2) * [ 1, 1, 0 ]; +rcv_end_pos = 3/sqrt(2) * [ 1, -1, 0 ]; +infw = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc); +finw = itaFiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc, len); +apex_dir = infw.aperture_direction; +% w.aperture_direction = apex_dir; +apex_point = infw.get_aperture_point(src, rcv_start_pos); +delta = 0.05; +c = 344; % Speed of sound + +% screen +n3 = [1, 0.00001, 0]; +n4 = [-1, 0.00001, 0]; +screen = itaFiniteWedge(n3/norm(n3), n4/norm(n4), loc, len); + + +wdg = finw; +%freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]; +f_sampling = 44100; +filter_length = 1024; +ref_face = wdg.point_facing_main_side( src ); +alpha_start = wdg.get_angle_from_point_to_wedge_face(rcv_start_pos, ref_face); +alpha_end = wdg.get_angle_from_point_to_wedge_face(rcv_end_pos, ref_face); +alpha_res = 200; + +alpha_d = linspace( alpha_start, alpha_end, alpha_res ); + + + +% Set Receiver Positions +rcv_positions = ita_align_points_around_aperture( wdg, rcv_start_pos, alpha_d, apex_point, ref_face ); + +% Direct field component for normalization of total field +f = linspace( 0, f_sampling/2, filter_length/2 + 1 ); +k = 2* pi * f / c; +R_dir = repmat( sqrt( sum( ( rcv_positions - src ).^2, 2 ) ), 1, numel(f) ); +E_dir = 1 ./ R_dir .* exp( -1i .* k .* R_dir ); + + +%% Calculations +in_shadow_zone = ita_diffraction_shadow_zone( wdg, src, rcv_positions ); + +%%% BTM infinite wedge %%%------------------------- +att_sum_btm_inf = ita_diffraction_btm_infinite_wedge( infw, src, rcv_positions, f_sampling, filter_length ); + +% f = att_sum_btm_inf.freqVector'; +% k = 2 * pi * f ./ c; + +% R_dir = repmat( sqrt( sum( ( rcv_positions - src ).^2, 2 ) ), 1, numel(f) ); +% E_dir = 1 ./ R_dir .* exp( -1i .* k .* R_dir ); + +% normalization +att_sum_btm_inf.freqData( :, ~in_shadow_zone ) = att_sum_btm_inf.freqData( :, ~in_shadow_zone ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_btm_inf.freqData = att_sum_btm_inf.freqData ./ E_dir'; +%--------------------------------------------------- + + +%%% BTM finite wedge %%%---------------------------- +R0 = norm( apex_point - src ) + norm( rcv_start_pos - apex_point ); +tau0 = R0 / c; +tau = tau0 : 1/f_sampling : tau0 + ( (filter_length - 1) * 1/f_sampling ); + +att = itaAudio(); +% att.timeVector = tau; +att.signalType = 'energy'; +att.samplingRate = f_sampling; +att.nSamples = filter_length; + +att_sum_btm_fin1 = itaAudio; +for j = 1 : numel(rcv_positions(:,1)) + att.timeData = ita_diffraction_btms_approx( wdg, src, rcv_positions(j, :), f_sampling, filter_length, c, true ); + att_sum_btm_fin1 = ita_merge( att_sum_btm_fin1, att ); +end +% normalization +att_sum_btm_fin1.freqData( :, ~in_shadow_zone ) = att_sum_btm_fin1.freqData( :, ~in_shadow_zone ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_btm_fin1.freqData = att_sum_btm_fin1.freqData ./ E_dir'; + + +att_sum_btm_fin2 = itaAudio; +for l = 1 : numel(rcv_positions(:,1)) + att.timeData = ita_diffraction_btms( wdg, src, rcv_positions(l, :), f_sampling, filter_length, c, true ); + att_sum_btm_fin2 = ita_merge( att_sum_btm_fin2, att ); +end +% normalization +att_sum_btm_fin2.freqData( :, ~in_shadow_zone ) = att_sum_btm_fin2.freqData( :, ~in_shadow_zone ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_btm_fin2.freqData = att_sum_btm_fin2.freqData ./ E_dir'; +%--------------------------------------------------- + + +%%% UTD %%%----------------------------------------- +att_sum_utd = itaResult; +att_sum_utd.freqVector = f; +att_sum_utd.freqData = ita_diffraction_utd( infw, src, rcv_positions, f, c ); +% normalization +att_sum_utd.freqData( :, ~in_shadow_zone' ) = att_sum_utd.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_utd.freqData = att_sum_utd.freqData ./ E_dir'; +%--------------------------------------------------- + + +%%% Maekawa %%%------------------------------------- +%--------------------------------------------------- + +%%% Svensson implementation %%%- + +% % fs = f_sampling; +% % closwedang = finw.wedge_angle; +% % rs = sqrt( src(1)^2 + src(2)^2 ); +% % thetas = infw.get_angle_from_point_to_wedge_face( src, ref_face ); +% % zs = 0; +% % zr = rcv_positions(1, 3) - src(3); +% % zw = [ loc(3) - src(3), len - src(3) ]; +% % Method = 'New'; +% % +% % att_sum_svensson = itaResult; +% % att_temp = itaResult; +% % for j = 1 : numel( rcv_positions(:, 1) ) +% % rr = sqrt( rcv_positions(j, 1)^2 + rcv_positions(j, 2)^2 ); +% % thetar = infw.get_angle_from_point_to_wedge_face( rcv_positions(j, :), ref_face ); +% % [ir,initdelay,singularterm] = EDB2wedge1st_int( fs, closwedang, rs, thetas, zs, rr, thetar, zr, zw, Method ); +% % att_temp.freqData = ir( ir ~= 0 ); +% % att_sum_svensson = ita_merge( att_sum_svensson, att_temp ); +% % end +% % +% % % normalization +% % % % att_sum_svensson.freqData( :, ~in_shadow_zone ) = att_sum_svensson.freqData( :, ~in_shadow_zone ) + ( E_dir( ~in_shadow_zone, : ) )'; +% % % % att_sum_svensson.freqData = att_sum_svensson.freqData ./ E_dir'; + +%-------------------------------- + + +%%% deviation between utd and btm %%%------------------- +% att_deviation = itaAudio; +% att_deviation.freqData = att_sum_utd.freqData ./ att_sum_btm_fin.freqData; +%------------------------------------------------------- + + +%% Plot +% f_plot = att_sum_btm_fin.freqVector( 6 : 50 : end ); +% +% figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); +% subplot( 2, 2, 1 ); +% res_utd = att_sum_utd.freqData_dB; +% plot( rad2deg( alpha_d ), (res_utd( 6 : 30 : end, : ) )' ); +% title( 'Tsingos et al.: UTD total wave field plot (Figure 6a)' ); +% legend( num2str( f_plot' ), 'Location', 'southwest' ) +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ); +% ylabel( 'dB SPL' ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% ylim( [-35, 10] ); +% grid on; +% +% subplot( 2, 2, 2 ); +% res_btm_inf = att_sum_btm_inf.freqData_dB; +% plot( rad2deg(alpha_d), ( res_btm_inf( 6 : 30 : end, : ) )' ); +% title( 'BTM diffraction for infinite wedge' ) +% legend( num2str( f_plot' ), 'Location', 'southwest' ) +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +% ylabel( 'dB SPL' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; +% +% subplot( 2, 2, 3 ); +% res_btm_fin = att_sum_btm_fin.freqData_dB; +% plot( rad2deg(alpha_d), ( res_btm_fin( 6 : 30 : end, : ) )' ); +% title( 'BTM diffraction for finite wedge' ); +% legend( num2str( f_plot' ), 'Location', 'southwest' ); +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ); +% ylabel( 'dB SPL' ); +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; +% +% subplot( 2, 2, 4 ); +% res_deviation = att_deviation.freqData_dB; +% plot( rad2deg(alpha_d), ( res_deviation( 6 : 30 : end, : ) )' ); +% title( 'Deviation between UTD and BTM' ) +% legend( num2str( f_plot' ), 'Location', 'southwest' ) +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +% ylabel( 'dB SPL' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; + +%% plots for thesis +f_plot = att_sum_btm_fin1.freqVector( 6 : 50 : end ); +str_hz = repmat( ' Hz', numel( f_plot ), 1 ); +res_btm_fin_approx = att_sum_btm_fin2.freqData_dB; +res_btm_fin = att_sum_btm_fin1.freqData_dB; + +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); + +subplot( 2, 2, 1 ); + +plot( rad2deg(alpha_d), ( res_btm_fin_approx( 6 : 50 : end, : ) )' ); +title( 'BTMS' ); +legend( [num2str( round( f_plot' ) ), str_hz], 'Location', 'southwest' ); +xlabel( 'theta_R [°]' ); +ylabel( 'p_{total} [dB]' ); +ylim( [-35, 10] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; + +subplot( 2, 2, 2 ); + +plot( rad2deg(alpha_d), ( res_btm_fin( 6 : 50 : end, : ) )' ); +title( 'BTMS with approx' ); +legend( [num2str( round( f_plot' ) ), str_hz], 'Location', 'southwest' ); +xlabel( 'theta_R [°]' ); +ylabel( 'p_{total} [dB]' ); +ylim( [-35, 10] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; + + + +% % subplot( 2, 2, 4 ); +% % plot( rad2deg(alpha_d), att_sum_svensson.freqData_dB ); +% % title( 'BTM diffraction with svensson implementation' ) +% % % legend( num2str( f_plot' ), 'Location', 'southwest' ) +% % xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +% % ylabel( 'dB SPL' ) +% % % ylim( [-35, 10] ); +% % xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% % grid on; diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_btms.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_btms.m new file mode 100644 index 0000000..3839c52 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_btms.m @@ -0,0 +1,367 @@ +function [res, offs] = ita_diffraction_btms( wedge, source_pos, receiver_pos, f_s, N, speed_of_sound, approx ) +%ITA_DIFFRACTION_BTM_FINITE_WEDGE Summary of this function goes here +% Detailed explanation goes here + +% if nargin < 4 +% sampling_rate = 44100; +% end +% if nargin < 5 +% filter_length_samples = 1024; +% end +% if nargin < 6 +% boundary_condition = 'hard'; +% end +% if ~( isequal(boundary_condition, 'hard') || isequal(boundary_condition, 'soft') ) +% error('invalid boundary condition. Use hard or soft.'); +% end +if nargin < 7 + approx = true; +end +if ~isa( wedge, 'itaFiniteWedge') + error( 'Wedge must be object from class itaFiniteWedge' ); +end +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +if dim_src(2) ~= 3 + if sim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3' ) + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Reciever point(s) must be of dimension 3' ) + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end + +%% Variables +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +Apex_Dir = wedge.aperture_direction; +ref_face = wedge.point_facing_main_side( source_pos ); +theta_S = wedge.get_angle_from_point_to_wedge_face( source_pos, ref_face ); +theta_R = wedge.get_angle_from_point_to_wedge_face( receiver_pos, ref_face ); +theta_i = wedge.get_angle_from_point_to_aperture( source_pos, Apex_Point ); +theta_w = wedge.opening_angle; +SA = norm( Apex_Point - source_pos ); +AR = norm( receiver_pos - Apex_Point ); +r_S = SA .* sin( theta_i ); +r_R = AR .* sin( theta_i ); +z_S = dot( source_pos - wedge.aperture_start_point, Apex_Dir ); +z_R = dot( receiver_pos - wedge.aperture_start_point, Apex_Dir ); +% z_apex = dot( Apex_Point - wedge.aperture_start_point, Apex_Dir ); + +ny = pi / theta_w; +R_0 = sqrt( (r_S + r_R).^2 + (z_R - z_S).^2 ); +c = speed_of_sound; + +% res = itaAudio(); +% res.signalType = 'energy'; +% res.samplingRate = sampling_rate; +% res.nSamples = filter_length_samples; + +T = 1 / f_s; +tau_0 = R_0 / c; +% tau = ( tau_0 ) : T : ( tau_0 + T * (N - 1) ); + +tau = 0 : T : T*(N-1); + +tau_offset = tau_0 / T; + +% +% tau = time_vec; +% T = tau(2) - tau(1); + +%% Calculate positions of secondary sources on aperture +Z = secondary_source_coordinates_on_aperture( r_S, z_S, r_R, z_R, tau, c ); +z_n_u = Z(1, :); +z_n_l = Z(2, :); + +dZ = secondary_source_coordinates_on_aperture( r_S, z_S, r_R, z_R, tau + 0.5*T, c ) -... + secondary_source_coordinates_on_aperture( r_S, z_S, r_R, z_R, tau - 0.5*T, c ); + + +tau_end_upper_branch = ( norm( wedge.aperture_end_point - source_pos ) + norm( receiver_pos - wedge.aperture_end_point) ) / c; +tau_end_lower_branch = ( norm( wedge.aperture_start_point - source_pos ) + norm( receiver_pos - wedge.aperture_start_point) ) / c; + +mask_u = ( tau <= tau_end_upper_branch ); +mask_l = ( tau <= tau_end_lower_branch ); + +% Choose longer branch of the edge from aperture point +if sum(mask_u) >= sum(mask_l) % upper branch + z_n = ( z_n_u(mask_u) )'; + mask_n = mask_u; + dz = dZ(1, :); + dz = ( dz(mask_u) )'; +else % lower branch + z_n = ( z_n_l(mask_l) )'; + mask_n = mask_l; + dz = dZ(2, :); + dz = ( dz(mask_l) )'; +end + +Apex_Start = repmat( wedge.aperture_start_point, numel( z_n ), 1 ); +Z_n = repmat( z_n, 1, 3 ); +Apex_Dir = repmat( wedge.aperture_direction, numel( z_n ), 1 ); +Sec_Source_Pos_On_Apex_n = Apex_Start + ( Z_n .* Apex_Dir ); + +%% Filter variables +% alpha_n = pi/2 - wedge.get_angle_from_point_to_aperture( source_pos, Sec_Source_Pos_On_Apex_n ); +% gamma_n = pi/2 - wedge.get_angle_from_point_to_aperture( receiver_pos, Sec_Source_Pos_On_Apex_n ); + +m_n = Norm( Sec_Source_Pos_On_Apex_n - source_pos ); +l_n = Norm( Sec_Source_Pos_On_Apex_n - receiver_pos ); + +% Calculation of dz for numerical integration +% dz_n = z_n - [ z_apex; z_n(1 : end-1) ]; +% dz = ( [dz_n(2 : end); dz_n(end)] ./ 2 ) + ( dz_n ./ 2 ); + +%% Filter Calculation +y = ( m_n .* l_n + (z_n - z_S) .* (z_n - z_R) ) ./ ( r_S * r_R ); +A = y + sqrt( y.^2 - 1 ); +cosh_approx = 0.5 * ( A.^(ny) + A.^(-ny) ); + +% eta = acosh( ( 1 + sin( alpha_n ) .* sin( gamma_n ) ) ./ ( cos( alpha_n ) .* cos( gamma_n ) ) ); + + +phi1 = pi + theta_S + theta_R; +phi2 = pi + theta_S - theta_R; +phi3 = pi - theta_S + theta_R; +phi4 = pi - theta_S - theta_R; + +beta_pp = sin( ny * phi1 ) ./ ( cosh_approx - cos( ny * phi1 ) ); +beta_pm = sin( ny * phi2 ) ./ ( cosh_approx - cos( ny * phi2 ) ); +beta_mp = sin( ny * phi3 ) ./ ( cosh_approx - cos( ny * phi3 ) ); +beta_mm = sin( ny * phi4 ) ./ ( cosh_approx - cos( ny * phi4 ) ); + +%% Seperate consideration of potential singularities at reflection or shadow boundary +z_rel = dz(1); +psi = theta_i; +rho = r_R / r_S; + +% theta_RB = pi - theta_S; +% theta_SB = pi + theta_S; + +symmetrical = (psi == pi/2) || (rho == 1); +if symmetrical + B0_pp = B0( R_0, rho, psi, ny, phi1 ); + B0_pm = B0( R_0, rho, psi, ny, phi2 ); + B0_mp = B0( R_0, rho, psi, ny, phi3 ); + B0_mm = B0( R_0, rho, psi, ny, phi4 ); + + B1_pp = B1( R_0, rho, ny, phi1 ); + B1_pm = B1( R_0, rho, ny, phi2 ); + B1_mp = B1( R_0, rho, ny, phi3 ); + B1_mm = B1( R_0, rho, ny, phi4 ); + + B3_pp = B3( R_0, rho, psi ); + B3_pm = B3( R_0, rho, psi ); + B3_mp = B3( R_0, rho, psi ); + B3_mm = B3( R_0, rho, psi ); + + if (psi == pi/2) && (rho == 1) + B4_pp = B4( ny, phi1 ); + B4_pm = B4( ny, phi2 ); + B4_mp = B4( ny, phi3 ); + B4_mm = B4( ny, phi4 ); + + h1_0 = - ( ny / (2*pi) ) * ( B4_pp / sqrt(B1_pp) ) * atan( z_rel / sqrt(B1_pp) ); + h2_0 = - ( ny / (2*pi) ) * ( B4_pm / sqrt(B1_pm) ) * atan( z_rel / sqrt(B1_pm) ); + h3_0 = - ( ny / (2*pi) ) * ( B4_mp / sqrt(B1_mp) ) * atan( z_rel / sqrt(B1_mp) ); + h4_0 = - ( ny / (2*pi) ) * ( B4_mm / sqrt(B1_mm) ) * atan( z_rel / sqrt(B1_mm) ); + else + h1_0 = - ( ny / (2*pi) ) * ( B0_pp / (B3_pp - B1_pp) ) * ( ( 1 / sqrt(B1_pp) ) * atan( z_rel / sqrt(B1_pp) ) - ( 1 / sqrt(B3_pp) ) * atan( z_rel / sqrt(B3_pp) ) ); + h2_0 = - ( ny / (2*pi) ) * ( B0_pm / (B3_pm - B1_pm) ) * ( ( 1 / sqrt(B1_pm) ) * atan( z_rel / sqrt(B1_pm) ) - ( 1 / sqrt(B3_pm) ) * atan( z_rel / sqrt(B3_pm) ) ); + h3_0 = - ( ny / (2*pi) ) * ( B0_mp / (B3_mp - B1_mp) ) * ( ( 1 / sqrt(B1_mp) ) * atan( z_rel / sqrt(B1_mp) ) - ( 1 / sqrt(B3_mp) ) * atan( z_rel / sqrt(B3_mp) ) ); + h4_0 = - ( ny / (2*pi) ) * ( B0_mm / (B3_mm - B1_mm) ) * ( ( 1 / sqrt(B1_mm) ) * atan( z_rel / sqrt(B1_mm) ) - ( 1 / sqrt(B3_mm) ) * atan( z_rel / sqrt(B3_mm) ) ); + end +else + B0_pp = B0( R_0, rho, psi, ny, phi1 ); + B0_pm = B0( R_0, rho, psi, ny, phi2 ); + B0_mp = B0( R_0, rho, psi, ny, phi3 ); + B0_mm = B0( R_0, rho, psi, ny, phi4 ); + + B1_pp = B1( R_0, rho, ny, phi1 ); + B1_pm = B1( R_0, rho, ny, phi2 ); + B1_mp = B1( R_0, rho, ny, phi3 ); + B1_mm = B1( R_0, rho, ny, phi4 ); + + B2_pp = B2( R_0, rho, psi ); + B2_pm = B2( R_0, rho, psi ); + B2_mp = B2( R_0, rho, psi ); + B2_mm = B2( R_0, rho, psi ); + + B3_pp = B3( R_0, rho, psi ); + B3_pm = B3( R_0, rho, psi ); + B3_mp = B3( R_0, rho, psi ); + B3_mm = B3( R_0, rho, psi ); + + q1 = 4 * B3_pp - B2_pp^2; + q2 = 4 * B3_pm - B2_pm^2; + q3 = 4 * B3_mp - B2_mp^2; + q4 = 4 * B3_mm - B2_mm^2; + + + C1 = ( B0_pp * B2_pp ) / ( B1_pp * B2_pp^2 + ( B1_pp - B3_pp )^2 ); + C2 = ( B0_pm * B2_pm ) / ( B1_pm * B2_pm^2 + ( B1_pm - B3_pm )^2 ); + C3 = ( B0_mp * B2_mp ) / ( B1_mp * B2_mp^2 + ( B1_mp - B3_mp )^2 ); + C4 = ( B0_mm * B2_mm ) / ( B1_mm * B2_mm^2 + ( B1_mm - B3_mm )^2 ); + + D1 = 0.5 * log( abs( ( B3_pp * ( z_rel^2 + B1_pp ) ) / ( B1_pp * ( z_rel^2 + B2_pp * z_rel + B3_pp ) ) ) ); + D2 = 0.5 * log( abs( ( B3_pm * ( z_rel^2 + B1_pm ) ) / ( B1_pm * ( z_rel^2 + B2_pm * z_rel + B3_pm ) ) ) ); + D3 = 0.5 * log( abs( ( B3_mp * ( z_rel^2 + B1_mp ) ) / ( B1_mp * ( z_rel^2 + B2_mp * z_rel + B3_mp ) ) ) ); + D4 = 0.5 * log( abs( ( B3_mm * ( z_rel^2 + B1_mm ) ) / ( B1_mm * ( z_rel^2 + B2_mm * z_rel + B3_mm ) ) ) ); + + E1 = ( ( B1_pp - B3_pp ) / ( sqrt(B1_pp) * B2_pp ) ) * atan( z_rel / sqrt(B1_pp) ); + E2 = ( ( B1_pp - B3_pm ) / ( sqrt(B1_pm) * B2_pm ) ) * atan( z_rel / sqrt(B1_pm) ); + E3 = ( ( B1_pp - B3_mp ) / ( sqrt(B1_mp) * B2_mp ) ) * atan( z_rel / sqrt(B1_mp) ); + E4 = ( ( B1_pp - B3_mm ) / ( sqrt(B1_mm) * B2_mm ) ) * atan( z_rel / sqrt(B1_mm) ); + + G1 = ( ( 2 * ( B3_pp - B1_pp ) - B2_pp^2 ) / ( 2 * B2_pp ) ) * F( q1, z_rel, B2_pp, rho, psi ); + G2 = ( ( 2 * ( B3_pm - B1_pm ) - B2_pm^2 ) / ( 2 * B2_pm ) ) * F( q2, z_rel, B2_pm, rho, psi ); + G3 = ( ( 2 * ( B3_mp - B1_mp ) - B2_mp^2 ) / ( 2 * B2_mp ) ) * F( q3, z_rel, B2_mp, rho, psi ); + G4 = ( ( 2 * ( B3_mm - B1_mm ) - B2_mm^2 ) / ( 2 * B2_mm ) ) * F( q4, z_rel, B2_mm, rho, psi ); + + h1_0 = ( ny / (2*pi) ) * C1 * ( D1 + E1 + G1 ); + h2_0 = ( ny / (2*pi) ) * C2 * ( D2 + E2 + G2 ); + h3_0 = ( ny / (2*pi) ) * C3 * ( D3 + E3 + G3 ); + h4_0 = ( ny / (2*pi) ) * C4 * ( D4 + E4 + G4 ); +end + +%% Filter Results +switch wedge.is_boundary_condition_hard + case true + beta = ( beta_pp + beta_pm + beta_mp + beta_mm ); % factor 2 from Svensson - An analytic secondary source model... eq (17) + h_0 = h1_0 + h2_0 + h3_0 + h4_0; + case false + beta = ( -beta_pp + beta_pm + beta_mp - beta_mm ); + h_0 = - h1_0 + h2_0 + h3_0 - h4_0; +end + +scaling = ( mask_u + mask_l )'; % factor 2 from Svensson - An analytic secondary source model... eq (17) +integrand = ( (scaling(mask_n) .* beta) ./ (m_n .* l_n) ); + +h_diffr = zeros( numel(tau), 1 ); +h_diffr(mask_n) = - ( ny / (4*pi) ) * integrand .* dz; +if approx + h_diffr(1) = h_0; +end +% res.timeData = h_diffr; + +if nargout > 1 + res = h_diffr; + offs = tau_0; +else + res = h_diffr; +end + +%% test output + +% res = [h1_0 , h2_0, h3_0, h4_0]; + + +end + + +% function Z = secondary_source_coordinates_on_aperture( r_S, z_S, r_R, z_R, tau, c ) +% z_RS = z_R - z_S; +% +% p_n = ( c * tau ).^2; +% +% a = z_RS.^2 + r_S.^2 + r_R.^2; +% b = 2 * z_RS; +% q = (z_RS.^2 + r_R.^2 - r_S.^2).^2; +% +% u_n = b.^2 - 4 * p_n; +% v_n = 4 * z_RS .* p_n - 4 * z_RS .* (z_RS.^2 + r_R.^2 - r_S.^2); +% w_n = p_n.^2 - 2 * a .* p_n + q; +% +% z_n_u = (-v_n ./ (2 * u_n)) + sqrt( 0.25 * (v_n ./ u_n).^2 - (w_n ./ u_n) ) + z_S; +% z_n_l = (-v_n ./ (2 * u_n)) - sqrt( 0.25 * (v_n ./ u_n).^2 - (w_n ./ u_n) ) + z_S; +% +% z_n_u = real(z_n_u); +% z_n_l = real(z_n_l); +% +% Z = [ z_n_u; z_n_l ]; +% end + +function Z = secondary_source_coordinates_on_aperture( r_S, z_S, r_R, z_R, tau, c ) + a = z_S^2 + z_R^2 + r_S^2 + r_R^2; + b = 2 * (z_S + z_R); + d = 2 * ( z_S * (z_R^2 + r_R^2) + z_R * (z_S^2 + r_S^2) ); + e = (z_S^2 + r_S^2) * (z_R^2 + r_R^2); + f = z_S^2 + z_R^2 + r_S^2 + r_R^2 + 4*z_R*z_S; + + p_n = (c * tau).^2; + u_n = b.^2 + 4*a - 4.*p_n - 4*f; + v_n = 2 .* p_n .* b - 2 .* a .* b + 4*d; + w_n = p_n.^2 - 2 .* a .* p_n + a.^2 - 4*e; + + det = 0.25 * (v_n ./ u_n).^2 - (w_n ./ u_n); + det(det < 0) = 0; + + z_n_u = (-v_n ./ (2 * u_n)) + sqrt( det ); + z_n_l = (-v_n ./ (2 * u_n)) - sqrt( det ); + + Z = [z_n_u; z_n_l]; +end +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end + +function res = B0( R_0, rho, psi, ny, phi ) + res = ( 4 * R_0^2 * rho^3 * sin( ny * phi ) ) / ( ny^2 * (1 + rho)^4 * ( (1 + rho)^2 * sin( psi )^2 - 2 * rho ) ); +end + +function res = B1( R_0, rho, ny, phi ) + res = ( 4 * R_0^2 * rho^2 * sin( (ny * phi) / 2 )^2 ) / ( ny^2 * (1 + rho)^4 ); +end + +function res = B2( R_0, rho, psi ) + res = - ( 2 * R_0 * (1 - rho) * rho * cos( psi ) ) / ( (1 + rho) * ( (1 + rho)^2 * sin( psi )^2 - 2 * rho ) ); +end + +function res = B3( R_0, rho, psi ) + res = ( 2 * R_0^2 * rho^2 ) / ( (1 + rho)^2 * ( (1 + rho)^2 * sin( psi )^2 - 2 * rho ) ); +end + +function res = B4( ny, phi ) + res = sin( ny * phi ) / ( 2 * ny^2 ); +end + +function res = F( q, z_rel, B2, rho, psi ) + if q < 0 + res = ( 1 / sqrt(-q) ) * log( abs( ( 2*z_rel + B2 - sqrt(-q) * B2 + sqrt(-q) ) / ( 2*z_rel + B2 - sqrt(-q) * B2 - sqrt(-q) ) ) ); + elseif q > 0 + res = ( 2 / sqrt(q) ) * ( atan( (2*z_rel + B2) / sqrt(q) ) - atan( B2 / sqrt(q) ) ); + else + res = ( 4 * z_rel ) / ( B2 * ( 2*z_rel + B2 ) ); + end + if abs( sin(psi)^2 - ( 2*rho ) / ( ( 1 + rho )^2 ) ) <= 10^(-6) + res = 0; + end +end + +% function Z = test_sec_source_position_on_aperture(r_S, r_R, z_S, z_R, z_A, tau) +% c = 343.21; +% +% z_S_rel = z_S - z_A; +% z_R_rel = z_R - z_A; +% +% D = c .* tau; +% m_0 = sqrt( r_S.^2 + z_S_rel.^2 ); +% l_0 = sqrt( r_R.^2 + z_R_rel.^2 ); +% K = m_0.^2 - l_0.^2 - D.^2; +% B = ( 2 * D.^2 .* z_R_rel - K .* ( z_S_rel - z_R_rel ) ) ./ ( ( z_S_rel - z_R_rel ).^2 - D.^2 ); +% C = ( 0.25 .* K.^2 - D.^2 * l_0.^2 ) ./ ( ( z_S_rel - z_R_rel ).^2 - D.^2 ); +% +% z_rel1 = ( -B + sqrt( B.^2 - 4 .* C ) ) ./ 2; +% z_rel2 = ( -B - sqrt( B.^2 - 4 .* C ) ) ./ 2; +% +% Z = [ z_rel1; z_rel2 ]'; +% end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_btms_approx.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_btms_approx.m new file mode 100644 index 0000000..7567c05 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_btms_approx.m @@ -0,0 +1,149 @@ +function [res, offs] = ita_diffraction_btms_approx( wedge, source, receiver, f_s, filter_length, c, integrand_approx, theta_0 ) +%ITA_DIFFRACTION_BTMS_APPROX Summary of this function goes here +% Detailed explanation goes here + +%% Assertions +if ~isa( wedge, 'itaFiniteWedge') + error( 'Wedge must be object from class itaFiniteWedge' ); +end +dim_src = size( source ); +dim_rcv = size( receiver ); +% dim_tau = size( time_vec ); +if dim_src(2) ~= 3 + if sim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3' ) + end + source = source'; + dim_src = size( source ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Reciever point(s) must be of dimension 3' ) + end + receiver = receiver'; + dim_rcv = size( receiver ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +% if dim_tau(1) ~= 1 +% if dim_tau(2) ~= 1 +% error( 'Invalid time vector. Use (1 x n) or (n x 1) for n-dimensional time vector!' ); +% end +% time_vec = time_vec'; +% dim_tau = size( time_vec ); +% end +if nargin < 8 + theta_0 = 0.1; +end +if nargin < 7 + integrand_approx = true; +end + + +%% +% if ~in_shadow_zone +% h_diffr_approx = ( zeros( dim_tau ) )'; +% else +% apex_point = wedge.get_aperture_point( source, receiver ); +% src_ap_dir = (apex_point - source) / norm(apex_point - source); +% rho_S = norm( apex_point - source ); +% rho_R = norm( receiver - apex_point ); +% receiver_SB = ( rho_S + rho_R ) * src_ap_dir + source; +% +% ref_face = wedge.point_facing_main_side( source ); +% theta_tolerance = deg2rad(5); +% theta_R = wedge.get_angle_from_point_to_wedge_face( receiver, ref_face ); +% theta_SB = wedge.get_angle_from_point_to_wedge_face( receiver_SB, ref_face ) + theta_tolerance; +% +% receiver_SB = ita_align_points_around_aperture( wedge, receiver_SB, theta_tolerance, apex_point, ref_face ); +% theta_rel = theta_R - theta_SB; +% if theta_rel > pi/4 +% theta_rel = pi/4; +% end +% +% R_dir = norm( receiver - source ); +% +% % determine nearest time sample at direct impact +% tau_eps = min( abs( time_vec - R_dir/c ) ); +% tau_dir = R_dir/c + tau_eps; +% +% h_dir_SB = ( zeros( dim_tau ) )'; +% h_dir_SB( time_vec == tau_dir ) = 1 / R_dir; +% h_diffr_SB = ita_diffraction_btm_finite_wedge( wedge, source, receiver_SB, time_vec, c, integrand_approx ); +% +% c_norm = h_dir_SB ./ h_diffr_SB; +% c_approx = 1 + ( c_norm - 1 ) .* exp( - theta_rel / theta_0 ); +% +% h_diffr = ita_diffraction_btm_finite_wedge( wedge, source, receiver, time_vec, c, integrand_approx ); +% h_diffr_approx = c_approx .* h_diffr; +% end + +%% Variables +in_shadow_zone = ita_diffraction_shadow_zone( wedge, source, receiver ); +ref_face = wedge.point_facing_main_side( source ); +apex_point = wedge.get_aperture_point( source, receiver ); +src_ap_dir = (apex_point - source) / norm(apex_point - source); +rho_S = norm( apex_point - source ); +rho_R = norm( receiver - apex_point ); +receiver_SB = ( rho_S + rho_R ) * src_ap_dir + source; +theta_SB = wedge.get_angle_from_point_to_wedge_face( receiver_SB, ref_face ); + +tau_0 = (rho_S + rho_R)/c; + +if ~in_shadow_zone + if nargout > 1 + res = zeros( filter_length, 1 ); + offs = tau_0 / (1/f_s); + else + res = zeros( filter_length, 1 ); + end +else + + ref_face = wedge.point_facing_main_side( source ); + theta_tolerance = deg2rad(0.000001); + theta_R = wedge.get_angle_from_point_to_wedge_face( receiver, ref_face ); + theta_SB = wedge.get_angle_from_point_to_wedge_face( receiver_SB, ref_face ) + theta_tolerance; + + receiver_SB = ita_align_points_around_aperture( wedge, receiver_SB, theta_tolerance + theta_SB, apex_point, ref_face ); + theta_rel = theta_R - theta_SB; + if theta_rel > pi/4 + theta_rel = pi/4; + end + + f1 = linspace( 0, f_s/2, filter_length/2 + 1 ); + f2 = linspace( -f_s/2, 0, filter_length/2 + 1 ); + f = [f1, f2(2 : end-1)]; + k = 2* pi * f / c; + R_dir = norm( receiver_SB - source ); + + H_i_SB = 1 / R_dir;% * exp( -1i * k * R_dir ); + + if nargout > 1 + [h_diffr_SB, tau_offs_SB] = ita_diffraction_btms( wedge, source, receiver_SB, f_s, filter_length, c, integrand_approx ); + [h_diffr, tau_offs] = ita_diffraction_btms( wedge, source, receiver, f_s, filter_length, c, integrand_approx ); + else + h_diffr_SB = ita_diffraction_btms( wedge, source, receiver_SB, f_s, filter_length, c, integrand_approx ); + h_diffr = ita_diffraction_btms( wedge, source, receiver, f_s, filter_length, c, integrand_approx ); + end + H_diffr_SB_double_sided = fft( h_diffr_SB ); + H_diffr_double_sided = fft( h_diffr ); + + H_diffr_SB = H_diffr_SB_double_sided; + H_diffr = H_diffr_double_sided; +% H_diffr_SB = H_diffr_SB_double_sided( 1 : filter_length/2 + 1 ); +% H_diffr = H_diffr_double_sided( 1 : filter_length/2 + 1 ); + + c_norm = H_i_SB' ./ H_diffr_SB; + c_approx = 1 + ( c_norm - 1 ) .* exp( - theta_rel / theta_0 ); + + H_diffr_approx = c_approx .* H_diffr; + h_diffr_approx = ifft( H_diffr_approx ); + if nargout > 1 + res = h_diffr_approx; + offs = tau_offs; + else + res = h_diffr_approx; + end +end + diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel.m new file mode 100644 index 0000000..962b2ac --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel.m @@ -0,0 +1,22 @@ +function Y = ita_diffraction_kawai_approx_fresnel( X ) +%ita_diffraction_kawai_approx_fresnel Evaluates the approximation of the Fresnel integral equation by Kawai et al. +% +% Example: +% +% X = logspace( -3, 1 ) +% Y = ita_diffraction_kawai_approx_fresnel( X ) +% + +if any( X < 0 ) + error( 'No negative values for Kawai approximation of Fresnel integral allowed' ) +end + +X_l_idx = X < 0.8; +X_l = X( X_l_idx ); +X_geq_idx = ( X >= 0.8 ); +X_geq = X( X_geq_idx ); + +Y( X_l_idx ) = sqrt( pi * X_l ) .* ( 1 - sqrt( X_l ) ./ ( 0.7 * sqrt( X_l ) + 1.2 ) ) .* exp( 1i * pi/4 * ( 1 - sqrt( X_l ./ ( X_l + 1.4 ) ) ) ); +Y( X_geq_idx ) = ( 1 - 0.8 ./ ( X_geq + 1.25 ) .^ 2 ) .* exp( 1i * pi/4 * ( 1 - sqrt( X_geq ./ ( X_geq + 1.4 ) ) ) ); + +end diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel_plot.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel_plot.m new file mode 100644 index 0000000..a13f271 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_kawai_approx_fresnel_plot.m @@ -0,0 +1,17 @@ +%% Plot the Kawai approximation of Fresnel integral in range 0.001 to 10.0 + +% Evaluation +X = logspace( -3, 1 ); +Y = ita_diffraction_kawai_approx_fresnel( X ); + +% Plot +title( 'Kawai approximation of Fresnel integral' ) +yyaxis left +semilogx( X, abs( Y ) ) +xlabel( 'X' ); +ylabel( 'magnitude' ) +yyaxis right +semilogx( X, rad2deg( angle( Y ) ) ) +ylabel( 'phase in degree' ); +ylim( [0, 50] ); +legend( 'magnitude', 'phase', 'location', 'north' ) diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa.m new file mode 100644 index 0000000..aabf257 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa.m @@ -0,0 +1,67 @@ +function H_diffr = ita_diffraction_maekawa( wedge, source_pos, receiver_pos, frequencies, speed_of_sound ) +% Calculates the attenuation filter(s) at a diffraction wedge for source +% and receiver location(s) at given frequencies +% +% wedge: diffracting wedge (itaInfiniteWedge or derived class +% itaFiniteWedge) +% source_pos: position of the source +% receiver_pos: position of the receiver + +%% Assertions +assert( isa( wedge, 'itaInfiniteWedge' ) ) +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +dim_f = size( frequencies ); +if dim_src(2) ~= 3 + if dim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3' ) + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Receiver point(s) must be of dimension 3' ) + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +if dim_f(1) ~= 1 + if dim_f(2) ~= 1 + error( 'Invalid frequency. Use row or column vector' ); + end + frequencies = frequencies'; +end + +%% Variables +% att = itaResult(); +% att.freqVector = frequencies'; + +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +r_dir = Norm( receiver_pos - source_pos ); +detour = Norm( Apex_Point - source_pos ) + Norm( receiver_pos - Apex_Point ) - Norm( receiver_pos - source_pos ); +c = speed_of_sound; % Speed of sound in air with a temprature of 20°C +lambda = c ./ frequencies; + +N = 2 * detour ./ lambda; % Fresnel number N + +in_shadow_zone = ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ); + +% N( ~in_shadow_zone, : ) = - N( ~in_shadow_zone, : ); + +%% Transfer function +H_dir = repmat( 1 ./ r_dir, 1, numel(frequencies) ); +H_diffr( :, ~in_shadow_zone' ) = zeros( numel( frequencies ), sum( ~in_shadow_zone ) ); +% From Handbook of Acoustics page 117 eq. 4.13 +H_diffr( :, in_shadow_zone' ) = ( ( 10^(5/20) * sqrt( 2 * pi * N(in_shadow_zone, :) ) ./ tanh( sqrt( 2*pi*N(in_shadow_zone, :) ) ) ).^(-1) .* H_dir(in_shadow_zone, :) )'; + +% att.freqData = 1 ./ ( 10^(5/20) * sqrt( 2*pi * N ) ./ tanh( sqrt( 2*pi * N ) ) )'; + +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_approx.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_approx.m new file mode 100644 index 0000000..1fdfb17 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_approx.m @@ -0,0 +1,76 @@ +function H_diffr = ita_diffraction_maekawa_approx( wedge, source_pos, receiver_pos, frequencies, speed_of_sound, transition_const ) +% Calculates the attenuation filter(s) at a diffraction wedge for source +% and receiver location(s) at given frequencies with. For purpose of +% continuity at shadow boundary an interpolation with exponential function +% is used. +% +% wedge: diffracting wedge (itaInfiniteWedge or derived class +% itaFiniteWedge) +% source_pos: position of the source +% receiver_pos: position of the receiver + + +%% Assertions +assert( isa( wedge, 'itaInfiniteWedge' ) ) +if nargin < 6 + transition_const = 0.2; +end +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +dim_f = size( frequencies ); +if dim_src(2) ~= 3 + if dim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3') + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Receiver point(s) must be of dimension 3') + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +if dim_f(1) ~= 1 + if dim_f(2) ~= 1 + error( 'Invalid frequency. Use row or column vector' ); + end + frequencies = frequencies'; +end + + +%% Variables +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +Src_Apex_Dir = ( Apex_Point - source_pos ) ./ Norm( Apex_Point - source_pos ); +Apex_Rcv_Dir = ( receiver_pos - Apex_Point ) ./ Norm( receiver_pos - Apex_Point ); +detour = Norm( Apex_Point - source_pos ) + Norm( receiver_pos - Apex_Point ) - Norm( receiver_pos - source_pos ); +c = speed_of_sound; +lambda = c ./ frequencies; +N = 2 * detour ./ lambda; % Fresnel number N +r_dir = Norm( receiver_pos - source_pos ); + +in_shadow_zone = ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ); + +phi = acos( dot( Apex_Rcv_Dir( in_shadow_zone, : ), Src_Apex_Dir( in_shadow_zone, : ), 2 ) ); % angle between receiver and shadow boundary +phi( phi > pi/4 ) = pi/4; +phi_0 = transition_const; + + +c_norm = 10^(5/20); +c_approx = repmat( 1 + ( c_norm - 1 ) .* exp( -phi/phi_0 ), 1, numel(frequencies) ); + +%% Transfer function +% From Handbook of Acoustics page 117 eq. 4.13 +H_dir = repmat( 1 ./ r_dir, 1, numel(frequencies) ); +H_diffr( :, ~in_shadow_zone ) = zeros( numel( frequencies ), sum( ~in_shadow_zone ) ); +H_diffr( :, in_shadow_zone ) = c_approx' .* ( ( 10^(5/20) * sqrt( 2 * pi * N(in_shadow_zone, :) ) ./ tanh( sqrt( 2*pi*N(in_shadow_zone, :) ) ) ).^(-1) .* H_dir(in_shadow_zone, :) )'; + +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_total_field_plots.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_total_field_plots.m new file mode 100644 index 0000000..8920150 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_maekawa_total_field_plots.m @@ -0,0 +1,205 @@ +%% Config +% screen +n1 = [1 0 0]; +n2 = [-1 0 0]; + +% % wedge +% n1 = [1, 1, 0]; +% n2 = [-1, 1, 0]; + +apex_dir = [0 0 1]; +loc = [0 0 2]; +len = 10; + +% screen +w = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc); +w.aperture_direction = apex_dir; + +n3 = [ 1, 0.00001, 0 ]; +n4 = [-1, 0.00001, 0 ]; +finw = itaFiniteWedge( n3 / norm(n3), n4 / norm(n4), loc, len ); +% finw.aperture_direction = apex_dir; +% finw.aperture_end_point = (len); + +% source and receiver +src = 3/sqrt(1) * [-1 0 0]; +rcv_start_pos = 3/sqrt(2) * [ 1 1 0]; +rcv_end_pos = 3/sqrt(2) * [ 1 -1 0]; +apex_point = w.get_aperture_point(src, rcv_start_pos); + +% parameters +delta = 0.05; +c = 344; % Speed of sound +angle_resolution = 100; +freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]; + +f_s = 44100; +n = 1024; +R_0 = norm( apex_point - src ) + norm( rcv_start_pos - apex_point ); +tau_0 = R_0 / c; +tau = tau_0 : 1/f_s : tau_0 + (n - 1) * 1/f_s; + +%% Variables + +ref_face = w.point_facing_main_side( src ); +alpha_start = w.get_angle_from_point_to_wedge_face(rcv_start_pos, ref_face); +alpha_end = w.get_angle_from_point_to_wedge_face(rcv_end_pos, ref_face); +alpha_d = linspace( alpha_start, alpha_end, angle_resolution ); + + + +% Set different receiver positions rotated around the aperture +rcv_positions = ita_align_points_around_aperture( w, rcv_start_pos, alpha_d, apex_point, ref_face ); + +detour_vec = norm( apex_point - src ) + Norm( rcv_positions - apex_point ) - Norm( rcv_positions - src ); + +% Fresnel Number +lambda = c ./ freq; +N_vec = 2 * detour_vec ./ lambda; + +%% Calculations +% total wave field +k = 2 * pi * freq ./ c; +R_dir = repmat( Norm( rcv_positions - src ), 1, numel(freq) ); +E_dir = 1 ./ R_dir .* exp( -1i .* k .* R_dir ); + +in_shadow_zone = ita_diffraction_shadow_zone( w, src, rcv_positions ); + + +% Maekawa diffraction attenuation +att_sum_maekawa = itaResult; +att_sum_maekawa.freqVector = freq; +att_sum_maekawa.freqData = ita_diffraction_maekawa( w, src, rcv_positions, freq, c ); +att_sum_maekawa.freqData( :, ~in_shadow_zone' ) = att_sum_maekawa.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_maekawa.freqData = att_sum_maekawa.freqData ./ E_dir'; + + +% Maekawa approximated diffraction attenuation +att_sum_maekawa_approx = itaResult; +att_sum_maekawa_approx.freqVector = freq; +att_sum_maekawa_approx.freqData = ita_diffraction_maekawa_approx( w, src, rcv_positions, freq, c ); +att_sum_maekawa_approx.freqData( :, ~in_shadow_zone ) = att_sum_maekawa_approx.freqData( :, ~in_shadow_zone ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_maekawa_approx.freqData = att_sum_maekawa_approx.freqData ./ E_dir'; + + +% UTD total wave field +att_sum_utd = itaResult; +att_sum_utd.freqVector = freq; +att_sum_utd.freqData = ita_diffraction_utd( w, src, rcv_positions, freq, c ); +att_sum_utd.freqData( :, ~in_shadow_zone' ) = att_sum_utd.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_utd.freqData = att_sum_utd.freqData ./ E_dir'; + + +% UTD approximated wave field +att_sum_utd_approx = itaResult; +att_sum_utd_approx.freqVector = freq; +att_sum_utd_approx.freqData = ita_diffraction_utd_approximated( w, src, rcv_positions, freq, c ); +att_sum_utd_approx.freqData( :, ~in_shadow_zone' ) = att_sum_utd_approx.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_utd_approx.freqData = att_sum_utd_approx.freqData ./ E_dir'; + + +% BTMS total field +att_btms = itaAudio; +att_btms.signalType = 'energy'; +att_btms.samplingRate = f_s; +att_btms.nSamples = n; + +att_sum_btms = itaAudio; +for j = 1 : numel( alpha_d ) + att_btms.timeData = ita_diffraction_btm_finite_wedge( finw, src, rcv_positions( j, : ), tau, c, true ); + att_sum_btms = ita_merge( att_sum_btms, att_btms ); +end +f_btms = att_sum_btms.freqVector'; +f_plot = f_btms( 6 : 50 : end ); +k2 = 2 * pi * f_btms ./ c; +R_dir2 = repmat( Norm( rcv_positions - src ), 1, numel(f_btms) ); +E_dir2 = 1 ./ R_dir2 .* exp( -1i .* k2 .* R_dir2 ); + +att_sum_btms.freqData( :, ~in_shadow_zone ) = att_sum_btms.freqData( :, ~in_shadow_zone ) + ( E_dir2( ~in_shadow_zone, : ) )'; +att_sum_btms.freqData = att_sum_btms.freqData ./ E_dir2'; + + +% maekawa with btms freq +att_sum_maekawa_comp = itaResult; +att_sum_maekawa_comp.freqData = ita_diffraction_maekawa( w, src, rcv_positions, f_btms, c ); +att_sum_maekawa_comp.freqData( :, ~in_shadow_zone' ) = att_sum_maekawa_comp.freqData( :, ~in_shadow_zone' ) + ( E_dir2( ~in_shadow_zone, : ) )'; +att_sum_maekawa_comp.freqData = att_sum_maekawa_comp.freqData ./ E_dir2'; + +% error calculation +att_sum_error = itaResult; +att_sum_error.freqVector = freq; +att_sum_error.freqData = ( att_sum_maekawa_comp.freqData ./ att_sum_btms.freqData ).^(-1); + + +% error +att_sum_error_approx = itaResult; +att_sum_error_approx.freqVector = freq; +att_sum_error_approx.freqData = att_sum_maekawa_approx.freqData ./ att_sum_utd_approx.freqData; + +%% Comparison to UTD plot +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); +% subplot( 2, 2, 1 ); +% plot( rad2deg( alpha_d ), att_sum_maekawa.freqData_dB' ) +% title( 'Maekawa diffraction' ) +% legend( [num2str( freq' ), repmat(' Hz', numel(freq), 1)], 'Location', 'southwest' ) +% xlabel( 'theta_R [°]' ) +% ylabel( 'p_{total} [dB]' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; +% +% subplot( 2, 2, 2 ); +% plot( rad2deg( alpha_d ), att_sum_maekawa_approx.freqData_dB' ) +% title( 'Maekawa approximation' ) +% legend( [num2str( freq' ), repmat(' Hz', numel(freq), 1)], 'Location', 'southwest' ) +% xlabel( 'theta_R [°]' ) +% ylabel( 'p_{total} [dB]' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; + +res_btms = att_sum_btms.freqData_dB'; +subplot( 2, 2, 1 ); +plot( rad2deg( alpha_d ), res_btms( :, 6 : 50 : end ) ) +title( 'diffraction by screen modeled with BTMS method' ) +legend( [num2str( round( f_plot' ) ), repmat(' Hz', numel(f_plot), 1)], 'Location', 'southwest' ) +xlabel( 'theta_R [°]' ) +ylabel( 'p_{total} [dB]' ) +ylim( [-35, 10] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; + +res_error = att_sum_error.freqData_dB'; +subplot( 2, 2, 2 ); +plot( rad2deg( alpha_d ), res_error( :, 6 : 50 : end ) ) +title( 'deviation BTMS from Maekawa method' ) +legend( [num2str( round( f_plot' ) ), repmat(' Hz', numel(f_plot), 1)], 'Location', 'southwest' ) +xlabel( 'theta_R [°]' ) +ylabel( 'p_{total} [dB]' ) +ylim( [-35, 10] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; +% +% subplot( 3, 2, 5 ); +% plot( rad2deg( alpha_d ), att_sum_utd_approx.freqData_dB' ) +% title( 'UTD total wave field plot' ) +% legend( num2str( freq' ), 'Location', 'southwest' ) +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +% ylabel( 'dB SPL' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; +% +% subplot( 3, 2, 6 ); +% plot( rad2deg( alpha_d ), att_sum_error_approx.freqData_dB' ) +% title( 'UTD total wave field plot' ) +% legend( num2str( freq' ), 'Location', 'southwest' ) +% xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +% ylabel( 'dB SPL' ) +% ylim( [-35, 10] ); +% xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% grid on; + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_shadow_zone.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_shadow_zone.m new file mode 100644 index 0000000..abe65c2 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_shadow_zone.m @@ -0,0 +1,85 @@ +function in_shadow = ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ) +%ITA_DIFFRACTION_SHADOW_ZONE Returns true if receiver is in shadow zone of +%source + +%% assertions +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +if dim_src(2) ~= 3 + if dim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3') + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Receiver point(s) must be of dimension 3') + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +switch dim_src(1) >= dim_rcv(1) + case true + dim_n = dim_src(1); + S = source_pos; + R = repmat( receiver_pos, dim_n, 1 ); + case false + dim_n = dim_rcv(1); + S = repmat( source_pos, dim_n, 1 ); + R = receiver_pos; +end + +%% function begin +in_shadow = false( dim_n, 1 ); +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +eps = wedge.set_get_geo_eps; + +ref_face = wedge.point_facing_main_side( S ); +N1 = zeros( dim_n, 3 ); +N2 = zeros( dim_n, 3 ); + +N1( ref_face, : ) = repmat( wedge.main_face_normal, sum(ref_face), 1 ); +N2( ref_face, : ) = repmat( wedge.opposite_face_normal, sum(ref_face), 1 ); +Apex_Dir( ref_face, : ) = repmat( wedge.aperture_direction, sum(ref_face), 1 ); + +N1( ~ref_face, : ) = repmat( wedge.opposite_face_normal, sum(~ref_face), 1 ); +N2( ~ref_face, : ) = repmat( wedge.main_face_normal, sum(~ref_face), 1 ); +Apex_Dir( ~ref_face, : ) = repmat( -wedge.aperture_direction, sum(~ref_face), 1 ); + + +% Distances of source from wedge faces +dist1_src = dot( S - Apex_Point, N1, 2 ); +dist2_src = dot( S - Apex_Point, N2, 2 ); + +if any( ( dist1_src < -eps ) & ( dist2_src < -eps ) ) + error('Invalid source location!') +end + +% Use of a auxiliary shadow plane which describes the border of the +% shadow zone +source_apex_direction = ( Apex_Point - S ) ./ Norm( Apex_Point - S ) ; + +% Check right orientation of normal vector of shadow plane by consideration +% of source facing wedge side +shadow_plane_normal = zeros( dim_n, 3 ); +mask = dist1_src >= -eps; +shadow_plane_normal(mask, :) = cross( source_apex_direction(mask, :), Apex_Dir(mask, :), 2 ) ./ Norm( cross( source_apex_direction(mask, :), Apex_Dir(mask, :), 2 ) ); % ./ sqrt( sum( (cross( source_apex_direction(mask, :), Apex_Dir(mask, :), 2 )).^2, 2 ) ); +shadow_plane_normal(~mask, :) = cross( Apex_Dir(~mask, :), source_apex_direction(~mask, :), 2 ) ./ Norm( cross( Apex_Dir(~mask, :), source_apex_direction(~mask, :), 2 ) ); % ./ sqrt( sum( (cross( Apex_Dir(~mask, :), source_apex_direction(~mask, :), 2 )).^2, 2 ) ); + +% Distance of receiver from shadow plane +dist_rcv = dot( R - Apex_Point, shadow_plane_normal, 2 ); + +% Checking for receiver position in shadow zone and considering special +% cases +case1 = ( ( dist1_src >= -eps ) & ( dist2_src >= -eps ) ); +case2 = dot( R - Apex_Point, N1, 2 ) >= -eps; +in_shadow( ~(case1 | case2) ) = dist_rcv( ~(case1 | case2) ) < -eps; +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_svensson_plots.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_svensson_plots.m new file mode 100644 index 0000000..1c34fec --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_svensson_plots.m @@ -0,0 +1,96 @@ +%% Config +n1 = [1 1 0]; +n2 = [-1 1 0]; +% apex_dir = [0 0 1]; +loc = [0 0 -3]; +len = 8; +src = 0.2/sqrt(1) * [-1 0 0]; +rcv_start_pos = 1/sqrt(2) * [ -1, -1, 0]; +infw = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc); +finw = itaFiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc, len); +apex_dir = infw.aperture_direction; +% w.aperture_direction = apex_dir; +apex_point = infw.get_aperture_point(src, rcv_start_pos); +delta = 0.05; +c = 343.21; % Speed of sound + +%freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]; +f_sampling = 384000; +filter_length = 1024; +angle_resolution = 2000; + +%% Calculations +ref_face = infw.point_facing_main_side( src ); +alpha_d = linspace( infw.get_angle_from_point_to_wedge_face(rcv_start_pos, ref_face), infw.opening_angle, angle_resolution ); +rcv_positions = ita_align_points_around_aperture( infw, rcv_start_pos, alpha_d, apex_point, ref_face ); + +h1_0 = zeros( 1, numel(alpha_d) ); +h2_0 = zeros( 1, numel(alpha_d) ); +h3_0 = zeros( 1, numel(alpha_d) ); +h4_0 = zeros( 1, numel(alpha_d) ); +R_0 = norm( src ) + norm( rcv_start_pos ); + +for j = 1 : numel( alpha_d ) + H = ita_diffraction_btm_finite_wedge( finw, src, rcv_positions( j, : ), f_sampling, filter_length ); + h1_0(j) = H( 1 ); + h2_0(j) = H( 2 ); + h3_0(j) = H( 3 ); + h4_0(j) = H( 4 ); +end + +% % %% EDB +% % % variables +% % fs = sampling_rate; +% % closwedang = finw.wedge_angle; +% % rs = sqrt( src(1)^2 + src(2)^2 ); +% % thetas = infw.get_angle_from_point_to_wedge_face( src, ref ); +% % zs = 0; +% % rr = sqrt( rcv(1)^2 + rcv(2)^2 ); +% % thetar = infw.get_angle_from_point_to_wedge_face( rcv, ref ); +% % zr = rcv(3) - src(3); +% % zw = [ loc(3) - src(3), len - src(3) ]; +% % Method = 'New'; +% % +% % % filter +% % [ir,initdelay,singularterm] = EDB2wedge1st_int( fs, closwedang, rs, thetas, zs, rr, thetar, zr, zw, Method ); +% % res1 = ir(ir ~= 0); + +%% Plots +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); + +subplot( 2, 2, 1 ); +plot( rad2deg( alpha_d(2:end) ), h1_0(2:end) .* R_0 ); +title( 'h1(n_0)' ); +xlabel( 'theta_R in degree' ); +ylabel( 'h1(n_0) rel. 1/R_0' ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +% ylim( [-0.002, 0.01] ); +grid on; + +subplot( 2, 2, 2 ); +plot( rad2deg( alpha_d(2:end) ), h2_0(2:end) .* R_0 ); +title( 'h2(n_0)' ); +xlabel( 'theta_R in degree' ); +xlabel( 'theta_R in degree' ) +ylabel( 'h2(n_0) rel. 1/R_0' ); +ylim( [-0.6, 0.6] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; + +subplot( 2, 2, 3 ); +plot( rad2deg( alpha_d(2:end) ), h3_0(2:end) .* R_0 ); +title( 'h3(n_0)' ); +xlabel( 'theta_R in degree' ); +ylabel( 'h3(n_0) rel. 1/R_0' ); +% ylim( [-0.003, 0.003] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; + +subplot( 2, 2, 4 ); +plot( rad2deg( alpha_d(2:end) ), h4_0(2:end) .* R_0 ); +title( 'h4(n_0)' ); +xlabel( 'theta_R in degree' ); +ylabel( 'h4(n_0) rel. 1/R_0' ); +ylim( [-0.6, 0.6] ); +xlim( [rad2deg(alpha_d(1)), rad2deg(alpha_d(end))] ); +grid on; \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_approx_plot.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_approx_plot.m new file mode 100644 index 0000000..38e4412 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_approx_plot.m @@ -0,0 +1,64 @@ +%% Config +n1 = [-1 1 0]; +n2 = [1 1 0]; +loc = [0 0 2]; +src = 3 * [-1 0 0]; +rcv = 3 * [1/sqrt(2) 1/sqrt(2) 0]; +w = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc); +apex_point = w.get_aperture_point(src, rcv); +apex_dir = w.aperture_direction; +c = 344; % Speed of sound + +freq = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]; +alpha_d = linspace( pi, 3 * pi / 2, 500 ); + +% Set different receiver positions rotated around the aperture +rcv_positions = zeros(numel(alpha_d), 3); +rcv_positions(1, :) = rcv; + +% Coordinate transformation +n3 = apex_dir; +n2 = w.main_face_normal; +n1 = cross( n2, n3 ); + +rho = norm( rcv - apex_point ); +z = rcv(3); + +rcv_pos_cylindrical = zeros(numel(alpha_d), 3); +for i = 1 : numel(alpha_d) + rcv_pos_cylindrical(i, :) = [rho, alpha_d(i), z]; +end + +for j = 2 : numel(rcv_positions(:, 1)) + rcv_positions(j, 1) = rcv_pos_cylindrical(j, 1) * cos( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 2) = rcv_pos_cylindrical(j, 1) * sin( pi*5/4 - rcv_pos_cylindrical(j, 2) ); + rcv_positions(j, 3) = - rcv_pos_cylindrical(j, 3); +end + +%% Calculation +att_sum = itaResult; +k = 2 * pi * freq ./ c; + +for j = 1 : numel(rcv_positions(:, 1)) + r_dir = norm( rcv_positions(j, :) - src ); + E_dir = 1 / r_dir * exp( -1i .* k * r_dir ); + if ~ita_diffraction_shadow_zone( w, src, rcv_positions(j, :) ) + E_approx = itaResult; + E_approx.freqVector = freq; + E_approx.freqData = ones( numel( freq ), 1 ); + else + E_approx = ita_diffraction_normalized_utd(w, src, rcv_positions(j, :), freq ) ./ E_dir; + end + att_sum = ita_merge( att_sum, E_approx ); +end + + +%% Tsingos paper plot +figure +plot( rad2deg( alpha_d ), att_sum.freqData_dB' ) +title( 'Tsingos et al.: UTD total wave field plot (Figure 6b)' ) +legend( num2str( freq' ) ) +xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +ylabel( 'dB SPL' ) +ylim( [-35, 10] ); +grid on; \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_total_field_plots.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_total_field_plots.m new file mode 100644 index 0000000..9f0ca1c --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_tsingos_utd_total_field_plots.m @@ -0,0 +1,91 @@ +%% 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]; +w = itaInfiniteWedge(n1 / norm( n1 ), n2 / norm( n2 ), loc); +apex_point = w.get_aperture_point(src, rcv_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]; +ref_face = w.point_facing_main_side( src ); +alpha_d = linspace( w.get_angle_from_point_to_wedge_face(rcv_start_pos, ref_face), w.opening_angle, 1200 ); + +% Set different receiver positions rotated around the aperture +rcv_positions = ita_align_points_around_aperture( w, rcv_start_pos, alpha_d, apex_point, ref_face ); + +%% Calculations +k = 2 * pi * freq ./ c; +R_dir = repmat( sqrt( sum( ( rcv_positions - src ).^2, 2 ) ), 1, numel(freq) ); +E_dir = 1 ./ R_dir .* exp( -1i .* k .* R_dir ); + +in_shadow_zone = ita_diffraction_shadow_zone( w, src, rcv_positions ); + +% UTD total wave field +att_sum = itaResult; +att_sum.freqVector = freq; +att_sum.freqData = ita_diffraction_utd( w, src, rcv_positions, freq', c ); +att_sum.freqData( :, ~in_shadow_zone' ) = att_sum.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum.freqData = att_sum.freqData ./ E_dir'; + +% UTD combined with Approximation by Tsingos et. al. +att_sum_approx = itaResult; +att_sum_approx.freqVector = freq; +att_sum_approx.freqData = ita_diffraction_utd_approximated( w, src, rcv_positions, freq, c ); +att_sum_approx.freqData( :, ~in_shadow_zone' ) = att_sum_approx.freqData( :, ~in_shadow_zone' ) + ( E_dir( ~in_shadow_zone, : ) )'; +att_sum_approx.freqData = att_sum_approx.freqData ./ E_dir'; + +% Error caused by approximation +att_sum_error = itaResult; +att_sum_error.freqVector = freq; +att_sum_error.freqData = att_sum_approx.freqData ./ att_sum.freqData; + +att_sum_diffr = itaResult; +att_sum_diffr.freqVector = freq; +att_sum_diffr.freqData = ita_diffraction_utd( w, src, rcv_positions, freq, c ); +att_sum_diffr.freqData = att_sum_diffr.freqData ./ E_dir'; + + +%% Tsingos paper plot +figure( 'units', 'normalized', 'outerposition', [0 0 1 1] ); +subplot( 2, 2, 1 ); +plot( rad2deg( alpha_d ), att_sum.freqData_dB' ) +title( 'total wave field with UTD' ) +legend( num2str( freq' ), 'Location', 'southwest' ) +xlabel( 'theta_R in degree (shadow boundary at 225deg)' ) +ylabel( 'dB SPL' ) +xlim( [alpha_d(1), w.opening_angle_deg] ); +ylim( [-35, 10] ); +grid on; + +subplot( 2, 2, 2 ); +plot( rad2deg( alpha_d ), att_sum_approx.freqData_dB' ) +title( 'Tsingos et al.: UTD total wave field plot with Approximation (Figure 6b)' ) +legend( num2str( freq' ), 'Location', 'southwest' ) +xlabel( 'alpha_d in degree (shadow boundary at 225deg)' ) +ylabel( 'dB SPL' ) +xlim( [alpha_d(1), w.opening_angle_deg] ); +ylim( [-35, 10] ); +grid on; + +subplot( 2, 2, 3 ); +plot( rad2deg( alpha_d ), att_sum_error.freqData_dB' ) +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( [alpha_d(1), w.opening_angle_deg] ); +ylim( [-35, 10] ); +grid on; + +subplot( 2, 2, 4 ); +plot( rad2deg( alpha_d ), att_sum_diffr.freqData_dB' ) +title( 'Diffracted field by UTD' ) +legend( [num2str( freq' ), repmat(' Hz', numel(freq),1)], 'Location', 'southeast' ) +xlabel( 'theta_R [°]' ) +ylabel( 'dB SPL' ) +xlim( [alpha_d(1), w.opening_angle_deg] ); +grid on; \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_utd.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd.m new file mode 100644 index 0000000..d21f19a --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd.m @@ -0,0 +1,203 @@ +function diffr_field = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequency_vec, speed_of_sound ) +%ITA_DIFFRACTION_UTD Calculates the diffraction filter based on uniform +%theory of diffraction (with Kawai approximation) +% +% Literature: +% [1] Tsingos, Funkhouser et al. - Modeling Acoustics in Virtual Environments using the Uniform Theory of Diffraction +% [2] Kouyoumjian and Pathak - A Uniform Geometrical Theory of Diffraction for an Edge in a Perfectly Conducting Surface +% +% Example: +% att = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequenc_vec ) +% + +%% Assertions +dim_freq = size( frequency_vec ); +dim_src = size( source_pos ); +dim_rcv = size( receiver_pos ); +if dim_freq(1) ~= 1 + if dim_freq(2) ~= 1 + error( 'Invalid frequency vector' ); + end + frequency_vec = frequency_vec'; +end +if dim_src(2) ~= 3 + if dim_src(1) ~= 3 + error( 'Source point(s) must be of dimension 3') + end + source_pos = source_pos'; + dim_src = size( source_pos ); +end +if dim_src(1) == 0 + error( 'no source found' ); +end +if dim_rcv(2) ~= 3 + if dim_rcv(1) ~= 3 + error( 'Receiver point(s) must be of dimension 3') + end + receiver_pos = receiver_pos'; + dim_rcv = size( receiver_pos ); +end +if dim_rcv(1) == 0 + error( 'no receiver found' ); +end +if dim_src(1) ~= 1 && dim_rcv(1) ~= 1 && dim_src(1) ~= dim_rcv(1) + error( 'Number of receiver and source positions do not match' ) +end +switch dim_src(1) >= dim_rcv(1) + case true + dim_n = dim_src(1); + case false + dim_n = dim_rcv(1); +end + +%% Variables +if dim_src(1) == dim_n + S = source_pos; + if dim_rcv(1) == dim_n + R = receiver_pos; + else + R = repmat( receiver_pos, dim_n, 1 ); + end +else + R = receiver_pos; + if dim_src(1) == dim_n + S = source_pos; + else + S = repmat( source_pos, dim_n, 1 ); + end +end +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +rho = Norm( Apex_Point - S ); % Distance of source to aperture point +r = Norm( R - Apex_Point ); % Distance of receiver to aperture point +c = speed_of_sound; + +face = wedge.point_facing_main_side( S ); +alpha_i = wedge.get_angle_from_point_to_wedge_face( S, face ); +alpha_d = wedge.get_angle_from_point_to_wedge_face( R, face ); +theta_i = wedge.get_angle_from_point_to_aperture( S, Apex_Point ); + +n = wedge.opening_angle / pi; % Variable dependend on opening angle of the wedge + +lambda = c ./ frequency_vec; % Wavelength +k = (2 * pi) ./ lambda; % Wavenumber + + +% Diffraction coefficient D +assert( all( rho + r ~= 0 ) && all( r ~= 0 ) ); +A = repmat( sqrt( rho ./ ( r .* ( rho + r ) ) ), 1, numel( frequency_vec ) ); +L = repmat( ( ( rho .* r ) ./ ( rho + r ) ) .* ( sin( theta_i ) ).^2, 1, numel( frequency_vec ) ); +H_i = 1 ./ rho .* exp( -1i * k .* rho ); % Consideration of transfer path from source to apex point + +D_factor = -exp( -1i * pi / 4 ) ./ ( 2 * n * sqrt( 2* pi * k ) .* sin( theta_i ) ); + +Cot1 = repmat( cot( ( pi + ( alpha_d - alpha_i ) ) ./ ( 2 * n ) ), 1, numel( frequency_vec ) ); +Cot2 = repmat( cot( ( pi - ( alpha_d - alpha_i ) ) ./ ( 2 * n ) ), 1, numel( frequency_vec ) ); +Cot3 = repmat( cot( ( pi + ( alpha_d + alpha_i ) ) ./ ( 2 * n ) ), 1, numel( frequency_vec ) ); +Cot4 = repmat( cot( ( pi - ( alpha_d + alpha_i ) ) ./ ( 2 * n ) ), 1, numel( frequency_vec ) ); + +a1 = 2 * ( cos( ( 2 * pi * n * N_p( n, alpha_d - alpha_i ) - ( alpha_d - alpha_i ) ) / 2 ) ).^2; +a2 = 2 * ( cos( ( 2 * pi * n * N_n( n, alpha_d - alpha_i ) - ( alpha_d - alpha_i ) ) / 2 ) ).^2; +a3 = 2 * ( cos( ( 2 * pi * n * N_p( n, alpha_d + alpha_i ) - ( alpha_d + alpha_i ) ) / 2 ) ).^2; +a4 = 2 * ( cos( ( 2 * pi * n * N_n( n, alpha_d + alpha_i ) - ( alpha_d + alpha_i ) ) / 2 ) ).^2; + +F1 = kawai_approx_fresnel( k .* L .* a1 ); +F2 = kawai_approx_fresnel( k .* L .* a2 ); +F3 = kawai_approx_fresnel( k .* L .* a3 ); +F4 = kawai_approx_fresnel( k .* L .* a4 ); + +% Avoid eventual singularities of the cot terms at the shadow or reflection boundary with a approximation by +% Kouyoumjian and Pathak +mask1 = ( alpha_d - alpha_i ) - 2 * pi * n * N_p( n, alpha_d - alpha_i ) + pi == 0; +mask2 = - ( alpha_d - alpha_i ) + 2 * pi * n * N_n( n, alpha_d - alpha_i ) + pi == 0; +mask3 = ( alpha_d + alpha_i ) - 2 * pi * n * N_p( n, alpha_d + alpha_i ) + pi == 0; +mask4 = - ( alpha_d + alpha_i ) + 2 * pi * n * N_n( n, alpha_d + alpha_i ) + pi == 0; + +singularities = [ any( mask1 ~= 0 ), any( mask2 ~= 0 ), any( mask3 ~= 0 ), any( mask4 ~= 0 ) ]; + + +if any( singularities ) + if singularities(1) + eps1 = ( alpha_d(mask1) - alpha_i(mask1) ) - 2 * pi * n * N_p( n, alpha_d(mask1) - alpha_i(mask1) ) + pi; + else + eps1 = 0; + end + if singularities(2) + eps2 = - ( alpha_d(mask2) - alpha_i(mask2) ) + 2 * pi * n * N_n( n, alpha_d(mask2) - alpha_i(mask2) ) + pi; + else + eps2 = 0; + end + if singularities(3) + eps3 = ( alpha_d(mask3) + alpha_i(mask3) ) - 2 * pi * n * N_p( n, alpha_d(mask3) + alpha_i(mask3) ) + pi; + else + eps3 = 0; + end + if singularities(4) + eps4 = - ( alpha_d(mask4) + alpha_i(mask4) ) + 2 * pi * n * N_n( n, alpha_d(mask4) + alpha_i(mask4) ) + pi; + else + eps4 = 0; + end + + term1(mask1, :) = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k .* L(mask1, :) ) .* sgn( eps1 ) - 2 .* k .* L(mask1, :) .* eps1 * exp( 1i * pi/4 ) ); + term2(mask2, :) = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k .* L(mask2, :) ) .* sgn( eps2 ) - 2 .* k .* L(mask2, :) .* eps2 * exp( 1i * pi/4 ) ); + term3(mask3, :) = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k .* L(mask3, :) ) .* sgn( eps3 ) - 2 .* k .* L(mask3, :) .* eps3 * exp( 1i * pi/4 ) ); + term4(mask4, :) = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k .* L(mask4, :) ) .* sgn( eps4 ) - 2 .* k .* L(mask4, :) .* eps4 * exp( 1i * pi/4 ) ); +end + +term1(~mask1, :) = Cot1(~mask1, :) .* F1(~mask1, :); +term2(~mask2, :) = Cot2(~mask2, :) .* F2(~mask2, :); +term3(~mask3, :) = Cot3(~mask3, :) .* F3(~mask3, :); +term4(~mask4, :) = Cot4(~mask4, :) .* F4(~mask4, :); + +if wedge.is_boundary_condition_hard + s = 1; +else + s = -1; +end + +D = D_factor .* ( term1 + term2 + s * ( term3 + term4 ) ); + +% Combined diffracted sound field filter at receiver +diffr_field = ( H_i .* D .* A .* exp( -1i .* k .* r ) )'; + +end + + +%% Auxiliary functions +% euclidean norm row wise +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end + +% N+ function +function N = N_p( n, beta ) + N = zeros( numel( beta ), 1 ); + N( beta > pi * ( 1 - n ) ) = 1; +end + +% N- function +function N = N_n( n, beta ) + N = zeros( numel( beta ), 1 ); + N( beta < pi * ( 1 - n ) ) = -1; + N( beta > pi * ( 1 + n ) ) = 1; +end + +% signum function +function res = sgn(x) + if all( size(x) == 0 ) + res = 1; + return; + end + res = ones( size(x) ); + res( x <= 0 ) = -1; +end + +function Y = kawai_approx_fresnel( X ) + if any( X < 0 ) + error( 'No negative values for Kawai approximation of Fresnel integral allowed' ) + end + X_s_idx = X < 0.8; + X_geq_idx = ( X >= 0.8 ); + Y = zeros( size(X) ); + Y( X_s_idx ) = sqrt( pi * X( X_s_idx ) ) .* ( 1 - sqrt( X( X_s_idx ) ) ./ ( 0.7 * sqrt( X( X_s_idx ) ) + 1.2 ) ) .* exp( 1i * pi/4 * ( 1 - sqrt( X( X_s_idx ) ./ ( X( X_s_idx ) + 1.4 ) ) ) ); + Y( X_geq_idx ) = ( 1 - 0.8 ./ ( X( X_geq_idx ) + 1.25 ) .^ 2 ) .* exp( 1i * pi/4 * ( 1 - sqrt( X( X_geq_idx ) ./ ( X( X_geq_idx ) + 1.4 ) ) ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_approximated.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_approximated.m new file mode 100644 index 0000000..06f234d --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_approximated.m @@ -0,0 +1,57 @@ +function diffr_field = ita_diffraction_utd_approximated( wedge, source_pos, receiver_pos, frequency_vec, speed_of_sound, transition_const ) +%ITA_DIFFRACTION_UTD Calculates the diffraction filter based on uniform +%theory off diffraction (with Kawai approximation) only in shadow regions. +%To preserve continuity of the total sound field in e.g. shadow boundaries, +%a normalization is taken to account (approximation by Tsingos et. al.) +% +% Literature: +% [1] Tsingos, Funkhouser et al. - Modeling Acoustics in Virtual Environments using the Uniform Theory of Diffraction +% +% Example: +% att = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequency_vec ) +% +%% Assertions +if nargin < 6 + transition_const = 0.1; +end + +%% Variables +Apex_Point = wedge.get_aperture_point( source_pos, receiver_pos ); +Src_Apex_Dir = ( Apex_Point - source_pos ) ./ Norm( Apex_Point - source_pos ); +Apex_Rcv_Dir = ( receiver_pos - Apex_Point ) ./ Norm( receiver_pos - Apex_Point ); +rho = Norm( Apex_Point - source_pos ); % Distance Source to Apex point +r = Norm( receiver_pos - Apex_Point ); % Distance Apex point to Receiver +receiver_SB = source_pos + Src_Apex_Dir .* ( r + rho ); % Virtual position of receiver at shadow boundary +c = speed_of_sound; +k_vec = 2 * pi * frequency_vec ./ c; % Wavenumber + +in_shadow_zone = ita_diffraction_shadow_zone( wedge, source_pos, receiver_pos ); + +phi = repmat( acos( dot( Apex_Rcv_Dir( in_shadow_zone, : ), Src_Apex_Dir( in_shadow_zone, : ), 2 ) ), 1, numel( frequency_vec ) ); % angle between receiver and shadow boundary +phi( phi > pi/4 ) = pi/4; +phi_0 = transition_const; + +% Incident field at shadow boundary +E_incident_SB = ( 1 ./ Norm( receiver_SB - source_pos ) .* exp( -1i .* k_vec .* ( r + rho ) ) )'; +% Diffracted field at shadow boundary +E_diff_SB = ita_diffraction_utd( wedge, source_pos, receiver_SB, frequency_vec, c ); + +%% Filter Calculation +% Normalization factor +C_norm = E_incident_SB( :, in_shadow_zone ) ./ E_diff_SB( :, in_shadow_zone ); +% Considering Interpolation to standard UTD form +C_total = 1 + ( C_norm - 1 ) .* ( exp( -phi/phi_0 ) )'; + +% if any( any( receiver_pos( in_shadow_zone, : ) ) ) + att_shadow_zone = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequency_vec, c ); + att_shadow_zone = C_total .* att_shadow_zone; % Normalization and interpolation + diffr_field( :, in_shadow_zone ) = att_shadow_zone; +% end + +diffr_field( :, ~in_shadow_zone ) = zeros( numel( frequency_vec ), sum( ~in_shadow_zone ) ); + +end + +function res = Norm( A ) + res = sqrt( sum( A.^2, 2 ) ); +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_cot_singularity_approx_test.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_cot_singularity_approx_test.m new file mode 100644 index 0000000..41e8061 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_cot_singularity_approx_test.m @@ -0,0 +1,147 @@ +% Config +n1 = [-1 1 0]; +n2 = [1 1 0]; +loc = [0 0 0]; +wedge = itaInfiniteWedge(n1, n2, loc); +source_pos = [-1 0 0]; +receiver_pos = [1 -0.0001 0]; +frequency_vec = ita_ANSI_center_frequencies; + +% Variables +apex_point = wedge.get_aperture_point( source_pos, receiver_pos ); +apex_dir = wedge.aperture_direction; +source_apex_direction = ( apex_point - source_pos ) / norm( apex_point - source_pos ); +rho = norm( apex_point - source_pos ); % Distance of source to aperture point +r = norm( receiver_pos - apex_point ); % Distance of receiver to aperture point +attenuation = itaResult; +attenuation.freqVector = frequency_vec; +attenuation.freqData = ones( numel( frequency_vec ), 1 ); +c = 344; % Speed of Sound + +if size( frequency_vec, 1 ) ~= 1 && size( frequency_vec, 2 ) ~= 1 + error( 'Invalid frequency vector' ); +end + +% Coordinate transformation with aperture point as origin +z = apex_dir; +y = wedge.main_face_normal; +x = cross( y, z ); % direction of main face of the wedge + +% Calculate angle between incedent ray from source to aperture point and source facing wedge +% side +u_i = dot( source_pos - apex_point, x ); % coordinates in new coordinate system +v_i = dot( source_pos - apex_point, y ); +alpha_i = atan2( v_i, u_i ); + +if alpha_i < 0 + alpha_i = alpha_i + pi; +end + + +% Calculate angle between ray from aperture point to receiver and receiver facing wedge +% side +u_d = dot( receiver_pos - apex_point, x ); % coordinates in new coordinate system +v_d = dot( receiver_pos - apex_point, y ); +alpha_d = atan2( v_d, u_d ); + +if alpha_d < 0 + alpha_d = alpha_d + 2*pi; +end + +% Angle between incedent ray from source to aperture point and aperture +% direction +theta_i = acos( dot( source_apex_direction, apex_dir ) ); +if theta_i > pi/2 + theta_i = pi - theta_i; +end + +n = wedge.opening_angle / pi; % Variable dependend on opening angle of the wedge + +lambda = c ./ frequency_vec; % Wavelength +k = (2 * pi) ./ lambda; % Wavenumber + + +% Diffraction coefficient D +assert( rho + r ~= 0 ) +A = sqrt( rho / ( r * ( rho + r ) ) ); +L = rho * r / ( rho + r ) * ( sin( theta_i ) )^2; + +D_factor = -exp( -1i * pi / 4 ) ./ ( 2 * n * sqrt( 2* pi * k ) * sin( theta_i ) ); + +cot1 = cot( ( pi + ( alpha_d - alpha_i ) ) ./ ( 2 * n ) ); +cot2 = cot( ( pi - ( alpha_d - alpha_i ) ) ./ ( 2 * n ) ); +cot3 = cot( ( pi + ( alpha_d + alpha_i ) ) ./ ( 2 * n ) ); +cot4 = cot( ( pi - ( alpha_d + alpha_i ) ) ./ ( 2 * n ) ); + +a1 = 2 * ( cos( ( 2 * pi * n * N_p( n, alpha_d - alpha_i ) - ( alpha_d - alpha_i ) ) / 2 ) ).^2; +a2 = 2 * ( cos( ( 2 * pi * n * N_n( n, alpha_d - alpha_i ) - ( alpha_d - alpha_i ) ) / 2 ) ).^2; +a3 = 2 * ( cos( ( 2 * pi * n * N_p( n, alpha_d + alpha_i ) - ( alpha_d + alpha_i ) ) / 2 ) ).^2; +a4 = 2 * ( cos( ( 2 * pi * n * N_n( n, alpha_d + alpha_i ) - ( alpha_d + alpha_i ) ) / 2 ) ).^2; + +F1 = ita_diffraction_kawai_approx_fresnel( k * L * a1 ); +F2 = ita_diffraction_kawai_approx_fresnel( k * L * a2 ); +F3 = ita_diffraction_kawai_approx_fresnel( k * L * a3 ); +F4 = ita_diffraction_kawai_approx_fresnel( k * L * a4 ); + + +%%% ------ approximation ------ %%% +eps1 = ( alpha_d - alpha_i ) - 2 * pi * n * N_p( n, alpha_d - alpha_i ) + pi; +eps2 = - ( alpha_d - alpha_i ) + 2 * pi * n * N_n( n, alpha_d - alpha_i ) + pi; +eps3 = ( alpha_d + alpha_i ) - 2 * pi * n * N_p( n, alpha_d + alpha_i ) + pi; +eps4 = - ( alpha_d + alpha_i ) + 2 * pi * n * N_n( n, alpha_d + alpha_i ) + pi; + +approx1 = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k * L) * sgn( eps1 ) - 2 .* k * L * eps1 * exp( 1i * pi/4 ) ); +approx2 = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k * L) * sgn( eps2 ) - 2 .* k * L * eps2 * exp( 1i * pi/4 ) ); +approx3 = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k * L) * sgn( eps3 ) - 2 .* k * L * eps3 * exp( 1i * pi/4 ) ); +approx4 = n * exp( 1i * pi/4 ) * ( sqrt( 2 * pi .* k * L) * sgn( eps4 ) - 2 .* k * L * eps4 * exp( 1i * pi/4 ) ); +%%% --------------------------- %%% + +%%% ------- Test Plots ------- %%% +term1 = cot1 .* F1; +term2 = cot2 .* F2; +term3 = cot3 .* F3; +term4 = cot4 .* F4; + + +plot_data = [term2; approx2]; +loglog( frequency_vec, abs(plot_data) ); +legend('summand', 'approx', 'location', 'south' ); +%%% -------------------------- %%% + + +D = D_factor .* ( cot1 .* F1 + cot2 .* F2 + cot3 .* F3 + cot4 .* F4 ); + + +% Combined diffracted sound field propagation filter at receiver +attenuation.freqData = ( D .* A .* exp( -1i * k * r ) )'; + + + +%% Auxiliary functions + +% N+ function +function N = N_p( n, beta ) + N = 0; + if beta > pi * ( 1 - n ) + N = 1; + end +end + +% N- function +function N = N_n( n, beta ) + N = 0; + if beta < pi * ( 1 - n ) + N = -1; + elseif beta > pi * ( 1 + n ) + N = 1; + end +end + +% signum function +function res = sgn(x) + if x > 0 + res = 1; + else + res = -1; + end +end \ No newline at end of file diff --git a/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_test.m b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_test.m new file mode 100644 index 0000000..0de8b9f --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_diffraction_utd_test.m @@ -0,0 +1,37 @@ +%% Config +n1 = [ 1 1 0 ]; +n2 = [ -1 1 0 ]; +loc = [ 0 0 -2 ]; +src = [ -5 0 0 ]; +rcv = 3 * [ 1 -0.5 0 ]; +rcv2 = [ 1 -0.99 0 ]; +w = itaInfiniteWedge( n1, n2, loc ); +freq = linspace( 20, 20000, 1000); +r_dir = norm( rcv - src ); +c = 343; +k = 2 * pi * freq ./ c; + + +%% Filter +att = ita_diffraction_utd( w, src, rcv, freq ); + +% % E_dir = ( 1 / r_dir * exp( -1i .* k * r_dir ) )'; +% % if ~ita_diffraction_shadow_zone(w, src, rcv) +% % att.freqData = att.freqData + E_dir; +% % end +% % att.freqData = att.freqData ./ E_dir; + +% plot(freq, att.freqData_dB); +semilogx(freq, att.freqData_dB); +xlim([freq(1), freq(end)]); +grid on +%att2 = ita_diffraction_utd( w, src, rcv2, ita_ANSI_center_frequencies ); + +%% +% % % a = itaAudio; +% % % a.samplingRate = 44100; +% % % a.nSamples = 44100; +% % % a.timeData( 2 ) = 1; +% % % a.pt +% % % +% % % sl = ita_generate_sweep( 'mode', 'lin' ) diff --git a/applications/VirtualAcoustics/diffraction/ita_reflection_zone.m b/applications/VirtualAcoustics/diffraction/ita_reflection_zone.m new file mode 100644 index 0000000..51eb522 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_reflection_zone.m @@ -0,0 +1,153 @@ +function reflection_zone = ita_reflection_zone( wedge, source_pos, receiver_pos ) +%ITA_DIFFRACTION_REFLECTION_ZONE returns true if receiver is within range +%of reflective waves by a wedge face +reflection_zone = false; + +if ~wedge.point_outside_wedge( source_pos ) || ~wedge.point_outside_wedge( receiver_pos ) + error( 'invalid source or receiver location!' ); +end + +apex_point = wedge.get_aperture_point(source_pos, receiver_pos); +source_apex_direction = ( apex_point - source_pos ) / norm( apex_point - source_pos ); +eps = wedge.set_get_geo_eps; + +% Distances of source from wedge faces +d1 = dot( source_pos, wedge.main_face_normal ); +d2 = dot( source_pos, wedge.opposite_face_normal ); + +% Check for the source facing side(s) of the wedge +if d1 >= -eps && d2 >= -eps +%% Source facing both sides of the wedge + % Check if source position is below face normal from aperture point + check_normal_main_face = cross( wedge.main_face_normal, wedge.aperture_direction ); + check_normal_opposite_face = cross( -wedge.opposite_face_normal, wedge.aperture_direction ); + + dist_below_main_face_normal = dot( -source_apex_direction, check_normal_main_face ); + dist_below_opposite_face_normal = dot( -source_apex_direction, check_normal_opposite_face ); + + if dist_below_main_face_normal >= -eps && dist_below_opposite_face_normal >= -eps % Opening angle of wedge > 180° + reflection_zone = true; + elseif dist_below_main_face_normal >= -eps % Opening angle of wedge between 90° and 180° + % Use of an auxiliary plane which is spanned by source apex direction and + % aperture direction for evaluating if receiver is between source + % position and wedge face. + aux_norm = cross( -source_apex_direction, wedge.aperture_direction ); + dist = dot( receiver_pos - apex_point, aux_norm ); + if dist >= -eps + reflection_zone = true; + else + % Angle of incidence + beta = acos( dot( -source_apex_direction, wedge.opposite_face_normal ) ); + + % Use of auxiliary plane describing the boundary of the + % reflection zone + refl_boundary_opposite_face = ita_rotation_rodrigues( -source_apex_direction, wedge.aperture_direction, 2*beta ); + refl_plane_normal_opposite_face = cross( refl_boundary_opposite_face, wedge.aperture_direction ); + + % Check if receiver is in reflection zone + dist_refl_opposite_face = dot( receiver_pos, refl_plane_normal_opposite_face ); + if dist_refl_opposite_face >= -eps + reflection_zone = true; + end + end + elseif dist_below_opposite_face_normal >= -eps % Opening angle of wedge between 90° and 180° + % Use of an auxiliary plane which is spanned by source apex direction and + % aperture direction for evaluating if receiver is between source + % position and wedge face. + aux_norm = cross( -source_apex_direction, wedge.aperture_direction ); + dist = dot( receiver_pos - apex_point, aux_norm ); + if dist >= -eps + reflection_zone = true; + else + % Angle of incidence + alpha = acos( dot( -source_apex_direction, wedge.main_face_normal ) ); + + % Use of auxiliary plane describing the boundary of the + % reflection zone + refl_boundary_main_face = ita_rotation_rodrigues( -source_apex_direction, -wedge.aperture_direction, 2*alpha ); + refl_plane_normal_main_face = cross( -refl_boundary_main_face, wedge.aperture_direction ); + + % Check if receiver is in reflection zone + dist_refl_main_face = dot( receiver_pos, refl_plane_normal_main_face ); + if dist_refl_main_face >= -eps + reflection_zone = true; + end + end + else % Opening angle of wedge < 90° + % Angles of incidence + alpha = acos( dot( -source_apex_direction, wedge.main_face_normal ) ); + beta = acos( dot( -source_apex_direction, wedge.opposite_face_normal ) ); + + % Use of auxiliary planes describing the boundary of the + % reflection zone + refl_boundary_main_face = ita_rotation_rodrigues( -source_apex_direction, -wedge.aperture_direction, 2*alpha ); + refl_boundary_opposite_face = ita_rotation_rodrigues( -source_apex_direction, wedge.aperture_direction, 2*beta ); + refl_plane_normal_main_face = cross( -refl_boundary_main_face, wedge.aperture_direction ); + refl_plane_normal_opposite_face = cross( refl_boundary_opposite_face, wedge.aperture_direction ); + + % Check if receiver is in reflection zone + dist_refl_main_face = dot( receiver_pos, refl_plane_normal_main_face ); + dist_refl_opposite_face = dot( receiver_pos, refl_plane_normal_opposite_face ); + if dist_refl_main_face >= -eps || dist_refl_opposite_face >= -eps + reflection_zone = true; + end + end +elseif d1 >= -eps +%% Source facing main face of the wedge + % Check if source position is below face normal from aperture point + check_normal = cross( wedge.main_face_normal, wedge.aperture_direction ); + if dot( -source_apex_direction, check_normal ) >= - eps + % Angle of incidence + alpha = acos( dot( -source_apex_direction, wedge.main_face_normal ) ); % Angle of incidence + + % Use of an auxiliary plane describing the boundary of the + % reflection zone + refl_boundary = ita_rotation_rodrigues( -source_apex_direction, wedge.aperture_direction, 2*alpha ); + refl_plane_normal = cross( refl_boundary, wedge.aperture_direction ); + + % Check if receiver is in reflection zone + dist = dot( receiver_pos - apex_point, refl_plane_normal ); + if dist >= -eps + reflection_zone = true; + end + else + % Use of an auxiliary plane which is spanned by source apex direction and + % aperture direction for evaluating if receiver is between source + % position and wedge face. + aux_normal = cross( -source_apex_direction, wedge.aperture_direction ); + dist = dot( receiver_pos - apex_point, aux_normal ); + if dist >= -eps + reflection_zone = true; + end + end +else +%% Source facing opposite face of the wedge + % Check if source position is below face normal from aperture point + check_normal = cross( wedge.opposite_face_normal, -wedge.aperture_direction ); + if dot( -source_apex_direction, check_normal ) >= - eps + % Angle of incidence + alpha = acos( dot( -source_apex_direction, wedge.opposite_face_normal ) ); % Angle of incidence + + % Use of an auxiliary plane describing the boundary of the + % reflection zone + refl_boundary = ita_rotation_rodrigues( -source_apex_direction, -wedge.aperture_direction, 2*alpha ); + refl_plane_normal = cross( wedge.aperture_direction, refl_boundary ); + + % Check if receiver is in reflection zone + dist = dot( receiver_pos - apex_point, refl_plane_normal ); + if dist >= -eps + reflection_zone = true; + end + else + % Use of an auxiliary plane which is spanned by source apex direction and + % aperture direction. + aux_normal = cross( -source_apex_direction, -wedge.aperture_direction ); + dist = dot( receiver_pos - apex_point, aux_normal ); + if dist >= -eps + reflection_zone = true; + end + end +end + +end + diff --git a/applications/VirtualAcoustics/diffraction/ita_rotation_rodrigues.m b/applications/VirtualAcoustics/diffraction/ita_rotation_rodrigues.m new file mode 100644 index 0000000..8a7a351 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/ita_rotation_rodrigues.m @@ -0,0 +1,32 @@ +function v_rot = ita_rotation_rodrigues( v, k, theta ) +%ITA_ROTATION_RODRIGUES rotates array of vectors by the angle theta +% around the axis k. Direction of the rotation is determined by a right +% handed system. +% v: Array can be composed of N rows of 3D row vectors or N columns of 3D column vectors. +% If v is 3x3 array, it is assumed that it is 3 rows of 3 3D row vectors. +% k: Rotation axis doesn't have to be normalized. +% theta: Rotation angle in radiant. + +[m, n] = size(v); +if m ~= 3 && n ~= 3 + error('input vector is/are not three dimensional') +end +if size(v) ~= size(k) + error('rotation vector v and axis k have different dimensions') +end + +k = k / norm( k ); +N = numel( v ) / 3; +v_rot = v; % Initialize result vector array + +if n == 3 % Row vectors + for i = 1 : N + v_rot(i, :) = v(i, :) * cos(theta) + cross( k, v ) * sin(theta) + k * dot( k, v ) * ( 1 - cos(theta) ); + end +else % Column vectors + for i = 1 : N + v_rot(:, i) = v(:, i) * cos(theta) + cross( k, v ) * sin(theta) + k * dot( k, v ) * ( 1 - cos(theta) ); + end +end +end + diff --git a/applications/VirtualAcoustics/diffraction/maekawa.m b/applications/VirtualAcoustics/diffraction/maekawa.m new file mode 100644 index 0000000..61af800 --- /dev/null +++ b/applications/VirtualAcoustics/diffraction/maekawa.m @@ -0,0 +1,88 @@ +% Test Maekawa diffraction algorithm + +%% Config +S = [ 0 -12 0.5 ]; +R = [ 0 2 0.5 ]; + +mainNormal = [0 -1 1]; +oppositeNormal = [0 1 1]; +pos = [4 0 1]; +c = 344; +A = [mainNormal; oppositeNormal; pos]'; +freq = ita_ANSI_center_frequencies; +w = itaInfiniteWedge(mainNormal, oppositeNormal, pos); + +%% Test operation for itaWedge + +% Wedge rotation around z axis +% rotAngle = pi/5; +% rotMatrix = [ cos(rotAngle), sin(rotAngle), 0; ... +% -sin(rotAngle), cos(rotAngle), 0; ... +% 0, 0, 1 ]; +% B = rotMatrix * A; + +% % Wedge rotation around aperture axis +% A(:, 3) = [0; 0; 0]; +% rotAngle = -pi/6; +% rotMatrix = [ 1, 0, 0; ... +% 0, cos(rotAngle), sin(rotAngle); ... +% 0, -sin(rotAngle), cos(rotAngle) ]; +% B = rotMatrix * A; +% B(:, 3) = B(:, 3) + pos'; +% +% itaWedge.set_get_geo_eps( 1e-6 ); % Micro meter precision for geo calcs +% +% w = itaInfiniteWedge( B(:, 1)', B(:, 2)', B(:, 3)' ); +% +% assert( w.point_outside_wedge( S ) ) +% assert( w.point_outside_wedge( R ) ) +% +% apex_point = w.get_aperture_point( S, R ) +% +% %% Maekawa Diffraction +% att = ita_diffraction_maekawa( w, S, R, ita_ANSI_center_frequencies ); + +% N = 0 : 0.5 : 100; +% att = zeros(2, numel(N)); +% att(1, :) = - ( 5 + 20 * log10( sqrt( 2*pi.*N ) ./ tanh( sqrt( 2*pi.*N ) ) ) ); % ( 10^(5/20) * sqrt(2*pi*N) ./ tanh( sqrt( 2*pi*N ) ) ); +% att(2, :) = - (10 * log10(3 + 20.*N) ); +% +% figure; +% semilogx( N, att ); +% % ylim( [0 , 33] ); +% grid on; + +att2 = itaResult; +att2.freqVector = freq; +att2.freqData = ita_diffraction_maekawa(w, S, R, freq, c); + +att2.pf +ylim([-40, 0]) + + +%% maekawa curve against fresnel number N +n1 = [1, 1, 0]; +n2 = [-1, 1, 0]; +loc = [0, 2, 0]; +src = [-4, 0, 0]; +rcv = [4, 0, 0]; +w = itaInfiniteWedge(n1, n2, loc); +ap = w.get_aperture_point(src, rcv); + +N = 0 : 0.001 : 100; +d = norm(ap - src) + norm(rcv - ap) - norm(rcv - src); +c = 344; +f = c * N / (2 * d); + +att = itaResult; +att.freqVector = freq; +att.freqData = ita_diffraction_maekawa(w, src, rcv, f, c); +res = -att.freqData_dB; +att_res = res(2:end); + +figure(); +semilogx(N(2:end), att_res); +xlabel('Fresnel number N'); +ylabel('Attenuation [dB]'); +ylim([0; 35]); +grid on; -- 2.26.2