% TjanModel.m
% Bosco's model is to assume an equivalent input noise (like the human's)
% to account for the measured contrast sensitivity function with gratings.
% Except for that the model is ideal. My implementation adds the extra
% assumption that human threshold corresponds to an arbitrary E/N for the
% ideal, ENThreshold, which is 2 (i think) for the ideal observer, and
% in practice is around 100 for the human. Changing that factor, changes
% Neq proportionally. I think we need that because
% if we measure threshold w and w/o noise then we can estimate Neq independently of
% the threshold value of E/N and Bosco's assumption of 2 would fail.

% 8/20/02 dgp Wrote it.

f=[0 0.3 1 3 10 30];

% CSF: assume grating area of 10 square cycles.
c0=[1 0.03 .01 .003 .02 .2];
EGrating=c0.*c0;
EGrating(2:end)=EGrating(2:end)./(10*f(2:end).^2);
EGrating(1)=EGrating(1)/(10*f(2).^2);
loglog(f,EGrating)
xlabel('Spatial frequency (c/deg)');
ylabel('EGrating (s deg^2)');

% deduce Neq from CSF assuming E/N=ENThreshold.
ENThreshold=2;
Neq=EGrating./ENThreshold;
loglog(f,Neq);
xlabel('Spatial frequency (c/deg)');
ylabel('Neq (s deg^2)');

% specify power spectrum of a letter
% f=[0 0.3 1 3 10 30];
ELetter1=[.5 .07 .1 .15 .1 .07];
ELetter1=[0 .07 .1 .15 .1 .07].^4;
% ELetter1=[.8 .1 .2 .2 .1 .03];
ELetter1=ELetter1./sum(ELetter1);
loglog(f,ELetter1);
xlabel('Spatial frequency (c/deg)');
ylabel('ELetter1 (s deg^2)');
title('Letter spectrum');

for fi=1:length(f)
	fCutoff=f(fi);
	% for any noise spectrum N, add in Neq to get Nef.
	N0=4e-5; % s deg^2
	NLo=(f>fCutoff).*N0;
	NHi=(f<=fCutoff).*N0;
% 	semilogx(f,N)
% 	axis([f(2) f(end) 0 N0])
	NefLo=NLo+Neq;
	NefHi=NHi+Neq;
	figure(1)
	loglog(f,NefLo,'-',f,NefHi,'--');
	hold on
	xlabel('Spatial frequency (c/deg)');
	ylabel('Nef (s deg^2)');
	title('Effective noise spectrum');
	
	EN1Lo=ELetter1./NefLo;
	EN1Hi=ELetter1./NefHi;
	figure(2);
	loglog(f,EN1Lo,'-',f,EN1Hi,'--');
	hold on
	xlabel('Spatial frequency (c/deg)');
	ylabel('E/N at contrast 1.');
	title('Stimulus snr spectrum');
	
	cSquaredLo(fi)=ENThreshold/sum(EN1Lo);
	cSquaredHi(fi)=ENThreshold/sum(EN1Hi);
	cSquared0(fi)=ENThreshold/sum(ELetter1./Neq);
end
figure(3)
loglog(f,cSquaredLo-cSquared0,'-',f,cSquaredHi-cSquared0,'--');
hold on
loglog(f,cSquared0);
xlabel('Cut-off frequency (c/deg)');
ylabel('c^2');
title('Contrast power threshold vs noise cut-off frequency');
