function [B,BINT,R,RINT,STATS, d, tcf] = mri_GLM(tc, EVs, TR, filter);

% function [B,BINT,R,RINT,STATS, d, tcf] = mri_GLM(tc, EVs, TR, filter);
%
% Aug 2003, Farshad Moradi
% California Institute of Technology
%
% Example:
%
% %         EVs(1).name = 'Stimulus 1';
% %         EVs(1).t = t1; 
% %         EVs(1).dur = 0.200; % 200 ms
% %         EVs(2).name = 'Stimulus 2';
% %         EVs(7).t = t2; 
% %         EVs(7).dur = 0.300; % 300 ms
% %         [B,BINT,R,RINT,STATS, d, tcf] = mri_GLM(tc, EVs, TR, 5);
%


if nargin<4,
    filter = 0;
end;

resolution = 10;
margin = 10; % secs
N = length(tc);
n = (N*TR+margin)*resolution;
m = margin*resolution;
hrf = gampdf(0:1/resolution:30, 5, 1)-gampdf(0:1/resolution:30, 16, 1)/6;
hrf = hrf/sum(hrf(:))*resolution;

design = [];
for i=1:length(EVs),
    t1 = round((EVs(i).t+margin)*resolution+1);
    t2 = round((EVs(i).t+EVs(i).dur+margin)*resolution+1);
    t = [];
    for j=1:length(t1),
        t = [t, t1(j):t2(j)];
    end;
    f = find(t>0 & t<=n);
    temp = zeros(n,1);
    temp(t(f)) = 1;
    temp = conv(temp, hrf');
    temp = temp(m+1:TR*resolution:end);
    temp = temp(1:N);
    design = [design, temp];
end;

if filter>0,
    lp = exp(-(-4:1/filter:4).^2/2); lp = lp/sum(lp);
    w = (length(lp)-1)/2;
    tcf = tcfilter(tc, lp, w);
    for i=1:size(design,2),
        design(:,i) = tcfilter(design(:,i), lp, w);
    end;
else,
    tcf = tc;
end;
d = design;

design = [ones(N,1) design];
[B,BINT,R,RINT,STATS] = REGRESS(tcf', design, 0.05);
% plot(1:length(tc), tcf, 1:length(tc), (design*B)', 1:length(tc), sum(design(:,2:end)'));

function y = tcfilter(x, lp, w);
	baseline = conv(x-mean(x), lp);
	y = x - baseline(w: length(x)+w-1);