function [IRF,t] = NExT(y,dt,Ts,method)
%
% [IRF] = NExT(y,ys,T,dt) implements the Natural Excitation Technique to
% retrieve the Impulse Response FUnction (IRF) from the cross-correlation
% of the measured output y.
%
% [IRF] = NExT(y,dt,Ts,1) calculate the IRF with cross-correlation
% calculated by using the inverse fast fourier transform of the
% cross-spectral power densities  (method = 1).
%
% [IRF] = NExT(y,dt,Ts,2) calculate the IRF with cross-correlation
% calculated by using the unbiased cross-covariance function (method = 2)
%
%
% y: time series of ambient vibrations: vector of size [1xN]
% dt : Time step
% method: 1 or 2 for the computation of cross-correlation functions
% T: Duration of subsegments (T<dt*(numel(y)-1))
% IRF: impusle response function
% t: time vector asociated to IRF
%%
if nargin<4, method = 2; end % the fastest method is the default method
if ~ismatrix(y), error('Error: y must be a vector or a matrix'),end


[Nyy,N]=size(y);
if Nyy>N,
    y=y';
    [Nyy,N]=size(y);
end

% get the maximal segment length fixed by T
M = round(Ts/dt);
switch method
    case 1
        clear IRF
        for ii=1:Nyy,
            for jj=1:Nyy,
                y1 = fft(y(ii,:));
                y2 = fft(y(jj,:));
                h0 = ifft(y1.*conj(y2));
                IRF(ii,jj,:) = h0(1:M);
            end
        end
        % get time vector t associated to the IRF
        t = linspace(0,dt.*(size(IRF,3)-1),size(IRF,3));
        if Nyy==1,
            IRF = squeeze(IRF)'; % if Nyy=1
        end
    case 2
        IRF = zeros(Nyy,Nyy,M+1);
        for ii=1:Nyy,
            for jj=1:Nyy,
                [dummy,lag]=xcov(y(ii,:),y(jj,:),M,'unbiased');
                IRF(ii,jj,:) = dummy(end-round(numel(dummy)/2)+1:end);
            end
        end
        if Nyy==1,
            IRF = squeeze(IRF)'; % if Nyy=1
        end
        % get time vector t associated to the IRF
        t = dt.*lag(end-round(numel(lag)/2)+1:end);
end
% normalize the IRF
if Nyy==1,
IRF = IRF./IRF(1);
else
end