function NoiseDemo2
% NoiseDemo2. Make and show, side by side, two movies consisting of dynamic white noise.
% For each movie you can adjust check size, check duration, and choose any of several amplitude 
% distributions: binary, uniform, or gaussian. 
%
% See also NoiseDemo, MovieDemo, GratingDemo, and DriftDemo.

% 2/26/02	dgp	 Wrote it.

fprintf('NoiseDemo2\n');

% Open a window 
screenNumber=0;
noiseKind1='binary';
noiseKind2='binary';
imageSize=256; % width & height of displayed image, in pixels
checkSize1=4; % width & height of each noise check, in pixels
checkSize2=2; % width & height of each noise check, in pixels
checkFrames1=1; % duration of each noise check, in frames
checkFrames2=1; % duration of each noise check, in frames
contrast1=1;
contrast2=1;
frames=2*FrameRate(screenNumber); % 2 second's worth
imageSize=round(imageSize/checkSize1)*checkSize1;
frames=round(frames/checkFrames1)*checkFrames1;
black=BlackIndex(screenNumber);
white=WhiteIndex(screenNumber);
gray=(white+black)/2;
% gray=round((white+black)/2);
window=Screen(screenNumber,'OpenWindow',gray,[],8);
black=BlackIndex(screenNumber);
white=WhiteIndex(screenNumber);
rand('state',getsecs); % random seed

% noiseList
noiseList=MakeNoiseList(noiseKind1);
noiseListPower2=std(noiseList,1)^2;

% Make movie 1.
inc=white-gray;
checks=imageSize/checkSize1;
rect=round([0 0 imageSize imageSize]);
noise=zeros(checks,checks);
for i=1:frames
	if checkFrames1==1 | mod(i,checkFrames1)==1
		noise=gray+contrast1*inc*randsample(noiseList,size(noise)); % one sample per check
	end
    w1(i)=Screen(screenNumber,'OpenOffscreenWindow',gray,rect);
    Screen(w1(i),'PutImage',noise,rect); % expand each sample to make a check
end;
w1(frames+1)=Screen(screenNumber,'OpenOffscreenWindow',gray,rect);

unique(noise(:));

% noiseList
noiseList=MakeNoiseList(noiseKind2);
noiseListPower1=std(noiseList,1)^2;

% Make movie 2.
inc=white-gray;
checks=imageSize/checkSize2;
rect=round([0 0 imageSize imageSize]);
noise=zeros(checks,checks);
for i=1:frames
	if checkFrames2==1 | mod(i,checkFrames2)==1
		noise=gray+contrast2*inc*randsample(noiseList,size(noise)); % one sample per check
	end
    w2(i)=Screen(screenNumber,'OpenOffscreenWindow',gray,rect);
    Screen(w2(i),'PutImage',noise,rect); % expand each sample to make a check
end;
w2(frames+1)=Screen(screenNumber,'OpenOffscreenWindow',gray,rect);

% Show the movies, side by side, using Rush to minimize interruptions.
rect1=CenterRect(rect,Screen(window,'Rect'));
rect1=OffsetRect(rect1,-round(RectWidth(rect)/2),0);
rect2=OffsetRect(rect1,RectWidth(rect),0);
loop={
'for j=1:3;'
	'for i=1:frames;'
		'Screen(window,''WaitBlanking'');'
		'Screen(''CopyWindow'',w1(i),window,rect,rect1);'
		'Screen(''CopyWindow'',w2(i),window,rect,rect2);'
	'end;'
'end;'
'Screen(window,''WaitBlanking'');'
'Screen(''CopyWindow'',w1(frames+1),window,rect,rect1);'
'Screen(''CopyWindow'',w2(frames+1),window,rect,rect2);'
};
HideCursor;
priorityLevel=MaxPriority(screenNumber,'WaitBlanking');
Rush(loop,priorityLevel);
ShowCursor;
Screen('CloseAll');

% Report the power spectal density of the noise.
fprintf('Dynamic white noise: two versions, left (and right), drawn from a %s\n',noiseKind1);
fprintf('(vs %s) distribution. ',noiseKind2);
fprintf('The rms contrast is\n');
fprintf('crms = %g (vs %g)\n', contrast1*noiseListPower1^0.5*(inc/gray), contrast2*noiseListPower2^0.5*(inc/gray));
fprintf('However, the rms contrast (and its squared value, the contrast power) is\n');
fprintf('of little relevance in specifying this as a visual stimulus. Usually\n');
fprintf('your noise will be white, i.e. have a flat spectrum out to a bandwidth\n');
fprintf('that exceeds that of the mechanism (visual channel) under study. The\n');
fprintf('total power reflects the noise bandwidth, most of which is invisible to\n');
fprintf('the channel. The proper statistic to specify the noise as a visual\n');
fprintf('stimulus is the power spectral density N. It is the contrast power per\n');
fprintf('unit frequency interval (i.e. 1 c/deg c/deg Hz). For white noise made\n');
fprintf('up of independent checks, N equals the product of the contrast power\n');
fprintf('and the area*duration of one check.\n');
fprintf('N = %g sec pixel^2 ',noiseListPower1*(inc/gray)^2*checkSize1*checkSize1*checkFrames1/FrameRate(screenNumber));
fprintf('(vs %g sec pixel^2)\n',noiseListPower2*(inc/gray)^2*checkSize2*checkSize2*checkFrames2/FrameRate(screenNumber));
fprintf('Once you know the viewing distance you should convert "pixel" to "deg"\n');
fprintf('of visual angle.\n\n');
fprintf('See Pelli, D.G. and Farell, B. (1999) Why use noise? J Opt Soc Am A Opt Image Sci Vis, 16 (3), 647-53.\n');
fprintf('PDF file available at:\n');
fprintf('web http://psych.nyu.edu/people/faculty/pelli/ ;\n');

function noiseList=MakeNoiseList(noiseKind)
switch noiseKind
case 'gaussian', %clipped at +-2 standard deviations.
	noiseBound=2;
	temp=randn([1,20000]);
	noiseList=find(sign(temp.^2-noiseBound^2)-1);
	noiseList=temp(noiseList);
	clear temp;
	noiseList=noiseList/noiseBound;
case 'uniform',
	noiseList=-1:1/1024:1;
case 'binary',
	noiseList=[-1 1];
otherwise,
	fprintf('Unknown noiseKind "%s".\n',noiseKind);
	return;
end
