Improving UTD diffraction calculation and adding more propagation functions

parent 54b14e3f
......@@ -137,9 +137,9 @@ classdef itaInfiniteWedge
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;
b = false;
if ( norm( obj.main_face_normal ) - 1 ) < eps && ( norm( obj.opposite_face_normal ) -1 ) < eps
b = true;
end
end
......
function [ diffr_field, D ] = ita_diffraction_utd( wedge, source_pos, receiver_pos, frequency_vec, speed_of_sound )
function [ diffr_field, D, A ] = 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 )
%
......@@ -71,9 +71,9 @@ 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 );
main_face = wedge.point_facing_main_side( S );
alpha_i = wedge.get_angle_from_point_to_wedge_face( S, main_face );
alpha_d = wedge.get_angle_from_point_to_wedge_face( R, main_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
......@@ -84,25 +84,26 @@ k = (2 * pi) ./ lambda; % Wavenumber
% Diffraction coefficient D
assert( all( rho + r ~= 0 ) && all( r ~= 0 ) );
L = repmat( ( ( rho .* r ) ./ ( rho + r ) ) .* ( sin( theta_i ) ).^2, 1, numel( frequency_vec ) );
L = ( ( rho .* r ) ./ ( rho + r ) ) .* ( sin( theta_i ) ).^2; % -> distance dependency
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 ) );
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 = kawai_approx_fresnel( k .* L .* a1 );
F1 = kawai_approx_fresnel( k .* L .* a1 ); % -> frequency dependent
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;
......@@ -141,7 +142,7 @@ if any( singularities )
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, :);
term1(~mask1, :) = Cot1(~mask1, :) .* F1(~mask1, :); % -> frequency dependent
term2(~mask2, :) = Cot2(~mask2, :) .* F2(~mask2, :);
term3(~mask3, :) = Cot3(~mask3, :) .* F3(~mask3, :);
term4(~mask4, :) = Cot4(~mask4, :) .* F4(~mask4, :);
......@@ -152,15 +153,15 @@ else
s = -1;
end
D = D_factor .* ( term1 + term2 + s * ( term3 + term4 ) );
D = ( D_factor .* ( term1 + term2 + s * ( term3 + term4 ) ) )'; % -> frequency dependent
%% Combined diffracted sound field filter at receiver
H_i = 1 ./ rho .* exp( -1i * k .* rho ); % Consideration of transfer path from source to apex point
A = repmat( sqrt( rho ./ ( r .* ( rho + r ) ) ), 1, numel( frequency_vec ) ); % Amplitude
H_i = 1 ./ rho .* exp( -1i * k .* rho )'; % Consideration of transfer path from source to apex point
A = repmat( sqrt( rho ./ ( r .* ( rho + r ) ) ), 1, numel( frequency_vec ) )'; % Amplitude
diffr_field = ( H_i .* D .* A .* exp( -1i .* k .* r ) )';
diffr_field = ( H_i .* D .* A .* exp( -1i .* k' .* r ) );
end
......@@ -168,20 +169,20 @@ end
%% Auxiliary functions
% N+ function
% N+ function (plus)
function N = N_p( n, beta )
N = zeros( numel( beta ), 1 );
N( beta > pi * ( 1 - n ) ) = 1;
end
% N- function
% N- function (minus)
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
% Signum function
function res = sgn(x)
if all( size(x) == 0 )
res = 1;
......@@ -191,6 +192,7 @@ function res = sgn(x)
res( x <= 0 ) = -1;
end
% Approximation of the Fresnel integral by Kawaii et al.
function Y = kawai_approx_fresnel( X )
if any( X < 0 )
error( 'No negative values for Kawai approximation of Fresnel integral allowed' )
......@@ -200,4 +202,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
%% 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, c );
% % 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' )
%% Simple example
n1 = [ 1 0 1 ] ./ sqrt( 2 );
n2 = [ -1 0 1 ] ./ sqrt( 2 );
loc = [ 0 0 0 ]; % edge along Y axis
src = [ -1 0 -0.5 ]; % meter
rcv = [ 1 0 -0.75 ]; % meter
w = itaInfiniteWedge( n1, n2, loc ); % wedge object
freq = linspace( 20, 20000, 1000); % Hz
c = 343; % m/s
[ d_tf, D, A ] = ita_diffraction_utd( w, src, rcv, freq, c );
figure
semilogx( freq, 10*log10( abs( d_tf ) ) )
hold on
semilogx( freq, 10*log10( abs( D ) ) )
semilogx( freq, 10*log10( abs( A ) ) )
legend( 'transmission', 'diffraction coefficient', 'amplitude coefficient' )
title( 'UTD diffraction transfer function (simple example scene)' )
xlabel( 'Frequency / Hz' )
ylabel( 'Amplitude / dB' )
%% Distance dependency
far_distance_scaling = 10;
[ d_tf_far, D_far, A_far ] = ita_diffraction_utd( w, far_distance_scaling * src, far_distance_scaling * rcv, freq, c );
figure
semilogx( freq, [ 10*log10( abs( d_tf ) ), 10*log10( abs( d_tf_far ) ) ] )
hold on
semilogx( freq, [ 10*log10( abs( D ) ), 10*log10( abs( D_far ) ) ], 'lineWidth', 2 )
semilogx( freq, [ 10*log10( abs( A ) ), 10*log10( abs( A_far ) ) ] )
legend( 'transmission (near)', 'transmission (far)', ...
'diffraction coefficient (near)', 'diffraction coefficient (far)', ...
'amplitude coefficient (near)', 'amplitude coefficient (far)' )
title( 'UTD diffraction transfer function (distance dependency)' )
xlabel( 'Frequency / Hz' )
ylabel( 'Amplitude / dB' )
function phase_by_delay = ita_propagation_delay( distance, c, fs, fft_degree )
%ITA_RPOPAGATION_SPREADING_LOSS Calculates spreading loss for different
%wave types for a given distance (straight line between emitter center
%point and sensor point)
if nargin ~= 4
error 'Expecting 4 arguments'
end
if distance <= 0
error 'Distance cannot be zero or negative'
end
delay_tf = itaAudio();
delay_tf.samplingRate = fs;
delay_tf.fftDegree = fft_degree;
f = delay_tf.freqVector;
lambda = c ./ f; % Wavelength
k = 2 * pi ./ lambda; % Wavenumber
phase_by_delay = exp( -1i .* k .* distance ); % Note: DC value set to 1
end
......@@ -46,7 +46,7 @@ w = itaInfiniteWedge( n1, n2, loc );
switch( diffraction_model )
case 'utd'
[ H, D ] = ita_diffraction_utd( w, effective_source_position( 1:3 )', effective_receiver_position( 1:3 )', specrefl_tf.freqVector( 2:end ), c );
[ ~, D, ~ ] = ita_diffraction_utd( w, effective_source_position( 1:3 )', effective_receiver_position( 1:3 )', specrefl_tf.freqVector( 2:end ), c );
specrefl_tf.freqData = [ 0 D ]';
case 'maekawa'
......
function spreading_loss_factor = ita_propagation_spreading_loss( distance, wave_type )
%ITA_RPOPAGATION_SPREADING_LOSS Calculates spreading loss for different
%wave types for a given distance (straight line between emitter center
%point and sensor point)
if nargin == 1
wave_type = 'spherical';
end
if distance <= 0
error 'Distance cannot be zero or negative'
end
switch( wave_type )
case 'plain'
spreading_loss_factor = 1;
case 'line'
spreading_loss_factor = 1 / sqrt( distance );
otherwise
spreading_loss_factor = 1 / distance;
end
end
\ No newline at end of file
......@@ -24,13 +24,15 @@ prop_tfs.freqData = ones( prop_tfs.nBins, N );
lambda = ita_speed_of_sound ./ prop_tfs.freqVector( 2:end ); % Wavelength
k = (2 * pi) ./ lambda; % Wavenumber
r = ita_propagation_path_length( propagation_path );
if r / ita_speed_of_sound > prop_tfs.trackLength
distance = ita_propagation_path_length( propagation_path );
if distance / ita_speed_of_sound > prop_tfs.trackLength
error 'Propagation path length too long, increase fft degree to generate transfer function for this propagation path'
end
freq_data_linear = [ 0 exp( 1i .* k .* r ) ]' ./ r; % 1/r spherical spreading / distance law & propagation delay
phase_by_delay = ita_propagation_delay( distance, ita_speed_of_sound, fs, fft_degree );
spreading_loss = ita_propagation_spreading_loss( distance, 'spherical' );
freq_data_linear = phase_by_delay .* spreading_loss; % oder quadratisch?!
for m = 1 : N
......
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