[h,h_dist]=sweepIR(rec,Nimp,T,f1,f2,offset,fs)

function [h,h_dist]=sweepIR(rec,Nimp,T,f1,f2,offset,fs)

% calculate impulse response from exponential sweep recording
% inputs:
% rec = mono sweep recording
% T = sweep duration in seconds
% f1 = start frequency in Hz
% f2 = end frequency in Hz
% offset = offset length before impulse response in samples
% outputs:
% h = impulse response
% h_dist = impulse responses of the first 6 nonlinear distortion harmonics
%          (for detailed information see AES paper "Simultaneous Measurement of Impulse 
%           Response and Distortion with a Swept-Sine Technique" by Angelo Farina)

warning('off');
if nargin < 7
    fs=44100;
end
if nargin < 6
    offset=10;
end

if size(rec,2) > size(rec,1)
    rec=rec';
end

if size(rec,2) > 1
    error('only mono input supported');
end

[temp,filt]=expsweep(T,f1,f2,0,fs);

t_end=find(20*log10(abs(rec(end:-1:1)))>-20,1,'first');
t=length(rec)-t_end-T*fs-0.5*fs;
rec=rec(t:end,:);
filt=[filt;zeros(round((length(filt)+length(rec))/2),1)];
h_full=fftfilt(rec,filt);
h_full=h_full./max(max(abs(h_full)));

t_start=find(20*log10(abs(h_full))>-20,1,'first')-offset;
h=h_full(t_start+1:t_start+Nimp,:);

if nargout > 1
    td1=T*log(2)/log(f2/f1);
    td2=T*log(3)/log(f2/f1);
    td3=T*log(4)/log(f2/f1);
    td4=T*log(5)/log(f2/f1);
    td5=T*log(6)/log(f2/f1);
    td6=T*log(7)/log(f2/f1);
    h1=h_full(t_start-fix(fs*td1):t_start-fix(fs*td1)+Nimp-1);
    h2=h_full(t_start-fix(fs*td2):t_start-fix(fs*td2)+Nimp-1);
    h3=h_full(t_start-fix(fs*td3):t_start-fix(fs*td3)+Nimp-1);
    h4=h_full(t_start-fix(fs*td4):t_start-fix(fs*td4)+Nimp-1);
    h5=h_full(t_start-fix(fs*td5):t_start-fix(fs*td5)+Nimp-1);
    h6=h_full(t_start-fix(fs*td6):t_start-fix(fs*td6)+Nimp-1); 
    
    h_dist=[h1 h2 h3 h4 h5 h6];
end