tp@8: function simulateBinauralSignals(inputstruct)
tp@8: %Creates .wav file of 2sec long binaural signals for specific board dimensions and gaussian white noise stimulus
tp@8: %
tp@8: % INPUTS:
tp@8: % The input must be a structure with fields:
tp@8: % .dist (1st input argument) is the distance in meters
tp@8: % .azim (2nd input argument) is the azimuth in degrees (0 means straight ahead, positive angles are left and negative
tp@8: % are right)
tp@8: % .orient (3rd input argument) must be either 'horz' or 'angled' corresponding to flat and angled descriptions in
tp@8: % Papadopoulos et al. BSPC 2011.
Chris@10: % .dirweight (4th input argument) must be a nonnegative real scalar determining what is the relative weight of the
tp@8: % emission path to the echo path (i.e. due to directivity focus in the frontal direction of the source, the emission
tp@8: % which is directed upwards and backwards in our specific geometry is significantly attenuated, typically by factor in
tp@8: % the vicinity of 0.2)
Chris@12: % .outdir (5th input argument) must be the directory in which the wave file
Chris@10: % is to be written. It may be an absolute path, or relative to the current
Chris@12: % directory. Temporary files will be written in the same directory.
Chris@15: % .outname (6th input argument) must be the name of the wave file to write,
Chris@15: % within the outdir, without the .wav suffix (which will be added).
Chris@15: % .indir (7th input argument) must be the directory from which input data
Chris@15: % such as IR files are to be read. It may be an absolute path, or relative
Chris@15: % to the current directory.
tp@8: % 
tp@8: 
tp@8: %% Internal workings description
tp@8: %
tp@8: % Center of spherical coordinates system is taken to be the center of head. The azimuth, elevation and distance
tp@8: % coordinates are as in matlabs sph2cart function (azimuth and elevation are angular displacements from the positive
tp@8: % x-axis and from the x-y plane, respectively) with positive x axis taken to be extending to the right of the head from
tp@8: % top view, positive y axis to be extending forward of the head in top view and positive z axis to be extending upwards
Chris@14: % in top view.http://getpermanent.com/index
tp@8: % 
tp@8: % Board dimensions are defined as
tp@8: % params.board_size_x (width in meters for the 'center' orientation of  Papadopoulos et al. BSPC 2011)
tp@8: % params.board_size_y (depth in meters for the 'center' orientation of  Papadopoulos et al. BSPC 2011)
tp@8: % params.board_size_z (height in meters for the 'center' orientation of  Papadopoulos et al. BSPC 2011)
tp@8: % 
tp@8: % Board center position is defined as
tp@8: % params.board_distance (distance in meters from coords origin to center of board following coordinate system described
tp@8: %                       above), can be row vector.
tp@8: % params.board_azim     (azimuth in radians of the center of the board following coordinate system described above),
tp@8: %                       can be row vector.
tp@8: % params.board_elev     (elevation in radians of the center of the board following coordinate system described above),
tp@8: %                       must be scalar
tp@8: % 
tp@8: % params.board_orientation  scalar or 1x2 cell with elements 'horz' or/and 'angled'
tp@8: % Board is taken to always be vertically positioned (i.e. with its width-height plane vertical to the y coordinate) and
tp@8: % two cases of orientation are considered: horizontal and angled corresponding to flat and angled descriptions in
tp@8: % Papadopoulos et al. BSPC 2011.
tp@8: % 
tp@8: % params.source_down    non-negative scalar in meters
tp@8: % params.source_front   non-negative scalar in meters
tp@8: % The source is assumed to always be in front of the chest and below the head. (params.source_down and
tp@8: % params.source_front cannot both be 0)
tp@8: %
tp@8: %
tp@8: % params.horz_disamb    string which can be either 'BSPC2011' or 'collocated'. If equal to 'BSPC2011' then board distance
tp@8: %                       for the horizontal (flat) case is as described in Papadopoulos et al. BSPC 2011, i.e. the board
tp@8: %                       distance is the distance between the centre of the head and the PLANE OF THE BOARD.
tp@8: %                       Alternatively, in the 'collocated' case, the board distance in the horizontal (flat) case is taken
tp@8: %                       as the distance between the centre of the head and the centre of the board
tp@8: % 
tp@8: % NOTE: in the 'BSPC2011' case of params.horz_disamb, it is easy to see that the geometry is ill-specified for board
tp@8: %       positions far away from the median plane. 
tp@8: % 
tp@8: % OUTPUTS:
tp@8: % 
tp@8: % --- simulation_data is a structure with the following fields: 
tp@8: % - simulation_data.echo is a cell of dimensions 
tp@8: %       length(params.board_dist) x length(params.board_azim) x length(params.board_orientation)
tp@8: % each element of which is a 2-row matrix with the left (top row) and right (bottom row) ear responses corresponding the
tp@8: % echo only.
tp@8: % 
tp@8: % - simulation_data.emission is a cell of dimensions 
tp@8: %       length(params.board_dist) x length(params.board_azim) x length(params.board_orientation)
tp@8: % each element of which is a 2-row matrix with the left (top row) and right (bottom row) ear responses corresponding the
tp@8: % direct source-to-receiver path only. (NOTE: source term (binaural_irs_emission) does not need to be computed for each
tp@8: % different board distance/azimuth/orientation as it is the same irs but containing a different number of trailing zero
tp@8: % samples. I include this redundancy because it simplified the computation a bit and also because this way the emission
tp@8: % part is always equal in length with the corresponding echo part and can be added easier)
tp@8: % 
tp@8: % Adding any given
tp@8: %       directivity_weighting * simulation_data.emission{ii,jj,kk}(1,:) + simulation_data.echo{ii,jj,kk}(1,:)
tp@8: % for any given ii, jj, kk will give the total left ear IR and the same for right ear by
tp@8: %       directivity_weighting * simulation_data.emission{ii,jj,kk}(2,:) + simulation_data.echo{ii,jj,kk}(2,:)
tp@8: % 
tp@8: % - simulation_data.params is the stucture params described above
tp@8: % - simulation_data.board_coords is an array of dimensions
tp@8: %       8 x 3 x length(params.board_dist) x length(params.board_azim) x length(params.board_orientation)
tp@8: % containing the x y z coordinates (2nd dimension) of the 8 edges (1st dimension) of the board geometries for different
tp@8: % board distances, azimuths and orientations.
tp@8: % - simulation_data.azerr and simulation_data.elerr are arrays of dimensions
tp@8: %       2 x 1+length(params.board_azim))
tp@8: % which contain the error (in degrees) between the azimuth and elevation respectively of the HRTFs that are used (as
tp@8: % existing in the CIPIC spherical grid) and the actual board center azimuth and elevation.
tp@8: % 
tp@8: % 
tp@8: % 
tp@8: % DEVELOPMENT NOTES
tp@8: % 
tp@8: % TODO: include case of multiple boards and of non-vertically oriented boards (i.e. all cases of pitch, roll, yaw for
tp@8: %       board orientation)
tp@8: % 
tp@8: % TODO: include geometrical description for all possible positions of the source
tp@8: % 
tp@8: % TODO: include validateattributes for all parameters
tp@8: % 
tp@8: % TODO: include parameter controls for specorder (now set to 1), difforder (now set to 1), elemsize (now set to [1]) and
tp@8: % nedgesubs (now set to 2)
tp@8: %
tp@8: % 
tp@8: 
tp@8: %% Take inputs args from input structure
tp@8: dist = inputstruct.dist;
tp@8: azim = inputstruct.azim;
tp@8: orient = inputstruct.orient;
tp@8: dirweight = inputstruct.dirweight;
Chris@12: outdir = inputstruct.outdir;
Chris@12: outname = inputstruct.outname;
Chris@15: indir = inputstruct.indir;
tp@8: 
tp@8: %% Validate attributes
tp@8: validateattributes(dist,{'double'},{'scalar','>',0})
tp@8: validateattributes(azim,{'double'},{'scalar','<=',90,'>=',-90})
tp@8: validateattributes(orient,{'char'},{'nonempty'})
tp@8: validateattributes(dirweight,{'double'},{'scalar','>',0})
Chris@12: validateattributes(outdir,{'char'},{'nonempty'})
Chris@12: validateattributes(outname,{'char'},{'nonempty'})
Chris@15: validateattributes(indir,{'char'},{'nonempty'})
tp@8: 
tp@8: %% Computation parameters (Some fixed, some taken from input arguments)
tp@8: params.board_size_x = .55;
tp@8: params.board_size_y = .02;
tp@8: params.board_size_z = .55;
tp@8: 
tp@8: params.board_azim = mod(90+azim,360)*pi/180;
tp@8: params.board_elev = 0*pi/180;
tp@8: params.board_dist = dist;
tp@8: %
tp@8: % params.board_azim = mod(90+[-17 17],360)*pi/180; % corresponds to a board 17degress to the right and a board 17
tp@8: % degrees to the left. The elements of the vector ([-17 17] in this example) must be in the range (-180,180] and are
tp@8: % converted by the line in this example to the coordinate system described in the help preample.
tp@8: %
tp@8: % params.board_elev = 10*pi/180; % corresponds to a board 10 degrees above the azimuthal plane. The value (0 in this
tp@8: % example) must be in the range [-90,90] and is converted by the line in this example to the coordinate system described
tp@8: % in the help preample.
tp@8: % 
tp@8: % params.board_dist is in meters and it can be a vector of (strictly positive) distances as needed
tp@8: 
tp@8: params.board_orientation = {validatestring(orient,{'horz','angled'})};
tp@8: 
tp@8: params.source_down = 0.25;   % Take the source to always be directly below the chin.
tp@8: params.source_front = 0.05;  % Take the source to always be directly in front of chest.
tp@8:                   
tp@8: params.Fs = 44100;
tp@8: params.Cair = 344;
tp@8: params.Rhoair = 1.21;
tp@8: params.WavDurationSec = 2;
tp@8: params.WavScaling = 0.1;
tp@8: params.dirweight = dirweight;
tp@8: 
tp@8: params.horz_disamb = 'BSPC2011';
Chris@12: params.wavfilename = [outdir filesep outname '.wav']
tp@8: 
tp@8: %% Compute free field (no head) IRs
tp@8: 
tp@8: temp_edges = [
tp@8:     +params.board_size_x/2 +params.board_size_y/2 +params.board_size_z/2
tp@8:     -params.board_size_x/2 +params.board_size_y/2 +params.board_size_z/2
tp@8:     -params.board_size_x/2 -params.board_size_y/2 +params.board_size_z/2
tp@8:     +params.board_size_x/2 -params.board_size_y/2 +params.board_size_z/2
tp@8:     +params.board_size_x/2 +params.board_size_y/2 -params.board_size_z/2
tp@8:     -params.board_size_x/2 +params.board_size_y/2 -params.board_size_z/2
tp@8:     -params.board_size_x/2 -params.board_size_y/2 -params.board_size_z/2
tp@8:     +params.board_size_x/2 -params.board_size_y/2 -params.board_size_z/2
tp@8:      ];
tp@8: 
tp@8: FFresp{length(params.board_dist),length(params.board_azim),length(params.board_orientation)} = [];
tp@8: binaural_irs_echo{length(params.board_dist),length(params.board_azim),length(params.board_orientation)} = [];
tp@8: binaural_irs_emission{length(params.board_dist),length(params.board_azim),length(params.board_orientation)} = [];
tp@8: % NOTE: in the above initialisation the source term (binaural_irs_emission) does not need to be computed for each
tp@8: % different board distance/azimuth/orientation as it is the same. I include this redundancy so that the emission
tp@8: % binaural irs that I compute each time in in a few lines are correct, i.e. I can check that all the binaural responses
tp@8: % in the binaural_irs_emission cell are equal.
tp@8: board_coords(8,3,length(params.board_dist),length(params.board_azim),length(params.board_orientation)) = 0; 
tp@8: 
tp@8: 
tp@8: geom.Fs = params.Fs;
tp@8: geom.Cair = params.Cair;
tp@8: geom.Rhoair = params.Rhoair;
tp@8: 
tp@8: geom.source_position = [0 params.source_front -params.source_down];
tp@8: 
tp@8: for ii = 1:length(params.board_dist)
tp@8:     for jj = 1:length(params.board_azim)
tp@8:         for kk = 1:length(params.board_orientation)
tp@8:             
tp@8:             if strcmp(params.board_orientation{kk},'angled')
tp@8:                 theta_angled = params.board_azim(jj) - pi/2;
tp@8:                 Rot_mat = [
tp@8:                     cos(theta_angled) -sin(theta_angled) 0
tp@8:                     sin(theta_angled) cos(theta_angled)  0
tp@8:                     0                 0                  1
tp@8:                     ];
tp@8:                 temp2_edges = ( Rot_mat * temp_edges.' ) .' ;
tp@8:                 % for the rotation matrix see http://en.wikipedia.org/wiki/Rotation_matrix
tp@8:                 [trans_x, trans_y, trans_z] = sph2cart(params.board_azim(jj),params.board_elev,params.board_dist(ii));
tp@8:                 geom.scatterer_edges = ...
tp@8:                     temp2_edges + repmat([trans_x, trans_y, trans_z],8,1);
tp@8:                 board_coords(:,:,ii,jj,kk) = geom.scatterer_edges;
tp@8:             elseif strcmp(params.board_orientation{kk},'horz')
tp@8:                 temp2_edges = temp_edges;
tp@8:                 if strcmp(params.horz_disamb,'collocated')
tp@8:                 [trans_x, trans_y, trans_z] = sph2cart(params.board_azim(jj),params.board_elev,params.board_dist(ii));
tp@8:                 elseif strcmp(params.horz_disamb,'BSPC2011')
tp@8:                 [trans_x, trans_y, trans_z] = sph2cart(...
tp@8:                     params.board_azim(jj),...
tp@8:                     params.board_elev,...
tp@8:                     params.board_dist(ii)/abs(cos(params.board_azim(jj)-pi/2)));
tp@8:                 else
tp@8:                     error('params.horz_disamb must be set to either ''BSPC2011'' or ''collocated''')
tp@8:                 end
tp@8:                 geom.scatterer_edges = ...
tp@8:                     temp2_edges + repmat([trans_x, trans_y, trans_z],8,1);
tp@8:                 board_coords(:,:,ii,jj,kk) = geom.scatterer_edges;
tp@8:             end
tp@8:             
Chris@14:             FFresp(ii,jj,kk) = edhrir_single_cuboid(outdir, geom);
tp@8:             
tp@8:         end
tp@8:     end 
tp@8: end
tp@8: 
tp@8: %% Compute binaural IRs
tp@8: 
Chris@15: temp = load([indir filesep 'hrir_final.mat']); %This is subject_165 CIPIC data
tp@8: h3D_lear = temp.hrir_l;
tp@8: h3D_rear = temp.hrir_r;
tp@8: clear temp
tp@8: 
tp@8: cipic_source_az = 0*180/pi;
tp@8: cipic_source_el = atan(-params.source_down/params.source_front)*180/pi;
tp@8: [source_lear_ir, azerr(1,1+length(params.board_azim)), elerr(1,1+length(params.board_azim))] = ...
tp@8:     getNearestUCDpulse(cipic_source_az,cipic_source_el,h3D_lear);
tp@8: [source_rear_ir, azerr(2,1+length(params.board_azim)), elerr(2,1+length(params.board_azim))] = ...
tp@8:     getNearestUCDpulse(cipic_source_az,cipic_source_el,h3D_rear);
tp@8: 
tp@8: % NOTE: in the following computation the source term (binaural_irs_emission) does not need to be computed for each
tp@8: % different board distance/azimuth/orientation as it is the same. I include the redundancy as a check that the geometry
tp@8: % modelling and computations are correct, i.e. I can check that all the binaural responses in the binaural_irs_emission
tp@8: % cell are equal.
tp@8: 
tp@8: for ii = 1:length(params.board_azim)
tp@8:     
tp@8:     cipic_board_az = 90-params.board_azim(ii)*180/pi;
tp@8:     cipic_board_el = params.board_elev*180/pi;
tp@8:     [lear_ir, azerr(1,ii), elerr(1,ii)] = getNearestUCDpulse(cipic_board_az,cipic_board_el,h3D_lear);
tp@8:     [rear_ir, azerr(2,ii), elerr(2,ii)] = getNearestUCDpulse(cipic_board_az,cipic_board_el,h3D_rear);
tp@8:     
tp@8:     for jj = 1:length(params.board_dist)
tp@8:         for kk = 1:length(params.board_orientation)
tp@8:             
tp@8:             binaural_irs_emission{jj,ii,kk}(1,:) = fftconv(source_lear_ir,FFresp{jj,ii,kk}(2,:)); 	%source term
tp@8:             binaural_irs_emission{jj,ii,kk}(2,:) = fftconv(source_rear_ir,FFresp{jj,ii,kk}(2,:)); 	%source term
tp@8:             
tp@8:             binaural_irs_echo{jj,ii,kk}(1,:) = fftconv(lear_ir,FFresp{jj,ii,kk}(1,:)+FFresp{jj,ii,kk}(3,:)); 
tp@8:                                                                                                     %board echo term
tp@8:             binaural_irs_echo{jj,ii,kk}(2,:) = fftconv(rear_ir,FFresp{jj,ii,kk}(1,:)+FFresp{jj,ii,kk}(3,:));
tp@8:                                                                                                     %board echo term
tp@8: 
tp@8:         end
tp@8:     end
tp@8: end
tp@8: 
tp@8: simulation_data.echo = binaural_irs_echo;
tp@8: simulation_data.emission =binaural_irs_emission;
tp@8: simulation_data.params = params;
tp@8: simulation_data.board_coords = board_coords;
tp@8: simulation_data.azerr = azerr;
tp@8: simulation_data.elerr = elerr;
tp@8: 
tp@8: %% Compute binaural signals and save binaural wav file
tp@8: 
tp@8: binaural_ir = params.dirweight*simulation_data.emission{1,1,1} + simulation_data.echo{1,1,1};
tp@8: stim = params.WavScaling*randn(1,params.WavDurationSec*params.Fs);
tp@8: temp = [fftfilt(binaural_ir(1,:),stim);fftfilt(binaural_ir(2,:),stim)].';
tp@8: if max(max(abs(temp))) > 1
tp@8:     warning('wavfile clipped')
tp@8: end
tp@8: audiowrite(params.wavfilename,temp,params.Fs);
tp@8: 
tp@8: %% Subfunctions
tp@8: 
tp@8: function [pulse, azerr, elerr] = getNearestUCDpulse(azimuth,elevation,h3D)
tp@8: azimuth = pvaldeg(azimuth);
tp@8: elevation = pvaldeg(elevation);
tp@8: 
tp@8: elmax = 50;
tp@8: elindices = 1:elmax;
tp@8: elevations = -45 + 5.625*(elindices - 1);
tp@8: el = round((elevation + 45)/5.625 + 1);
tp@8: el = max(el,1);
tp@8: el = min(el,elmax);
tp@8: elerr = pvaldeg(elevation - elevations(el));
tp@8: 
tp@8: azimuths = [-80 -65 -55 -45:5:45 55 65 80];
tp@8: [azerr, az] = min(abs(pvaldeg(abs(azimuths - azimuth))));
tp@8: 
tp@8: pulse = squeeze(h3D(az,el,:));
tp@8: 
tp@8: function angle = pvaldeg(angle)
tp@8: dtr = pi/180;
tp@8: angle = atan2(sin(angle*dtr),cos(angle*dtr))/dtr;
tp@8: if angle < - 90
tp@8:     angle = angle + 360;
tp@8: end
tp@8: 
tp@8: function out=fftconv(in1,in2)
tp@8: % 
tp@8: % Computes the convolution output out of the input vectors in1 and in2
tp@8: % either with time-domain convolution (as implemented by conv) or by
tp@8: % frequency-domain filtering (as implemented by zero-padded fftfilt).
tp@8: %
tp@8: % out=fftconv(in1,in2)
tp@8: % 
tp@8: in1=in1(:).';in2=in2(:).';
tp@8: 
tp@8: if (length(in1)>100 && length(in2)>100)
tp@8:     if length(in1)>=length(in2)
tp@8:         out=fftfilt(in2,[in1 zeros(1,length(in2)-1)]);
tp@8:     else
tp@8:         out=fftfilt(in1,[in2 zeros(1,length(in1)-1)]);
tp@8:     end
tp@8: else
tp@8:     out=conv(in1,in2);
tp@8: end
tp@8: 
Chris@14: function ir_matrix=edhrir_single_cuboid(outdir, geom)
tp@8: %
tp@8: % Matrix of irs for a single cuboid scatterer using EDTB
tp@8: %
tp@8: % Works for external geometry. Definition of corners and planes in the created CAD file should be pointing outwards (in
tp@8: % cartesian coordinates as defined below).
tp@8: %
tp@8: % x coordinate  (increasing to right in top view, right hand thumb)
tp@8: % y coordinate  (increasing in front in top view, right hand index)
tp@8: % z coordinate  (increasing upwards in top view, right hand middle)
tp@8: % 
tp@8: % For more details about geometry description see the main help preample.
tp@8: %
tp@8: % geom input parameter should be a structure with the following fields:
tp@8: %
tp@8: % geom.scatterer_edges is a 8x3 matrix with elements:
tp@8: % x y z coordinates of back  right top    edge (point 1)
tp@8: % x y z coordinates of back  left  top    edge (point 2)
tp@8: % x y z coordinates of front left  top    edge (point 3)
tp@8: % x y z coordinates of front right top    edge (point 4)
tp@8: % x y z coordinates of back  right bottom edge (point 1)
tp@8: % x y z coordinates of back  left  bottom edge (point 2)
tp@8: % x y z coordinates of front left  bottom edge (point 3)
tp@8: % x y z coordinates of front right bottom edge (point 4)
tp@8: %
tp@8: % (or any sequence of rotations of a cuboid with edges as above)
tp@8: % 
tp@8: % geom.source_position  (Single) source coordinates in cartesian system defined 1x3 positive real vector. (Single)
tp@8: % receiver always taken to be at coordinates origin (0,0,0)
tp@8: % 
tp@8: % geom.Fs
tp@8: % geom.Cair
tp@8: % geom.Rhoair
tp@8: %
tp@8: 
tp@8: %% 1
Chris@13: temp_filename=['a' num2str(now*1e12,'%-24.0f')];
Chris@12: cad_filename=[temp_filename '.cad']
tp@8: 
Chris@13: fid=fopen([outdir filesep cad_filename],'w');
tp@8: % TODO add onCleanup to fclose this file
tp@8: fprintf(fid,'%%CORNERS\n\n');
tp@8: fprintf(fid,'%.0f %9.6f %9.6f %9.6f\n',...
tp@8:     1, geom.scatterer_edges(1,1), geom.scatterer_edges(1,2), geom.scatterer_edges(1,3),...
tp@8:     2, geom.scatterer_edges(2,1), geom.scatterer_edges(2,2), geom.scatterer_edges(2,3),...
tp@8:     3, geom.scatterer_edges(3,1), geom.scatterer_edges(3,2), geom.scatterer_edges(3,3),...
tp@8:     4, geom.scatterer_edges(4,1), geom.scatterer_edges(4,2), geom.scatterer_edges(4,3),...
tp@8:     5, geom.scatterer_edges(5,1), geom.scatterer_edges(5,2), geom.scatterer_edges(5,3),...
tp@8:     6, geom.scatterer_edges(6,1), geom.scatterer_edges(6,2), geom.scatterer_edges(6,3),...
tp@8:     7, geom.scatterer_edges(7,1), geom.scatterer_edges(7,2), geom.scatterer_edges(7,3),...
tp@8:     8, geom.scatterer_edges(8,1), geom.scatterer_edges(8,2), geom.scatterer_edges(8,3));
tp@8: fprintf(fid,'\n%%PLANES\n');
tp@8: fprintf(fid,'\n1 / /RIGID\n1 2 3 4\n');
tp@8: fprintf(fid,'\n2 / /RIGID\n5 8 7 6\n');
tp@8: fprintf(fid,'\n3 / /RIGID\n1 4 8 5\n');
tp@8: fprintf(fid,'\n4 / /RIGID\n1 5 6 2\n');
tp@8: fprintf(fid,'\n5 / /RIGID\n2 6 7 3\n');
tp@8: fprintf(fid,'\n6 / /RIGID\n7 8 4 3\n');
tp@8: fprintf(fid,'\n%%EOF');
tp@8: fclose(fid);
tp@8: 
Chris@12: setup_filename = [temp_filename '_setup.m']
luis@9: 
Chris@13: fid=fopen([outdir filesep setup_filename],'wt');
tp@8: % TODO add onCleanup to fclose this file
tp@8: fprintf(fid,'\n global FSAMP CAIR RHOAIR SHOWTEXT');
tp@8: fprintf(fid,['\n FSAMP = ' num2str(geom.Fs) ';']);
tp@8: fprintf(fid,['\n CAIR = ' num2str(geom.Cair) ';']);
tp@8: fprintf(fid,['\n RHOAIR = ' num2str(geom.Rhoair) ';']);
Chris@14: fprintf(fid,'\n SHOWTEXT = 2;');
tp@8: fprintf(fid,'\n SUPPRESSFILES = 0;');
Chris@18: temp = strrep([outdir filesep],'\','\\');
Chris@18: fprintf(fid,['\n Filepath=''' temp ''';']);
tp@8: fprintf(fid,['\n Filestem=''' temp_filename ''';']);
Chris@18: fprintf(fid,['\n CADfile=''' temp cad_filename ''';']);
tp@8: fprintf(fid,'\n open_or_closed_model = ''open'';');
tp@8: fprintf(fid,'\n int_or_ext_model = ''ext'';');
tp@8: fprintf(fid,'\n EDcalcmethod = ''n'';');
tp@8: fprintf(fid,'\n directsound = 1;');
tp@8: fprintf(fid,'\n specorder = 1;');
tp@8: fprintf(fid,'\n difforder = 1;');
tp@8: fprintf(fid,'\n elemsize = [1];');
tp@8: fprintf(fid,'\n nedgesubs = 2;');
tp@8: fprintf(fid,'\n calcpaths = 1;');
tp@8: fprintf(fid,'\n calcirs = 1;');
tp@8: fprintf(fid,'\n sources=[');
tp@8: for n=1:1
tp@8:     fprintf(fid,['\n' num2str(geom.source_position(1)) ' ' ...
tp@8:         num2str(geom.source_position(2)) ' ' ...
tp@8:         num2str(geom.source_position(3))]);
tp@8: end
tp@8: fprintf(fid,'\n ];');
tp@8: fprintf(fid,'\n receivers=[');
tp@8: for n=1:1
tp@8:     fprintf(fid,'\n 0 0 0');
tp@8: end
tp@8: fprintf(fid,'\n ];');
tp@8: fprintf(fid,'\n skipcorners = 1000000;');
tp@8: fprintf(fid,'\n Rstart = 0;');
tp@8: 
tp@8: fclose(fid);
tp@8: 
Chris@16: %pause(1);
tp@8: 
Chris@14: [~, ir_matrix]=myver_edtb(outdir, [outdir filesep setup_filename]);
tp@8: 
Chris@13: delete([outdir filesep cad_filename]);
Chris@13: delete([outdir filesep setup_filename]);
tp@8: 
Chris@14: function [ir,varargout]=myver_edtb(outdir,EDsetupfile)
tp@8: % implements the computation of EDToolbox (by Peter Svensson) with a front-end
tp@8: % designed for my needs.
tp@8: % 
tp@8: % Takes one input which a setup file as specified by Svensson and gives
tp@8: % either one vector output (which is the irtot output of Svensson's implementation for
tp@8: % the first source and first receiver in the input "EDsetupfile" setup file) or
tp@8: % two outputs, the first as above and the second being a {NxM} cell (N the number
tp@8: % of sources and M the number of receivers) each element of which is a 4xL matrix
tp@8: % with its 4 rows being the irdiff, irdirect, irgeom and irtot outputs of Svensson's
tp@8: % computation for the corresponding source and receiver.
tp@8: %
tp@8: % The code deletes all other .mat files created by Svensson's
tp@8: % implementation.
tp@8: % 
tp@8: % [ir ir_all_data]=myver_edtb(EDsetupfile);
tp@8: 
tp@8: %keep record of files in the current directory before computation
Chris@14: dir_bef=dir(outdir);
tp@8: ir=1;
tp@8: 
Chris@14: disp(['Will call ' EDsetupfile])
luis@9: 
tp@8: %run computation
tp@8: EDB1main(EDsetupfile);
tp@8: 
tp@8: %get files list in the current directory after computation
Chris@14: dir_aft=dir(outdir);
tp@8: 
tp@8: %get all genuine files (excluding other directories) in current directory
tp@8: %after computation
tp@8: k=1;
tp@8: for n=1:size(dir_aft,1)
tp@8:     if(~dir_aft(n).isdir)
tp@8:         names_aft{k}=dir_aft(n).name;k=k+1; %#ok<AGROW>
tp@8:     end
tp@8: end
tp@8: 
tp@8: %get all genuine files (excluding other directories) in current directory
tp@8: %before computation
tp@8: k=1;
tp@8: for n=1:size(dir_bef,1)
tp@8:     if(~dir_bef(n).isdir)
tp@8:         names_bef{k}=dir_bef(n).name;k=k+1; %#ok<AGROW>
tp@8:     end
tp@8: end
tp@8: 
tp@8: %get irtot for 1st source 1st receiver and delete files created by
tp@8: %computation.
tp@8: % for n=1:length(names_aft)
tp@8: %     if ~strcmp(names_aft(n),names_bef)
tp@8: %         if ~isempty(strfind(names_aft{n},'_1_1_ir.mat'))
tp@8: %             temp=load(names_aft{n});
tp@8: %             ir=full(temp.irtot);
tp@8: %             clear temp
tp@8: %         end
tp@8: %         delete(names_aft{n})
tp@8: %     end
tp@8: % end
tp@8: for n=1:length(names_aft)
tp@8:     if ~strcmp(names_aft(n),names_bef)
tp@8:         if ~isempty(strfind(names_aft{n},'_1_1_ir.mat'))
Chris@14:             temp=load([outdir filesep names_aft{n}]);
tp@8:             ir=full(temp.irtot);
tp@8:             clear temp
tp@8:         end
tp@8:         if nargout>1
tp@8:             if ~isempty(strfind(names_aft{n},'_ir.mat'))
Chris@14:                 temp=load([outdir filesep names_aft{n}]);
tp@8:                 temp2=regexp(regexp(names_aft{n},'_\d*_\d*_','match'),'\d*','match');
tp@8:                 temp3{str2double(temp2{1}(1)),str2double(temp2{1}(2))}(4,:)=full(temp.irtot); %#ok<AGROW>
tp@8:                 temp3{str2double(temp2{1}(1)),str2double(temp2{1}(2))}(1,:)=full(temp.irdiff); %#ok<AGROW>
tp@8:                 temp3{str2double(temp2{1}(1)),str2double(temp2{1}(2))}(2,:)=full(temp.irdirect); %#ok<AGROW>
tp@8:                 temp3{str2double(temp2{1}(1)),str2double(temp2{1}(2))}(3,:)=full(temp.irgeom); %#ok<AGROW>
tp@8:                 clear temp temp2
tp@8:             end
tp@8:         end
Chris@14:         delete([outdir filesep names_aft{n}])
tp@8:     end
tp@8: end
tp@8: if nargout>1
tp@8:     varargout(1)={temp3};
tp@8: end
tp@8: