Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Institute of Technical Acoustics (ITA)
toolbox
Commits
684e8595
Commit
684e8595
authored
Dec 02, 2019
by
Dipl.-Ing. Jonas Stienen
Browse files
Great clean-up before merging
parent
153229b6
Changes
21
Hide whitespace changes
Inline
Side-by-side
applications/SoundFieldSimulation/Diffraction/@itaFiniteWedge/point_on_aperture.m
View file @
684e8595
...
...
@@ -18,4 +18,3 @@ else
end
end
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/
get
_aperture_point
2
.m
→
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/
approx
_aperture_point.m
View file @
684e8595
function ap = get_aperture_point2( obj, source_pos, receiver_pos )
start = obj.location;
dir = obj.aperture_direction;
function ap = approx_aperture_point( obj, source_pos, receiver_pos, spatial_precision )
if nargin
<
4
spatial_precision =
1e-3;
%
mm
end
ap_start =
obj.location;
ap_dir =
obj.aperture_direction;
S_on_ap = orthogonal_projection( start, dir, source_pos ); %project the source to the aperture
S_t = (S_on_ap - start) / dir; S_t = S_t(1,1); %find the parametric distance along the aperture
R_on_ap = orthogonal_projection( start, dir, receiver_pos ); %same as above but for receiver
R_t = (R_on_ap - start) / dir; R_t = R_t(1,1);
S_on_ap =
orthogonal_projection(
ap_start
,
ap_dir
,
source_pos
);
%
project
the
source
to
the
aperture
S_t =
(S_on_ap
-
ap_start
)
/
ap_dir
;
S_t =
S_t(1,1);
%
find
the
parametric
distance
along
the
aperture
R_on_ap =
orthogonal_projection(
ap_start
,
ap_dir
,
receiver_pos
);
%
same
as
above
but
for
receiver
R_t =
(R_on_ap
-
ap_start
)
/
ap_dir
;
R_t =
R_t(1,1);
start_t =
min(
S_t
,
R_t
);
%
start
the
optimisation
at
whichever
of
the
projected
source
/
receiver
comes
first
on
the
aperture
end_t =
max(
S_t
,
R_t
);
%
finish
at
the
other
t = fminbnd(@(t)total_path_distance(t, source_pos, receiver_pos, start, dir), start_t, end_t );
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 = start + t
*
dir; %using the optimised parameter, find the aperture position
ap =
ap_
start
+
t
.
*
ap_
dir
;
%
using
the
optimised
parameter
,
find
the
aperture
position
end
...
...
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/get_aperture_point.m
→
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/get_aperture_point
_far_field
.m
View file @
684e8595
function
aperture_point
=
get_aperture_point
(
obj
,
source_pos
,
receiver_pos
)
% get_aperture_point Returns aperture point on wedge (closest point on wedge
% between source and receiver)
aperture_point
=
obj
.
get_aperture_point2
(
source_pos
,
receiver_pos
);
return
% ... the rest is faulty!
function
aperture_point
=
get_aperture_point_far_field
(
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)
assert
(
numel
(
source_pos
)
==
3
)
assert
(
numel
(
receiver_pos
)
==
3
)
...
...
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/itaInfiniteWedge.m
View file @
684e8595
...
...
@@ -17,7 +17,6 @@ classdef itaInfiniteWedge
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
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
...
...
@@ -140,6 +139,13 @@ classdef itaInfiniteWedge
l
=
obj
.
l
;
end
function
obj
=
set
.
location
(
obj
,
location
)
if
numel
(
location
)
~=
3
error
(
'Location must be of dimension 3'
)
end
obj
.
l
=
location
;
end
function
b
=
validate_normals
(
obj
)
% Returns true, if the normals of the faces are both normalized
b
=
false
;
...
...
@@ -148,7 +154,7 @@ classdef itaInfiniteWedge
end
end
function
et
=
get
.
wedge_type
(
obj
)
function
et
=
wedge_type
(
obj
)
et
=
obj
.
edge_type
;
warning
'Function
''
wedge_type
''
is deprecated, use
''
edge_type
''
instead.'
end
...
...
@@ -165,6 +171,12 @@ classdef itaInfiniteWedge
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
);
end
function
obj
=
set
.
boundary_condition
(
obj
,
bc
)
if
ischar
(
bc
)
if
strcmpi
(
'hard'
,
bc
)
...
...
applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/point_on_aperture.m
View file @
684e8595
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
if
numel
(
point
)
~=
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
;
d
=
point
-
obj
.
location
;
if
norm
(
d
)
<
obj
.
set_get_geo_eps
b
=
true
;
% Point is indescriminably close to wedge location
else
dir1
=
d
/
norm
(
d
);
d2
=
bj
.
aperture_direction
-
obj
.
location
;
dir2
=
d2
/
norm
(
d2
);
% Point should have same (or opposite) direction as aperture
if
1
-
abs
(
dot
(
dir1
,
dir2
)
)
<
obj
.
set_get_geo_eps
b
=
true
;
end
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
applications/SoundFieldSimulation/Diffraction/ita_diffraction_biot_tolstoy_jst.m
deleted
100644 → 0
View file @
153229b6
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
applications/SoundFieldSimulation/Diffraction/ita_diffraction_btm_infinite_wedge.m
View file @
684e8595
function
att
=
ita_diffraction_btm_infinite_wedge
(
wedge
,
source_pos
,
receiver_pos
,
sampling_rate
,
filter_length_samples
,
boundary_condition
)
function
att
_ir
=
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
...
...
@@ -18,49 +18,17 @@ 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
A
pex_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
);
A
=
wedge
.
get_aperture_point
(
source_pos
,
receiver_pos
);
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
,
A
);
theta_w
=
wedge
.
opening_angle
;
SA
=
N
orm
(
A
pex_Point
-
S
);
AR
=
N
orm
(
R
-
A
pex_Point
);
SA
=
n
orm
(
A
-
S
);
AR
=
n
orm
(
R
-
A
);
r_S
=
SA
.*
sin
(
theta_i
);
r_R
=
AR
.*
sin
(
theta_i
);
z_R
=
(
SA
+
AR
)
.*
cos
(
theta_i
);
...
...
@@ -68,24 +36,17 @@ ny = pi / theta_w;
L_0
=
sqrt
(
(
r_S
+
r_R
)
.^
2
+
z_R
.^
2
);
c
=
343.21
;
% Speed of sound at 20C
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
)
);
eta
=
acosh
(
(
(
c
.*
tau
)
.^
2
-
(
r_S
.^
2
+
r_R
.^
2
+
z_R
.^
2
)
)
.
/
(
2
*
r_S
.*
r_R
)
);
phi1
=
pi
+
T
heta_S
+
T
heta_R
;
phi2
=
pi
+
T
heta_S
-
T
heta_R
;
phi3
=
pi
-
T
heta_S
+
T
heta_R
;
phi4
=
pi
-
T
heta_S
-
T
heta_R
;
phi1
=
pi
+
t
heta_S
+
t
heta_R
;
phi2
=
pi
+
t
heta_S
-
t
heta_R
;
phi3
=
pi
-
t
heta_S
+
t
heta_R
;
phi4
=
pi
-
t
heta_S
-
t
heta_R
;
beta_pp
=
sin
(
ny
*
phi1
)
.
/
(
cosh
(
ny
.*
eta
)
-
cos
(
ny
*
phi1
)
);
beta_pm
=
sin
(
ny
*
phi2
)
.
/
(
cosh
(
ny
.*
eta
)
-
cos
(
ny
*
phi2
)
);
...
...
@@ -99,10 +60,6 @@ switch boundary_condition
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
att_ir
=
(
(
-
((
c
*
ny
)
/
(
2
*
pi
))
*
beta
)
.
/
(
r_S
.*
r_R
.*
sinh
(
eta
)
)
)
'
*
T
;
function
res
=
Norm
(
A
)
res
=
sqrt
(
sum
(
A
.^
2
,
2
)
);
end
applications/SoundFieldSimulation/Diffraction/ita_diffraction_maekawa.m
View file @
684e8595
...
...
@@ -25,9 +25,9 @@ end
%% Calculation
apex_
P
oint
=
wedge
.
get_aperture_point
(
source_pos
,
receiver_pos
);
apex_
p
oint
=
wedge
.
get_aperture_point
(
source_pos
,
receiver_pos
);
r_dir
=
norm
(
receiver_pos
-
source_pos
);
detour
=
norm
(
apex_
P
oint
-
source_pos
)
+
norm
(
receiver_pos
-
apex_
P
oint
)
-
norm
(
receiver_pos
-
source_pos
);
detour
=
norm
(
apex_
p
oint
-
source_pos
)
+
norm
(
receiver_pos
-
apex_
p
oint
)
-
norm
(
receiver_pos
-
source_pos
);
lambda
=
speed_of_sound
.
/
frequencies
;
N
=
2
*
detour
.
/
lambda
;
% Fresnel number N
...
...
applications/SoundFieldSimulation/Diffraction/ita_diffraction_maekawa_approx.m
View file @
684e8595
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.
% Calculates the attenuation filter at a diffraction wedge for source
% and receiver location(s) at given frequencies. For purpose of
% continuity at shadow boundary, an interpolation with exponential function
% towards target response in shadow region is used.
%
% Usage: see ita_diffraction_maekawa
% Parameter: transition_const = 0.2 rad (default) determines the angle into
% shadow region where transition reaches target filter of Maekawa's
% formula.
%
% wedge: diffracting wedge (itaInfiniteWedge or derived class
% itaFiniteWedge)
% source_pos: position of the source
% receiver_pos: position of the receiver
%% Assertions
...
...
@@ -15,62 +16,34 @@ 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
A
pex_
P
oint
=
wedge
.
get_aperture_point
(
source_pos
,
receiver_pos
);
S
rc_Apex_Dir
=
(
A
pex_
P
oint
-
source_pos
)
.
/
N
orm
(
A
pex_
P
oint
-
source_pos
);
A
pex_Rcv_Dir
=
(
receiver_pos
-
A
pex_
P
oint
)
.
/
N
orm
(
receiver_pos
-
A
pex_
P
oint
);
detour
=
N
orm
(
A
pex_
P
oint
-
source_pos
)
+
N
orm
(
receiver_pos
-
A
pex_
P
oint
)
-
N
orm
(
receiver_pos
-
source_pos
);
a
pex_
p
oint
=
wedge
.
get_aperture_point
(
source_pos
,
receiver_pos
);
S
A
=
(
a
pex_
p
oint
-
source_pos
)
.
/
n
orm
(
a
pex_
p
oint
-
source_pos
);
A
R
=
(
receiver_pos
-
a
pex_
p
oint
)
.
/
n
orm
(
receiver_pos
-
a
pex_
p
oint
);
detour
=
n
orm
(
a
pex_
p
oint
-
source_pos
)
+
n
orm
(
receiver_pos
-
a
pex_
p
oint
)
-
n
orm
(
receiver_pos
-
source_pos
);
c
=
speed_of_sound
;
lambda
=
c
.
/
frequencies
;
N
=
2
*
detour
.
/
lambda
;
% Fresnel number N
r_dir
=
N
orm
(
receiver_pos
-
source_pos
);
r_dir
=
n
orm
(
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
=
acos
(
AR
,
SA
);
% angle between receiver and shadow boundary
if
phi
>
pi
/
4
phi
=
pi
/
4
;
end
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, :) )'
;
if
in_shadow_zone
H_di
ff
r
=
zeros
(
numel
(
frequencies
),
1
);
else
% From Handbook of Acoustics page 117 eq. 4.13
H_diffr
=
c_approx
.*
(
c_norm
*
sqrt
(
2
*
pi
*
N
.
/
tanh
(
sqrt
(
2
*
pi
*
N
)
)
)
.^
(
-
1
)
.
/
r_dir
);
end
function
res
=
Norm
(
A
)
res
=
sqrt
(
sum
(
A
.^
2
,
2
)
);
end
\ No newline at end of file
end
applications/SoundFieldSimulation/Diffraction/ita_reflection_zone.m
→
applications/SoundFieldSimulation/Diffraction/ita_
diffraction_
reflection_zone
_main
.m
View file @
684e8595
function
reflection_zone
=
ita_reflection_zone
(
wedge
,
source_pos
,
receiver_pos
)
function
reflection_zone
=
ita_
diffraction_
reflection_zone
_main
(
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
;
...
...
@@ -7,7 +7,7 @@ if ~wedge.point_outside_wedge( source_pos ) || ~wedge.point_outside_wedge( recei
error
(
'invalid source or receiver location!'
);
end
apex_point
=
wedge
.
get
_aperture_point
(
source_pos
,
receiver_pos
);
apex_point
=
wedge
.
approx
_aperture_point
(
source_pos
,
receiver_pos
);
source_apex_direction
=
(
apex_point
-
source_pos
)
/
norm
(
apex_point
-
source_pos
);
eps
=
wedge
.
set_get_geo_eps
;
...
...
@@ -80,8 +80,8 @@ if d1 >= -eps && d2 >= -eps
% Use of auxiliary planes describing the boundary of the
% reflection zone
refl_boundary_main_face
=
ita_
rotation_rodrigue
s
(
-
source_apex_direction
,
-
wedge
.
aperture_direction
,
2
*
alpha
);
refl_boundary_opposite_face
=
ita_
rotation_rodrigue
s
(
-
source_apex_direction
,
wedge
.
aperture_direction
,
2
*
beta
);
refl_boundary_main_face
=
ita_
diffraction_rotate_vectors_around_axi
s
(
-
source_apex_direction
,
-
wedge
.
aperture_direction
,
2
*
alpha
);
refl_boundary_opposite_face
=
ita_
diffraction_rotate_vectors_around_axi
s
(
-
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
);
...
...
@@ -102,7 +102,7 @@ elseif d1 >= -eps
% Use of an auxiliary plane describing the boundary of the
% reflection zone
refl_boundary
=
ita_
rotation_rodrigue
s
(
-
source_apex_direction
,
wedge
.
aperture_direction
,
2
*
alpha
);
refl_boundary
=
ita_
diffraction_rotate_vectors_around_axi
s
(
-
source_apex_direction
,
wedge
.
aperture_direction
,
2
*
alpha
);
refl_plane_normal
=
cross
(
refl_boundary
,
wedge
.
aperture_direction
);
% Check if receiver is in reflection zone
...
...
@@ -130,7 +130,7 @@ else
% Use of an auxiliary plane describing the boundary of the
% reflection zone
refl_boundary
=
ita_
rotation_rodrigue
s
(
-
source_apex_direction
,
-
wedge
.
aperture_direction
,
2
*
alpha
);
refl_boundary
=
ita_
diffraction_rotate_vectors_around_axi
s
(
-
source_apex_direction
,
-
wedge
.
aperture_direction
,
2
*
alpha
);
refl_plane_normal
=
cross
(
wedge
.
aperture_direction
,
refl_boundary
);
% Check if receiver is in reflection zone
...
...
applications/SoundFieldSimulation/Diffraction/ita_diffraction_rotate_vectors_around_axis.m
View file @
684e8595
function
v_rot
=
ita_diffraction_rotate_vectors_around_axis
(
v
,
k
,
theta
)
%
ITA_ROTATION_RODRIGUES
rotates array of vectors by the angle theta
%
ita_diffraction_rotate_vectors_around_axis
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.
...
...
applications/SoundFieldSimulation/Diffraction/ita_diffraction_shadow_zone.m
View file @
684e8595
...
...
@@ -3,10 +3,10 @@ function in_shadow = ita_diffraction_shadow_zone( wedge, source_pos, receiver_po
%point of view, covered by the wedge and therefor inside the shadow region
%% assertions
if
~
numel
(
source_pos
)
=
=
3
if
numel
(
source_pos
)
~
=
3
error
(
'Source point must be of dimension 3'
)
end
if
~
numel
(
receiver_pos
)
if
numel
(
receiver_pos
)
~=
3
error
(
'Receiver point must be of dimension 3'
)
end
...
...
applications/SoundFieldSimulation/Diffraction/ita_diffraction_utd.m
View file @
684e8595
...
...
@@ -20,7 +20,7 @@ if numel( receiver_pos ) ~= 3
end
if
nargin
<
6
apex_point
=
wedge
.
get
_aperture_point
2
(
source_pos
,
receiver_pos
);
apex_point
=
wedge
.
approx
_aperture_point
(
source_pos
,
receiver_pos
);
end
%% Variables
...
...
applications/SoundFieldSimulation/Diffraction/ita_diffraction_visualize_scene.m
View file @
684e8595
...
...
@@ -8,7 +8,7 @@ if nargin < 4
is_opengl_cs
=
false
;
end
apx
=
w
.
get
_aperture_point
2
(
source_pos
,
target_pos
);
apx
=
w
.
approx
_aperture_point
(
source_pos
,
target_pos
);
d1
=
norm
(
source_pos
-
apx
);
d2
=
norm
(
target_pos
-
apx
);
d3
=
norm
(
w
.
location
-
source_pos
);
...
...
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_align_point_around_aperture.m
View file @
684e8595
...
...
@@ -15,16 +15,12 @@ point_pos_end_angle = inf_wdg.get_angle_from_point_to_wedge_face(point_end_pos,
point_angles
=
linspace
(
point_pos_start_angle
,
point_pos_end_angle
,
num_of_positions
);
%% Set different receiver positions rotated around the aperture
aligned_positions
=
ita_align_points_around_aperture
(
inf_wdg
,
point_start_pos
,
point_angles
,
apex_point
,
false
);
aligned_positions
=
ita_
diffraction_
align_points_around_aperture
(
inf_wdg
,
point_start_pos
,
point_angles
,
apex_point
,
false
);
orig_dist
=
norm
(
point_start_pos
);
for
i
=
1
:
size
(
aligned_positions
,
1
)
assert
(
point_keeps_same_dist
(
aligned_positions
(
i
,
:),
orig_dist
,
inf_wdg
.
set_get_geo_eps
*
10
)
)
diff
=
norm
(
aligned_positions
(
i
,
:)
)
-
orig_dist
;
assert
(
abs
(
diff
)
<
inf_wdg
.
set_get_geo_eps
*
10
)
end
disp
(
'Test passed'
);
%% aux function
function
res
=
point_keeps_same_dist
(
point
,
orig_dist_from_aperture
,
tolerance
)
diff
=
norm
(
point
)
-
orig_dist_from_aperture
;
res
=
abs
(
diff
)
<
tolerance
;
end
\ No newline at end of file
disp
(
'Test passed'
)
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_aperture_point_algorithm.m
→
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_aperture_point_
far_field_
algorithm.m
View file @
684e8595
...
...
@@ -16,31 +16,33 @@ w = itaInfiniteWedge( n_main / norm( n_main ), n_opposite / norm( n_opposite ),
% Initial setup
source_pos
=
[
1
-
1
1
];
receiver_pos
=
[
1
1
-
1
];
apx
=
w
.
get_aperture_point
(
source_pos
,
receiver_pos
);
apx
=
w
.
get_aperture_point
_far_field
(
source_pos
,
receiver_pos
);
% Plain scaling is valid
lambda
=
10
;
assert
(
all
(
apx
==
w
.
get_aperture_point
(
source_pos
*
lambda
,
receiver_pos
*
lambda
)
)
)
w_scaled
=
w
;
w_scaled
.
location
=
w_scaled
.
location
.*
lambda
;
assert
(
all
(
apx
==
w
.
get_aperture_point_far_field
(
source_pos
*
lambda
,
receiver_pos
*
lambda
)
)
)
% Same movement along source-receiver-direction is valid
source_receiver_vec
=
receiver_pos
-
source_pos
;
assert
(
all
(
apx
==
w
.
get_aperture_point
(
source_pos
-
lambda
*
source_receiver_vec
,
receiver_pos
+
lambda
*
source_receiver_vec
)
)
)
assert
(
all
(
apx
==
w
.
get_aperture_point
_far_field
(
source_pos
-
lambda
*
source_receiver_vec
,
receiver_pos
+
lambda
*
source_receiver_vec
)
)
)
% Arbitrary uneven movement along source-receiver-direction is valid
lambda_1
=
5
;
lambda_2
=
13
;
assert
(
all
(
apx
==
w
.
get_aperture_point
(
source_pos
-
lambda_1
*
source_receiver_vec
,
receiver_pos
+
lambda_2
*
source_receiver_vec
)
)
)
assert
(
all
(
apx
==
w
.
get_aperture_point
_far_field
(
source_pos
-
lambda_1
*
source_receiver_vec
,
receiver_pos
+
lambda_2
*
source_receiver_vec
)
)
)
% Same movement along source-apex-direction and receiver-apex-direction is valid
source_pos_scaled
=
source_pos
-
lambda
*
(
apx
-
source_pos
);
receiver_pos_scaled
=
receiver_pos
-
lambda
*
(
apx
-
receiver_pos
);
apx_scaled
=
w
.
get_aperture_point
(
source_pos_scaled
,
receiver_pos_scaled
);
apx_scaled
=
w
.
get_aperture_point
_far_field
(
source_pos_scaled
,
receiver_pos_scaled
);
assert
(
all
(
apx
==
apx_scaled
)
)
% Aribtrary uneven movement along source-apex-direction and receiver-apex-direction is NOT valid
source_pos_tilted
=
source_pos
-
lambda_1
*
(
apx
-
source_pos
);
receiver_pos_tilted
=
receiver_pos
-
lambda_2
*
(
apx
-
receiver_pos
);
apx_tilted
=
w
.
get_aperture_point
(
source_pos_tilted
,
receiver_pos_tilted
);
apx_tilted
=
w
.
get_aperture_point
_far_field
(
source_pos_tilted
,
receiver_pos_tilted
);
assert
(
all
(
apx
~=
apx_tilted
)
&&
all
(
apx_scaled
~=
apx_tilted
)
)
disp
(
'all good'
)
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_utd_shadow_zone_transition.m
View file @
684e8595
...
...
@@ -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
.
get
_aperture_point
(
source_pos
,
receiver_start_pos
);
apex_point
=
w
.
approx
_aperture_point
(
source_pos
,
receiver_start_pos
);
apex_dir
=
w
.
aperture_direction
;
freq
=
[
20
,
50
,
100
,
200
,
400
,
800
,
1600
,
3200
,
6400
,
12800
,
24000
]
'
;
...
...
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_visualization.m
View file @
684e8595
...
...
@@ -8,10 +8,8 @@ w = itaInfiniteWedge( n1 / norm( n1 ), n2 / norm( n2 ), loc );
s
=
[
3.1
0
0
];
r
=
[
-
3.1
0
-
0.1
];
figure
ita_diffraction_visualize_scene
(
w
,
s
,
r
)
w_inner
=
itaInfiniteWedge
(
n1
/
norm
(
n1
),
n2
/
norm
(
n2
),
loc
,
'inner_edge'
);
figure
ita_diffraction_visualize_scene
(
w_inner
,
s
,
r
)
applications/SoundFieldSimulation/Diffraction/tests/ita_diffraction_test_zones.m
View file @
684e8595