destination.m

% Optimize parameters C and D to match
% the exponential function our empirical
% shape within a large a large shift range
% P. G. Sämann 2016, E-Mail: saemann@psych.mpg.de

% The file 'destination.mat' can be created by the script hrf_generation.m
clear all
load ~/destination.mat
x=0; y=0;
resultmatrix (1:41,1:81)=1000;
figure, plot(d,delay,'r.'), hold on;

% The following two correction values were added
% to correct for shifts of the entire sigmoid
% transformation function over the total range
% in the x and y direction. See formula below.

Ycorr =  0.25; 
Xcorr = -0.50;

% Loop around combinations of C and D
for C= -3:0.25:-2 % C range for optimization
    x = x + 1;
    y = 0;
    for D = -4:0.05:0 % D range for optimization
        y = y +1;
        s = 1;
        for ratio = -20:0.025:20
            testfunction(s) = 2 * C/(1 + exp(D*(ratio-Xcorr))) - C + Ycorr;
            s = s + 1;
        end
        plot(d,testfunction,'b'); drawnow;
        % now subtract both functions from eachother
        diff = abs(delay - testfunction);
        diffsum = sum(diff);
        tmpmin = min(resultmatrix(:));
        resultmatrix(x,y)=diffsum;
        if diffsum < tmpmin
            storeC = C;
            storeD = D;
            fprintf(['New optimum for C: ',num2str(C),' and D: ', num2str(D),' diff: ', num2str(diffsum), '\n']);
        end
    end
    plot(d,delay,'r.'),
    fprintf(['X = ',num2str(x),'\n']);
end

% Set C and D to achieved optima:
C = storeC;
D = storeD;

% Calculate resulting 'testfunction_final'
% and compare it with gold standard
s = 1;
for ratio = -20:0.025:20
    testfunction_final(s) = 2*C /(1+exp(D*(ratio-Xcorr))) - C + Ycorr;
    d_testfunction(s) = ratio;
    s = s + 1; 
end

figure, plot(d,delay,'.'), hold on,
plot(d_testfunction,testfunction_final,'r');
title ('In blue: goldstandard (destination.mat), in red: achieved by C, D, xcorr, ycorr');