% program to estimate the optimal high-pass cut-off period
% fitted for an analysis dir on 64 highpass filter analysis

% SPM creates the following high-pass counfounds in this order:
% K1=1440 sec (or a period twice of the scan duration = 720/0.5)
% K2= 720 sec (720/1.0)
% K3= 480 sec (720/1.5)
% K4= 360 sec (720/2.0)
% K5= 288 sec (720/2.5)
% K6= 240 sec (720/3.0)
% K7= 206 sec (720/3.5)
% K8= 180 sec (720/4.0)
% K9= 160 sec (720/4.5)
% K10=144 sec (720/5.0)
% K11=131 sec (720/5.5)
% K12=120 sec (720/6.0)
% K13=111 sec (720/6.5)
% K14=103 sec (720/7.0)
% K15= 96 sec (720/7.5)
% K16= 90 sec (720/8.0)
% K17= 85 sec (720/8.5)
% K18= 80 sec (720/9.0)
% K19= 76 sec (720/9.5)
% K20= 72 sec (720/10.0)
% K21= 69 sec (720/10.5)
% K22= 65 sec (720/11.0)
% K23= 63 sec (720/11.5)
% K24= 60 sec (720/12.0) 
% K25= 58 sec (720/12.5)
% K26= 55 sec (720/13.0)
% K27= 53 sec (720/13.5)
% K28= 51 sec (720/14.0)
% K29= 50 sec (720/14.5)
% K30= 48 sec (720/15.0)
% .... and then it gives an error .....

load SPM

% the following variables are found in SPM
% C - (#frames)x(#contrasts) matrix, containing X's for non-0 contrast value
% B - (#frames)x(#scans) matrix indicating the scan (ith 1s)
% CONTRAST - (#(stimul_types*scans+confounds) x #(countrasts +?zeroes)


C1=CONTRAST(:,1:9);
size_G=min(size(G));
%for i=1:size_G/3
for i=1:30
GG=G(:,1:i);
GGG=[GG';GG';GG']';

CC=[C1';zeros(3,3+min(size(GGG)))']';

D=[C';B';GGG']';
v=CC*pinv(D'*D)*CC';
var1(i)=v(1,1);
var2(i)=v(2,2);
var3(i)=v(3,3);
end

D=[C';B']';
CC=[C1';zeros(3,3)']';
v=CC*pinv(D'*D)*CC';
var01=v(1,1)
var02=v(2,2)
var03=v(3,3)


% "cut_off" - a vector of the possible highpass values

for j = 1:30

cut_off(j)=(720 / (j*0.5));

end

% To plot the variance of each of the 3 contrasts vs. the
% possible high-pass cut_off periods:
 
plot (cut_off(3:30),var1(3:30),'w',cut_off(3:30),var2(3:30),'g',cut_off(3:30),var3(3:30),'r')
title('Contrast variance vs. high-pass cut-off periods')

% To plot the SNR of each of the 3 contrasts vs. the
% possible high-pass cut_off periods [SNR = 1/sqrt(var)]

for i = 1:30
SNR1(i) = 1/sqrt(var1(i));
SNR2(i) = 1/sqrt(var2(i));
SNR3(i) = 1/sqrt(var3(i));
end

figure

plot (cut_off(3:30),SNR1(3:30),'w',cut_off(3:30),SNR2(3:30),'g',cut_off(3:30),SNR3(3:30),'r')
title('SNR for each contrast vs. high-pass cut-off periods')






