Commit ad3f6de7 authored by Michael Kohnen's avatar Michael Kohnen

added metadata handling to itaHRTF

parents 9b2df9f3 2a77cc7c
......@@ -19,6 +19,11 @@ function cThis = interp(this,varargin)
% set to 1 if no range extrapolation is required
% order ... order of spherical harmonics matrix (default: 50)
% epsilon ... regularization coefficient (default: 1e-8)
% shiftToEar to a shift to approximate ear position to
% improve sh transformation (see [2])
% shiftAxis shift along this axis ('x','y' (default),'z')
% shiftOffset shift ears (L - R) by these values
% (default: [-0.0725 0.0725])
%
% OUTPUT:
% itaHRTF object
......@@ -26,28 +31,34 @@ function cThis = interp(this,varargin)
% .timeData: interpolated / range-extrapolated HRIRs for defined field points
% .dirCoord: itaCoordinates object
%
% Required: SphericalHarmonics functions of ITA Toolbox
%
%
% [1] Pollow, Martin et al., "Calculation of Head-Related Transfer Functions
% for Arbitrary Field Points Using Spherical Harmonics Decomposition",
% Acta Acustica united with Acustica, Volume 98, Number 1, January/February 2012,
% pp. 72-82(11)
%
% [2] Richter, Jan-Gerrit et al. "Spherical harmonics based hrtf datasets:
% Implementation and evaluation for real-time auralization",
% 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
% TODO: check why this is still not working (coordinate assignment???)
sArgs = struct('order',50,'eps',1e-5);
sArgs = struct('order',50,'eps',1e-5,'shiftToEar',false,'shiftAxis','y','shiftOffset',[-0.0725 0.0725]);
sArgs = ita_parse_arguments(sArgs,varargin,2);
if ~isa(varargin{1},'itaCoordinates'),error('itaHRTF:interp', ' An itaCoordinate object is needed!')
if isempty(varargin) || ~isa(varargin{1},'itaCoordinates')
error('itaHRTF:interp', ' An itaCoordinate object is needed!')
end
field_in = varargin{1};
% only take unique direction coordinates (round to 0.01deg resolution)
tempfield = unique(round([field_in.phi_deg*100 field_in.theta_deg*100]),'rows'); % may cause problems with older Matlab versions (<=R2013)!
tempfield = tempfield./100;
temp_r = this.dirCoord.r(1);
temp_r = field_in.r(1);
field = itaCoordinates(size(tempfield,1));
field.r = repmat(temp_r,size(tempfield,1),1);
field.phi_deg = tempfield(:,1);
......@@ -67,46 +78,34 @@ 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_r0 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr0);
hankel_r1 = ita_sph_besselh(ita_sph_linear2degreeorder(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
nSH = (Nmax+1).^2;
I = sparse(eye(nSH));
n = ita_sph_linear2degreeorder(1:round(nSH)).';
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)];
%% move data from earcenter
copyData = this;
if sArgs.shiftToEar
ear_d = sArgs.shiftOffset;
for ear=1:2
movedData = moveFullDataSet(this.ch(ear:2:this.nChannels),sArgs,ear_d(ear));
copyData.freqData(:,ear:2:copyData.nChannels) = movedData;
end
counter = counter + nTimes;
end
%% move data from earcenter
ear_d = [-0.07 0.07];
%% Weights
[~,w]= this.dirCoord.spherical_voronoi; % calculate weighting coefficients (Voronoi surfaces <-> measurement points)
% sG_full = ita_sph_sampling_equiangular(73,144,'theta_type','[]','phi_type','[)');
% % sG_full = ita_sph_sampling_gaussian(71);
% sG = sG_full;
% isOnArc = sG.theta < 2.75;
% sG.cart = sG.cart(isOnArc,:);
% sG.weights = sG.weights(isOnArc);
% % rescale the weights for the sum to be 4*pi again
% w = sG.weights .* 4.*pi ./ sum(sG.weights);
W = sparse(diag(w)); % diagonal matrix containing weights
D = sparse(diag(dweights_rep)); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'real'); % 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
......@@ -117,73 +116,58 @@ end
hrtf_arbi = zeros(this.nBins,2*field.nPoints); % columns: LRLRLR...
for ear=1:2
freqData_temp = this.freqData(:,ear:2:end);
% sh transformation
freqData_temp = copyData.freqData(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
% newCoords = this.dirCoord;
% newCoords.y = newCoords.y + ear_d(ear);
%
% data.channelCoordinates = newCoords;
% data.freqData = freqData_temp;
% data.freqVector = this.freqVector;
% data.c_meas = 344;
% data = process_result_delay_correction(data,this.dirCoord.r(1));
% % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
%
% freqData_temp = data.freqData;
% Y = ita_sph_base(data.channelCoordinates,Nmax,'real');
% %% test the sh transformation results by plotting both spatial and sh data with surf
% s = itaSamplingSph(field);
% s.nmax = Nmax;
% figure
% surf(s,a0(:,10))
% figure
% surf(s,freqData_temp(10,:))
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
% range extrapolation
if ~isequal(this.dirCoord.r(1),field.r(1))
% calculate range-extrapolated HRTFs
a1 = a0 .* hankel_rep.';
Yest = ita_sph_base(field,N,'real'); % use real-valued SH's
%%% test here to see extrapolation results in spatial domain
% surf(s,a1(:,10))
% 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
else
Yest = ita_sph_base(field,Nmax,'real'); % use real-valued SH's
% 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
end
end
%% move back to head center
% todo
% for ear=1:2
%
% newCoords = field;
% newCoords.y = newCoords.y - ear_d(ear);
%
% freqData = hrtf_arbi(:,ear:2:end);
%
%
% data.channelCoordinates = newCoords;
% data.freqData = freqData;
% data.freqVector = this.freqVector;
% data.c_meas = 344;
% data = process_result_delay_correction(data,this.dirCoord.r(1));
%
% hrtf_arbi(:,ear:2:end) = data.freqData;
%
% 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 = this;
cThis.freqData = hrtf_arbi;
cThis.channelCoordinates.sph= sph;
cThis = itaHRTF(cAudio);
cThis.freqData = hrtf_arbi;
%% move back to head center
if sArgs.shiftToEar
ear_d_back = -ear_d;
% movedData = zeros(size(hrtf_arbi));
for ear=1:2
movedData = moveFullDataSet(cThis.ch(ear:2:cThis.nChannels),sArgs,ear_d_back(ear));
cThis.freqData(:,ear:2:cThis.nChannels) = movedData;
end
end
if ~isequal(cThis.dirCoord.r(1),field.r(1))%???
cThis.dirCoord.r = field.r;
......@@ -195,15 +179,44 @@ end
end
function result = process_result_delay_correction(result, target_d)
% shift every measurement point to target_d by applying a
% phase shift to the channel: (simplified!)
freq_L = result.freqVector;
add_phase_L = (result.channelCoordinates.r - target_d)* (freq_L./result.c_meas)' .* 2.*pi;
function [ data ] = moveFullDataSet(data,options,offsetShift)
fullCoords = data.channelCoordinates;
freqVector = data.freqVector;
shiftedData = zeros(size(data.freqData));
axis = options.shiftAxis;
for index = 1:length(freqVector)
shiftedData(index,:) = moveHRTF(fullCoords,data.freqData(index,:),freqVector(index),axis,offsetShift);
end
data = shiftedData;
end
function [data] = moveHRTF(s, data, frequency, axis, offset)
% the offset is given in m
result.freqData = result.freqData .* exp(1i.*add_phase_L');
origAxis = s.r;
if (size(data,2) > size(data,1))
data = data.';
end
offset = real(offset); % ??
switch axis
case 'x'
s.x = s.x + offset;
case 'y'
s.y = s.y + offset;
case 'z'
s.z = s.z + offset;
end
result.channelCoordinates.r = target_d;
newAxis = s.r;
k = 2*pi*frequency/340;
% the phase is moved by the difference of the axis points
data = data .* exp(1i*k*(newAxis - origAxis));
% amplitude manipulation did not yield better results
% data = data .* newAxis ./ origAxis;
end
......@@ -117,10 +117,21 @@ classdef itaHRTF < itaAudio
methods % Special functions that implement operations that are usually performed only on instances of the class
%% Input
function this = itaHRTF(varargin)
this = this@itaAudio();
iniAudio = [];
% initialize itaHRTF with itaAudio properties (only for nargin == 1)
if nargin > 1 || (nargin == 1 && (ischar(varargin{1}) || isa(varargin{1},'itaAudio')))
iniAudio = [];
elseif nargin == 1 && isstruct(varargin{1})
fNames = {'domain','data','signalType','samplingRate'};
for idxFN = 1:numel(fNames)
iniAudio.(fNames{idxFN}) = varargin{1}.(fNames{idxFN});
end
end
this = this@itaAudio(iniAudio);
if nargin >1
% itaAudio input
TF_types = this.propertiesTF_type;
for iTF = 1:numel(TF_types)
......@@ -148,6 +159,7 @@ classdef itaHRTF < itaAudio
end
elseif nargin == 1
if isa(varargin{1},'itaHRTF')
this = varargin{1};
......@@ -170,13 +182,17 @@ classdef itaHRTF < itaAudio
this.dirCoord.sph = this.mCoordSave;
% saving channelNames in itaHRTF does not work at the
% moment
for iCh = 1:this.dimensions
this.channelNames{iCh} = this.mChNames(iCh,:);
end
this.channelNames = cellstr(this.mChNames);
elseif isa(varargin{1},'itaAudio')
this.itaAudio2itaHRTF = varargin{1};
end
elseif ischar(varargin{1}) % openDaff/ sofa/ hdf5 input
if strfind(lower(varargin{1}),'.daff'), this.openDAFF2itaHRTF = varargin{1};
elseif strfind(lower(varargin{1}),'.hdf5'), this.hdf2itaHRTF = varargin{1};
elseif strfind(lower(varargin{1}),'.sofa'), this.sofa2itaHRTF = varargin{1};
end
end
end
end
......@@ -373,7 +389,7 @@ classdef itaHRTF < itaAudio
counter= counter+2;
end
tempMetadata=DAFFv17('getMetadata', handleDaff);
metadata = DAFFv17('getMetadata', handleDaff);
catch
disp( 'Could not read DAFF file right away, falling back to old version and retrying ...' );
......@@ -388,7 +404,27 @@ classdef itaHRTF < itaAudio
counter = 1;
data = zeros(props.filterLength,props.numRecords*2,'double' ) ;
coordDaff = zeros(props.numRecords,2) ;
tempMetadata=DAFFv15('getMetadata', handleDaff);
tempMetadata = DAFFv15('getMetadata', handleDaff);
% Convert old-style metadata format to v17.
names = fieldnames( tempMetadata );
for k = 1:numel( tempMetadata )
switch class(tempMetadata.(names{k}))
case 'logical'
datatype='bool';
case 'char'
datatype='string';
case 'double'
if rem(tempMetadata.(names{k}),1)==0
datatype='int';
else
datatype='float';
end
end
metadata = daffv17_add_metadata( metadata,cell2mat(names(k)),datatype,tempMetadata.(names{k}) );
end
for iDir = 1:props.numRecords
data(:,[counter counter+1]) = DAFFv15( 'getRecordByIndex', handleDaff,iDir )';
coordDaff(iDir,:) = DAFFv15( 'getRecordCoords', handleDaff, 'data', iDir )';
......@@ -396,6 +432,7 @@ classdef itaHRTF < itaAudio
end
end
<<<<<<< HEAD
% Proceed (version independent)
names=fieldnames(tempMetadata);
for k=1:(numel(names))
......@@ -415,36 +452,36 @@ classdef itaHRTF < itaAudio
end
metadata=daffv17_add_metadata(metadata,cell2mat(names(k)),datatype,tempMetadata.(names{k}));
end
phiM = coordDaff(:,1)*pi/180;
%phiM = mod(coordDaff(:,1),360)*pi/180;
%if ~isempty(find(0<coordDaff(:,2),1,'first'))
thetaM = coordDaff(:,2)*pi/180;
%thetaM = mod(180-(coordDaff(:,2)+90),180)*pi/180;
%else
% thetaM = coordDaff(:,2)*pi/180;
%end
radius = ones(props.numRecords,1);
chCoord = itaCoordinates;
chCoord.sph = ones(size(data,2),3);
chCoord.phi(1:2:2*props.numRecords) = phiM;
chCoord.phi(2:2:2*props.numRecords) = phiM;
chCoord.theta(1:2:2*props.numRecords) = thetaM;
chCoord.theta(2:2:2*props.numRecords) = thetaM;
this.mMetadata = metadata;
this.data = data;
this.mDirCoord = itaCoordinates([radius thetaM phiM],'sph');
this.channelCoordinates = chCoord;
this.mEarSide = repmat(['L'; 'R'],props.numRecords, 1);
this.signalType = 'energy';
% channelnames coordinates
this.channelNames = ita_sprintf('%s ( %2.0f, \\theta= %2.0f)',...
this.mEarSide , this.channelCoordinates.theta_deg, this.channelCoordinates.phi_deg);
end
phiM = coordDaff(:,1)*pi/180;
%phiM = mod(coordDaff(:,1),360)*pi/180;
%if ~isempty(find(0<coordDaff(:,2),1,'first'))
thetaM = coordDaff(:,2)*pi/180;
%thetaM = mod(180-(coordDaff(:,2)+90),180)*pi/180;
%else
% thetaM = coordDaff(:,2)*pi/180;
%end
radius = ones(props.numRecords,1);
chCoord = itaCoordinates;
chCoord.sph = ones(size(data,2),3);
chCoord.phi(1:2:2*props.numRecords) = phiM;
chCoord.phi(2:2:2*props.numRecords) = phiM;
chCoord.theta(1:2:2*props.numRecords) = thetaM;
chCoord.theta(2:2:2*props.numRecords) = thetaM;
this.mMetadata = metadata;
this.data = data;
this.mDirCoord = itaCoordinates([radius thetaM phiM],'sph');
this.channelCoordinates = chCoord;
this.mEarSide = repmat(['L'; 'R'],props.numRecords, 1);
this.signalType = 'energy';
% channelnames coordinates
this.channelNames = ita_sprintf('%s ( %2.0f, \\theta= %2.0f)',...
this.mEarSide , this.channelCoordinates.theta_deg, this.channelCoordinates.phi_deg);
end
function this = set.init(this,var)
% TO DO !!!!!!!!!!!!!!!!!!!!!!!!!!!
......@@ -514,6 +551,7 @@ classdef itaHRTF < itaAudio
% data
% the data is saved as positions x channel x filterdata
this.samplingRate = handleSofa.Data.SamplingRate;
data = zeros(size(handleSofa.Data.IR,3),numPositions*2);
data(:,1:2:numPositions*2) = squeeze(handleSofa.Data.IR(:,1,:)).';
......@@ -645,6 +683,10 @@ classdef itaHRTF < itaAudio
%HRTFout = this.direction(idxCoord);
end
function this = buildsearchdatabase(this)
this.dirCoord = this.dirCoord.build_search_database;
end
function obj = direction(this, idxCoord)
idxDir = zeros(numel(idxCoord)*2,1);
idxDir(1:2:numel(idxCoord)*2,:) = 2*idxCoord-1;
......@@ -658,14 +700,14 @@ classdef itaHRTF < itaAudio
end
function thetaUni = theta_Unique(this,varargin)
thetaUni = unique(this.dirCoord.theta);
thetaUni = uniquetol(this.dirCoord.theta,eps);
if nargin == 2
thetaUni = unique(this.dirCoord.theta,'stable');
end
end
function phiUni = phi_Unique(this,varargin)
phiUni = unique(this.dirCoord.phi);
phiUni = uniquetol(this.dirCoord.phi,eps);
if nargin == 2
phiUni = unique(this.dirCoord.phi,'stable');
end
......@@ -679,18 +721,44 @@ classdef itaHRTF < itaAudio
phiUni = rad2deg(phi_Unique(this,varargin));
end
function slice = sphericalSlice(this,dirID,dir_deg)
function slice = sphericalSlice(this,dirID,dir_deg,exactSearch)
% dir in degree
% dirID [phi, theta]
if ~exist('exactSearch','var')
exactSearch = 0;
end
phiU = rad2deg(this.phi_Unique);
thetaU = rad2deg(this.theta_Unique);
switch dirID
case {'phi_deg', 'p'}
slice = this.findnearestHRTF(thetaU,dir_deg);
case {'theta_deg', 't'}
slice = this.findnearestHRTF(dir_deg,phiU);
if ~exactSearch
phiU = rad2deg(this.phi_Unique);
thetaU = rad2deg(this.theta_Unique);
switch dirID
case {'phi_deg', 'p'}
slice = this.findnearestHRTF(thetaU,dir_deg);
case {'theta_deg', 't'}
slice = this.findnearestHRTF(dir_deg,phiU);
end
else
earCoords = this.getEar('L').channelCoordinates;
switch dirID
case {'phi_deg', 'p'}
phiValues = uniquetol(earCoords.phi_deg);
[~,index] = min(abs(phiValues - dir_deg));
exactPhiValue = phiValues(index);
tmp = earCoords.n(earCoords.phi_deg == exactPhiValue);
thetaU = tmp.theta_deg;
slice = this.findnearestHRTF(thetaU,dir_deg);
case {'theta_deg', 't'}
thetaValues = uniquetol(earCoords.theta_deg);
[~,index] = min(abs(thetaValues - dir_deg));
exactThetaValue = thetaValues(index);
tmp = earCoords.n(earCoords.theta_deg == exactThetaValue);
phiU = tmp.phi_deg;
slice = this.findnearestHRTF(dir_deg,phiU);
end
end
end
function slice = ss(this,dirID,dir_deg)
......@@ -762,13 +830,15 @@ classdef itaHRTF < itaAudio
end
function surf(varargin)
sArgs = struct('pos1_data','itaHRTF', 'earSide', 'L', 'freq' , 5000,'type','directivity','log',1);
[this,sArgs] = ita_parse_arguments(sArgs,varargin);
sArgs = struct('pos1_data','itaHRTF', 'earSide', 'L', 'freq' , 1000,'type','directivity','log',0);
[this,sArgs,unused] = ita_parse_arguments(sArgs,varargin);
idxF = this.freq2index(sArgs.freq);
position = get(0,'ScreenSize');
figure('Position',[10 50 position(3:4)*0.85]);
if ~(~isempty(unused) && find(strcmp(unused,'parent')))
position = get(0,'ScreenSize');
figure('units','normalized','outerposition',[0 0 1 1]);
end
if sArgs.log
freqData_dB = this.getEar(sArgs.earSide).freqData_dB;
else
......@@ -776,14 +846,14 @@ classdef itaHRTF < itaAudio
end
switch sArgs.type
case 'directivity'
surf(this.dirCoord,freqData_dB(idxF,:));
surf(this.dirCoord,freqData_dB(idxF,:),unused{:});
c = colorbar; ylabel(c,'Magnitude in dB')
case 'sphere'
surf(this.dirCoord,this.dirCoord.r,freqData_dB(idxF,:));
surf(this.dirCoord,this.dirCoord.r,freqData_dB(idxF,:),unused{:});
c = colorbar;ylabel(c,'Magnitude in dB')
case 'phase'
phase = unwrap(angle(this.getEar(sArgs.earSide).freqData(idxF,:)));
surf(this.dirCoord,freqData_dB(idxF,:),phase);
surf(this.dirCoord,freqData_dB(idxF,:),phase,unused{:});
c = colorbar;ylabel(c,'Phase in rad')
end
title([sArgs.earSide ', f = ' num2str(round(this.freqVector(idxF)/100)/10) ' kHz'])
......@@ -869,10 +939,12 @@ classdef itaHRTF < itaAudio
ITD = phasenDiff./(2*pi*repmat(thisC.freqVector,1,size(phase1,2)));