function [CVR,ts]=CVR(time_etco2,etco2,time_fnirs,yfnirs)

%% time_etco2: time points issue from capnograph
%% etco2 : etco2 value at the time points issue from capnograph (same lenght as time_etco2)
%% time_fnirs : one column time points where fnirs data were acquired.
%% yfnirs : fnirs variable (HBO, HbR or HbO-HbR). dimension: 16 columns x lenghts of time_fnirs (16 channels)

%% Mean value of fNIRS signal for delay measurement.
mean_yfnirs=mean(yfnirs,2);

%% interpolate yfnirs values on etco2 timepoint.
YFNIRS=interp1(time_fnirs,mean_yfnris,time_etco2);

%% time delay
%% shift etco2 for finding the maimum correlation with YFNIRS
deltatime=size(etco2,2)-size(YFNIRS,2);
for i=1:deltatime-1
    tempcorr=corrcoef(YFNIRS,etco2(1,1+i:size(YFNIRS,2)+i)');
    c(i)=tempcorr(2,1);
end

[maxCC,j]=max(c);
Etco2_shifted=etco2(1,j+1:size(YFNIRS,2)+j)';
ts=time_etco2(j);

%% CVR values. Correlation between the shifted etco2 and fNIRS on each channel (16 CVR values)

for i=1:16
    X = [ones(size(Etco2_shifted)) yfnirs(:,i)];
    s_coeffs(:,i) = X\Etco2_shifted;
end

CVR=s_coeffs(2,:);