interp.m 7.79 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
function cThis = interp(this,varargin)
% 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)
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
22 23 24 25 26
%     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])
27 28 29 30 31 32 33
%
% OUTPUT:
%     itaHRTF object
%     .freqData: interpolated / range-extrapolated HRTFs for defined field points
%     .timeData: interpolated / range-extrapolated HRIRs for defined field points
%     .dirCoord: itaCoordinates object
%
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
34
% 
35 36 37 38 39 40
%
% [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)
%
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
41 42 43 44 45 46 47
% [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)
%
%
%
48 49 50
% Author:  Florian Pausch <fpa@akustik.rwth-aachen.de>
% Version: 2016-02-05

Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
51
sArgs           = struct('order',50,'eps',1e-5,'shiftToEar',false,'shiftAxis','y','shiftOffset',[-0.0725 0.0725]);
52
sArgs           = ita_parse_arguments(sArgs,varargin,2);
53 54
if isempty(varargin) || ~isa(varargin{1},'itaCoordinates')
    error('itaHRTF:interp', ' An itaCoordinate object is needed!')
55 56 57 58 59 60
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;
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
61
temp_r          = field_in.r(1);
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
field           = itaCoordinates(size(tempfield,1));
field.r         = repmat(temp_r,size(tempfield,1),1);
field.phi_deg   = tempfield(:,1);
field.theta_deg = tempfield(:,2);

N               = sArgs.order;
epsilon         = sArgs.eps;                                   % regularization parameter
k               = this.wavenumber;                             % wave number
k(1)            = eps;
% add eps to avoid NaN's

Nmax            = floor(sqrt(this.nDirections/4)-1);

% construct vector of length (N+1)^2 regularization weights and,
% if needed, spherical hankel functions of second kind (for r0 and r1)
if ~isequal(this.dirCoord.r(1),field.r(1))
    kr0 = k*this.dirCoord.r(1);                         % measurement radius
    kr1 = k*field.r(1);                                 % extrapolation radius
    
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
81 82
    hankel_r0 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr0);
    hankel_r1 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr1);
83 84 85 86 87
    hankel_div = hankel_r1 ./ hankel_r0;
    
    hankel_rep = hankel_div(:,1);
end

Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
88 89 90
nSH = (Nmax+1).^2;
I = sparse(eye(nSH));
n = ita_sph_linear2degreeorder(1:round(nSH)).';
91

Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
92 93 94 95 96 97 98
%% 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;
99 100 101 102 103 104 105 106
    end
end

%% Weights
[~,w]= this.dirCoord.spherical_voronoi;         % calculate weighting coefficients (Voronoi surfaces <-> measurement points)


W = sparse(diag(w));                                      % diagonal matrix containing weights
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
107 108
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
109 110 111 112 113 114 115 116 117 118

%% Calculate HRTF data for field points
if Nmax > 25
    ita_disp('[itaHRTF.interp] Be patient...')
end
    
% init.
hrtf_arbi = zeros(this.nBins,2*field.nPoints); % columns: LRLRLR...
for ear=1:2
    
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
119 120 121
    % sh transformation
    freqData_temp   = copyData.freqData(:,ear:2:end);
    a0              = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
122

Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
123 124 125 126 127 128 129
%     %% 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,:))
130 131
    
    
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
132
    % range extrapolation
133 134 135 136
    if ~isequal(this.dirCoord.r(1),field.r(1))
        % calculate range-extrapolated HRTFs
        a1 = a0 .* hankel_rep.';
        
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
137 138 139 140 141
        %%% 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
142 143
        hrtf_arbi(:,ear:2:end) = (Yest*a1).';           % interpolated + range-extrapolated HRTFs
    else
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
144 145
        % reconstruction to spatial data
        Yest = ita_sph_base(field,Nmax,'orthonormal',true);  % use real-valued SH's
146 147 148 149 150 151 152 153 154 155 156
        hrtf_arbi(:,ear:2:end) = (Yest*a0).';           % interpolated HRTFs
    end
end


% set new direction coordinates
sph                         = zeros(field.nPoints*2 ,3);
sph(1:2:end,:)              = field.sph;
sph(2:2:end,:)              = field.sph;

% write new HRTF data set
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
157 158 159
cThis                       = this;
cThis.freqData = hrtf_arbi;
cThis.channelCoordinates.sph= sph;
160

Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
161 162 163 164 165 166 167 168 169 170

%% 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
171 172 173 174 175 176 177 178 179 180 181 182

if ~isequal(cThis.dirCoord.r(1),field.r(1))%???
    cThis.dirCoord.r = field.r;
end

if N > 25
    ita_disp('[itaHRTF.interp] ...calculation finished!')
end
end



Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
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
199
    
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213
    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
214
    
Jan-Gerrit Richter's avatar
Jan-Gerrit Richter committed
215 216 217 218 219 220
    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;
221 222

end