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

Improving code in get path data / DSP coeff getter

parent 419afea6
function alpha = get_angle_from_point_to_aperture( obj, field_point, point_on_aperture )
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.
% 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 ~obj.point_outside_wedge( field_point )
error( 'Field point must be outside wedge' );
end
......@@ -32,13 +11,10 @@ if ~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);
dir_vec = ( point_on_aperture - field_point ) / norm( point_on_aperture - field_point );
alpha_rad =dot( dir_vec, obj.aperture_direction );
if alpha_rad > pi/2
alpha_rad = pi - alpha_rad;
end
end
......
function [frequency_mags, gain, delay] = ita_propagation_path_get_data( path_struct, f, c )
path_data = path_struct.propagation_anchors;
gain = 1;
frequency_mags = ones(1,length(f));
total_distance = 0;
is_diff = 0; %set to 1 whenever a diffraction is encountered
%{
figure
for i = 1:numel(path_data)-1
plot3([path_data{i}.interaction_point(1),path_data{i+1}.interaction_point(1)],[path_data{i}.interaction_point(2),path_data{i+1}.interaction_point(2)],[path_data{i}.interaction_point(3),path_data{i+1}.interaction_point(3)])
hold on
function [ frequency_mags, gain, delay, valid ] = ita_propagation_path_get_data( pp, f, c )
% ita_propagation_path_get_data Returns frequency magnitudes, gain and
% delay for DSP processing. 4th return value is a validity flag.
%
% Example: [ frequency_mags, gain, delay, valid ] = ita_propagation_path_get_data( path_struct, f, c )
%
if ~isfield( pp, 'propagation_anchors' )
error 'Did not find a field ''propagation_anchors'' in this propagation path.'
end
gain = 1;
frequency_mags = ones( 1, numel( f ) );
total_distance = 0;
pd = pp.propagation_anchors;
N = numel( pd );
if N < 1
% No path constructable
valid = false;
return;
elseif N == 2
source = pd( 1 );
receiver = pd( 2 );
total_distance = norm( receiver.interaction_point - source.interaction_point );
end
is_diff = 0; %set to 1 whenever a diffraction is encountered
valid = true;
for i = 2:N-1 %start from 2, first entry is always source, -1 as receiver always the last
if isa( pd, 'struct' )
a_prev = pd( i-1 );
a_curr = pd( i );
a_next = pd( i+1 );
else
a_prev = pd{ i-1 };
a_curr = pd{ i };
a_next = pd{ i+1 };
end
%}
for i = 2:numel(path_data)-1 %start from 2, first entry is always source, -1 as receiver always the last
anchor_type = path_data{i}.anchor_type;
segment_distance = norm(path_data{ i }.interaction_point(1:3) - path_data{ i-1 }.interaction_point(1:3));
total_distance = total_distance + segment_distance;
switch anchor_type
anchor_type = a_curr.anchor_type;
segment_distance = norm( a_curr.interaction_point - a_prev.interaction_point );
total_distance = total_distance + segment_distance;
switch anchor_type
case 'outer_edge_diffraction' %case for diffraction
w = ita_propagation_wedge_from_diffraction_anchor( a_curr );
w.set_get_geo_eps( 1e-6 );
source_pos = a_prev.interaction_point(1:3);
target_pos = a_next.interaction_point(1:3);
case 'outer_edge_diffraction' %case for diffraction
main_face_normal(1,:) = path_data{i}.main_wedge_face_normal(1:3);
opposite_face_normal(1,:) = path_data{i}.opposite_wedge_face_normal(1:3);
aperture_start(1,:) = path_data{i}.vertex_start(1:3); %aperture point
vertex_length(1,:) = norm( path_data{i}.vertex_start(1:3) - path_data{i}.vertex_end(1:3) );
%wedge_type = path_struct{i}.anchor_type; %FOR NOW ALWAYS USE THE DEFAULT WEDGE TYPE
w = itaFiniteWedge( main_face_normal, opposite_face_normal, aperture_start, vertex_length, 'outer_edge' );
w.set_get_geo_eps( 1e-6 );
source_pos(1,:) = path_data{i-1}.interaction_point(1:3);
receiver_pos(1,:) = path_data{i+1}.interaction_point(1:3);
rho = ita_propagation_effective_source_distance( path_struct, i ); %effective distance from aperture point to source
last_pos_dirn(1,:) = path_data{i-1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %direction to the last source
eff_source_pos(1,:) = ( last_pos_dirn .* rho ./ norm(last_pos_dirn) ) + path_data{i}.interaction_point(1:3)';
r = ita_propagation_effective_target_distance( path_struct, i ); %effective distance from aperture point to receiver
next_pos_dirn(1,:) = path_data{i+1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %"receiver"
eff_receiver_pos(1,:) = ( next_pos_dirn .* r ./ norm(next_pos_dirn) ) + path_data{i}.interaction_point(1:3)';
if( w.point_outside_wedge( eff_source_pos ) == 0 ) %catch error if source is inside wedge
delay = -1;
return
end
%{
plot3([w.aperture_start_point(1),w.aperture_end_point(1)],[w.aperture_start_point(2),w.aperture_end_point(2)],[w.aperture_start_point(3),w.aperture_end_point(3)])
plot3([eff_source_pos(1),path_data{i}.interaction_point(1)],[eff_source_pos(2),path_data{i}.interaction_point(2)],[eff_source_pos(3),path_data{i}.interaction_point(3)])
plot3([eff_receiver_pos(1),path_data{i}.interaction_point(1)],[eff_receiver_pos(2),path_data{i}.interaction_point(2)],[eff_receiver_pos(3),path_data{i}.interaction_point(3)])
%}
%{
smallest_dist = 1000000;
n0 = 10000; %number of iterations
its = norm(w.aperture_start_point - w.aperture_end_point) / norm(w.aperture_direction); %number of aperture directions along aperture
for it = 0:n0
point_on_ap = w.aperture_start_point + (it*its/n0)*w.aperture_direction;
source_length = norm( point_on_ap - source_pos );
receiver_length = norm( eff_receiver_pos - point_on_ap );
total_dist = source_length + receiver_length;
if( total_dist < smallest_dist )
smallest_dist = total_dist;
ap_point = point_on_ap;
end
end
plot3(ap_point(1),ap_point(2),ap_point(3),'.k') %true closest point
ap_point2 = w.get_aperture_point( source_pos, receiver_pos );
ap_point3 = w.get_aperture_point( eff_source_pos, eff_receiver_pos );
ap_point4 = w.get_aperture_point2( source_pos, receiver_pos );
plot3(ap_point2(1),ap_point2(2),ap_point2(3),'.r') %point which should be predicted from aperture point method
plot3(ap_point3(1),ap_point3(2),ap_point3(3),'.b') %point which should be predicted from aperture point method
plot3(ap_point4(1),ap_point4(2),ap_point4(3),'.g') %point which should be predicted from aperture point method
%}
aperture_point = w.get_aperture_point2( source_pos, receiver_pos );
if( w.point_on_aperture( aperture_point ) == 0 )
warning('Skipping path, aperture point calculated not on the aperture');
continue
end
[~, D, A] = ita_diffraction_utd( w, eff_source_pos, eff_receiver_pos, f, c, aperture_point );
if( is_diff == 0 )
gain = gain * (A / rho);
is_diff = 1;
else
gain = gain * A;
end
frequency_mags = frequency_mags .* abs(D);
case 'inner_edge_diffraction'
source_pos(1,:) = path_data{i-1}.interaction_point(1:3);
receiver_pos(1,:) = path_data{i+1}.interaction_point(1:3);
opposite_face_normal(1,:) = path_data{i}.main_wedge_face_normal(1:3);
main_face_normal(1,:) = path_data{i}.opposite_wedge_face_normal(1:3);
aperture_start(1,:) = path_data{i}.vertex_start(1:3); %aperture point
vertex_length(1,:) = norm( path_data{i}.vertex_start(1:3) - path_data{i}.vertex_end(1:3) );
%wedge_type = path_struct{i}.anchor_type; %FOR NOW ALWAYS USE THE DEFAULT WEDGE TYPE
w = itaFiniteWedge( main_face_normal, opposite_face_normal, aperture_start, vertex_length, 'inner_edge' );
w.set_get_geo_eps( 1e-6 );
rho = ita_propagation_effective_source_distance( path_struct, i ); %effective distance from aperture point to source
last_pos_dirn(1,:) = path_data{i-1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %direction to the last source
eff_source_pos(1,:) = ( last_pos_dirn .* rho ./ norm(last_pos_dirn) ) + path_data{i}.interaction_point(1:3)';
r = ita_propagation_effective_target_distance( path_struct, i ); %effective distance from aperture point to receiver
next_pos_dirn(1,:) = path_data{i+1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %"receiver"
eff_receiver_pos(1,:) = ( next_pos_dirn .* r ./ norm(next_pos_dirn) ) + path_data{i}.interaction_point(1:3)';
if( w.point_outside_wedge( eff_source_pos ) == 0 ) %catch error if source is inside wedge
delay = -1;
return
end
aperture_point = w.get_aperture_point2( source_pos, receiver_pos );
if( w.point_on_aperture( aperture_point ) == 0 )
warning('Skipping path, aperture point calculated not on the aperture');
continue
end
[~, D, A] = ita_diffraction_utd( w, eff_source_pos, eff_receiver_pos, f, c, aperture_point );
if( is_diff == 0 )
gain = gain * (A / rho);
is_diff = 1;
else
gain = gain * A;
rho = ita_propagation_effective_source_distance( pp, i ); %effective distance from aperture point to source
last_pos_dirn = a_prev.interaction_point(1:3) - a_curr.interaction_point(1:3); %direction to the last source
eff_source_pos = ( last_pos_dirn .* rho ./ norm(last_pos_dirn) ) +a_curr.interaction_point(1:3);
r = ita_propagation_effective_target_distance( pp, i ); %effective distance from aperture point to receiver
next_pos_dirn = a_next.interaction_point(1:3) - a_curr.interaction_point(1:3); %"receiver"
eff_receiver_pos = ( next_pos_dirn .* r ./ norm(next_pos_dirn) ) + a_curr.interaction_point(1:3);
if ~w.point_outside_wedge( eff_source_pos ) %catch error if source is inside wedge
delay = -1; % @todo remove
valid = false;
return
end
%{
plot3([w.aperture_start_point(1),w.aperture_end_point(1)],[w.aperture_start_point(2),w.aperture_end_point(2)],[w.aperture_start_point(3),w.aperture_end_point(3)])
plot3([eff_source_pos(1),path_data{i}.interaction_point(1)],[eff_source_pos(2),path_data{i}.interaction_point(2)],[eff_source_pos(3),path_data{i}.interaction_point(3)])
plot3([eff_receiver_pos(1),path_data{i}.interaction_point(1)],[eff_receiver_pos(2),path_data{i}.interaction_point(2)],[eff_receiver_pos(3),path_data{i}.interaction_point(3)])
%}
%{
smallest_dist = 1000000;
n0 = 10000; %number of iterations
its = norm(w.aperture_start_point - w.aperture_end_point) / norm(w.aperture_direction); %number of aperture directions along aperture
for it = 0:n0
point_on_ap = w.aperture_start_point + (it*its/n0)*w.aperture_direction;
source_length = norm( point_on_ap - source_pos );
receiver_length = norm( eff_receiver_pos - point_on_ap );
total_dist = source_length + receiver_length;
if( total_dist < smallest_dist )
smallest_dist = total_dist;
ap_point = point_on_ap;
end
frequency_mags = frequency_mags .* abs(D);
case 'specular_reflection' %case for specular reflection
%path = 'C:\ITASoftware\Raven\RavenDatabase\MaterialDatabase';
%data = load(fullfile(path,'brickwall'));
%frequency_mags = frequency_mags .* FREQ_DATA_FOR_REFLECTION_SURFACE; %INSERT LOOKUP FIR FREQ DATA BASED ON VERTEX NUMBER
otherwise
error('Unrecognised anchor type');
end
end
frequency_mags = frequency_mags .* ita_atmospheric_absorption_factor( f, total_distance ); %flter contribution from atmospheric absorption
drawnow
if( is_diff == 0 ) %if there was no diffraction in path, apply 1/r distance law for gain
gain = 1/total_distance;
end
plot3(ap_point(1),ap_point(2),ap_point(3),'.k') %true closest point
ap_point2 = w.get_aperture_point( source_pos, receiver_pos );
ap_point3 = w.get_aperture_point( eff_source_pos, eff_receiver_pos );
ap_point4 = w.get_aperture_point2( source_pos, receiver_pos );
plot3(ap_point2(1),ap_point2(2),ap_point2(3),'.r') %point which should be predicted from aperture point method
plot3(ap_point3(1),ap_point3(2),ap_point3(3),'.b') %point which should be predicted from aperture point method
plot3(ap_point4(1),ap_point4(2),ap_point4(3),'.g') %point which should be predicted from aperture point method
%}
aperture_point = w.get_aperture_point2( source_pos, target_pos );
if ~w.point_on_aperture( aperture_point )
warning('Skipping path, aperture point calculated not on the aperture');
continue
end
[~, D, A] = ita_diffraction_utd( w, eff_source_pos, eff_receiver_pos, f, c, aperture_point );
if( is_diff == 0 )
gain = gain * (A / rho);
is_diff = 1;
else
gain = gain * A;
end
frequency_mags = frequency_mags .* abs(D);
case 'inner_edge_diffraction'
source_pos = a_prev.interaction_point(1:3);
target_pos = a_next.interaction_point(1:3);
w = ita_propagation_wedge_from_diffraction_anchor( a_curr );
w.set_get_geo_eps( 1e-6 );
rho = ita_propagation_effective_source_distance( pp, i ); %effective distance from aperture point to source
last_pos_dirn = a_prev.interaction_point(1:3) - a_curr.interaction_point(1:3); %direction to the last source
eff_source_pos = ( last_pos_dirn .* rho ./ norm(last_pos_dirn) ) + a_curr.interaction_point(1:3);
r = ita_propagation_effective_target_distance( pp, i ); %effective distance from aperture point to receiver
next_pos_dirn = a_next.interaction_point(1:3) - a_curr.interaction_point(1:3); %"receiver"
eff_receiver_pos = ( next_pos_dirn .* r ./ norm(next_pos_dirn) ) + a_curr.interaction_point(1:3);
if ~w.point_outside_wedge( eff_source_pos ) %catch error if source is inside wedge
delay = -1; % @todo remove
valid = false;
return
end
aperture_point = w.get_aperture_point2( source_pos, target_pos );
if ~w.point_on_aperture( aperture_point )
warning('Skipping path, aperture point calculated not on the aperture');
continue
end
if norm( eff_receiver_pos' - aperture_point ) < w.set_get_geo_eps
warning('Skipping a path where aperture point and receiver point coincide');
continue
end
[~, D, A] = ita_diffraction_utd( w, eff_source_pos, eff_receiver_pos, f, c, aperture_point );
if( is_diff == 0 )
gain = gain * (A / rho);
is_diff = 1;
else
gain = gain * A;
end
frequency_mags = frequency_mags .* abs(D);
case 'specular_reflection' %case for specular reflection
%path = 'C:\ITASoftware\Raven\RavenDatabase\MaterialDatabase';
%data = load(fullfile(path,'brickwall'));
%frequency_mags = frequency_mags .* FREQ_DATA_FOR_REFLECTION_SURFACE; %INSERT LOOKUP FIR FREQ DATA BASED ON VERTEX NUMBER
otherwise
error('Unrecognised anchor type');
end
delay = total_distance / c;
end
%% Determine DSP coefficients / path data
frequency_mags = frequency_mags .* ita_atmospheric_absorption_factor( f, total_distance ); %flter contribution from atmospheric absorption
if( is_diff == 0 ) %if there was no diffraction in path, apply 1/r distance law for gain
gain = 1 / total_distance;
end
delay = total_distance / c;
drawnow % can be removed?
end
\ No newline at end of file
function [ w ] = ita_propagation_wedge_from_diffraction_anchor( a )
% ita_propagation_wedge_from_diffraction_anchor Generates a wedge based on the data from a
% diffraction anchor point
%
% Example: [ wedge ] = ita_propagation_wedge_from_diffraction_anchor( anchor )
%
if ~isfield( a, 'anchor_type' )
error 'Anchor does not contain a type description'
end
if strcmpi( a.anchor_type, 'outer_edge_diffraction' )
edge_type = 'outer_edge';
elseif strcmpi( a.anchor_type, 'inner_edge_diffraction' )
edge_type = 'inner_edge';
else
error 'Invalid anchor type'
end
main_face_normal = a.main_wedge_face_normal( 1:3 );
opposite_face_normal = a.opposite_wedge_face_normal( 1:3 );
aperture_start = a.vertex_start( 1:3 );
vertex_length = norm( a.vertex_start - a.vertex_end );
w = itaFiniteWedge( main_face_normal, opposite_face_normal, aperture_start, vertex_length, edge_type );
end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment