Skip to content

GitLab

Commit b73ca0a5 authored by jonasstienen's avatar jonasstienen
Browse files

Improving interface of wedges and modifications of tests for utd

parent 3fc57de9
Hide whitespace changes
Inline Side-by-side
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/angle_main_face.m
View file @ b73ca0a5
......@@ -2,7 +2,7 @@ function alpha = angle_main_face( obj, point )
%angle_main_face Returns angle (radiant) between given point and a wedge face.
% input: point arbitrary field point outside the wedge.
% output: theta azimuth angle (radiant) in cylinder coordinates with
% the aperture as z axis. Between [0, 2*pi]
% the edge as z axis. Between [0, 2*pi]
if numel( point ) ~= 3
error( 'Point has to be of dimension 3' );
......@@ -16,7 +16,7 @@ end
% and use cylinder coordinates
e_y = obj.main_face_normal;
e_z = obj.aperture_direction;
e_z = obj.edge_direction;
e_x = cross( e_y, e_z );
x_cylinder = dot( point - obj.location, e_x );
......
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/apex_point.m
View file @ b73ca0a5
function aperture_point = get_aperture_point_far_field( obj, source_pos, receiver_pos )
function ap = apex_point( obj, source_pos, receiver_pos )
% GET_APERTURE_POINT_FAR_FIELD Returns aperture point on wedge (closest point on wedge
% between source and receiver if both are in the far field)
......@@ -12,7 +12,7 @@ assert( numel( receiver_pos ) == 3 )
source_receiver_vec = receiver_pos - source_pos;
source_receiver_dir = source_receiver_vec / norm( source_receiver_vec );
aperture_dir = obj.aperture_direction / norm( obj.aperture_direction );
edge_dir = obj.edge_direction / norm( obj.edge_direction );
if norm( source_receiver_vec ) == 0
warning( '@todo auxiliar plane must be created differently if source and receiver positions are equal. Trying to continue.' )
......@@ -20,7 +20,7 @@ end
% Auxilary plane spanned by source_receiver_dir and aux_plane_dir (closest
% line between aperture vector and source-receiver-vector
aux_plane_vec = cross( source_receiver_dir, aperture_dir );
aux_plane_vec = cross( source_receiver_dir, edge_dir );
if norm( aux_plane_vec ) ~= 0
......@@ -29,11 +29,11 @@ if norm( aux_plane_vec ) ~= 0
aux_plane_normal = cross( source_receiver_dir, aux_plane_dir );
% Determine intersection of line (aperture) and auxiliary plane
lambda_divisor = dot( aux_plane_normal, aperture_dir );
lambda_divisor = dot( aux_plane_normal, edge_dir );
assert( lambda_divisor ~= 0 )
d = dot( aux_plane_normal, obj.location ); % Distance to origin
lambda = ( d - dot( aux_plane_normal, source_pos ) );% / lambda_divisor; % .. or receiver_pos
aperture_point = obj.location + lambda * aperture_dir;
ap = obj.location + lambda * edge_dir;
else
......@@ -41,11 +41,11 @@ else
% onto aperture
% Project middle point onto aperture
lambda = dot( ( source_pos + source_receiver_vec ./ 2 )' - obj.location, aperture_dir );
aperture_point = obj.location + lambda * aperture_dir;
lambda = dot( ( source_pos + source_receiver_vec ./ 2 )' - obj.location, edge_dir );
ap = obj.location + lambda * edge_dir;
end
assert( any( ~isnan( aperture_point ) ) )
assert( any( ~isnan( ap ) ) )
end
\ No newline at end of file
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/apex_point_approx.m
View file @ b73ca0a5
function ap = approx_aperture_point( obj, source_pos, receiver_pos, spatial_precision )
function apx = apex_point_approx( obj, source_pos, receiver_pos, spatial_precision )
if nargin < 4
spatial_precision = 1e-3; % mm
......@@ -19,7 +19,7 @@ function ap = approx_aperture_point( obj, source_pos, receiver_pos, spatial_prec
opts_t = optimset( 'TolX', spatial_precision, 'TolFun', spatial_precision, 'FunValCheck', 'on' );
t = fminbnd(@(t)total_path_distance(t, source_pos, receiver_pos, ap_start, ap_dir), start_t, end_t, opts_t );
ap = ap_start + t .* ap_dir; %using the optimised parameter, find the aperture position
apx = ap_start + t .* ap_dir; %using the optimised parameter, find the aperture position
end
......
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/get_angle_from_point_to_apex.m
View file @ b73ca0a5
function alpha_rad = 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.
function alpha_rad = get_angle_from_point_to_apex( obj, field_point, point_on_edge )
%Returns angle (radiant) between the ray from field point to apex point and the
%edge of the wedge.
% output angle alpha: 0 <= alpha <= pi/2
if ~obj.point_outside_wedge( field_point )
error( 'Field point must be outside wedge' );
end
if ~obj.point_on_aperture( point_on_aperture )
error( 'No point on aperture found' )
if ~obj.point_on_edge( point_on_edge )
error( 'No point on edge found' )
end
dir_vec = ( point_on_aperture - field_point ) / norm( point_on_aperture - field_point );
alpha_rad = acos( dot( dir_vec, obj.aperture_direction ) );
dir_vec = ( point_on_edge - field_point ) / norm( point_on_edge - field_point );
alpha_rad = acos( dot( dir_vec, obj.edge_direction ) );
if alpha_rad > pi/2
alpha_rad = pi - alpha_rad;
......
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/itaInfiniteWedge.m
View file @ b73ca0a5
......@@ -3,7 +3,7 @@ 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)
ed % 3-dim edge direction vector (internal)
l % Internal location variable
et % type of edge (internal)
bc_hard % Internal boundary condition (hard = true)
......@@ -12,8 +12,9 @@ classdef itaInfiniteWedge
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)
aperture_direction % 3-dim normal vector of aperture direction (normalized) LEGACY
edge_direction % 3-dim normal vector of edge direction (normalized)
location % Location of wedge (somewhere along edge)
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)
edge_type % 'wedge' for opening angles > pi or 'corner' for opening angles < pi
......@@ -21,15 +22,16 @@ classdef itaInfiniteWedge
end
methods
function obj = itaInfiniteWedge( main_face_normal, opposite_face_normal, location, edge_type, aperture_direction )
function obj = itaInfiniteWedge( main_face_normal, opposite_face_normal, location, edge_type, edge_direction )
% 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: Point on edge which defines
% location of the wedge in 3_dim sapce
% edge_type: use 'outer_edge' for opening angles > pi (default) and
% 'inner_edge' for opening angles < pi
% edge_direction Edge direction vector (3-dim)
% Note: 3-dim direction vectors will be normalized automatically
%
if nargin < 4
......@@ -63,12 +65,12 @@ classdef itaInfiniteWedge
if nargin < 5
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.'
warning 'Normals are linear dependent and edge direction could not be determined. Please set edge direction manually.'
else
obj.ad = n_scaled ./ norm( n_scaled );
obj.ed = n_scaled ./ norm( n_scaled );
end
else
obj.ad = aperture_direction;
obj.ed = edge_direction;
end
end
......@@ -80,32 +82,32 @@ classdef itaInfiniteWedge
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
function n = get.edge_direction( obj )
% Returns normalized direction of edge. Vectors main face normal, opposite face normal and edge direction
% form a clockwise system.
if isempty( obj.ad )
error 'Invalid wedge, aperture direction not set and face normals are linear dependent'
if isempty( obj.ed )
error 'Invalid wedge, edge direction not set and face normals are linear dependent'
end
n = obj.ad;
n = obj.ed;
end
function obj = set.aperture_direction( obj, aperture_direction )
% Sets aperture direction manually (in case of linear
function obj = set.edge_direction( obj, edge_direction )
% Sets edge 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'
error 'Edge of linear independent normals is fixed can not be modified'
end
if ~norm( aperture_direction )
error 'Aperture vector must be a valid direction'
if ~norm( edge_direction )
error 'Edge vector must be a valid direction'
end
if norm( aperture_direction ) ~= 1
warning ' Normalizing aperture direction'
aperture_direction = aperture_direction / norm( aperture_direction );
if norm( edge_direction ) ~= 1
warning ' Normalizing edge direction'
edge_direction = edge_direction / norm( edge_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'
if ~( dot( edge_direction, obj.n1 ) == 0 && dot( edge_direction, obj.n2 ) == 0 )
error 'Invalid edge direction, vector must be perpendicular to face normals'
end
obj.ad = aperture_direction;
obj.ed = edge_direction;
end
function beta = get.wedge_angle( obj )
......@@ -170,14 +172,7 @@ classdef itaInfiniteWedge
bc = 'soft';
end
end
function aperture_point = get_aperture_point( obj, source_pos, receiver_pos )
% GET_APERTURE_POINT Returns approximated aperture point on wedge
%warning 'Aperture cannot be analytically determined, using approximation. See also get_aperture_point_far_field.'
%aperture_point = obj.approx_aperture_point( source_pos, receiver_pos );
aperture_point = obj.get_aperture_point_far_field( source_pos, receiver_pos );
end
function obj = set.boundary_condition( obj, bc )
if ischar(bc)
if strcmpi('hard', bc)
......@@ -202,7 +197,38 @@ classdef itaInfiniteWedge
function bc = is_boundary_condition_soft( obj )
bc = ~obj.bc_hard;
end
end
% Legacy support (before renaming aperture to apex)
function apx = approx_aperture_point( obj, source_pos, receiver_pos, spatial_precision )
apx = obj.apex_point_approx( obj, source_pos, receiver_pos, spatial_precision );
end
function ap = get_aperture_point_far_field( obj, source_pos, receiver_pos )
ap = obj.apex_point( obj, source_pos, receiver_pos );
end
function b = point_on_aperture( obj, point )
b = obj.point_on_edge( obj, point );
end
function alpha_rad = get_angle_from_point_to_aperture( obj, field_point, point_on_edge )
alpha_rad = obj.get_angle_from_point_to_apex( obj, field_point, point_on_edge );
end
function ap = get_aperture_point( obj, source_pos, receiver_pos )
ap = obj.apex_point( obj, source_pos, receiver_pos );
end
function n = get.aperture_direction( obj )
n = obj.edge_direction;
end
function obj = set.aperture_direction( obj, edge_direction_ )
obj.edge_direction = edge_direction_;
end
end
......
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/point_on_edge.m
View file @ b73ca0a5
function b = point_on_aperture( obj, point )
function b = point_on_edge( obj, point )
% Returns true if point is on aperture of the wedge
if numel( point ) ~= 3
......
applications/SoundFieldSimulation/Diffraction/ita_diffraction_shadow_zone.m
View file @ b73ca0a5
......@@ -49,7 +49,7 @@ end
% opposite face plane
if source_pos_above_main_plane
source_shadow_boundary_plane_dir = cross( wedge.aperture_direction, source_pos - wedge.location );
source_shadow_boundary_plane_dir = cross( wedge.edge_direction, source_pos - wedge.location );
source_shadow_boundary_plane_normal = source_shadow_boundary_plane_dir / norm( source_shadow_boundary_plane_dir );
if dot( source_shadow_boundary_plane_normal, receiver_pos - source_pos ) >= 0
in_shadow = false;
......@@ -60,7 +60,7 @@ if source_pos_above_main_plane
end
if source_pos_above_opposite_plane
source_shadow_boundary_plane_dir = cross( wedge.aperture_direction, source_pos - wedge.location );
source_shadow_boundary_plane_dir = cross( wedge.edge_direction, source_pos - wedge.location );
source_shadow_boundary_plane_normal = source_shadow_boundary_plane_dir / norm( source_shadow_boundary_plane_dir );
if dot( source_shadow_boundary_plane_normal, receiver_pos - source_pos ) >= 0
in_shadow = true;
......
applications/SoundFieldSimulation/Diffraction/ita_diffraction_utd.m
View file @ b73ca0a5
function [ diffr_field, D, A ] = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequency_vec, speed_of_sound, apex_point )
%ITA_DIFFRACTION_UTD Calculates the diffraction filter based on uniform
%theory of diffraction (with Kawai approximation). Apex point is optional
%and can
%theory of diffraction (UTD, with Kawai approximation). Apex point is optional.
%
% Literature:
% [1] Tsingos, Funkhouser et al. - Modeling Acoustics in Virtual Environments using the Uniform Theory of Diffraction
......@@ -20,26 +19,26 @@ if numel( receiver_pos ) ~= 3
end
if nargin < 6
apex_point = wedge.approx_aperture_point( source_pos, receiver_pos );
apex_point = wedge.apex_point( source_pos, receiver_pos );
end
%% Variables
rho = norm( apex_point - source_pos ); % Distance of source to aperture point
r = norm( receiver_pos - apex_point ); % Distance of receiver to aperture point
rho = norm( apex_point - source_pos ); % Distance of source to edge point
r = norm( receiver_pos - apex_point ); % Distance of receiver to edge point
assert( rho + r ~= 0 && r ~= 0 );
alpha_i = wedge.angle_main_face( source_pos );
alpha_d = wedge.angle_main_face( receiver_pos );
theta_i = wedge.get_angle_from_point_to_aperture( source_pos, apex_point );
theta_i = wedge.get_angle_from_point_to_apex( source_pos, apex_point );
n = wedge.opening_angle / pi; % Variable dependend on opening angle of the wedge
k = (2 * pi * frequency_vec) ./ speed_of_sound; % Wavenumber
k = ( 2 * pi * frequency_vec ) ./ speed_of_sound; % Wavenumber
%% Calculations
H_i = 1 ./ rho .* exp( -1i * k .* rho ); % direct path from source to apex point
A = sqrt( rho ./ ( r .* ( rho + r ) ) ); % attenuation factor at receiver
D = get_diffr_coeff( wedge, k, alpha_d, alpha_i, rho, r, theta_i, n );
H_i = 1 ./ rho .* exp( -1i * k .* rho ); % direct path from source to apex point, ideal point source
A = sqrt( rho ./ ( r .* ( rho + r ) ) ); % attenuation factor at receiver after diffraction at edge
D = get_diffr_coeff( wedge, k, alpha_d, alpha_i, rho, r, theta_i, n ); % diffraction coefficient
% Combined diffracted sound field filter at receiver
diffr_field = H_i .* D .* A .* exp( -1i .* k .* r );
......@@ -135,7 +134,7 @@ end
%% Auxiliary functions
% N+ function
function N = N_p( n, beta )
if beta > pi * (n - 1)
if beta > pi * ( n - 1 )
N = 1;
else
N = 0;
......@@ -144,9 +143,9 @@ end
% N- function
function N = N_n( n, beta )
if beta > pi * (1 + n)
if beta > pi * ( 1 + n )
N = 1;
elseif beta < pi * (n - 1)
elseif beta < pi * ( n - 1 )
N = -1;
else
N = 0;
......@@ -172,4 +171,4 @@ function Y = kawai_approx_fresnel( X )
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
end
applications/SoundFieldSimulation/Diffraction/ita_diffraction_utd_approx.m
View file @ b73ca0a5
......@@ -22,7 +22,7 @@ if ~in_shadow_zone
end
%% Variables
apex_point = wedge.approx_aperture_point( source_pos, receiver_pos );
apex_point = wedge.apex_point_approx( source_pos, receiver_pos );
dir_src_2_apex_pt = ( apex_point - source_pos ) ./ norm( apex_point - source_pos );
dir_apex_pt_2_rcv = ( receiver_pos - apex_point ) ./ norm( receiver_pos - apex_point );
r = norm( apex_point - source_pos ); % Distance Source to Apex point
......
applications/SoundFieldSimulation/Diffraction/plots/ita_diffraction_plot_utd_total_field.m
View file @ b73ca0a5
......@@ -16,7 +16,7 @@ speed_of_sound = 344;
%% Calculations
% Set receiver positions alligned around the aperture
%apex_point = inf_wdg.get_aperture_point(src_pos, rcv_start_pos);
%src_is_facing_main_face = ~inf_wdg.point_facing_main_side( src_pos );
src_is_facing_main_face = ~inf_wdg.point_facing_main_side( src_pos );
alpha_start = inf_wdg.get_angle_from_point_to_wedge_face(rcv_start_pos, src_is_facing_main_face); % implement without "facing" anything
alpha_end = inf_wdg.get_angle_from_point_to_wedge_face(rcv_end_pos, src_is_facing_main_face);
alpha_d = linspace( alpha_start, alpha_end, angle_resolution );
......
applications/SoundFieldSimulation/Diffraction/plots/ita_diffraction_plot_utd_tsingos_total_field.m
View file @ b73ca0a5
......@@ -6,7 +6,7 @@ source_pos = 5 * [ -1 0 0];
receiver_start_pos = 5 * [ -1 -1 0 ] / sqrt( 2 );
w = itaInfiniteWedge( n1, n2, loc );
apex_point = w.approx_aperture_point( source_pos, receiver_start_pos );
apex_point = w.apex_point( source_pos, receiver_start_pos );
apex_dir = w.aperture_direction;
c = 344; % Speed of sound
......
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_utd_shadow_zone_transition.m
View file @ b73ca0a5
......@@ -42,7 +42,7 @@ utd_tf.channelNames = { 'Diffracted field (shadow)', 'Diffracted field (illumina
receiver_start_pos = 5 * [ -1 -1 0 ] / sqrt( 2 );
apex_point = w.approx_aperture_point( source_pos, receiver_start_pos );
apex_point = w.apex_point_approx( source_pos, receiver_start_pos );
apex_dir = w.aperture_direction;
freq = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment