% Lecture_042512.m

dt = 1;  %or TR for fMRI
maxth = 15; 
th = 0:dt:(maxth-dt);  %time vector for the hdr
k= 1.25;
delay = 2; %seconds
nCascades = 3;

noise = 0;
amp =10;
h = amp*Gamma(nCascades,k,th-delay)';

%h = randn(size(h));

figure(1)
clf
plot(th,h,'ro-','MarkerFaceColor','w');
set(gca,'XLim',[0,maxth]);

%stimulus
maxt = 255;
t = [0:dt:(maxt-dt)]';

%boring 'slow event-related 'stimulus
s = zeros(size(t));
s(1:20:end) = 1;

%rapid event-related stimulus
prob = .25; 
s = floor(rand(size(t))+prob);


y = conv(s,h)*dt;
y = y(1:length(t));

y = y+noise*randn(size(y));
%dc
y = y+100;
%ramp
 y = y+linspace(0,20,length(t))';



figure(2)
clf
subplot(2,1,1)
stairs(t,s,'r-');
xlabel('Time (s)');
set(gca,'XLim',[0,maxt]);
set(gca,'YLim',[0,1.2]);

subplot(2,1,2)
plot(t,y,'b-');
xlabel('Time (s)');
set(gca,'XLim',[0,maxt]);

%% design matrix
n = length(h);
m = length(s);

X = zeros(m,n);
tmp = s;
for i=1:n
    X(:,i) = tmp;
    tmp = [0;tmp(1:end-1)];
end




figure(3)
clf
imagesc(X);
colormap(gray)
axis equal
axis off

yhat = X*h;
figure(2)
subplot(2,1,2)
hold on
plot(t,yhat,'r.');


%% deconvolution!


%dc
XX = [X,ones(m,1)];
XX = [XX,linspace(0,1,m)'];

hhat = inv(XX'*XX)*XX'*y;

%or:
hhat = pinv(XX)*y;
hhat = hhat(1:end-2);

figure(1)
hold on
plot(th,hhat,'b-');

%%

h_unscaled = Gamma(nCascades,k,th-delay)';

Xh = X*h_unscaled;

XXh = [Xh,ones(m,1),linspace(0,1,m)'];

hhath = pinv(XXh)*y;
    
    
    
    




