diff simulateBinauralSignals.m @ 8:aeb703465220

added simulateBinauralSignals.m

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