VA_example_ppa.m 14.1 KB
Newer Older
1
2
3
4
5
6
7
8
%% VA offline simulation/auralization example for outputs of the propagation path algorithm

% Requires VA to run with a virtual audio device that can be triggered by
% the user. Also the generic path prototype rendering module(s) has to record the output
% to hard drive.

% Requires ITA-Toolbox, obtain from http://www.ita-toolbox.org

9
ppa_folder = '../../../dist/win32-x64.vc12/bin/UrbanTrajectory';
10
ppa_diffraction_model = 'utd';
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

%% Preps

file_listing = dir( strcat( ppa_folder, '/*.json' ) );

c = ita_propagation_speed_of_sound_air_approx( 20 );
f = ita_ANSI_center_frequencies;


%% VA

va = VA( 'localhost' );
va.add_search_path( 'D:/Users/andrew/dev/ITASuite/VA/VACore/data' );
va.add_search_path( 'D:/Users/andrew/dev/ITASuite/VA/VAMatlab/matlab' );
va.add_search_path( 'C:/dev/VA/VACore/data' );
va.add_search_path( 'D:/Users/stienen/dev/VA/VACore/data' );
va.add_search_path( 'C:\ITASoftware\Raven\RavenDatabase\SoundDatabase' );

c = va.get_homogeneous_medium_sound_speed();
f = ita_ANSI_center_frequencies;

R = va.create_sound_receiver( 'PPA_sensor' );
H = va.create_directivity_from_file( 'ITA_Artificial_Head_5x5_44kHz_128.v17.ir.daff' );
va.set_sound_receiver_directivity( R, H );

S = va.create_sound_source( 'PPA_emitter' );
%X = va.create_signal_source_buffer_from_file( 'WelcomeToVA.wav' );
X = va.create_signal_source_buffer_from_file( 'Full song - Vulfpeck - Dean Town.wav' );

va.set_signal_source_buffer_playback_action( X, 'play' )
va.set_signal_source_buffer_looping( X, true );
va.set_sound_source_signal_source( S, X )

44
45
46
path_identifiers = containers.Map('KeyType','char','ValueType','uint16'); %a map linking the path number to name
path_status = cell(1); %cell array indexed by path number containing a struct with onfo on path
path_count = 0; %The total number of paths, NOTE -NOT a count of paths in current frame, if paths are deleted, they are not removed from this count
47
48
49
50
51
52
53
   
%% Run synchronized scene update & audio processing simulation/auralization

frame_rate = 128 / 44100; % here: depends on block size and sample rate
manual_clock = 0;
va.set_core_clock( 0 );

Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
54
N = numel( file_listing );
55
56
57
58
59
disp( [ 'Simulation result duration: ' num2str( N * frame_rate ) ' s' ] )

h = waitbar( 0, 'Hold on, running auralization' );

% Iterate over frames
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
60
for n = 1:N
61
62
63

    % Load propagation paths for current frame
    ppa_file_path = fullfile( file_listing( n ).folder, file_listing( n ).name );
64
    pps = ita_propagation_load_paths( ppa_file_path ); %pps = struct containing the current time frame
65
66
67
    
    % Update source (first anchor)
    source_pos = pps( n ).propagation_anchors{ 1 }.interaction_point; % OpenGL coordinates?! -> transform using 
68
    source_pos_OpenGL = ita_matlab2openGL( source_pos(1:3)' );
69
70
71
72
    va.set_sound_source_position( S, source_pos_OpenGL );
    
    % Update receiver (last anchor)
    receiver_pos = pps( n ).propagation_anchors{ end }.interaction_point; % OpenGL coordinates?! -> transform using 
73
    receiver_pos_OpenGL = ita_matlab2openGL( receiver_pos(1:3)' );    
74
75
76
77
78
79
    va.set_sound_receiver_position( R, source_pos_OpenGL );
    
    paths_update = struct();
    
    for p = 1:numel( pps )
        
80
        pp = pps( p ); % Propagation path
81
82
83
84
85
86
        
        pu = struct(); % Path update
        
        % Assemble DSP settings (gain, delay & filter coefficients)
        pu.source = S;
        pu.receiver = R;
87
88
        path_name = get_path(pp);
        pu.path = path_name;
89
        
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        exists = isKey( path_identifiers, path_name );
        if( ~exists ) %if path does not already exist, add a new path with a unique int ID
            path_count = path_count + 1;
            ID = path_count;
            new_ID = struct;
            new_ID.ID = ID;
            new_ID.updated = true;
            new_ID.deleted = false;
            new_ID.name = path_name;
            path_status{ID} = new_ID;
            path_identifiers(path_name) = ID;
        else
            ID = path_identifiers(path_name);
            path_status{ID}.updated = true; %if path already exists, mark that it has been updated
        end
105
        
106
107
        pu.ID = ID;
        %pu.path = strcat( 'direct_', num2str(S), '_to_', num2str(R) ); %sound path ID
108
        
109
110
111
112
113
114
        [frequency_mags, gain, delay] = ita_get_propagation_path_data( pp, f, c);
        
        pu.gain = gain;
        pu.delay = delay;
        pu.frequency_magnitudes = frequency_mags;
        pu.frequencies = f;
115
        pu.position = pp.propagation_anchors{ end-1 }.interaction_point; % next to last anchor
116
117
       
        pu.delete = false;
118
119
        pu.audible = true;
        
120
        paths_upate.(strcat('prop_path_',num2str(p))) = pu; %create a new field in the update struct containing the new path data
121
122
123
        
    end
    
124
125
126
127
128
129
130
131
132
133
134
135
    for i = 1:numel(path_status) %add any paths which need to be deleted to the update struct
        if( (path_status{i}.updated == false) && (path_status{i}.deleted == false) )
            p = p + 1;
            dpu = struct();
            dpu.ID = i;
            dpu.delete = true;          
            paths_upate.(strcat('prop_path_',num2str(p))) = dpu; %set a new field in the path update struct telling the renderer to delete the path which no longer appears
            path_status{i}.deleted = true;
        end
        path_status{i}.updated = false;
    end
    
136
137
138
139
140
141
142
143
144
145
    % Update all propagation paths
    va.set_rendering_module_parameters( 'MyBinauralOutdoorNoise', paths_upate );
    
    % Increment core clock
    manual_clock = manual_clock + frame_rate;
    va.call_module( 'manualclock', struct( 'time', manual_clock ) );
    
    % Process audio chain by incrementing one block
    va.call_module( 'virtualaudiodevice', struct( 'trigger', true ) );
    
146
    waitbar( n / N )
147
148
149
150
151
152
153
     
end
close( h )

va.disconnect

disp( 'Stop VA to export simulation results from rendering module(s)' )
154
155
156
157
158



%%
function path_name = get_path( path_struct )
159
160
    path_data = path_struct.propagation_anchors;
    num_interactions = numel( path_data );
161
162
    path_name = '';
    for i = 1:num_interactions
163
        switch path_data{i}.anchor_type
164
            case 'source'
165
                ID = path_data{i}.name;
166
167
                path_name = strcat( path_name, 'SRC', num2str(ID), '_' );
            case 'receiver'
168
                ID = path_data{i}.name;
169
170
                path_name = strcat( path_name, 'RCVR', num2str(ID) );
            case 'specular_reflection'
171
                ID = path_data{i}.polygon_id;
172
173
                path_name = strcat( path_name, 'SR', num2str(ID), '_' );
            case 'outer_edge_diffraction'
174
175
                ID1 = path_data{i}.main_wedge_face_id;
                ID2 = path_data{i}.opposite_wedge_face_id;
176
177
178
179
180
181
                path_name = strcat( path_name, 'EDMF', num2str(ID1), 'OF', num2str(ID2), '_' );
        end
    end
end

function [frequency_mags, gain, delay] = ita_get_propagation_path_data( path_struct, f, c )
182
    path_data = path_struct.propagation_anchors;
183
184
185
    gain = 1;
    frequency_mags = ones(1,length(f));
    total_distance = 0;
186
    is_diff = 0; %set to 1 whenever a diffraction is encountered
187
    
188
189
190
191
192
    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
    end
193
    
194
195
196
    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));
197
198
199
200
        total_distance = total_distance + segment_distance;
        
        switch anchor_type
            case 'outer_edge_diffraction' %case for diffraction
201
202
203
204
                main_face_normal = path_data{i}.main_wedge_face_normal(1:3);
                opposite_face_normal = path_data{i}.opposite_wedge_face_normal(1:3);
                aperture_start = path_data{i}.vertex_start(1:3); %aperture point    
                vertex_length = norm( path_data{i}.vertex_start(1:3) - path_data{i}.vertex_end(1:3) );
205
206
                %wedge_type = path_struct{i}.anchor_type; %FOR NOW ALWAYS USE THE DEFAULT WEDGE TYPE
                
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
                w = itaFiniteWedge( main_face_normal, opposite_face_normal, aperture_start, vertex_length );
                
                source_pos = path_data{i-1}.interaction_point(1:3);
                receiver_pos = path_data{i+1}.interaction_point(1:3);
                rho = ita_propagation_effective_source_distance( path_struct, i ); %effective distance from aperture point to source
                source_dist = rho;
                last_pos_dirn = path_data{i-1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %direction to the last source
                eff_source_pos = ( last_pos_dirn .* source_dist ./ 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
                target_dist= r;
                next_pos_dirn = path_data{i+1}.interaction_point(1:3) - path_data{i}.interaction_point(1:3); %"receiver"
                eff_receiver_pos = ( next_pos_dirn .* target_dist ./ norm(next_pos_dirn) ) + path_data{i}.interaction_point(1:3);
                
                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
246
247

                
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
                [~, D, ~] = ita_diffraction_utd( w, source_pos, receiver_pos, f, c );    
                
                A = sqrt( rho ./ ( r .* ( rho + r ) ) ); % attenuation factor at receiver
                
                if( is_diff == 0 )
                    gain = gain * (A / source_dist);
                    is_diff = 1;
                else
                    gain = gain * A;
                end
                frequency_mags = frequency_mags .* abs(D); 
            case 'inner_edge_diffraction'
                opposite_face_normal = path_data{i}.main_wedge_face_normal(1:3);
                main_face_normal = path_data{i}.opposite_wedge_face_normal(1:3);
                aperture_start = path_data{i}.vertex_start(1:3); %aperture point    
                vertex_length = 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 );     
                
                rho = ita_propagation_effective_source_distance( path_struct, i ); %effective distance from aperture point to source
                source_dist = rho;
                r = ita_propagation_effective_target_distance( path_struct, i ); %effective distance from aperture point to receiver
                target_dist= r;
                
                [~, D, ~] = ita_diffraction_utd( w, source_pos, receiver_pos, f, c );    
                
                A = sqrt( rho ./ ( r .* ( rho + r ) ) ); % attenuation factor at receiver
                
                if( is_diff == 0 )
                    gain = gain * (A / source_dist);
                    is_diff = 1;
                else
                    gain = gain * A;
                end
                frequency_mags = frequency_mags .* abs(D);               
284
            case 'specular_reflection' %case for specular reflection
285
286
                %path = 'C:\ITASoftware\Raven\RavenDatabase\MaterialDatabase';
                %data = load(fullfile(path,'brickwall'));
287
288
289
290
291
292
                %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
293
294
295
    drawnow
    if( is_diff == 0 ) %if there was no diffraction in path, apply 1/r distance law for gain
        gain = 1/total_distance;
296
297
298
    end
    delay = total_distance / c;
end
299
300
301
302
303
304
305

%{
function absorption = readMaterial(name)
    path = 'C:\ITASoftware\Raven\RavenDatabase\MaterialDatabase';
    data = dlmread( fullfile(path,name), ',' );
end
%}