From 4e1548d861ffae13bf56c9305f89662b1e17739e Mon Sep 17 00:00:00 2001 From: rbo Date: Fri, 3 Feb 2017 17:04:30 +0100 Subject: [PATCH] some fixes... coming soon --- .../HRTF_class/@itaHRTF/itaHRTF.m | 1467 +++++++++++++++++ 1 file changed, 1467 insertions(+) diff --git a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m index d92b218..e7b06ee 100644 --- a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m +++ b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m @@ -1,3 +1,1470 @@ +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,freqData_dB(idxF,:),'radius',1); + 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) + if numel(this) == 1 + 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']); + else + disp([ num2str(numel(this)) ' itaHRTFs']) + end + 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:numel(idxL) + 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 + + + classdef itaHRTF < itaAudio %ITAHRTF - class to deal with HRTFs -- GitLab