diff --git a/NEWS b/NEWS index 31d9ec95..1523d172 100644 --- a/NEWS +++ b/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) diff --git a/SFS_binaural_synthesis/ir_generic.m b/SFS_binaural_synthesis/ir_generic.m index c6e3d815..9515144f 100644 --- a/SFS_binaural_synthesis/ir_generic.m +++ b/SFS_binaural_synthesis/ir_generic.m @@ -75,17 +75,32 @@ %% ===== 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 @@ -93,8 +108,10 @@ 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'); diff --git a/SFS_ir/get_ir_loop.m b/SFS_ir/get_ir_loop.m new file mode 100644 index 00000000..1c8f30bf --- /dev/null +++ b/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 . * +% * +% 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 diff --git a/SFS_ir/sofa_get_data_fir.m b/SFS_ir/sofa_get_data_fir.m index 43baf1f0..c5faee54 100644 --- a/SFS_ir/sofa_get_data_fir.m +++ b/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]); diff --git a/SFS_ir/sofa_get_data_fire.m b/SFS_ir/sofa_get_data_fire.m index f3e63bc2..26c5bc13 100644 --- a/SFS_ir/sofa_get_data_fire.m +++ b/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,6 +14,9 @@ % 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 @@ -21,9 +24,9 @@ % 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,11 +68,15 @@ %% ===== 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 @@ -77,7 +84,9 @@ %% ===== 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) diff --git a/SFS_ir/sofa_is_file.m b/SFS_ir/sofa_is_file.m index 24cb679c..d86f9be7 100644 --- a/SFS_ir/sofa_is_file.m +++ b/SFS_ir/sofa_is_file.m @@ -1,20 +1,19 @@ function boolean = sofa_is_file(sofa) -%SOFA_CHECK returns 1 for a sofa file, 0 for a sofa struct or an error otherwise +%SOFA_IS_FILE returns 1 for a sofa file, 0 for a sofa struct or an error otherwise % -% Usage: number = sofa_check(sofa) +% Usage: boolean = sofa_is_file(sofa) % % Input parameters: % sofa - sofa struct or file name % % Output parameters: -% number - 1: sofa is a file -% 0: sofa is a struct +% boolean - true: sofa is a file +% false: sofa is a struct % -% SOFA_CHECK(sofa) checks if the given sofa is a file or a struct. In the -% first case a 1 is returned, in the second case a 0. If none of both is true -% an error will be thrown. +% SOFA_IS_FILE(sofa) checks if the given sofa variable is a file or a struct. +% If both are false an error will be thrown. % -% See also: sofa_get_header, sofa_get_data +% See also: sofa_get_header, sofa_get_data_fir %***************************************************************************** % Copyright (c) 2010-2015 Quality & Usability Lab, together with * @@ -50,17 +49,18 @@ %% ===== Checking of input parameters ================================== -nargmin = 1; -nargmax = 1; -narginchk(nargmin,nargmax) +% No checking for speedup +%nargmin = 1; +%nargmax = 1; +%narginchk(nargmin,nargmax) %% ===== Main =========================================================== -if ~isstruct(sofa) && exist(sofa,'file') - boolean = true; -elseif isstruct(sofa) && isfield(sofa,'GLOBAL_Conventions') && ... - strcmp('SOFA',sofa.GLOBAL_Conventions) +if isstruct(sofa) && isfield(sofa,'GLOBAL_Conventions') && ... + strcmp('SOFA',sofa.GLOBAL_Conventions) boolean = false; +elseif ~isstruct(sofa) && exist(sofa,'file') + boolean = true; else error('%s: sofa has to be a file or a SOFA struct.',upper(mfilename)); end