function [dtf,ctf]=SOFAhrtf2dtf(hrtf,varargin)
%SOFAHRTF2DTF converts HRTFs to DTFs (and CTFs)
%   Usage:      [dtf,ctf]=SOFAhrtf2dtf(hrtf)
%               [dtf,ctf]=SOFAhrtf2dtf(hrtf,f1,f2)
%
%   Input parameters:
%     hrtf:     SOFA object with SimpleFreeFieldHRIR convention
%
%   Output arguments:
%     dtf:      SOFA object with the directional transfer functions
%     ctf:      SOFA object with the common transfer functions
%
%   `SOFAhrtf2dtf(...)` calculates directional transfer functions (DTFs) using 
%   the method from Majdak et al. (2010), which is similar to the procedure 
%   of Middlebrooks (1999a). 
%   The magnitude of the common transfer function (CTF) is calculated by   
%   averaging the log-magnitude spectra across all HRTFs for each ear.
%   The phase spectrum of the CTF is the minimum phase 
%   corresponding to the amplitude spectrum of the CTF. 
%   A DTF results from filtering the HRTF with the inverse complex CTF. 
%
%   `SOFAhrtf2dtf` accepts the following optional parameters:
%
%     'f1',f1     start frequency of the filtering (in Hz; default: 50 Hz)
%     'f2',f2     end frequency of the filtering (in Hz; default: 18 kHz)
%     'atten',a   broadband attenuation in order to avoid clipping (in dB; default: 20 dB)
%
%   References: majdak2010methods middlebrooks1999scaling

% Author: Robert Baumgartner, 16.01.2014

%% Check Input
definput.keyvals.f1 = 50;     % Hz
definput.keyvals.f2 = 18000;  % Hz
definput.keyvals.atten = 20;  % dB

[flags,kv]=SOFAarghelper({'f1','f2','atten'},definput,varargin);

%% Settings
kv.fs = hrtf.Data.SamplingRate;
hrtfmtx = shiftdim(hrtf.Data.IR,2); % dim 1: time, dim 2: source position, dim 3: receiver/ear
N = size(hrtfmtx,1);
Nfft = N; %2^nextpow2(N);

%% Frequency bounds
df = kv.fs/Nfft;
f = 0:df:kv.fs-df;
idx = f >= kv.f1 & f <= kv.f2;
idx(Nfft/2+2:end) = fliplr(idx(2:Nfft/2));

%% CTF calculation
hrtff=fft(hrtfmtx,Nfft);
ctfflog=mean(log(abs(hrtff)+eps),2);

% Force minimum phase 
ctfcep = ifft(ctfflog,Nfft);
ctfcep(Nfft/2+2:Nfft,:,:) = 0;    % flip acausal part to causal part or simply multiply
ctfcep(2:Nfft/2,:,:) = 2*ctfcep(2:Nfft/2,:,:);    % causal part by 2 (due to symmetry)
ctfflog = fft(ctfcep,Nfft);
ctff = exp(ctfflog);
ctfmtx = ifft(ctff,Nfft);

%% DTF calculation
dtff = hrtff;
dtff(idx,:,:) = hrtff(idx,:,:)./repmat(ctff(idx,:,:),[1 size(hrtff,2) 1]);
dtfmtx = ifft(dtff,Nfft);

%% Attenuate to avoid clipping
ctfmtx = ctfmtx / 10^(kv.atten/20);
dtfmtx = dtfmtx / 10^(kv.atten/20);

%% Output Objects
dtf = hrtf;
dtf.Data.IR = shiftdim(dtfmtx,1);
dtf.GLOBAL_Comment = [dtf.GLOBAL_Comment '. Directional transfer functions (DTFs) were generated by removing from the HRTFs the direction-independent log-amplitude spectrum for each ear.'];

ctf = hrtf;
ctf.Data.IR = shiftdim(ctfmtx,1);
ctf.GLOBAL_Comment = [dtf.GLOBAL_Comment '. Common transfer functions (CTFs) were extracted from the HRTFs in terms of averaging the log-amplitude spectrum for each ear. The phase of a CTF is determined as minimum phase.'];
ctf.SourcePosition = [0 0 0];
if isfield(ctf,'MeasurementSourceAudioChannel')
  ctf.MeasurementSourceAudioChannel = 0;
end
if isfield(ctf,'MeasurementAudioLatency')
  ctf.MeasurementAudioLatency = zeros(1,size(ctf.MeasurementAudioLatency,2));
end
ctf = SOFAupdateDimensions(ctf);

end
