Skip to content

GitLab

toolbox
toolbox
Commit 5d1cc16c authored by Michael Kohnen's avatar Michael Kohnen
Browse files
Options

Merge branch 'master' of http://git.rwth-aachen.de/ita/toolbox

parents 4663b01d 8357bc22
Hide whitespace changes
Inline Side-by-side
Showing with 1181 additions and 953 deletions
applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m
View file @ 5d1cc16c
......@@ -43,12 +43,10 @@ function cThis = interp(this,varargin)
% Acta Acustica united with Acustica, Volume 100, Number 4, July/August 2014,
% pp. 667-675(9)
%
%
%
% Author: Florian Pausch <fpa@akustik.rwth-aachen.de>
% Version: 2016-02-05
sArgs = struct('order',50,'eps',1e-5,'shiftToEar',false,'shiftAxis','y','shiftOffset',[-0.0725 0.0725]);
sArgs = struct('order',50,'epsilon',1e-8,'shiftToEar',false,'shiftAxis','y','shiftOffset',[-0.0725 0.0725]);
sArgs = ita_parse_arguments(sArgs,varargin,2);
if isempty(varargin) || ~isa(varargin{1},'itaCoordinates')
error('itaHRTF:interp', ' An itaCoordinate object is needed!')
......@@ -65,18 +63,22 @@ field.phi_deg = tempfield(:,1);
field.theta_deg = tempfield(:,2);
N = sArgs.order;
epsilon = sArgs.eps; % regularization parameter
epsilon = sArgs.epsilon; % regularization parameter
k = this.wavenumber; % wave number
k(1) = eps;
% add eps to avoid NaN's
Nmax = floor(sqrt(this.nDirections/4)-1);
if N>Nmax
N=Nmax;
warning(['Order of SH reconstruction matrix was set to N = Nmax = ',num2str(Nmax),'.'])
end
% construct vector of length (N+1)^2 regularization weights and,
% if needed, spherical hankel functions of second kind (for r0 and r1)
if ~isequal(this.dirCoord.r(1),field.r(1))
kr0 = k*this.dirCoord.r(1); % measurement radius
kr1 = k*field.r(1); % extrapolation radius
kr0 = k*this.dirCoord.r(1); % measurement radius [m]
kr1 = k*field.r(1); % extrapolation radius [m]
hankel_r0 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr0);
hankel_r1 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr1);
......@@ -102,10 +104,9 @@ end
%% Weights
[~,w]= this.dirCoord.spherical_voronoi; % calculate weighting coefficients (Voronoi surfaces <-> measurement points)
W = sparse(diag(w)); % diagonal matrix containing weights
D = I .* diag(1 + n.*(n+1)); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',true); % calculate real-valued SHs using the measurement grid
D = I .* diag(1 + n.*(n+1)); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',true); % calculate real-valued SHs using the measurement grid
%% Calculate HRTF data for field points
if Nmax > 25
......@@ -116,7 +117,7 @@ end
hrtf_arbi = zeros(this.nBins,2*field.nPoints); % columns: LRLRLR...
for ear=1:2
% sh transformation
% SH transformation
freqData_temp = copyData.freqData(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
......@@ -128,7 +129,6 @@ for ear=1:2
% figure
% surf(s,freqData_temp(10,:))
% range extrapolation
if ~isequal(this.dirCoord.r(1),field.r(1))
% calculate range-extrapolated HRTFs
......@@ -139,11 +139,11 @@ for ear=1:2
% reconstruction to spatial data
Yest = ita_sph_base(field,N,'orthonormal',true); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a1).'; % interpolated + range-extrapolated HRTFs
hrtf_arbi(:,ear:2:end) = (Yest*a1(1:(N+1)^2,:)).'; % interpolated + range-extrapolated HRTFs
else
% reconstruction to spatial data
Yest = ita_sph_base(field,Nmax,'orthonormal',true); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a0).'; % interpolated HRTFs
Yest = ita_sph_base(field,N,'orthonormal',true); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a0(1:(N+1)^2,:)).'; % interpolated HRTFs
end
end
......@@ -158,6 +158,11 @@ cThis = this;
cThis.freqData = hrtf_arbi;
cThis.channelCoordinates.sph= sph;
% channelnames coordinates
cThis.channelNames = ita_sprintf('%s (\\theta=%2.0f, \\phi=%2.0f)',...
repmat(['L'; 'R'],cThis.dirCoord.nPoints, 1),...
cThis.channelCoordinates.theta_deg,...
cThis.channelCoordinates.phi_deg);
%% move back to head center
if sArgs.shiftToEar
......
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_coord_flute_mics.m deleted 100644 → 0
View file @ 4663b01d
function [ coord_mics ] = ita_HRTFarc_coord_flute_mics( phi, r_z, d, a, h, gamma, beta, alpha,xOff,yOff)
%TEST_ZILLEKENS_COORD_FLUTE_MICS returns an itaCoordinates containing the position of the
%flutes mics
%
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a, h)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a, h, gamma, beta, alpha)
%
% phi - vector of the angles phi
% r_z - distance of the mics to center of flute
% d - length of boom arm
%
% optional:
% a - distance between joint and boom arm
% h - height of stand to joint in center of turntable (default 0)
%
% optional uncertainties:
% gamma - angle between mic holder and flute (default 0)
% beta - angle between mic stand and mic boom arm (default 0)
% alpha - angle offset of reference point of turntable (default 0)
%
% See also TEST_ZILLEKENS_HT_DHTRAFO TEST_ZILLEKENS_ERROR_DISTANCE_LS
% Author: Stefan Zillekens
% Created: 2013-09-02
if(nargin < 6)
gamma = 0; % angle beetween mic holder and flute
beta = 0; % angle beetween mic stand and mic boom
alpha = 0; % angle offset of reference point of motor
end
if(nargin < 5)
h = 0;
end
if(nargin < 4)
a = 0;
end
% init
npos = numel(phi);
nmic = numel(r_z);
coord_mics = itaCoordinates(npos*nmic);
XYZ = zeros(npos*nmic,3);
for idp = 1:npos;
for idm = 1:nmic;
% homogeneous transformation matrices:
T1 = ita_HRTFarc_ht_DHtrafo(-pi/2+phi(idp)+alpha, h, 0, beta);
T2 = ita_HRTFarc_ht_DHtrafo(+pi/2,a,d,gamma);
T3 = ita_HRTFarc_ht_DHtrafo(0, r_z(idm), 0, 0);
T = T1 * T2 * T3;
XYZ((idp-1)*nmic+idm, :) = T(1:3,4)';
end
end
coord_mics.cart = XYZ;
coord_mics.x = coord_mics.x+xOff;
coord_mics.y = coord_mics.y+yOff;
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_cropMeasurementData.m deleted 100644 → 0
View file @ 4663b01d
function [ output_args ] = ita_HRTFarc_cropMeasurementData( iMS,dataPath )
%ITA_HRTFARC_CROPMEASUREMENTDATA crop of raw measurement data
% This function is intended to be used after a measurement with the HRTF
% arc. During the measurement, the cropping is not done to save time.
% The function assumes a data structure from itaEimar
% It will read all files from dataPath/data, crop them
% It will rename the folder to data_raw. Cropped data will be in
% dataPath/data to maintain backwards compability.
savePath = sprintf('%s/cropped/',dataPath);
dataPathInternal = sprintf('%s/data/',dataPath);
mkdir(savePath)
% load all data, crop and save again
dataDir = dir(dataPathInternal);
wb = itaWaitbar(length(dataDir)-2);
for index = 3:length(dataDir)
data = ita_read(sprintf('%s/%s',dataPathInternal,dataDir(index).name));
data_cropped = iMS.crop(data);
ita_write(data_cropped,sprintf('%s/%s',savePath,dataDir(index).name));
wb.inc;
end
wb.close;
rawPath = sprintf('%s/data_raw/',dataPath);
movefile(dataPathInternal,rawPath);
movefile(savePath,dataPathInternal);
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_error_distance_LS.m deleted 100644 → 0
View file @ 4663b01d
function [ err ] = ita_HRTFarc_error_distance_LS( x, toa, phi, r_z, SR )
%TEST_ZILLEKENS_ERROR_DISTANCE_LS returns the result of a least-mean square
%error
%
% err = TEST_ZILLEKENS_ERROR_DISTANCE_LS( x, toa, phi, r_z )
% err = TEST_ZILLEKENS_ERROR_DISTANCE_LS( x, toa, phi, r_z, SR )
%
% x - containg a vector of LS coordinates, latency sample between DA
% converters, alpha, beta, gamma, d, height, speed of sound
% toa - vector of time of arrivals
% phi - vector of phi position of flute mics
% r_z - vector with distances of mics to center of flute
% SR - sampling rate
%
% See also TEST_ZILLEKENS_COORD_FLUTE_MICS TEST_ZILLEKENS_OPTIMIZE_COORD_ARCLS
% Author: Stefan Zillekens
% Created: 2013-09-09
if nargin < 5
SR = 44100;
end
nopt = 8; % number of values to optimize
npos = numel(phi);
nmic = numel(r_z);
nLS = (numel(x)-nopt)/3;
if numel(toa) ~= nLS*npos*nmic;
error('Dimensions do not fit')
end
%% uncertainties:
a = x(end-7); % distance from joint to boom arm ~4cm
d = x(end-6); % length of boom arm
% height = x(end-4);
c_meas = x(end-5); % Speed of Sound
alpha = x(end-4); % angle offset of reference point of TT ~0
% alpha = 0;
beta = x(end-3); % angle between mic stand and mic boom ~pi/2
gamma = x(end-2); % angle between mic holder and flute ~0
xOff = x(end-1);
yOff = x(end);
% latency = x(end-4:end); % relative changes in latency of D/A converter
latency = zeros(1,8);
%% calculate mic coordinations via homogeneous transforms
height = 0;
coord_mics = ita_HRTFarc_coord_flute_mics(phi,r_z,d,a,height,gamma,beta,alpha,xOff,yOff);
c_meas = 344;
%% error calculation
micsCart = coord_mics.cart;
% tao MxLSxP
toa_reshape = reshape(toa,nmic,nLS,npos);
% only for the measurements from 28/29.05
% toa_reshape(:,8:8:32,:) = 0;
% toa_reshape(:,19:20,:) = 0;
toa_corrected = reshape(toa_reshape,1,nmic*nLS*npos);
%p1 3xMxLSxP
micPos = reshape(micsCart.',3,nmic,1,npos);
micPos = repmat(micPos,[1,1,nLS,1]);
micPos = reshape(micPos,3,nLS*npos*nmic);
%p2 3xMxLSxP
lsPos = reshape(x(1:end-nopt),3,nLS);
lsPos = repmat(lsPos,[1,1,nmic,npos]);
lsPos_reshape = permute(lsPos,[1 3 2 4]);
lsPos = reshape(lsPos_reshape,3,nLS*npos*nmic);
diffMatrix = lsPos-micPos;
dist = sqrt(diffMatrix(1,:).^2+diffMatrix(2,:).^2+diffMatrix(3,:).^2);
dist_reshape = reshape(dist,nmic,nLS,npos);
% dist_reshape(:,8:8:32,:) = 0;
% dist_reshape(:,19:20,:) = 0;
dist_corrected = reshape(dist_reshape,1,nmic*nLS*npos);
err = (dist_corrected-toa_corrected/SR*c_meas).';
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_DHtrafo.m deleted 100644 → 0
View file @ 4663b01d
function [ T ] = ita_HRTFarc_ht_DHtrafo( theta, d, a, alpha )
%TEST_ZILLEKENS_HT_DHTRAFO Homogeneous transformation matrix via
%Denavit-Hartenberg (DH) parameter
%
% T = TEST_ZILLEKENS_HT_DHTRAFO( theta, d, a, alpha )
%
% Rot('z_-1',theta) * Trans([0,0,d]) * Trans([a,0,0]) * Rot('x',aplpha)
%
% Rot(axis,angle) rotation around axis of coordinate system with angle
% Trans([x,y,z]) translation in x, y, z direction.
%
% angles in radiant
%
% EXAMPLE:
% T = test_zillekens_ht_DHtrafo( -pi, 1, 0, pi/8 )
%
% See also TEST_ZILLEKENS_HT_TRANS TEST_ZILLEKENS_HT_ROT
% Author: Stefan Zillekens
% Created: 2013-07-30
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
T = ita_HRTFarc_ht_rot('z',theta) ...
* ita_HRTFarc_ht_trans([0,0,d]) ...
* ita_HRTFarc_ht_trans([a,0,0]) ...
* ita_HRTFarc_ht_rot('x',alpha);
%%
% T = zeros(4);
% T(1:2,1) = [ cos(theta) sin(theta) ];
% T(1:3,2) = [ -cos(alpha)*sin(theta) cos(alpha)*cos(theta) sin(alpha) ];
% T(1:3,3) = [ sin(alpha)*sin(theta) -sin(alpha)*cos(theta) cos(alpha) ];
% T(:,4) = [ a*cos(theta) a*sin(theta) d 1 ];
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_rot.m deleted 100644 → 0
View file @ 4663b01d
function [ T_rot ] = ita_HRTFarc_ht_rot( rot_axis, phi )
%TEST_ZILLEKENS_HT_ROT Homogenous transformation matrix for a rotation
%
% T_rot = TEST_ZILLEKENS_HT_TRANS( 'rot_axis', angle )
%
% rot_axis - rotatation axis { 'x', 'y', 'z' }
% phi - rotation angle in radiant
%
% EXAMPLE:
% T_rot = test_zillekens_ht_rot( 'z', pi )
% % Rotation around z-axis with the angle pi
%
% See also TEST_ZILLEKENS_HT_TRANS TEST_ZILLEKENS_HT_DHTRAFO
% Author: Stefan Zillekens
% Created: 2013-07-30
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
if ~isscalar(phi)
error('Expecting a scalar')
end
T_rot = eye(4);
if rot_axis == 'x'
T_rot(2:3,2:3) = [ cos(phi), -sin(phi); ...
sin(phi), cos(phi) ];
elseif rot_axis == 'y'
T_rot(1:2:3,1:2:3) = [ cos(phi), sin(phi); ...
-sin(phi), cos(phi) ];
elseif rot_axis == 'z'
T_rot(1:2,1:2) = [ cos(phi), -sin(phi); ...
sin(phi), cos(phi) ];
else
error('Rotation axis should be x, y or z.')
end
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_trans.m deleted 100644 → 0
View file @ 4663b01d
function [ T_trans ] = test_rbo_zillekens_ht_trans( coord )
%TEST_ZILLEKENS_HT_TRANS Homogenous transformation matrix for a translation
%
% Ttrans = TEST_ZILLEKENS_HT_TRANS( coord )
% Ttrans = TEST_ZILLEKENS_HT_TRANS( itaCoordinates )
%
% EXAMPLE:
% Ttrans = test_zillekens_ht_trans( [1, 0, 0] )
% % Translation of 1 in x direction
% Ttrans = test_zillekens_ht_trans( itaCoordinates )
%
% See also TEST_ZILLEKENS_HT_ROT TEST_ZILLEKENS_HT_DHTRAFO
% Author: Stefan Zillekens
% Created: 2013-08-22
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
% convert from itaCoordinates
% if isa(coord, 'itaCoordinates')
% coord = coord.cart;
% end
%
% coord = coord(:);
if ~isvector(coord) || length(coord)~= 3;
error('Invalid size of input data. Expecting [x,y,z] or an itaCoordinate')
end
T_trans = eye(4);
T_trans(1:3,4) = coord;
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_measurementScript.m
View file @ 5d1cc16c
%% Init
% Before you start to measure the TFs of the HRTF, please make sure that:
% 0. set ita_preferences: playing/recording device: ASIO Hammerfall
% connect the first channel output to the first channel input
%% calibrate
ms_calibrate = itaMSTF;
ms_calibrate.inputChannels = 11;
ms_calibrate.outputChannels = 19;
ms_calibrate.run_latency;
latency = ms_calibrate.latencysamples;
%%
% MovtecCOM Port: COM3
% midi output device: multiface midi
% 1. turntable is runnig
% 2. you change the sensitivity of the microphones with ita_robocontrol
% 3. you created a itaEimar object. This object has to be reseted for each
% measurement
% 4. load interleaved sweep
%ita_robocontrol('-20dB','Norm','0dBu')
HRTF_Eimar = itaEimar;
outCh = 19:1:82;
inCh = 11:12;
% inCh = 11:14; %flute
% fluet 75mm and 200mm
outAmp = 35;
freqRange = [200 22050];
latency = 619;
iMS_HRTF = itaMSTFinterleaved; % please close it...
iMS_HRTF.inputChannels = inCh;
iMS_HRTF.outputChannels = outCh;
iMS_HRTF.freqRange = freqRange;
iMS_HRTF.latencysamples = latency;
iMS_HRTF.outputamplification = outAmp;
iMS_HRTF.optimize;
iMS_HRTF.skipCrop = 1;
% %% Old interleaved
phiRes = 2.5;
coord = ita_generateSampling_equiangular(phiRes,1); % create coordinates
coord_cut = coord.n(coord.theta_deg == 90);
coord_cut = coord_cut.n(coord_cut.phi_deg <= 180 );
% this prevents that the turntable move to every single direction%
%pause(30);
% coord_cut = coord_cut.n(1);
HRTF_Eimar.measurementPositions = coord_cut;
HRTF_Eimar.waitBeforeMeasurement = 1;
HRTF_Eimar.measurementSetup = iMS_HRTF;
HRTF_Eimar.doSorting = false;
%% Teller rotieren und Kamera
% Camera
% IP http://137.226.61.17/
% benutzername: itacam
HRTF_Eimar.reset
HRTF_Eimar.reference;
HRTF_Eimar.move_turntable(90)
%% Init/ Reset Turntable
% dataPath = [datestr(now,'yyyymmmmdd_HHMM') 'reference_' num2str(phiRes) 'degree_KE3'];
dataPath = [datestr(now,'yyyymmmmdd_HHMM') 'Jan' num2str(phiRes)];
HRTF_Eimar.dataPath = dataPath;
HRTF_Eimar.reset
HRTF_Eimar.reference;
disp('...................................................................')
disp('Reference position reached!')
disp('...................................................................')
%% Run measurement with Eimar
tic
HRTF_Eimar.run
t= toc;
disp('...................................................................')
disp(['Total time: ' num2str(round(t/60*100)/100)])
disp('...................................................................')
% run postprocessing
if (iMS_HRTF.skipCrop == 1)
ita_HRTFarc_cropMeasurementData(iMS_HRTF,dataPath);
function ita_HRTFarc_measurementScript(varargin)
%ita_HRTFarc_measurementScript - +++ example of continuous measurement with the new HRTFarc +++
%
% This is an example of a step-wise measurement. Only use this with the
% turntable and subject rotation
%
% Reference page in Help browser
% <a href="matlab:doc test">doc test</a>
% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% </ITA-Toolbox>
% Author: Jan-Gerrit Richter -- Email: jri@akustik.rwth-aachen.de
% Created: 04-Oct-2018
% set this to 1 if you want to measure the system latency
options.measureLatency = 0;
% init the motor object
iMS = test_itaEimarMotorControl;
% init the measurement object
ms = itaMSTFinterleaved;
if options.measureLatency == 0
ms.latencysamples = 629;
else
ms.inputChannels = 11;
ms.outputChannels = 20;
input('Connect the amp to ch 11')
ms.run_latency;
input('Disconnect the amp')
end
test_rbo_HRTF_meas_peakDetection(dataPath)
ms.outputChannels = (1:64)+18;
% for continuous measurement, three inputs are required, left right and the
% switch
ms.inputChannels = [11 12 13];
%define a freq range
ms.freqRange = [500 22050];
ms.optimize
ms.twait = 0.03;
coords = ita_generateSampling_equiangular(5,5);
coords_cut = coords.n(coords.theta_deg == 90);
iMS.measurementPositions = coords_cut;
saveName = 'test';
iMS.dataPath = saveName;
iMS.reference;
iMS.run;
% always leave the turntable in reference position
iMS.reference
end
\ No newline at end of file
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_measurementScript_continuous.m 0 → 100755
View file @ 5d1cc16c
function ita_HRTFarc_measurementScript_continuous(varargin)
%ita_HRTFarc_measurementScript_continuous - +++ example of continuous measurement with the new HRTFarc +++
%
% Reference page in Help browser
% <a href="matlab:doc test">doc test</a>
% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% </ITA-Toolbox>
% Author: Jan-Gerrit Richter -- Email: jri@akustik.rwth-aachen.de
% Created: 04-Oct-2018
% set this to 1 if you want to measure the system latency
options.measureLatency = 0;
% init the motor object
iMS = test_itaEimarMotorControl;
% init the measurement object
ms = itaMSTFinterleaved;
if options.measureLatency == 0
ms.latencysamples = 629;
else
ms.inputChannels = 11;
ms.outputChannels = 20;
input('Connect the amp to ch 11')
ms.run_latency;
input('Disconnect the amp')
end
ms.outputChannels = (1:64)+18;
% for continuous measurement, three inputs are required, left right and the
% switch
ms.inputChannels = [11 12 13];
%define a freq range
ms.freqRange = [500 22050];
ms.optimize
ms.twait = 0.03;
% the number of repetitions defines the rotation speed
% 64 is ~ 3 minutes
numRepetitions