Adding diffraction simulation scripts and wedge classes

parent bf723f4c
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
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
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
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
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
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
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
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
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
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
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
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
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
%% 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,