function [hrtf_out, hrtfArray] = extrapolate_hrtf_set(hrtf_name, Nout)

hrtf_sofa = SOFAload(strcat(hrtf_name,'.sofa'));
name = strcat(hrtf_sofa.GLOBAL_DatabaseName,'_',hrtf_sofa.GLOBAL_ListenerShortName);
hrtfMeasurementArray = spherical_design('Direct',hrtf_sofa.SourcePosition,name);

rls = mean(hrtf_sofa.SourcePosition(:,3));
if length(Nout) == 1
    hrtfArray = spherical_design('Design',Nout,'T',rls);
else
    rls = mean(hrtf_sofa.SourcePosition(:,3));
    hrtfArray = spherical_design('Direct',Nout*rls,name);
end
HRTF_spectrum = fft( hrtf_sofa.Data.IR,[],3 );
HRTF_amp = abs(HRTF_spectrum);
HRTF_phase = unwrap(angle(HRTF_spectrum),[],3);
hrtf_out = zeros(size(hrtfArray.x0,1),2,size(hrtf_sofa.Data.IR,3));
for n = 1 : size(hrtfArray.x0,1)
    ampOut = hrtfMeasurementArray.interpolate(hrtfArray.x0(n,:)', HRTF_amp,'linear');
    phaseOut = hrtfMeasurementArray.interpolate(hrtfArray.x0(n,:)', HRTF_phase,'linear');
    hrtf_out(n,:,:) = ifft(ampOut.*exp(1i*phaseOut),[],1,'symmetric').';
end

end

