% Look at how well use of the temporal derivative is
% in capturing delayed responses

TR=1;

% Get canonical
[h p]=spm_hrf(TR);

% Get temporal
% in SPM `proper' this is implemented in spm_get_bf lines 135-137
dp     = 1;
p6 = p(6);
p(6)   = p6 + dp;
d = (h - spm_hrf(TR,p))/dp;

X = [h d];

N=size(X,1);
dd=0.5*[1:1:9];
for i=1:length(dd),
    
    % True response starts later
    p(6)=p6+dd(i);
    y=spm_hrf(TR,p);
    vy=std(y)^2;
    
    % Estimated response
    beta=pinv(X)*y;
    hy=X*beta;
    evy=std(y-hy)^2;
    
    subplot(3,3,i);
    plot(y);
    hold on
    plot(hy,'r');
    title(sprintf('Delay=%1.2f, b=[%1.2f,%1.2f]',dd(i),beta(1),beta(2)));
    grid on

    r(i)=(vy-evy)/vy;
end


figure
plot(dd,r);
xlabel('Delay, s');
ylabel('R^2');
title('Proportion of variance explained using HRF + Deriv');
grid on