Commit ef7486f7 authored by Hark Braren's avatar Hark Braren
Browse files

new option in ITD phase_delay method and minor cleanup

parent 65de01bd
......@@ -198,6 +198,7 @@ classdef itaHRTF < itaAudio
end
end
end
%% ..............GET...............................................
function nDirections = get.nDirections(this)
[~,idxDim] = uniquetol([this.channelCoordinates.phi_deg this.channelCoordinates.theta_deg] ,'ByRows',true);
......@@ -602,9 +603,10 @@ classdef itaHRTF < itaAudio
this.mTF_type = TF_types{strcmpi(type, TF_types)};
end
end
end
methods %Class related functions
methods %Class related functions
%% Direction related functions
function HRTFout = findnearestHRTF(this,varargin)
if nargin ==2
......@@ -882,7 +884,7 @@ classdef itaHRTF < itaAudio
% end line
end
%% INTERAURAL PARAMETERS
%% Interaural Parameters
function varargout = ILD(this, varargin)
%Calculates the ILD either for each distinct frequency or
%by integrating between two frequencies.
......@@ -938,16 +940,40 @@ classdef itaHRTF < itaAudio
end
function varargout = ITD(this, varargin)
% -----------------------------------------------------------------
% See methods and options below
% -----------------------------------------------------------------
% Input
sArgs = struct('method', 'phase_delay', 'filter' , [200 2000] ,...
'thresh','10dB','energy',true,'centroid',false,'reshape',true);
% Calculate ITD of HRTF
% [ITD, azimuthAngle] = itaHRTF.ITD(options)
%
% Options (default):
% 'method' ('phase_delay'): phase_delay, xcorr, threshold
% 'filter' ([200 2000]) : [freqRange] or 'none' -> frequency
% dependened ITD with phase_delay
% 'thresh' ('10dB') : only with threshold method
% 'energy' (true) : xcorr on energtic (^2) not linear
% timeData
% 'centroid' (false) : only with xcorr
% 'reshape' (true) : Reshape the ITD in a matrix where
% the column defines the phi-
% direction and the row the theta-
% direction
% 'compensateMinPhase' (false) : only with phase_delay:
% only analyze linear phase comp
% -> phaseITD = phase-minPhase
%
% see code below for further details
% Parse Options
sArgs = struct('method', 'phase_delay',...
'filter' , [200 2000] ,...
'thresh', '10dB',...
'energy', true,...
'centroid', false,...
'reshape', true,...
'compensateMinPhase',false);
sArgs = ita_parse_arguments(sArgs,varargin);
if numel(this.theta_Unique)>1
ita_verbose_info(' More than one elevation in this object!', 0);
ita_verbose_info(['More than one elevation in this object!',...
' - behavior depends on chosen method'], 0);
%this = this.sphericalSlice('theta_deg',90);
end
......@@ -958,15 +984,14 @@ classdef itaHRTF < itaAudio
% study of interaural time delay estimation methods. In: The
% Journal of the Acoustical Society of America 135 (6), S.
% 3530-3540.
switch sArgs.method
case 'phase_delay'
% .....................................................
% options: filter
% options: filter, compensateMinPhase
% .....................................................
[~,tau] = ita_time_shift(this,'0dB');
[~,idxMin] = max(tau); % shift of trackLength/3 seems to be good for plotting - No idea
thisC = ita_time_shift(this,tau(idxMin)-this.trackLength/3,'time');
[~,tau] = ita_time_shift(this,'10dB');
[~,idxMin] = max(tau); %all tau are negative -> min(abs)== max
thisC = ita_time_shift(this,tau(idxMin),'time');
if ischar(sArgs.filter) % frequency dependent
pLeft = thisC.freqData(:,thisC.EarSide == 'L');
......@@ -974,8 +999,18 @@ classdef itaHRTF < itaAudio
phaseLeft = unwrap(angle(pLeft));
phaseRight = unwrap(angle(pRight));
phasenDiff = phaseLeft - phaseRight;
if sArgs.compensateMinPhase
thisCMin = ita_minimumphase(thisC);
pLeftMin = thisCMin.freqData(:,thisCMin.EarSide == 'L');
pRightMin = thisCMin.freqData(:,thisCMin.EarSide == 'R');
phaseLeftMin = unwrap(angle(pLeftMin));
phaseRightMin = unwrap(angle(pRightMin));
phaseLeft = phaseLeft-phaseLeftMin;
phaseRight = phaseRight-phaseRightMin;
end
phasenDiff = phaseLeft - phaseRight;
ITD = phasenDiff./(2*pi*repmat(thisC.freqVector,1,size(phaseLeft,2)));
else % averaged
......@@ -986,8 +1021,11 @@ classdef itaHRTF < itaAudio
ita_verbose_info('itaHRTF: Excluding invalid bin (f=0Hz) for ITD calculation (phase_delay method)', 2)
end
usedBins( usedBins <= 1 ) = [];
phase = unwrap(angle(thisC.freqData(usedBins,:)));
if sArgs.compensateMinPhase
minPhase = unwrap(angle(ita_minimumphase(thisC).freqData(usedBins,:)));
phase = unwrap(angle(thisC.freqData(usedBins,:)));
phase = phase-minPhase;
end
freqVector = thisC.freqVector(usedBins);
phaseDelay = bsxfun(@rdivide, phase,2*pi*freqVector);
......@@ -1294,7 +1332,8 @@ classdef itaHRTF < itaAudio
% calculate ITD
if ~isempty(sArgs.theta_deg)
thisS = this.sphericalSlice('theta_deg',sArgs.theta_deg);
else thisS = this;
else
thisS = this;
end
thetaC_deg = rad2deg(thisS.theta_Unique);
......@@ -1305,10 +1344,12 @@ classdef itaHRTF < itaAudio
[~, idxC] = sort(coord,2);
[~, idxCT] = uniquetol(thisS.dirCoord.theta_deg,eps);
ITD = thisS.ITD('method',...
sArgs.method, 'filter' , sArgs.filter , 'thresh',sArgs.thresh,...
'energy',sArgs.energy,'centroid',sArgs.centroid,'reshape',true);
ITD = thisS.ITD('method', sArgs.method,...
'filter', sArgs.filter ,...
'thresh', sArgs.thresh,...
'energy',sArgs.energy,...
'centroid',sArgs.centroid,...
'reshape',true);
ITD_S = ITD;
for idxT = 1:nTheta
ITD_S(idxT,:) = ITD(idxT,idxC(idxT,:));
......@@ -1325,9 +1366,10 @@ classdef itaHRTF < itaAudio
axes(sArgs.axes_handle);
hold on;
end
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)), '^\circ'))
shading flat
colorbar
......@@ -1342,10 +1384,12 @@ classdef itaHRTF < itaAudio
cMax = max(max(ITD(2:end,:)));
cMin = abs(min(min(ITD(2:end,:))));
if cMax>cMin,caxis([-cMax cMax]);
else caxis([-cMin cMin]);
if cMax>cMin
caxis([-cMax cMax]);
else
caxis([-cMin cMin]);
end
elseif strcmp(sArgs.plot_type,'color') && numel(sArgs.theta_deg)~= 1
elseif strcmp(sArgs.plot_type,'color') && numel(thetaC_deg)~= 1
% angle dependent ITD (theta & phi)
pcolor(thetaC_deg, phiC_deg,ITD_SS'*1000)
shading flat
......@@ -1358,7 +1402,7 @@ classdef itaHRTF < itaAudio
ylabel('Azimuth Angle in Degree');
set(gca,'xTick',0:15:360,'yTick',0:30:360)
title('ITD in Milliseconds')
elseif strcmp(sArgs.plot_type,'line') || numel(sArgs.theta_deg)== 1
elseif strcmp(sArgs.plot_type,'line') || numel(thetaC_deg)== 1
% angle dependent ITD (phi)
plot(phiC_deg,ITD_SS*1000)
yMax = max(abs(ITD_SS(:)));
......@@ -1553,7 +1597,7 @@ classdef itaHRTF < itaAudio
set(gca,'xTick',xticks,'xticklabel',xlabels)
set(gca,'yTick',yticks,'xticklabel',yticks)
caxis([cMin cMax]);
caxis([cMin cMax]);
title(strTitle)
shading interp
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment