Skip to content
GitLab
Commit 7455971c authored by Jan-Gerrit Richter's avatar Jan-Gerrit Richter
Browse files

merged last svn commit into the git toolbox

parent 692796ef
Hide whitespace changes
Inline Side-by-side
applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m
View file @ 7455971c
......@@ -35,7 +35,7 @@ classdef itaHRTF < itaAudio
% HRTF_left = getEar(earSide)
%
% itaHRTF Methods (play):
% play_gui(stimulus)
% play_gui(stimulus)
%
% itaHRTF Methods (store):
% audioHRTF = itaHRTF2itaAudio
......@@ -45,16 +45,14 @@ classdef itaHRTF < itaAudio
% ITD = ITD(varargin)
% t0 = meanTimeDelay(varargin)
% ILD = ILD(varargin)
% IACC = IACC(varargin)
%
% itaHRTF Methods (manipulation):
% DTF = calcDTF
% HRTF_int = interp(varargin)
%
% itaHRTF Methods (plot):
% plot_ILD
% plot_ITD(varargin)
% plot_freqSlice(varargin)
% plot_ITD(varargin)
% plot_freqSlice(varargin)
%
% See also:
......@@ -354,16 +352,16 @@ classdef itaHRTF < itaAudio
end
function this = set.openDaff2itaHRTF(this,pathDaff)
handleDaff = DAFFv17('open',pathDaff);
props = DAFFv17('getProperties', handleDaff);
handleDaff = DAFF('open',pathDaff);
props = DAFF('getProperties', handleDaff);
counter = 1;
data = zeros(props.filterLength,props.numRecords*2,'double' ) ;
coordDaff = zeros(props.numRecords,2) ;
for iDir = 1:props.numRecords
data(:,[counter counter+1]) = DAFFv17('getRecordByIndex', handleDaff,iDir)';
%coordDaff(iDir,:) = DAFFv17('getRecordCoords', handleDaff, 'object', iDir)';
coordDaff(iDir,:) = DAFFv17('getRecordCoords', handleDaff, 'data', iDir)';
data(:,[counter counter+1]) = DAFF('getRecordByIndex', handleDaff,iDir)';
%coordDaff(iDir,:) = DAFF('getRecordCoords', handleDaff, 'object', iDir)';
coordDaff(iDir,:) = DAFF('getRecordCoords', handleDaff, 'data', iDir)';
counter= counter+2;
end
......@@ -628,13 +626,17 @@ classdef itaHRTF < itaAudio
phiU = rad2deg(this.phi_Unique);
thetaU = rad2deg(this.theta_Unique);
switch dirID
case 'phi_deg'
case {'phi_deg', 'p'}
slice = this.findnearestHRTF(thetaU,dir_deg);
case 'theta_deg',
case {'theta_deg', 't'}
slice = this.findnearestHRTF(dir_deg,phiU);
end
end
function slice = ss(this,dirID,dir_deg)
slice = this.sphericalSlice(dirID,dir_deg);
end
function HRTFout = getEar(this,earSide)
switch earSide
case 'L',
......@@ -895,9 +897,6 @@ classdef itaHRTF < itaAudio
end
end
%% Florian's Functions
function cThis = interp(this,varargin)
% function this = interp(varargin)
%
% Function to calculate HRTFs for arbitrary field points using a N-th order
......@@ -935,112 +934,8 @@ classdef itaHRTF < itaAudio
% Author: Florian Pausch <fpa@akustik.rwth-aachen.de>
% Version: 2016-02-05
% TODO: check why this is still not working (coordinate assignment???)
sArgs = struct('order',50,'eps',1e-8);
sArgs = ita_parse_arguments(sArgs,varargin,2);
if ~isa(varargin{1},'itaCoordinates'),error('itaHRTF:interp', ' An itaCoordinate object is needed!')
end
field_in = varargin{1};
% only take unique direction coordinates (round to 1deg resolution)
tempfield = unique(round([field_in.phi_deg field_in.theta_deg]),'rows'); % may cause problems with older Matlab versions (<=R2013)!
temp_r = this.dirCoord.r(1);
field = itaCoordinates(size(tempfield,1));
field.r = repmat(temp_r,size(tempfield,1),1);
field.phi_deg = tempfield(:,1);
field.theta_deg = tempfield(:,2);
N = sArgs.order;
epsilon = sArgs.eps; % regularization parameter
k = this.wavenumber; % wave number
k(1) = eps;
% add eps to avoid NaN's
Nmax = floor(sqrt(this.nDirections/4)-1);
% 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
hankel_r0 = ita_sph_besselh(1:Nmax,2,kr0);
hankel_r1 = ita_sph_besselh(1:Nmax,2,kr1);
hankel_div = hankel_r1 ./ hankel_r0;
hankel_rep = hankel_div(:,1);
end
dweights = 1 + (0:Nmax).*((0:Nmax)+1); % calculate regularization weights
dweights_rep = zeros(sum(2*(0:Nmax)'+1),1);
dweights_rep(1)=dweights(1);
counter = 2;
for n=1:Nmax
nTimes = 2*n+1;
dweights_rep(counter:counter+nTimes-1)=dweights(n+1)*ones(nTimes,1);
if ~isequal(this.dirCoord.r(1),field.r(1))
hankel_rep=[hankel_rep, repmat(hankel_div(:,n),1,2*n+1)];
end
counter = counter + nTimes;
end
%% Weights
[~,vWeights]= this.dirCoord.spherical_voronoi; % calculate weighting coefficients (Voronoi surfaces <-> measurement points)
W = diag(vWeights); % diagonal matrix containing weights
D = diag(dweights_rep); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',false); % calculate real-valued SHs using the measurement grid
%% Calculate HRTF data for field points
if Nmax > 25
ita_disp('[itaHRTF.interp] Be patient...')
end
% init.
hrtf_arbi = zeros(this.nBins,2*field.nPoints); % columns: LRLRLR...
for ear=1:2
% calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
freqData_temp = this.freqData(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
%a0 = (Y.'*W*Y + epsilon*D) \ Y.' *
%this.freqData(:,ear:2:end).'; % fpa version
%a0 = pinv(Y)* this.freqData(:,ear:2:end).'; % jck version
if ~isequal(this.dirCoord.r(1),field.r(1))
% calculate range-extrapolated HRTFs
a1 = a0 .* hankel_rep.';
Yest = ita_sph_base(field,N,'orthonormal',false); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a1).'; % interpolated + range-extrapolated HRTFs
else
Yest = ita_sph_base(field,Nmax,'orthonormal',false); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a0).'; % interpolated HRTFs
end
end
% set new direction coordinates
sph = zeros(field.nPoints*2 ,3);
sph(1:2:end,:) = field.sph;
sph(2:2:end,:) = field.sph;
% write new HRTF data set
cAudio = itaAudio(hrtf_arbi, 44100, 'freq');
cAudio.channelCoordinates.sph= sph;
cThis = itaHRTF(cAudio);
cThis.freqData = hrtf_arbi;
if ~isequal(cThis.dirCoord.r(1),field.r(1))%???
cThis.dirCoord.r = field.r;
end
if N > 25
ita_disp('[itaHRTF.interp] ...calculation finished!')
end
end
function this = smooth_linphase(this,varargin)
% function this = smooth_linphase(varargin)
......@@ -1221,10 +1116,10 @@ classdef itaHRTF < itaAudio
end
thetaC_deg = rad2deg(thisS.theta_Unique);
phiC_deg = rad2deg(thisS.phi_Unique);
phiC_deg = sort(mod(round(rad2deg(thisS.phi_Unique)),360));
nTheta = numel(thetaC_deg);
nPhi = numel(phiC_deg);
coord = round(reshape(thisS.dirCoord.phi_deg,nTheta,nPhi));
coord = reshape(mod(round(thisS.dirCoord.phi_deg),360),nTheta,nPhi);
[~, idxC] = sort(coord,2);
[~, idxCT] = unique(thisS.dirCoord.theta_deg);
......@@ -1245,7 +1140,7 @@ classdef itaHRTF < itaAudio
set(gcf,'Position',[10 50 position(3:4)*0.85]);
if strcmp(sArgs.method,'phase_delay') && ischar(sArgs.filter) % frequency dependent ITD
pcolor(phiC_deg,this.freqVector,ITD)
title(strcat('\phi = ', num2str(round(thetaC_deg)), ''))
title(strcat('\phi = ', num2str(round(thetaC_deg)), '�'))
shading flat
colorbar
......@@ -1286,7 +1181,7 @@ classdef itaHRTF < itaAudio
xlabel('Azimuth Angle in Degree');
ylabel('ITD in Milliseconds');
set(gca,'xTick',0:30:360)
legend(ita_sprintf('%i', round(thetaC_deg)))
legend(ita_sprintf('%i�', round(thetaC_deg)))
end
end
......@@ -1337,11 +1232,11 @@ classdef itaHRTF < itaAudio
earSidePlot = sArgs.earSide;
if numel(phiC_deg)>1,
xData = phiC_deg;
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) ''];
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '�'];
strXlabel = '\phi in Degree';
else
xData = thetaC_deg;
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) ''];
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '�'];
strXlabel = '\theta in Degree';
end
......@@ -1382,6 +1277,165 @@ classdef itaHRTF < itaAudio
grid on;set(gca,'layer','top')
end
%% Jan's Functions
function writeDAFFFile(this, filePath)
% writes DAFF file to hard disc
%
% Input: filePath / fileName (string)
%
% Required: openDAFF matlab executables
%
% Output: none
if nargin == 2
fileName = filePath;
else
fileName = [ 'HRTF_Length' int2str(this.nSamples) '_' int2str(this.resAzimuth) 'x' int2str(this.resElevation) '.daff'];
end
%..............................................................
% Comment: Angles are not exact - improve it if you like to :)
precision = 3;
nThetaU = numel(unique(round(this.channelCoordinates.theta_deg*10)./10));
nPhiU = numel(unique(round(this.channelCoordinates.phi_deg*10)./10));
resElevation= round(median(diff(rad2deg(this.theta_Unique)))*10^precision)/10^precision;
rangeEl = round([min(rad2deg(this.theta_Unique)) max(rad2deg(this.theta_Unique))]*10^precision)/10^precision;
rangeElnew = [min(rangeEl) min(rangeEl)+(nThetaU-1)*resElevation];
resAzimuth= round(median(diff(rad2deg(this.phi_Unique)))*10^precision)/10^precision;
rangeAz = round([min(rad2deg(this.phi_Unique)) max(rad2deg(this.phi_Unique))]*10^precision)/10^precision;
rangeAznew = [min(rangeAz) min(rangeAz)+(nPhiU-1)*resAzimuth];
%..............................................................
%% config values
threshold_db = -20;
pre_taps = 12;
window_length = 128;
delay = 10;
gPeakL = 0;
gPeakR = 0;
gRangeStart = inf;
gRangeEnd = 0;
gRangeStartHit = [];
gRangeEndHit = [];
gPower = 0;
% Measurement distance [m], rounded on one digit
distance = delay / 44100 * 340;
distance = round(distance * 100)/100;
% this = ita_time_shift(this,delay);
% metadata.object = 'ITA Kunstkopf, an artificial head developed and designed at the Institute of Technical Acoustics (ITA), RWTH Aachen University';
% metadata.copyright = '(c) Copyright Institute of Technical Acoustics (ITA), RWTH Aachen University, Germany';
% metadata.contact = 'Frank Wefers (fwe@akustik.rwth-aachen.de)';
% metadata.environment = 'semi-anechoic chamber';
% metadata.session = 'Mess01: Farfield HRIRs, Tobias Lentz, 2001';
% metadata.processing = 'Loudspeaker deconvolved. Windowed using peak-oriented tukey window (0.1) to length 96 taps.';
%% find filter ranges
peak = max(abs(this.timeData));
threshold = peak * 10^(threshold_db/20);
timeData = abs(this.timeData);
for index = 1:this.nChannels
kTemp = find(timeData(:,index) >=threshold(index));
if (isempty(kTemp))
error('Impulse response is completely below the threshold')
end
r(index,:) = [ kTemp(1) kTemp(end) ];
k{index} = kTemp;
end
offset = max(r(:,1)-pre_taps,1);
for index = 1:this.nChannels
chObj = ita_time_window(this.ch(index),[offset(index) offset(index)+ window_length-1],'samples');
timeData(:,index) = chObj.timeData;
end
copyObj = this;
copyObj.timeData = timeData;
copyCoordinates = copyObj.channelCoordinates;
copyCoordinates = copyCoordinates.build_search_database;
%% global peaks
peak = max(abs(copyObj.timeData));
threshold = peak * 10^(threshold_db/20);
timeData = abs(copyObj.timeData);
for index = 1:copyObj.nChannels
kTemp = find(timeData(:,index) >=threshold(index));
if (isempty(kTemp))
error('Impulse response is completely below the threshold')
end
r(index,:) = [ kTemp(1) kTemp(end) ];
k{index} = kTemp;
end
gRangeStart = min(r(r(:,1) > delay,1));
gRangeEnd = max(r(r(:,2) < 1000,2));
nDegree = log2(window_length);
% create dataset
dataset = daff_create_dataset(...
'alphares', resAzimuth, ...
'alpharange',rangeAznew,...
'betares', resElevation, ...
'betarange', rangeElnew, ...
'channels', 2);
% set samplerate and metainfo
dataset.samplerate = this.samplingRate;
dataset.metadata.desc = 'Dummy HRTF';
dataset.metadata.delay_samples = int32(delay);
dataset.metadata.measurement_distance_meters = distance;
dataset.metadata.creation_date = datestr(now, 'yyyy-mm-dd HH:MM');
% assign data
for i=1:dataset.numrecords
%dataset.records{i}.data = this.data(:,(2*i-1:2*i))';
alpha = dataset.records{i}.alpha;
beta = dataset.records{i}.beta;
% get and save data
data = this.findnearestHRTF(180-beta,alpha);
% if ceil(nDegree) ~= floor(nDegree),data.nSamples = ceil(nDegree);
% end
dataset.records{i}.data = data.timeData(gRangeStart:gRangeStart+window_length-1,:).';
% Optionally you can supply individual metadata for the records
dataset.records{i}.metadata.filename = fileName;
end
% write file
daff_write( 'filename', fileName, ...
'content', 'IR', ...
'dataset', dataset, 'verbose');
end
end
methods(Hidden = true)
function sObj = saveobj(this)
......
applications/Binaural-HRTF/HRTF_class/@itaHRTF/writeDAFFFile.m
View file @ 7455971c
function writeDAFFFile( this, filePath, metadata_user )
function writeDAFFFile( this, file_path, metadata_user )
% Exports itaHRTF to a DAFF file
%
% Input: filePath / fileName (string) [optional]
% Input: file_path (string) [optional]
% user metadata (struct created with daff_add_metadata) [optional]
%
% Required: OpenDAFF matlab scripts, http://www.opendaff.org
......@@ -15,9 +15,9 @@ if nargin >= 3
end
hrtf_variable_name = inputname( 1 );
fileName = [ hrtf_variable_name '_' int2str( this.nSamples ) 'samples_' int2str( this.resAzimuth ) 'x' int2str(this.resElevation) '.daff'];
file_name = [ hrtf_variable_name '_' int2str( this.nSamples ) 'samples_' int2str( this.resAzimuth ) 'x' int2str( this.resElevation ) '.daff'];
if nargin >= 2
fileName = filePath;
file_name = file_path;
end
if nargin == 0
......@@ -32,12 +32,11 @@ if strcmp( this.domain, 'freq' )
ct_indicator = 'dft';
end
file_path = '';
file_path_base = strsplit( fileName, '.' );
if ~strcmp( file_path_base(end), 'daff' )
file_path = strjoin( [ file_path_base(:) ct_indicator 'daff' ], '.' );
[ file_path, file_base_name, file_suffix ] = fileparts( file_name );
if ~strcmp( file_suffix, '.daff' )
file_path = fullfile( file_path, strjoin( {file_base_name file_suffix 'v17' ct_indicator 'daff' }, '.' ) );
else
file_path = strjoin( [ file_path_base(1:end-1) ct_indicator 'daff' ], '.' );
file_path = fullfile( file_path, strjoin( {file_base_name 'v17' ct_indicator 'daff'}, '.' ) );
end
......@@ -47,8 +46,8 @@ theta_start_deg = rad2deg( min( this.channelCoordinates.theta ) );
theta_end_deg = rad2deg( max( this.channelCoordinates.theta ) );
theta_num_elements = size( unique( this.channelCoordinates.theta ), 1 );
phi_start_deg = rad2deg( min( this.channelCoordinates.phi ) );
phi_end_deg = rad2deg( max( this.channelCoordinates.phi ) );
phi_start_deg = rad2deg( min( mod( this.channelCoordinates.phi, 2*pi ) ) );
phi_end_deg = rad2deg( max( mod( this.channelCoordinates.phi, 2*pi ) ) );
phi_num_elements = size( unique( this.channelCoordinates.phi ), 1 );
assert( phi_num_elements ~= 0 );
......
applications/Binaural-HRTF/HpTF_class/itaHpTF_Audio.m
View file @ 7455971c
......@@ -60,6 +60,20 @@ classdef itaHpTF_Audio < itaHpTF
this = varargin;
elseif isa(varargin{1},'itaHpTF_MS')
this.init = varargin{1};
elseif nargin ==1 && isstruct(varargin{1}) % only for loading
obj = varargin{1};
this.data = obj.data;
this.signalType = 'energy';
% additional itaHRTF data
objFNsaved = this.propertiesSaved;
objFNload = this.propertiesLoad;
for i1 = 1:numel(objFNload)
this.(objFNload{i1}) = obj.(objFNsaved{i1});
end
end
end
end
......@@ -254,11 +268,13 @@ classdef itaHpTF_Audio < itaHpTF
end
function result = propertiesSaved
result = {'TF','fLower','fUpper','method','normalized','smoothing'};
result = {'nameHP','nameMic','nameSubj','repeat','mic','savePath',...
'TF','fLower','fUpper','method','normalized','smoothing'};
end
function result = propertiesLoad
result = {'mTF','m_fUpper','m_fLower','mMethod','mNormalized','mSmoothing'};
result = {'nameHP','nameMic','nameSubj','repeat','mic','savePath',...
'mTF','m_fUpper','m_fLower','mMethod','mNormalized','mSmoothing'};
end
function result = propertiesMethod
......
applications/Binaural-HRTF/HpTF_class/ita_itaHpTF_examples.m 0 → 100644
View file @ 7455971c
%% measure HpTF with GUI
HpTF_ms = itaHpTF_MS; % uses output and input channels 1,2
HpTF_subj = HpTF_ms.run %#ok<NOPTS> run measurement
HpTF_subj.TF.pf % show measured transfer functions
%% measure HpTF with a measurement setup
in = 1:2; % input channels
out = 1:2; % output channels
amp = -40; % output amplification
fftDeg = 16; % define fftDegree
MSTF = itaMSTF; % create measurement setup
MSTF.fftDegree = fftDeg;
MSTF.inputChannels = in;
MSTF.outputChannels = out;
MSTF.outputamplification = amp;
HpTF_ms = itaHpTF_MS(MSTF); % init HpTF measurement object
HpTF_ms.nameHP = 'HD 650'; % name of the headphones
HpTF_ms.nameMic = 'KE 3'; % name of the headphones
%HpTF_ms.mic % [transfer functions of the microphones]
HpTF_ms.nameSubj = name; % name of the subject
HpTF_ms.repeat = 4; % repeatitions [8]
HpTF_subj = HpTF_ms.run; % run measurement (press any button to continue
HpTF_subj.TF.pf % show measured transfer functions
%% calculate equalization curve (prepare HpTF)
HpTF_subj.method = 'mSTD'; % choose method [mSTD, mean, max]
HpTF_subj.smoothing = 1/6; % smoothing [1/6 octave bands]
HpTF_subj.fLower = 200; % [100 Hz] lowest frequency for the headphones (smooth spectrum below)
HpTF_subj.fUpper = 18000; % [18 kHz] highest freq for the headphones (regularization)
HpTF_eq = HpTF_subj.HP_equalization; % calculate equalization curve
HpTF_eq.pf
applications/HRTFarc/ita_HRTFarc_postprocessing.m 0 → 100644
View file @ 7455971c
function varargout = ita_HRTFarc_postprocessing(varargin)
%% parse arguments
sArgs = struct('dataPath',[],'savePath',[],'saveName',[],'path_ref',...
[],'ref_name',[],'tWin',[],'samples',[],'flute',[],'phiAdd',0,'eimar',true,'refIsCropped',0);
sArgs = ita_parse_arguments(sArgs,varargin);
dataPath = sArgs.dataPath; % Pfad des Ordners mit den HRTF Rohdaten
path_ref = sArgs.path_ref; % Pfad der Referenzmessung linkes Mikro
savePath = sArgs.savePath; % Pfad angeben, falls gefensterte HRTF gespeichert werden soll
saveName = sArgs.saveName; % Dateiname (ohne .ita Endung)
if isempty(sArgs.samples),
if numel(sArgs.tWin) ==1, t2 = sArgs.tWin(1); t1 = t2/1.25;
else t1 = sArgs.tWin(1); t2 = sArgs.tWin(2);
end
else
SR = 44100;
if numel(sArgs.samples) ==1, t2 = sArgs.samples(1)/SR; t1 = t2/1.25;
else t1 = sArgs.samples(1); t2 = sArgs.samples(2);
end
end
%% coord & userData
allFolders = dir(fullfile(dataPath,'*.ita'));
numAzAngle = size(allFolders,1);
elAngle_arc = sArgs.flute.theta;
coordHRTF = itaCoordinates;
coordHRTF.sph = zeros(numel(elAngle_arc)*numAzAngle ,3);
coordHRTF.r = ones(numel(elAngle_arc)*numAzAngle,1)*1.2;
%% Post processing Window
HRTF_TMP_L = itaAudio(numAzAngle,1);
HRTF_TMP_R = itaAudio(numAzAngle,1);
channelCounter = 1;
currentPath = dataPath;
%% Main processing Ref and Window
wb = itaWaitbar(numAzAngle, 'calculate HRTF', {'azimuth'});
% reference
% first determine how many channels the measurement has
currentDataEnding = num2str(1);
currentData = ita_merge(ita_read([currentPath filesep currentDataEnding '.ita']));
numChannels = currentData.nChannels;