Speedup sofa

NEWS

@@ -15,6 +15,7 @@ Changes in the Sound Field Synthesis-Toolbox. Recent changes on top.
       delayline()
     - add heuristic to find a good local wave field synthesis pre-filter
     - loudspeaker geometry can now be read from a SOFA file
+    - speedup ir_generic() by introducing get_ir_loop()
 
 1.2.0 (2. June 2015)
 

SFS_binaural_synthesis/ir_generic.m

@@ -75,26 +75,43 @@
 
 
 %% ===== BRIR ===========================================================
+warning('off','SFS:irs_intpol');
+warning('off','SOFA:upgrade');
+
 % Initial values
 ir_generic = zeros(N,2);
 
+% Get information about given SOFA data set.
+% Doing this outside the loop together with using get_ir_loop() speeds up
+% things by a large factor. If you use get_ir() instead you don't have to do
+% this.
+sofa_header = sofa_get_header(sofa);
+header.convention = sofa_header.GLOBAL_SOFAConventions;
+header.X = sofa_get_listener_position(sofa_header,'cartesian');
+header.x0 = sofa_get_secondary_sources(sofa_header,'spherical');
+[header.phi,header.theta] = sofa_get_head_orientations(sofa_header);
+header.API.N = sofa_header.API.N;
+header.API.R = sofa_header.API.N;
+header.API.E = sofa_header.API.N;
+
 % Create a BRIR for every single loudspeaker
-warning('off','SFS:irs_intpol');
 for ii=1:size(x0,1)
 
     % === Get the desired impulse response.
     % If needed interpolate the given impulse response set and weight, delay the
     % impulse for the correct distance
-    ir = get_ir(sofa,X,[phi 0],x0(ii,1:3),'cartesian',conf);
+    ir = get_ir_loop(sofa,header,X,[phi 0],x0(ii,1:3),'cartesian',conf);
 
     % === Sum up virtual loudspeakers/HRIRs and add loudspeaker time delay ===
     % Also applying the weights of the secondary sources including integration
     % weights or tapering windows etc.
     ir_generic = ir_generic + fix_length(convolution(ir,d(:,ii)),N).*x0(ii,7);
 
 end
-warning('on','SFS:irs_intpol');
 
 
 %% ===== Headphone compensation =========================================
 ir = compensate_headphone(ir_generic,conf);
+
+warning('on','SFS:irs_intpol');
+warning('on','SOFA:upgrade');

SFS_ir/get_ir_loop.m

@@ -0,0 +1,191 @@
+function [ir,x0] = get_ir_loop(sofa,header,X,head_orientation,xs,coordinate_system,conf)
+%GET_IR_LOOP returns an impulse response for the given apparent angle
+% This is an optimized version of get_ir() and should be called inside a loop.
+%
+%   Usage: ir = get_ir_loop(sofa,header,X,head_orientation,xs,coordinate_system,conf)
+%
+%   Input parameters:
+%       sofa              - impulse response data set (sofa struct/file)
+%       header            - struct, containing the following fields:
+%                           convention, X, x0, phi, theta, API.N, API.R, API.E.
+%                           See ir_generic() for more details.
+%       X                 - position of the listener, specified in cartesian
+%                           coordinates
+%       head_orientation  - orientation of the listener with [phi theta] /
+%                           (rad, rad)
+%       xs                - position of the desired source, specified in
+%                           cartesian coordinates. For SOFA convention
+%                           SimpleFreeFieldHRIR xs will be interpreted relative
+%                           to X, for MultiSpeakerBRIR as an absolute position.
+%       conf              - configuration struct (see SFS_config),
+%                           which will be passed to:
+%                             interpolate_ir
+%                             ir_correct_distance (only for SimpleFreeFieldHRIR)
+%
+%   Output parameters:
+%       ir      - impulse response for the given position (length of IR x 2)
+%       x0      - position corresponding to the returned impulse response
+%
+%   GET_IR_LOOP(sofa,header,X,head_orientation,xs,conf) returns a single impulse
+%   response from the given SOFA file or struct. The impulse response is
+%   determined by the position X and head orientation head_orientation of the
+%   listener, and the position xs of the desired point source.
+%   For the SOFA convention SimpleFreeFieldHRIR the desired distance between
+%   the point source and listener is achieved by delaying and weighting the
+%   impulse response. Distances larger than 10m are ignored and set constantly
+%   to 10m.
+%   If the desired angles are not present in the SOFA data set and
+%   conf.ir.useinterpolation is set to true an interpolation is applied to
+%   create the impulse response. For the SOFA convention MultiSpeakerBRIR the
+%   interpolation is only performed for the head orientations not the different
+%   loudspeaker positions.
+%   A further configuration setting that is considered is conf.ir.useoriglength,
+%   which indicates if additional zeros corresponding to the actual radius of
+%   the measured impulse responses should be added at the beginning of every
+%   impulse response (if set to false). If you know that the measured impulse
+%   responses include already the zeros from the measurement it can be safely
+%   set to true. This is important because the delaying of the impulse responses
+%   in order to achieve the correct distance require enough zeros at the
+%   beginning of every impulse response.
+%
+%   For a description of the SOFA file format see: http://sofaconventions.org
+%
+%   See also: get_ir, ir_generic, SOFAload, sofa_get_data_fir,
+%             sofa_get_data_fire, intpolate_ir, ir_correct_distance
+
+%*****************************************************************************
+% Copyright (c) 2010-2015 Quality & Usability Lab, together with             *
+%                         Assessment of IP-based Applications                *
+%                         Telekom Innovation Laboratories, TU Berlin         *
+%                         Ernst-Reuter-Platz 7, 10587 Berlin, Germany        *
+%                                                                            *
+% Copyright (c) 2013-2015 Institut fuer Nachrichtentechnik                   *
+%                         Universitaet Rostock                               *
+%                         Richard-Wagner-Strasse 31, 18119 Rostock           *
+%                                                                            *
+% This file is part of the Sound Field Synthesis-Toolbox (SFS).              *
+%                                                                            *
+% The SFS is free software:  you can redistribute it and/or modify it  under *
+% the terms of the  GNU  General  Public  License  as published by the  Free *
+% Software Foundation, either version 3 of the License,  or (at your option) *
+% any later version.                                                         *
+%                                                                            *
+% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY *
+% WARRANTY;  without even the implied warranty of MERCHANTABILITY or FITNESS *
+% FOR A PARTICULAR PURPOSE.                                                  *
+% See the GNU General Public License for more details.                       *
+%                                                                            *
+% You should  have received a copy  of the GNU General Public License  along *
+% with this program.  If not, see <http://www.gnu.org/licenses/>.            *
+%                                                                            *
+% The SFS is a toolbox for Matlab/Octave to  simulate and  investigate sound *
+% field  synthesis  methods  like  wave  field  synthesis  or  higher  order *
+% ambisonics.                                                                *
+%                                                                            *
+% http://github.com/sfstoolbox/sfs                      sfstoolbox@gmail.com *
+%*****************************************************************************
+
+
+%% ===== Computation ====================================================
+%
+% Get the listener position during the measurement
+X_sofa = header.X;
+x0 = header.x0;
+
+% === Get Impulse Response ===
+if strcmp('SimpleFreeFieldHRIR',header.convention)
+    %
+    % http://www.sofaconventions.org/mediawiki/index.php/SimpleFreeFieldHRIR
+    %
+    % Returns a single impulse response for the desired position. The impulse
+    % response is shifted in time and its amplitude is weighted according to the
+    % desired distance. The desired direction is done by returning the nearest
+    % neighbour or applying a linear interpolation.
+    % NOTE: for SimpleFreeFieldHRIR head orientation is always zero in the SOFA
+    % file and we handle a change in head orientation by changing the source
+    % position accordingly.
+    %
+    % For SimpleFreeFieldHRIR only the relative position between listener
+    % position and source position is of relevance.
+    xs = xs-X+X_sofa;
+    % Get measured loudspeaker positions and go to spherical coordinates
+    [xs(1),xs(2),xs(3)] = cart2sph(xs(1),xs(2),xs(3));
+    % Combine head orientation and desired direction of source (see note above)
+    xs(1) = correct_azimuth(xs(1)-head_orientation(1));
+    xs(2) = correct_elevation(xs(2)-head_orientation(2));
+    % Find nearest neighbours and interpolate if desired and needed
+    [neighbours,idx] = findnearestneighbour(x0(:,1:2)',xs(1:2),3);
+    ir = sofa_get_data_fir(sofa,idx,header);
+    ir = ir_correct_distance(ir,x0(idx,3),xs(3),conf);
+    [ir,x0] = interpolate_ir(ir,neighbours,xs(1:2)',conf);
+    [x0(1),x0(2),x0(3)] = sph2cart(x0(1),x0(2),xs(3));
+
+elseif strcmp('MultiSpeakerBRIR',header.convention)
+    %
+    % http://www.sofaconventions.org/mediawiki/index.php/MultiSpeakerBRIR
+    %
+    % This looks for the nearest loudspeaker corresponding to the desired source
+    % position. For that loudspeaker the impulse reponse for the specified head
+    % orientation is returned. If the head orientation could not perfectly
+    % matched, an interpolation is applied. If the specified head orientation is
+    % out of bounds, the nearest head orientation is returned.
+    %
+    % Check if we are requesting a listening position that is available
+    if norm(X-X_sofa)>0.01
+        warning('SFS:get_ir',['Your chosen listening position (%.2f,', ...
+                              '%.2f,%.2f)m is not available, using the ', ...
+                              'measured one (%.2f,%.2f,%.2f)m instead.'], ...
+                             X(1),X(2),X(3),X_sofa(1),X_sofa(2),X_sofa(3));
+    end
+    % Find nearest loudspeaker
+    [neighbours_emitter,idx_emitter] = findnearestneighbour(x0(:,1:3)',xs,1);
+    x0 = neighbours_emitter(:,1);
+    if norm(x0'-xs)>0.01
+        warning('SFS:get_ir',['Your chosen loudspeaker position (%.2f,', ...
+                              '%.2f,%.2f)m deviates from the measured ', ...
+                              'one (%.2f,%.2f,%.2f)m.'], ...
+                             xs(1),xs(2),xs(3),x0(1),x0(2),x0(3));
+    end
+    % Find nearest head orientation in the horizontal plane
+    [neighbours_head,idx_head] = ...
+        findnearestneighbour([header.phi header.theta]',head_orientation,3);
+    % Check if head orientation is out of bounds
+    if all(abs(head_orientation(1))>abs(neighbours_head(1,:)))
+        warning('SFS:get_ir',['Head azimuth %.1f deg out of bound, ', ...
+            'using %.1f deg instead.'], ...
+            deg(head_orientation(1)), ...
+            deg(neighbours_head(1,1)));
+        head_orientation(1) = neighbours_head(1,1);
+    end
+    if all(abs(head_orientation(2))>abs(neighbours_head(2,:)))
+        warning('SFS:get_ir',['Head elevation %.1f deg out of bound, ', ...
+            'using %.1f deg instead.'], ...
+            deg(head_orientation(2)), ...
+            deg(neighbours_head(2,1)));
+        head_orientation(2) = neighbours_head(2,1);
+    end
+    % Get the impulse responses, reshape and interpolate
+    ir = sofa_get_data_fire(sofa,idx_head,idx_emitter);
+    ir = reshape(ir,[size(ir,1) size(ir,2) size(ir,4)]); % [M R E N] => [M R N]
+    ir = interpolate_ir(ir,neighbours_head,head_orientation',conf);
+
+elseif strcmp('SingleRoomDRIR',header.convention)
+    %
+    % http://www.sofaconventions.org/mediawiki/index.php/SingleRoomDRIR
+    %
+    error(['%s: SingleRoomDRIR is not supported as it should handle ', ...
+           'microphone array recordings. If you used it for (multiple) ', ...
+           'loudspeakers in a room you should consider to use ', ...
+           'MultiSpeakerBRIR instead.'], upper(mfilename));
+
+else
+    error('%s: %s convention is currently not supported.', ...
+        upper(mfilename),header.convention);
+end
+
+% Reshape [1 2 N] to [N 2]
+ir = squeeze(ir)';
+% Convert x0 to the specified coordinate system
+if strcmp('spherical',coordinate_system)
+    [x0(1),x0(2),x0(3)] = cart2sph(x0(1),x0(2),x0(3));
+end

SFS_ir/sofa_get_data_fir.m

@@ -1,7 +1,7 @@
-function ir = sofa_get_data_fir(sofa,idx)
+function ir = sofa_get_data_fir(sofa,idx,header)
 %SOFA_GET_DATA_FIR returns impulse responses from a SOFA file or struct
 %
-%   Usage: ir = sofa_get_data_fir(sofa,[idx])
+%   Usage: ir = sofa_get_data_fir(sofa,[idx,[header]])
 %
 %   Input parameters:
 %       sofa    - impulse response data set (SOFA struct/file)
@@ -11,13 +11,16 @@
 %                 responses for the corresponding index positions will be
 %                 returned.
 %                 If no index is specified all data will be returned.
+%       header  - header of the sofa file. This will speed things up, as the
+%                 header has not to be extracted from the sofa file. This is
+%                 especially useful if you call this function inside a loop
 %
 %   Output parameters:
 %       ir      - impulse response (M,2,N), where
 %                   M ... number of impulse responses
 %                   N ... samples
 %
-%   SOFA_GET_DATA_FIR(sofa,idx) returns impulse response of the given
+%   SOFA_GET_DATA_FIR(sofa,idx,header) returns impulse response of the given
 %   SOFA file or struct, specified by idx. If no idx is specified all data
 %   contained in sofa is returned.
 %   For the struct the SOFA file has to loaded before with SOFAload().
@@ -60,12 +63,13 @@
 
 %% ===== Checking of input  parameters ==================================
 nargmin = 1;
-nargmax = 2;
+nargmax = 3;
 narginchk(nargmin,nargmax)
 if nargin==nargmax-1
+    header = [];
+elseif nargin==nargmax-2
+    header = [];
     idx = [];
-else
-    isargvector(idx);
 end
 
 
@@ -76,8 +80,10 @@
     end
     ir = sofa.Data.IR;
 else
-    header = sofa_get_header(sofa);
     if sofa_is_file(sofa)
+        if isempty(header)
+            header = sofa_get_header(sofa);
+        end
         ir = zeros(length(idx),2,header.API.N);
         for ii=1:length(idx)
             tmp = SOFAload(sofa,[idx(ii) 1]);

SFS_ir/sofa_get_data_fire.m

@@ -1,7 +1,7 @@
-function ir = sofa_get_data_fire(sofa,idxM,idxE)
+function ir = sofa_get_data_fire(sofa,idxM,idxE,header)
 %SOFA_GET_DATA_FIRE returns impulse responses from a SOFA file or struct
 %
-%   Usage: ir = sofa_get_data_fire(sofa,[idxM],[idxE])
+%   Usage: ir = sofa_get_data_fire(sofa,[idxM,[idxE,[header]]])
 %
 %   Input parameters:
 %       sofa    - impulse response data set (SOFA struct/file)
@@ -14,16 +14,19 @@
 %                 If no index is specified all measurements will be returned.
 %       idxE    - index of the emitter for which the single impulse
 %                 responses should be returned. The rest is identical to idxM.
+%       header  - header of the sofa file. This will speed things up, as the
+%                 header has not to be extracted from the sofa file. This is
+%                 especially useful if you call this function inside a loop
 %
 %   Output parameters:
 %       ir      - impulse response (M,2,E,N), where
 %                   M ... number of impulse responses
 %                   E ... number of emitters (loudspeakers)
 %                   N ... samples
 %
-%   SOFA_GET_DATA_FIRE(sofa,idxM,idxE) returns impulse response of the given
-%   SOFA file or struct, specified by idxM and idxE, where idxM defines the
-%   measurements and idxE the emitters for which impulse responses should be
+%   SOFA_GET_DATA_FIRE(sofa,idxM,idxE,header) returns impulse response of the
+%   given SOFA file or struct, specified by idxM and idxE, where idxM defines
+%   the measurements and idxE the emitters for which impulse responses should be
 %   returned. If no index is specified all data contained in sofa is returned.
 %   For the struct the SOFA file has to loaded before with SOFAload().
 %   For a description of the SOFA file format see: http://sofaconventions.org
@@ -65,19 +68,25 @@
 
 %% ===== Checking of input  parameters ==================================
 nargmin = 1;
-nargmax = 3;
+nargmax = 4;
 narginchk(nargmin,nargmax)
 if nargin==nargmax-1
+    header = [];
+if nargin==nargmax-2
+    header = [];
     idxE = ':';
-elseif nargin==nargmax-2
+elseif nargin==nargmax-3
+    header = [];
     idxE = ':';
     idxM = ':';
 end
 
 
 %% ===== Computation ====================================================
 if sofa_is_file(sofa)
-    header = sofa_get_header(sofa);
+    if isempty(header)
+        header = sofa_get_header(sofa);
+    end
     if isnumeric(idxE) && isnumeric(idxM)
         ir = zeros(length(idxM),header.API.R,length(idxE),header.API.N);
         for ii=1:length(idxM)