hrf_generation.m

% Simulation: Study on the relationship between different beta2/beta1
% ratios and empirical peak position resulting from the summed canonical
% HRF and its first derivative.
%
% P. Sämann 2016, E-Mail: saemann@psych.mpg.de

[hrf p]  = spm_hrf(0.0625, [6 16 1 1 6 0 32]);
[hrf2 p] = spm_hrf(0.0625, [6 16 1 1 6 1 32]);
hrfd     = 1*(hrf - hrf2); % order as in spm_get_bf.m

ortho = 1; % Set to 1 if an orthogonalization of the
% difference function hrfd on the canonical (hrf) is
% sought as described in Henson et al. 2002.

if ortho == 1
    m1  = mean(hrf);  hrf  = hrf - m1;
    m2  = mean(hrfd); hrfd = hrfd - m2;
    tmp = [hrf hrfd];
    tmp = spm_orth(tmp);
    hrf = tmp(:,1); hrfd=tmp(:,2);
    hrf = hrf + m1; hrfd = hrfd + m2;
end

figure, plot(hrf), hold on, plot(hrfd,'r'), title('HRF and 1 s-pseudoderivative');
figure, plot(hrf), hold on

s = 1;
for i=-20:0.025:20
    hrf_sum = hrf + i*hrfd;
    col='k';
    if i==0;
        plot(hrf_sum, 'r','LineWidth',2);
    else
        plot(hrf_sum, 'k');
    end
    tmp = find(hrf_sum==max(hrf_sum));
    m = tmp*0.0625;
    fprintf(['i set to ',num2str(i),' results in delay of ',num2str(m),'\n']);
    delay_tmp(s) = m;
    d(s) = i;
    s = s + 1;
end

s = 1;
for i=-20:0.025:20; delay(s) = delay_tmp(s) - delay_tmp(801); s = s + 1; end
figure, plot(d, delay, 'k.-'), title ('x-axis = factor d, y-axis = resulting delay in s'); hold on

% Comment: The resulting d- and delay-values were stored as destination.mat
% and used as gold standard to optimize the coefficients C, D as well as 
% correction coefficients xcorr and ycorr.