WIP

matlab/VA_example_outdoor_acoustics.m

@@ -4,13 +4,10 @@
 % the user. Also the generic path prototype rendering module(s) has to record the output
 % to hard drive.
 
-free_field_only = false;
-direct_field_only = false; % with shadowing of edge
-diffracted_field_only = false;
-
-assert( ~( free_field_only && diffracted_field_only ) )
-assert( ~( free_field_only && direct_field_only ) )
-assert( ~( diffracted_field_only && direct_field_only ) )
+free_field_only = false; % Now occlusion and no reflections
+direct_field_component = true;
+reflected_field_component = true;
+diffracted_field_component = true;
 
 receiver_pos = [ 3 0 3 ]; % OpenGL coordinates
 f = ita_ANSI_center_frequencies;
@@ -61,6 +58,8 @@ catch
     
 end
 
+k = 2 * pi * freq ./ c;
+
 %% Example for a synchronized scene update & audio processing simulation/auralization
 
 timestep = 128 / 44100; % here: depends on block size and sample rate
@@ -69,17 +68,18 @@ if ~dry_run
     va.set_core_clock( 0 );
 end
 
-spatialstep = 0.1;
+spatialstep = 0.1 / 2;
 disp( [ 'Resulting sound source speed: ' num2str( spatialstep / timestep ) ' m/s' ] )
 
-numsteps = 3400;
+numsteps = 3400 * 2;
 disp( [ 'Simulation result duration: ' num2str( numsteps * timestep ) ' s' ] )
 
-x = linspace( -1, 1, numsteps ) * 50; % motion from x = -50m to x = 50m
-
-n0_of_paths = 2;
+x = linspace( -1, 1, numsteps ) * 20; % motion from x = -50m to x = 50m
 
-H_log = []; % Diffracted field
+H_direct_log = [];
+H_diffracted_log = [];
+H_reflected_log = [];
+H_total_field_log = [];
 IL_log = []; % Insertion loss
 
 h = waitbar( 0, 'Hold on, running auralization' );
@@ -89,15 +89,14 @@ for n = 1:length( x )
     source_pos = [ x( n ) 0 -3 ]; % OpenGL coordinates    
     if ~dry_run
         va.set_sound_source_position( S, source_pos );
-    end
-        
+    end        
     
     %% Direct sound path
     
     % Manually create direct sound path and diffracted sound path
     shadow_zone = ita_diffraction_shadow_zone( w, source_pos, receiver_pos ); %is receiver in shadow zone?
     if n > 1 && last_shadow_zone ~= shadow_zone
-        fprintf( 'Shadow zone has been crossed on frame %i\n', n ) 
+        fprintf( 'Shadow zone boundary has been crossed on frame %i (time %.2f)\n', n, n * timestep ) 
     end
     last_shadow_zone = shadow_zone;
     
@@ -120,13 +119,15 @@ for n = 1:length( x )
         prop_path_direct.audible = ~shadow_zone; % Direct field: may be occluded
     end
     
-    H_direct = exp( -1i * 2 * pi * f / c / distance ) / distance;
+    H_direct = 1 ./ distance .* exp( -1i .* k .* distance );
+    
     
     %% Diffracted sound path
     
     apex = w.get_aperture_point( source_pos, receiver_pos );
     detour = norm( source_pos - apex ) + norm( apex - receiver_pos );
-    [ H, D, A ] = ita_diffraction_utd( w, source_pos, receiver_pos, f, c );    
+    H_diffracted = ita_diffraction_utd( w, source_pos, receiver_pos, f, c );
+    [ ~, D, A ] = ita_diffraction_utd( w, source_pos, receiver_pos, f, c );
     rho = norm( apex - source_pos );
         
     prop_path_diffracted = struct();    
@@ -135,28 +136,63 @@ for n = 1:length( x )
     prop_path_diffracted.identifier = 'diffracted_path';
     prop_path_diffracted.delete = false;
     prop_path_diffracted.position = apex;
-    prop_path_diffracted.audible = true;
     prop_path_diffracted.delay = detour / c;
     prop_path_diffracted.frequencies = f;
     prop_path_diffracted.gain = A / rho;
     prop_path_diffracted.frequency_magnitudes = abs( D );
+    prop_path_diffracted.audible = true;
+    
     
+    %% Reflected sound path
+    
+    reflection_zone = ita_diffraction_reflection_zone( w, source_pos, receiver_pos, true );
+    if n > 1 && last_reflection_zone ~= reflection_zone
+        fprintf( 'Reflection zone boundary has been crossed on frame %i (time %.2f)\n', n, n * timestep ) 
+    end
+    last_reflection_zone = reflection_zone;
+    
+    distance_source_main_face = dot( w.main_face_normal, source_pos - w.location );
+    assert( ~reflection_zone || distance_source_main_face >= 0 )
+    
+    source_pos_image = source_pos - 2 * w.main_face_normal * distance_source_main_face;
+    
+    distance_image = norm( source_pos_image - receiver_pos );
+    
+    prop_path_reflected = struct();
+    prop_path_reflected.identifier = 'reflected_path'; 
+    prop_path_reflected.gain = 1 / distance_image;
+    prop_path_reflected.delay = distance_image / c; %delay from sound emitting from source to being received at listener
+    prop_path_reflected.frequencies = f; %Frequencies corresponding to the mags
+    values = ones(1,31);
+    prop_path_reflected.frequency_magnitudes = values;
+    prop_path_reflected.source = S; %sound source ID
+    prop_path_reflected.receiver = L; %sound receiver ID
+    prop_path_reflected.delete = false; %set to true when an existing path should be deleted
+    prop_path_reflected.position = source_pos_image;
+    prop_path_reflected.audible = reflection_zone; % Reflected field: may be inaudible
+    
+    R = 1;
+    H_reflected = R ./ distance_image .* exp( -1i .* k ./ distance_image );
+        
     
     %% Update wave fronts in renderer
     paths_update = struct();
-    if ~diffracted_field_only
+    if direct_field_component || direct_field_only
         paths_update.direct_path = prop_path_direct;
     end
-    if ~free_field_only && ~direct_field_only
+    if ~direct_field_only && diffracted_field_component
         paths_update.diffracted_path = prop_path_diffracted;
     end
+    if ~direct_field_only && reflected_field_component
+        paths_update.reflected_path = prop_path_reflected;
+    end
     
     if ~dry_run
         
         rends = va.get_rendering_modules();
-        for rend = 1:numel( rends )
-            if strcmpi( rend.class, 'BinauralOutdoorNoise' )
-                va.set_rendering_module_parameters( rend.id, paths_update );
+        for rend_idx = 1:numel( rends )
+            if strcmpi( rends( rend_idx ).class, 'BinauralOutdoorNoise' )
+                va.set_rendering_module_parameters( rends( rend_idx ).id, paths_update );
             end
         end
 
@@ -171,12 +207,22 @@ for n = 1:length( x )
     
     waitbar( n / numsteps )
     
-    H_log = [ H_log; H ];
+    H_direct_log = [ H_direct_log; H_direct ];
+    H_diffracted_log = [ H_diffracted_log; H_diffracted ];
+    H_reflected_log = [ H_reflected_log; H_reflected ];
+    
     if shadow_zone
-        IL_log = [ IL_log; H ./ H_direct ];
+        H_total_field = H_diffracted;
     else
-        IL_log = [ IL_log; ( H + H_direct ) ./ H_direct ];
+        if reflection_zone
+            H_total_field = H_reflected + H_diffracted + H_direct;
+        else
+            H_total_field = H_diffracted + H_direct;
+        end
     end
+    
+    H_total_field_log = [ H_total_field_log; H_total_field ];
+    IL_log = [ IL_log; H_total_field ./ H_direct ];
      
 end
 close( h )
@@ -184,4 +230,7 @@ close( h )
 if ~dry_run    
     va.disconnect
     disp( 'Stop VA to export simulation results from rendering module(s)' )    
+else
+    figure
+    plot( db( IL_log ) )
 end