classdef itaHRTF < itaAudio
%ITAHRTF - class to deal with HRTFs
%
% Examples:
% hrtf = itaHRTF('sofa','TU-Berlin_QU_KEMAR_anechoic_radius_1m.sofa')
%
% These objects can be used like itaAudios and helps to find HRTF angles
% quickly. In addition different methods are implemented to evaluate
% binaural parameters and interpolate the data set.
%
% itaHRTF Properties:
% dirCoord Measured directions
% EarSide Ear side ('L' left or 'R' right) of each channel
% TF_type [HRTF DTF Recording]
% sphereType [ring cap sphere undefined]
%
% resAzimuth resolution in azimuth (only equiangular)
% resElevation resolution in elevation (only equiangular)
%
% rangeAzimuth min. and max. angle in azimuth
% rangeElevation min. and max. angle in elevation
%
% nPointsAzimuth number of directions in azimuth
% nPointsElevation number of directions in elevation
%
% nPoints total number of directions
%
% itaHRTF Methods (find & select directions):
% HRTFfind = findnearestHRTF(varargin)
% HRTFdir = direction(idxCoord)
% thetaUni = theta_Unique
% phiUni = phi_Unique
% slice = sphericalSlice(dirID,dir_deg)
% HRTF_left = getEar(earSide)
%
% itaHRTF Methods (play):
% play_gui(stimulus)
%
% itaHRTF Methods (store):
% audioHRTF = itaHRTF2itaAudio
% writeDAFFFile(filePath)
%
% itaHRTF Methods (binaural parameter):
% ITD = ITD(varargin)
% t0 = meanTimeDelay(varargin)
% ILD = ILD(varargin)
%
% itaHRTF Methods (manipulation):
% DTF = calcDTF
% HRTF_int = interp(varargin)
%
% itaHRTF Methods (plot):
% plot_ITD(varargin)
% plot_freqSlice(varargin)
%
% See also:
% itaAudio, test_rbo_postprocessing_HRTF_arc_CropDiv
%
% Reference page in Help browser
% doc itaHRTF
%
% This file is part of the application HRTF_class for the ITA-Toolbox. All rights reserved.
% You can find the license for this m-file in the application folder.
%
% Author: Ramona Bomhardt -- Email: rbo@akustik.rwth-aachen.de
% Created: 10-Jul-2014
properties (Access = private)
mCoordSave = [];
mChNames = [];
mDirCoord = itaCoordinates;
mEarSide = [];
mTF_type = 'HRTF';
mSphereType = 'undefined';
end
properties (Dependent = true, Hidden = false)
dirCoord = itaCoordinates;
EarSide = [];
TF_type = 'HRTF';
sphereType = 'undefined';
resAzimuth = 5;
resElevation = 5;
rangeAzimuth = [0 359];
rangeElevation = [0 180];
nPointsAzimuth = 72;
nPointsElevation= 37;
nPoints = [];
phi_Offset = zeros(37,1);
end
properties (Dependent = true, Hidden = true)
end
properties (Dependent = true, SetAccess = private)
openDaff2itaHRTF;
itaAudio2itaHRTF;
init;
hdf2itaHRTF;
sofa2itaHRTF;
nDirections = [];
end
methods % Special functions that implement operations that are usually performed only on instances of the class
%% Input
function this = itaHRTF(varargin)
this = this@itaAudio();
if nargin >1
% itaAudio input
TF_types = this.propertiesTF_type;
for iTF = 1:numel(TF_types)
if ~isempty(find(strcmpi(varargin, TF_types{iTF})==1, 1))
this.itaAudio2itaHRTF = varargin{find(strcmpi(varargin, TF_types{iTF})==1)-1};
this.TF_type = TF_types(iTF);
end
end
% init
if nargin == 4
this.init = varargin;
end
% openDaff input
if ~isempty(find(strcmpi(varargin,'Daff')==1, 1))
this.openDaff2itaHRTF = varargin{find(strcmpi(varargin,'Daff')==1)+1};
end
% hdf5 input
if ~isempty(find(strcmpi(varargin,'hdf5')==1, 1))
this.hdf2itaHRTF = varargin{find(strcmpi(varargin,'hdf5')==1)+1};
end
% sofa input
if ~isempty(find(strcmpi(varargin,'SOFA')==1, 1))
this.sofa2itaHRTF = varargin{find(strcmpi(varargin,'SOFA')==1)+1};
end
elseif nargin == 1
if isa(varargin{1},'itaHRTF')
this = varargin{1};
elseif nargin ==1 && isstruct(varargin{1}) % only for loading
obj = varargin{1};
this.data = obj.data;
this.signalType = 'energy';
% additional itaHRTF data
if datenum(2014,7,5)1,
if numel(HRTF)>1000 % takes a while
ita_verbose_info(' A lot of data ...please wait... don''t use itaAudio multi instances for the next time!', 0);
end
coordinates = HRTF(1).channelCoordinates;
if (coordinates.nPoints == 2) & (sum(isnan(coordinates.sph)) < numel(coordinates.sph))
ita_verbose_info('Found NaNs in the coordinates. I will copy existing coordinates');
for index = 1:length(HRTF)
coordinates = HRTF(index).channelCoordinates;
coordinates.sph = repmat(coordinates.sph(1,:),2,1);
HRTF(index).channelCoordinates = coordinates;
end
end
HRTFc = HRTF.merge;
else HRTFc = HRTF;
end
% coordinates available?
if isnan(HRTFc.channelCoordinates.cart)
error('itaHRTF:Def', ' No channelCoordinates available')
end
coord = HRTFc.channelCoordinates;
% find the corresponding left and right channel
pairs = zeros(coord.nPoints/2,2);
if coord.nPoints>10000 % takes a while
ita_verbose_info([num2str(coord.nPoints) ' Points has to be sorted ...please wait...'], 0);
end
counter = 1;
thetaPhi = round([coord.theta_deg coord.phi_deg]*10)/10;
deletedChannel = 0;
for i1 = 1:coord.nPoints
coordCurrent = thetaPhi(i1,:);
if isempty(find(pairs(:) == i1, 1)) % only if the corresponding channel is not found
% find corresponding channel
coordComp = thetaPhi([1:i1-1 i1+1:coord.nPoints],:);
diffCoord = bsxfun(@minus,coordCurrent,coordComp)== zeros(size(coordComp));
idxCoord = find(diffCoord(:,1).*diffCoord(:,2) ==1);
if length(idxCoord) > 1
% deletedChannel = deletedChannel + length(idxCoord) -1;
idxCoord = idxCoord(1);
end
% store the corresponding channel
pairs(counter,1) = i1;
if idxCoord 1,
thetaC = thetaC';
end
elseif numel(thetaC)==1 && numel(phiC)~=1,
thetaC = ones(numel(phiC),1)*thetaC;
if size(phiC,2)>1,
phiC = phiC';
end
end
coordC = itaCoordinates([r thetaC phiC],'sph');
end
idxCoord = this.dirCoord.findnearest(coordC);
[~, I] = unique(idxCoord);
idxCoordUnique = idxCoord(I);
% idxCoordUnique = unique(idxCoord,'stable');
if numel(idxCoord)~= numel(idxCoordUnique)
ita_verbose_info('Multiple coordinates are neglected!', 0);
end
if sum(this.EarSide == 'R') ~= sum(this.EarSide == 'L') % only one ear is available
ita_verbose_info('You use only one Ear! Conversion to itaAudio.', 0);
idxCoord = this.channelCoordinates.findnearest(coordC);
[~, I] = unique(idxCoord);
idxCoordUnique = idxCoord(I);
HRTFout = this.ch(idxCoordUnique).itaHRTF2itaAudio;
else
HRTFout = this.direction(idxCoordUnique);
end
%HRTFout = this.direction(idxCoord);
end
function obj = direction(this, idxCoord)
idxDir = zeros(numel(idxCoord)*2,1);
idxDir(1:2:numel(idxCoord)*2,:) = 2*idxCoord-1;
idxDir(idxDir==0)=1;
idxDir(2:2:numel(idxCoord)*2) = idxDir(1:2:numel(idxCoord)*2,:)+1;
hrtfTMP = this.ch(idxDir);
hrtfTMP.channelCoordinates = this.channelCoordinates.n(idxDir);
hrtfTMP.EarSide = this.EarSide(idxDir);
obj = itaHRTF(hrtfTMP);
end
function thetaUni = theta_Unique(this,varargin)
thetaUni = unique(this.dirCoord.theta);
if nargin == 2
thetaUni = unique(this.dirCoord.theta,'stable');
end
end
function phiUni = phi_Unique(this,varargin)
phiUni = unique(this.dirCoord.phi);
if nargin == 2
phiUni = unique(this.dirCoord.phi,'stable');
end
end
function thetaUni = theta_UniqueDeg(this,varargin)
thetaUni = rad2deg(theta_Unique(this,varargin));
end
function phiUni = phi_UniqueDeg(this,varargin)
phiUni = rad2deg(phi_Unique(this,varargin));
end
function slice = sphericalSlice(this,dirID,dir_deg)
% dir in degree
% dirID [phi, theta]
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
end
function slice = ss(this,dirID,dir_deg)
slice = this.sphericalSlice(dirID,dir_deg);
end
function HRTFout = getEar(this,earSide)
switch earSide
case 'L',
HRTFout = this.ch(this.EarSide == 'L');
HRTFout.mEarSide = repmat('L',HRTFout.nChannels,1);
case 'R',
HRTFout = this.ch(this.EarSide == 'R');
HRTFout.mEarSide = repmat('R',HRTFout.nChannels,1);
end
end
%% ITA Toolbox Functions
function stimuli = conv(this,stimulus)
if isa(stimulus, 'itaAudio')
stimuli = itaAudio(this.nDirections,1);
idxCh = 1:2:this.dimensions;
for idxDir = 1:this.nDirections
stimuli(idxDir) = ita_convolve(stimulus,this.ch([idxCh(idxDir) idxCh(idxDir)+1]));
end
end
end
function play_gui(this,stimulus)
if isa(stimulus, 'itaAudio')
% check size of input data
if this.nDirections>75,
thisTmp = this.direction(1:75);
ita_verbose_info(' A lot of data ... you cannot show everything in the GUI!', 0);
else thisTmp = this;
end
% convolve
stimuli = thisTmp.conv(stimulus);
% normalize level
stimuliAll = stimuli.merge;
maxLevel = max(abs(stimuliAll.timeData(:)))*1.05;
stimuliNorm = stimuli;
for idxDir = 1:thisTmp.nDirections
stimuliNorm(idxDir) = stimuli(idxDir)/maxLevel;
end
% play gui
ita_play_gui(stimuliNorm, thisTmp.channelNames(1:2:thisTmp.dimensions));
%ita_play_gui(stimuliNorm, ita_sprintf('phi= %2.0f� theta= %2.0f�',...
% thisTmp.dirCoord.phi_deg,thisTmp.dirCoord.theta_deg));
end
end
function audioHRTF = itaHRTF2itaAudio(this)
audioHRTF = itaAudio;
audioHRTF.samplingRate = this.samplingRate;
audioHRTF.timeData = this.timeData;
audioHRTF.channelNames = ita_sprintf('%s ( %2.0f, %2.0f)',...
this.mEarSide , this.channelCoordinates.theta_deg,this.channelCoordinates.phi_deg );
audioHRTF.channelCoordinates = this.channelCoordinates;
audioHRTF.signalType = 'energy';
end
function surf(varargin)
sArgs = struct('pos1_data','itaHRTF', 'earSide', 'L', 'freq' , 5000,'type','directivity');
[this,sArgs] = ita_parse_arguments(sArgs,varargin);
idxF = this.freq2index(sArgs.freq);
position = get(0,'ScreenSize');
figure('Position',[10 50 position(3:4)*0.85]);
freqData_dB = this.getEar(sArgs.earSide).freqData_dB;
switch sArgs.type
case 'directivity'
surf(this.dirCoord,freqData_dB(idxF,:));
c = colorbar; ylabel(c,'Magnitude in dB')
case 'sphere'
surf(this.dirCoord,this.dirCoord.r,freqData_dB(idxF,:));
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);
c = colorbar;ylabel(c,'Phase in rad')
end
title([sArgs.earSide ', f = ' num2str(round(this.freqVector(idxF)/100)/10) ' kHz'])
end
function display(this)
this.displayLineStart
this.disp
dir = num2str(this.nDirections,5);
stringD = [dir ' Directions (Type = ' this.mTF_type ')'];
middleLine = this.LINE_MIDDLE;
middleLine(3:(2+length(stringD))) = stringD;
fprintf([middleLine '\n']);
end
function disp(this)
disp@itaAudio(this)
sphType = [this.sphereType repmat(' ',1,9-length(this.sphereType))];
string = [' Sphere Type = ' sphType ];
% this block adds the class name
classnamestring = ['^--|' mfilename('class') '|'];
fullline = repmat(' ',1,this.LINE_LENGTH);
fullline(1:numel(string)) = string;
startvalue = length(classnamestring);
fullline(length(fullline)-startvalue+1:end) = classnamestring;
disp(fullline);
% end line
end
%% Ramonas' Functions
function varargout = ITD(varargin)
% -----------------------------------------------------------------
% See methods and options below
% -----------------------------------------------------------------
% Input
sArgs = struct('pos1_data','itaHRTF', 'method', 'phase_delay', 'filter' , [200 2000] ,...
'thresh','10dB','energy',true,'centroid',false,'reshape',true);
[this,sArgs] = ita_parse_arguments(sArgs,varargin);
if numel(this.theta_Unique)>1
ita_verbose_info(' More than one elevation in this object!', 0);
%this = this.sphericalSlice('theta_deg',90);
end
% -------------------------------------------------------------
% methods: phase_delay, xcorr, threshold
% -------------------------------------------------------------
% Katz, Brian F. G.; Noisternig, Markus (2014): A comparative
% 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
% .....................................................
[~,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');
if ischar(sArgs.filter) % frequency dependent
p1 = thisC.freqData(:,1:2:thisC.dimensions);
p2 = thisC.freqData(:,2:2:thisC.dimensions);
phase1 = unwrap(angle(p1));
phase2 = unwrap(angle(p2));
phasenDiff = phase1 - phase2;
ITD = phasenDiff./(2*pi*repmat(thisC.freqVector,1,size(phase1,2)));
else % averaged
phase = unwrap(angle(thisC.freqData));
t0_freq = bsxfun(@rdivide, phase,2*pi*thisC.freqVector);
t0_freq = t0_freq(~isnan(t0_freq(:,1)),:);
t0_mean = mean(t0_freq(unique(thisC.freq2index(sArgs.filter(1)):thisC.freq2index(sArgs.filter(2))),:)); %mean is smoother than max; lower freq smooths also the result
ITD = t0_mean(thisC.EarSide == 'L') - t0_mean(thisC.EarSide == 'R');
end
case 'xcorr'
% .....................................................
% options: energy, filter, centroid
% .....................................................
if ischar(sArgs.filter), thisF = this; % FILTER
else thisF = ita_mpb_filter(this,[sArgs.filter(1), sArgs.filter(2)]);
end
% Interpolation for smoother curves
xUpSample = 5;
SR = xUpSample*thisF.samplingRate;
tV_Interp = 0:1/SR:thisF.trackLength;
timeData_Interp = interp1(thisF.timeVector,thisF.timeData,tV_Interp,'spline');
% case: energy
if sArgs.energy ,timeData_Interp = timeData_Interp.^2;
end
idxL = find(thisF.EarSide== 'L'); idxR = find(thisF.EarSide == 'R');
corrIR = zeros(2*numel(tV_Interp)-1,this.nDirections);
for idxDir = 1:thisF.nDirections
corrIR(:,idxDir) = xcorr(timeData_Interp(:,idxL(idxDir)),timeData_Interp(:,idxR(idxDir)));
end
if ~sArgs.centroid % max
[~, idxMax] = max(abs(corrIR));
ITD = (numel(tV_Interp)- idxMax)/SR;
else % centroid
tV = 0:1/SR:(2*numel(tV_Interp)-2)/SR;
C = sum(bsxfun(@times,abs(corrIR),tV'))./sum(abs(corrIR));
ITD = thisF.trackLength-C;
end
case 'threshold'
% .....................................................
% options: filter
% .....................................................
if ischar(sArgs.filter), thisF = this; % FILTER
else thisF = ita_mpb_filter(this,[sArgs.filter(1), sArgs.filter(2)]);
end
[~,tau] = ita_time_shift(thisF,sArgs.thresh);
ITD = tau(thisF.EarSide== 'L')-tau(thisF.EarSide == 'R');
end
% Reshape the ITD in a matrix where the column defines the phi-
% direction and the row the theta-direction
if sArgs.reshape && ~ischar(sArgs.filter)
nPhi = numel(this.phi_Unique);
nTheta = numel(this.theta_Unique);
if nPhi*nTheta == this.nDirections
sITD = reshape(ITD,nTheta,nPhi);
else
ita_verbose_info(' ITD could not be reshape: nPhi*nTheta ~= nDir!', 0);
sITD = ITD;
end
else
sITD = ITD;
end
varargout{1} = sITD;
if nargout == 2, varargout{2} = rad2deg(this.phi_Unique('stable'));
end
end
function t0 = meanTimeDelay(this,varargin)
%-- OLD -------------------------------------------------------
[~,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*0.33,'time');
phase = unwrap(angle(thisC.freqData));
t0_freq = bsxfun(@rdivide, phase,2*pi*thisC.freqVector);
%t0_mean = t0_freq(thisC.freq2index(1000),:);
t0_mean = mean(t0_freq(thisC.freq2index(500):thisC.freq2index(2000),:)); %mean is smoother than max; lower freq smooths also the result
if nargin==2
if strcmpi(varargin{1},'L')
t0 = t0_mean(thisC.EarSide == 'L');
elseif strcmpi(varargin{1},'R')
t0 = t0_mean(thisC.EarSide == 'R');
end
else t0 = t0_mean;
end
end
function varargout = calcDTF(this)
if ~strcmpi(this.TF_type,'DTF')
[DTF,comm] = test_rbo_DTF_itaHRTF(this);
varargout{1} =DTF;
if nargout ==2,varargout{2} = comm;end
end
end
% function this = interp(varargin)
%
% Function to calculate HRTFs for arbitrary field points using a N-th order
% spherical harmonics (SH) interpolation / range extrapolation, as described in [1],
% SH expansion coefficients are calculated by means of a least-squares
% approach with Tikhonov regularization
%
% Function may also be used for spatial smoothing of HRTF using
% the method described in [2]. As field input use the original
% measurement grid and set the desired order of the SH matrix /
% truncation order.
%
% INPUT:
% varargin{1} ... itaCoordinates object (required)
% varargin{1}.phi: desired azimuth angles for HRTF interpolation [0 2*pi)
% varargin{1}.theta: desired zenith angles for HRTF interpolation [0 pi]
% varargin{1}.r: (optional) desired radius used for range extrapolation in [m],
% set to 1 if no range extrapolation is required
% order ... order of spherical harmonics matrix (default: 50)
% epsilon ... regularization coefficient (default: 1e-8)
%
% OUTPUT:
% itaHRTF object
% .freqData: interpolated / range-extrapolated HRTFs for defined field points
% .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)
%
% Author: Florian Pausch
% Version: 2016-02-05
function this = smooth_linphase(this,varargin)
% function this = smooth_linphase(varargin)
%
% Function to smooth HRTFs in the frequency domain based on the method proposed by Rasumov et al. in [3], complex smoothing
% is done via ita_smooth()
%
% Parameters:
% 'f_lin' ... frequency above which the phase is approximated by a linear phase term
% 'smoothtype' ... smoothing method, 'LinTimeSec', 'LinTimeSamp', 'LinFreqHertz', 'LinFreqBins',
% 'LogFreqOctave1' (default), 'LogFreqOctave2' or 'Gammatone'
% 'windowWidth' ... bandwidth of filter (depends on smoothtype - type help ita_smooth), e.g. 1/9 (default) in frequency domain
% 'dataTypes' ... defines on which data type smoothing is applied, 'Real', 'Complex', 'Abs' (default), 'GDelay', 'Abs+GDelay'
% or 'Abs+Phase' (type help ita_smooth)
%
% [2] Rasumow, Eugen et al, "Smoothing individual head-related transfer functions in the frequency and spatial domains"
% The Journal of the Acoustical Society of America, 135, 2012-2025 (2014), DOI:http://dx.doi.org/10.1121/1.4867372
%
% Author: Florian Pausch
% Version: 2015-11-04
sArgs = struct('f_lin',5000,'smoothtype','LogFreqOctave1','windowWidth',1/9,'dataTypes','Abs');
sArgs = ita_parse_arguments(sArgs,varargin,1);
f_lin = sArgs.f_lin; % frequency above which the phase is approximated by a linear phase term (f_lin=5000, default)
% parameters for ita_smooth()
smoothtype = sArgs.smoothtype; % smoothing method, 'LinTimeSec', 'LinTimeSamp', 'LinFreqHertz', 'LinFreqBins',
% 'LogFreqOctave1' (default), 'LogFreqOctave2' or 'Gammatone'
windowWidth = sArgs.windowWidth; % bandwidth of filter (depends on smoothtype - type help ita_smooth), e.g. 1/9 (default) in frequency domain
dataTypes = sArgs.dataTypes; % 'Real', 'Complex', 'Abs' (default), 'GDelay', 'Abs+GDelay' or 'Abs+Phase' (type help ita_smooth)
%% Step I: Estimation of the delay of the HRTF peak and the resulting linear phase
% HRTF_env = ita_envelope(this); % calculate the envelope of the HRIR
tau = ita_start_IR(ita_mpb_filter(this,[200,10000]),'threshold',0,'correlation',true);
tau = tau/this.samplingRate;
linphase = exp( -1i*2*pi .* repmat(this.freqVector(this.freq2index(f_lin)+1:end)',1,this.nChannels).*...
repmat(tau,length(this.freqVector(this.freq2index(f_lin)+1:end)),1) ); % linear phase of evaluated HRTF set
%% Step II: Linearize phase for f >= f_lin
this.freqData = abs(this.freqData) .* [exp( 1i*angle(this.freqData(1:this.freq2index(f_lin),:)) );...
linphase ] ;
%% Step III: Remove delay tau
this = ita_time_shift(this,-tau,'samples');
%% Step IV: Complex smoothing
this_smooth = ita_smooth(this,smoothtype,windowWidth,dataTypes);
this.timeData = this_smooth.timeData;
%% Step V: Reconstruct delay tau
this = ita_time_shift(this,tau,'samples');
end
function thisS = smooth_spatial(this, varargin)
% function this = smooth_spatial(varargin)
%
% Function to smooth HRTFs in the spatial domain as shown in [3]
%
% Parameters
% 'N' ... order of truncated spherical harmonics matrix (default: 4)
% a lower order results in less spatial detail/high-frequency detail
% in smoothed HRTF data set
% 'epsilon' ... regularization coefficient (default: 1e-8)
%
% Required: SphericalHarmonics functions of ITA Toolbox
%
% [3] Romigh, G.D.; Brungart, D.S.; Stern, R.M.; Simpson, B.D., "Efficient Real Spherical Harmonic Representation of Head-Related
% Transfer Functions," in Selected Topics in Signal Processing, IEEE Journal of , vol.9, no.5, pp.921-930, Aug. 2015
% doi: 10.1109/JSTSP.2015.2421876
%
% Author: Florian Pausch
% Version: 2016-02-12
tic;
sArgs = struct('N',4,'epsilon',1e-8,'type','min');
sArgs = ita_parse_arguments(sArgs,varargin);
N = sArgs.N;
epsilon = sArgs.epsilon;
Nmeas = floor(sqrt(this.nDirections/4)-1); % SH order of measurement grid (assuming equiangular grid)
if N>Nmeas
fprintf('[\b[itaHRTF.smooth_spatial] Chosen SH order is too high. Order is set to maximum SH order of measurement grid!]\b\n')
fprintf('[\b[itaHRTF.smooth_spatial] N = Nmeas = %s (assuming equiangular sampling)]\b\n',num2str(Nmeas))
N=Nmeas;
end
%% Weighting + regularization
regweights = ita_sph_degreeorder2linear(0:Nmeas,0); % construct vector of length (Nmeas+1) regularization weights
regweights_rep = zeros(sum(2*(0:Nmeas)'+1),1);
regweights_rep(1) = regweights(1);
cntr = 2;
for n=1:Nmeas % repeat regularization weights to get a (Nmeas+1)^2 x 1 vector (TODO: more elegant solution needed)
nTimes = 2*n+1;
regweights_rep(cntr:cntr+nTimes-1) = regweights(n+1)*ones(nTimes,1);
cntr = cntr + nTimes;
end
[~, vWeights] = this.dirCoord.spherical_voronoi; % calculate weighting coefficients (Voronoi surfaces <-> measurement points)
W = diag(vWeights); % diagonal matrix containing weights
D = diag(regweights_rep); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmeas,'orthonormal',false); % calculate real-valued SHs using the measurement grid (high SH-order)
%% Calculate spatially smoothed HRTF data set
hrtf_smoo_wo_ITD = zeros(this.nBins,2*this.dirCoord.nPoints); % init.: columns: LRLRLR...
for ear=1:2
% decompose logarithmic magnitude spectra of measured HRTF set into SH basis functions, as done in [3]
switch sArgs.type
case 'complex'
freqData_temp = this.freqData(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.'; % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
otherwise
freqData_dB = this.freqData_dB;
freqData_temp = freqData_dB(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.'; % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
end
Yest = Y(:,1:(N+1)^2); % eat first (N+1)^2 SH basis functions
a0_trunc = a0(1:(N+1)^2,:); % reduce number of coefficients
hrtf_smoo_wo_ITD(:,ear:2:end) = (Yest*a0_trunc).'; % spatially smoothed HRTF due to reduction of SH decomposition order
end
% % calculate magnitude spectrum and add original HRIR delays as linear phase component
% linphase = exp( -1i*2*pi * repmat(this.freqVector,1,this.nChannels).*...
% repmat(idxIRs_orig/this.samplingRate,this.nBins,1) );
% thisS = this;
% thisS.freqData = 10.^(hrtf_smoo_wo_ITD/20) .* linphase;
switch sArgs.type
case 'min'
this_minphase = ita_minimumphase(this);
idxIRs_orig = ita_start_IR(ita_mpb_filter(this,[200,2000]),'threshold',0,'correlation',true);
deltaT = idxIRs_orig./this_minphase.samplingRate*1.3;
if min(deltaT) < 0 % no negative shifts
deltaT = deltaT-min(deltaT);
end
thisMin = this; %smoothed HRTF
thisMin.freqData= 10.^(hrtf_smoo_wo_ITD/20);
thisS = test_rbo_FIR_lagrange_delay(deltaT,thisMin);
%thisS = ita_mpb_filter(thisS,[200 20000]);
case 'old'
oldPhase = angle(this.freqData);% rbo test
thisS = itaHRTF(this);
thisS.freqData = 10.^(hrtf_smoo_wo_ITD/20) .* exp(1i.*oldPhase); %rbo test
%thisS = ita_mpb_filter(thisS,[200 20000]);
case 'complex'
thisS = this;
thisS.freqData = hrtf_smoo_wo_ITD; %rbo test
end
t2 = toc;
fprintf(['[itaHRTF.smooth_spatial] Calculation finished after ',num2str(round(t2*100/60)/100),' min\n'])
end
%% Plot
function plot_ITD(varargin)
% init
sArgs = struct('pos1_data','itaHRTF', 'method', 'phase_delay', 'filter' , [200 2000] ,...
'thresh','10dB','energy',true,'centroid',false,'reshape',true,...
'theta_deg',[],'plot_type','color');
[this,sArgs] = ita_parse_arguments(sArgs,varargin);
% calculate ITD
if ~isempty(sArgs.theta_deg)
thisS = this.sphericalSlice('theta_deg',sArgs.theta_deg);
else thisS = this;
end
thetaC_deg = rad2deg(thisS.theta_Unique);
phiC_deg = sort(mod(round(rad2deg(thisS.phi_Unique)),360));
nTheta = numel(thetaC_deg);
nPhi = numel(phiC_deg);
coord = reshape(mod(round(thisS.dirCoord.phi_deg),360),nTheta,nPhi);
[~, idxC] = sort(coord,2);
[~, idxCT] = unique(thisS.dirCoord.theta_deg);
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,:));
end
ITD_SS = ITD_S(idxCT(1:nTheta),:);
%..............................................................
% create figure
position = get(0,'ScreenSize');
figure
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)), '�'))
shading flat
colorbar
ylabel('frequency');
ylim([this.freqVector(1) this.freqVector(end)])
xlabel('azimuth angle');
set(gca, 'YScale', 'log');
[xticks, xlabels] = ita_plottools_ticks('log');
set(gca,'yTick',xticks,'yticklabel',xlabels)
cMax = max(max(ITD(2:end,:)));
cMin = abs(min(min(ITD(2:end,:))));
if cMax>cMin,caxis([-cMax cMax]);
else caxis([-cMin cMin]);
end
elseif strcmp(sArgs.plot_type,'color') && numel(sArgs.theta_deg)~= 1
% angle dependent ITD (theta & phi)
pcolor(thetaC_deg, phiC_deg,ITD_SS'*1000)
shading flat
colorbar
cMax = max(abs(ITD_SS(:)));
caxis([-cMax cMax]*1100);
grid on
set(gca,'layer','top')
xlabel('Zenith Angle in Degree');
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
% angle dependent ITD (phi)
plot(phiC_deg,ITD_SS*1000)
yMax = max(abs(ITD_SS(:)));
ylim([-yMax yMax]*1100);
grid on
set(gca,'layer','top')
xlabel('Azimuth Angle in Degree');
ylabel('ITD in Milliseconds');
set(gca,'xTick',0:30:360)
legend(ita_sprintf('%i�', round(thetaC_deg)))
end
end
function plot_freqSlice(varargin)
% init
sArgs = struct('pos1_data','itaHRTF', 'earSide', 'L','plane','horizontal','axes_handle',gca);
[this,sArgs]= ita_parse_arguments(sArgs,varargin);
ah = sArgs.axes_handle;
phiC_deg = rad2deg(unique(round(this.phi_Unique*100)/100));
thetaC_deg = rad2deg(unique(round(this.theta_Unique*100)/100));
% create slice
if numel(thetaC_deg)>1 && numel( phiC_deg)>1
ita_verbose_info(' More than one elevation in this object!', 0);
if strcmp(sArgs.plane,'horizontal')
thetaC_deg = 90;
thisC = this.sphericalSlice('theta_deg', thetaC_deg);
elseif strcmp(sArgs.plane,'median')
phiC_deg = 0;
thisC = this.sphericalSlice('phi_deg', phiC_deg);
end
else thisC = this;
end
% multi defined coordinates
if numel(phiC_deg)1
[~,idxPhiS] = sort(thisC.dirCoord.phi_deg);
thisCs = thisC.direction(idxPhiS);
yticks = round(min(rad2deg(thisCs.phi_Unique))/10)*10:30:round(max(rad2deg(thisCs.phi_Unique))/10)*10;
else
[~,idxPhiS] = sort(thisC.dirCoord.theta_deg);
thisCs = thisC.direction(idxPhiS);
yticks = round(min(rad2deg(thisCs.theta_Unique))/10)*10:30: round(max(rad2deg(thisCs.theta_Unique))/10)*10;
end
% theta or phi slice
earSidePlot = sArgs.earSide;
if numel(phiC_deg)>1,
xData = phiC_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)) '°'];
strXlabel = '\theta in Degree';
end
% Plot properties
% position = get(0,'ScreenSize');
% figure
% set(gcf,'Position',[10 50 position(3:4)*0.85]);
idxfMax = find(this.freqVector>2e4,1,'first');
if isempty(idxfMax), idxfMax = this.nBins; end
fMax = thisCs.freqVector(idxfMax);
[tick, lab] = ita_plottools_ticks('log');
data_dB= thisCs.freqData_dB;
cMax = max(max(data_dB(2:idxfMax,:)));
cMin = min(min(data_dB(2:idxfMax,:)))*0.5;
pcolor(ah, thisCs.freqVector,xData,data_dB(:,thisCs.EarSide == earSidePlot)');
[xticks, xlabels] = ita_plottools_ticks('log');
set(ah,'xTick',xticks,'xticklabel',xlabels)
set(ah,'yTick',yticks,'xticklabel',yticks)
caxis([cMin cMax]);
set(ah, 'XScale', 'log')
title(strTitle)
shading interp
cb = colorbar;
zlab = get(cb,'ylabel');
set(zlab,'String','Level in [dB]');
set(ah,'xtick',tick,'xticklabel',lab)
xlabel('Frequency in Hertz');xlim([thisCs.freqVector(2) fMax ]);
ylabel(strXlabel);
grid on;set(ah,'layer','top')
end
function plot_timeSlice(varargin)
% init
sArgs = struct('pos1_data','itaHRTF', 'earSide', 'L','plane','horizontal');
[this,sArgs]= ita_parse_arguments(sArgs,varargin);
phiC_deg = rad2deg(unique(round(this.phi_Unique*100)/100));
thetaC_deg = rad2deg(unique(round(this.theta_Unique*100)/100));
% create slice
if numel(thetaC_deg)>1 && numel( phiC_deg)>1
ita_verbose_info(' More than one elevation in this object!', 0);
if strcmp(sArgs.plane,'horizontal')
thetaC_deg = 90;
thisC = this.sphericalSlice('theta_deg', thetaC_deg);
elseif strcmp(sArgs.plane,'median')
phiC_deg = 0;
thisC = this.sphericalSlice('phi_deg', phiC_deg);
end
else thisC = this;
end
% multi defined coordinates
if numel(phiC_deg)1
[~,idxPhiS] = sort(thisC.dirCoord.phi_deg);
thisCs = thisC.direction(idxPhiS);
yticks = round(min(rad2deg(thisCs.phi_Unique))/10)*10:30:round(max(rad2deg(thisCs.phi_Unique))/10)*10;
else
[~,idxPhiS] = sort(thisC.dirCoord.theta_deg);
thisCs = thisC.direction(idxPhiS);
yticks = round(min(rad2deg(thisCs.theta_Unique))/10)*10:30: round(max(rad2deg(thisCs.theta_Unique))/10)*10;
end
% theta or phi slice
earSidePlot = sArgs.earSide;
if numel(phiC_deg)>1,
xData = phiC_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)) '°'];
strXlabel = '\theta in Degree';
end
% Plot properties
position = get(0,'ScreenSize');
figure
set(gcf,'Position',[10 50 position(3:4)*0.85]);
idxfMax = find(this.freqVector>2e4,1,'first');
if isempty(idxfMax), idxfMax = this.nBins; end
fMax = thisCs.freqVector(idxfMax);
[tick, lab] = ita_plottools_ticks('log');
data_dB= thisCs.timeData;
cMax = max(max(data_dB(2:idxfMax,:)));
cMin = min(min(data_dB(2:idxfMax,:)))*0.5;
pcolor(thisCs.timeVector,xData,data_dB(:,thisCs.EarSide == earSidePlot)');
[xticks, xlabels] = ita_plottools_ticks('log');
set(gca,'xTick',xticks,'xticklabel',xlabels)
set(gca,'yTick',yticks,'xticklabel',yticks)
caxis([cMin cMax]);
set(gca, 'XScale', 'log')
title(strTitle)
shading interp
cb = colorbar;
zlab = get(cb,'ylabel');
set(zlab,'String','Level in [dB]');
set(gca,'xtick',tick,'xticklabel',lab)
xlabel('Frequency in Hertz');xlim([thisCs.freqVector(2) fMax ]);
ylabel(strXlabel);
grid on;set(gca,'layer','top')
end
end
methods(Hidden = true)
function sObj = saveobj(this)
% Called whenever an object is saved
% have to get save objects for both base classes
% Both options doesn't work at the moment...
this.mCoordSave = this.dirCoord.sph;
this.mChNames = char(this.channelNames);
sObj = saveobj@itaAudio(this);
% Copy all properties that were defined to be saved
propertylist = itaHRTF.propertiesSaved;
for idx = 1:numel(propertylist)
sObj.(propertylist{idx}) = this.(propertylist{idx});
end
end
end
methods(Static)
function this = loadobj(sObj)
this = itaHRTF(sObj);
end
function result = propertiesEarSide
result = {'L','R'};
end
function result = propertiesSaved
result = {'EarSide','sphereType','TF_type','mCoordSave','mChNames'};
end
function result = oldPropertiesSaved
result = {'EarSite','sphereType','TF_type','mCoordSave','mChNames'};
end
function result = propertiesLoad
result = {'mEarSide','mSphereType','mTF_type','mCoordSave','mChNames'};
end
function result = propertiesTF_type
result = {'HRTF', 'DTF','Recording', 'Common'};
end
function result = propertiesSphereType
result = {'cap', 'ring','full','undefined'};
end
function result = propertiesInit
result = {'channelCoordinates','domain','data'};
end
end
end