NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_CORE/datat.m
function [mixt, sm_sources, tm_sources] = datat(nr, nc, nt)

% Each col of mixt contains an image.
% nr=20; nc=20; [m U V]=datat(nr,nc,nr*nc);

if(nr ~= nc) error('nr=/=nc'); end
n = nr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKE DISC IMAGES.
[x, y] = meshgrid(linspace(-1, 1, n));


s1=rand(nr*nc,1)+0.5*randk(nr*nc,1,0.1); s1=s1(:);

s1=randk(nr*nc,1,0.1); s1=s1(:);
s2=randk(nr*nc,1,0.1);;s2=s2(:);
s3=randk(nr*nc,1,0.1);;s3=s3(:); 
s4=randk(nr*nc,1,0.1);s4=s4(:);
s_sources=[s1 s2 s3 s4];

y=randk(32,32,0.1,4);s_sources=reshape(y,32*32,4);

disc=0;
if disc
radius = 0.3;
X=-0.5;Y=-0.5;	im1 = make_disc_image(nr, nc, X,Y,radius);	im1=im1(:);
X=-0.5;Y=0.5;	im2 = make_disc_image(nr, nc, X,Y,radius);	im2=im2(:);
X=0.5;Y=-0.5;	im3 = make_disc_image(nr, nc, X,Y,radius);	im3=im3(:);
X=0.5;Y=0.5;	im4 = make_disc_image(nr, nc, X,Y,radius);	im4=im4(:);
ramp=linspace(-1,1,nt);   
im1=im1 +ramp';
s_sources = [im1 im2 im3 im4]; 
end;

noise_disc=1;
if noise_disc
nfac=0.5; % for neuroimage 2001 results!
radius = 0.15; % for neuroimage 2001 results!
nradius = 0.075; % for neuroimage 2001 results!
mradius = 0.05; % for neuroimage 2001 results!

X=-0.5;Y=-0.5;	im1 = make_disc_image(nr, nc, X,Y,radius);	
X=-0;Y=0;	im1b = make_disc_image(nr, nc, X,Y,mradius);	
im1=im1+im1b; 
X=-0.5;Y=0.5;
im1 = (1-nfac)*im1+nfac*make_disc_noise_image(nr, nc, X,Y,nradius);	%imagesc(im1);
im1=im1(:); im1 = im1-mean(im1);  im1 = im1/std(im1);
%im1=im1(:)+randn(size(im1(:)))*nfac;

X=-0.5;Y=0.5;	im2 = make_disc_image(nr, nc, X,Y,radius);	
X=-0;Y=0;	im2b = make_disc_image(nr, nc, X,Y,mradius);	
im2=im2+im2b; 
X=-0.5;Y=-0.5;	
im2 = (1-nfac)*im2+ nfac*make_disc_noise_image(nr, nc, X,Y,nradius); %imagesc(im2);
im2=im2(:); im2 = im2-mean(im2);  im2 = im2/std(im2);
%im2=im2(:)+randn(size(im2(:)))*nfac;

X=0.5;Y=-0.5;	im3 = make_disc_image(nr, nc, X,Y,radius);	
X=-0;Y=0;	im3b = make_disc_image(nr, nc, X,Y,mradius);	
im3=im3+im3b; 
X=0.5;Y=0.5;	
im3 = (1-nfac)*im3+nfac*make_disc_noise_image(nr, nc, X,Y,nradius);	%imagesc(im3);
im3=im3(:); im3 = im3-mean(im3);  im3 = im3/std(im3);
%im3=im3(:)+randn(size(im3(:)))*nfac;

X=0.5;Y=0.5;	im4 = make_disc_image(nr, nc, X,Y,radius);	
X=-0;Y=-0;	im4b = make_disc_image(nr, nc, X,Y,mradius);	
im4=im4+im4b;
X=0.5;Y=-0.5;	
im4 = (1-nfac)*im4+nfac*make_disc_noise_image(nr, nc, X,Y,nradius); %imagesc(im4);
im4=im4(:); im4 = im4-mean(im4);  im4 = im4/std(im4);
%im4=im4(:)+randn(size(im4(:)))*nfac;

%ramp=linspace(-1,1,nt);   
%im1=im1 +ramp';
s_sources = [im1 im2 im3 im4]; 
end;

sins=0;
if sins
[x, y] = meshgrid( linspace(0, 1, nr), linspace(0, 1, nc) );
% (2) 3 5 7 9 11 13
im1 = sin(2*pi*2*x);im1=im1(:);
im2 = sin(2*pi*3*y);im2=im2(:);
im3 = sin(2*pi*5*y);im3=im3(:);
im4 = sin(2*pi*4*x) .* sin(2*pi*5*y);im4=im4(:);
end;

%ramp = x*y';
% jminmax(ramp(:)) 0 16

%ramp=ramp'./5;ramp=ramp(:);
% jsize(ramp); jsize(im3);
%im3=im3; % + ramp(:);
s_sources = [im1 im2 im3 im4]; 

% jsize(s_sources,'s_sources')
s_sources = zm_uvar_cols(s_sources);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKE TIME COURSES.

t = linspace(0, 1, nt);
t_sources = s_sources;

s1=randk(nr*nc,1,0.1); s1=s1(:);
s2=randk(nr*nc,1,0.1); s2=s2(:);
s3=randk(nr*nc,1,0.1); s3=s3(:); 
s4=randk(nr*nc,1,0.1); s4=s4(:);
t_sources=[s4 s4 s1 s2];

ss= sin(2*pi*3*t);
cc= cos(2*pi*3*t);
ramp=linspace(-1,1,nt);   
t1 = sin(2*pi*2*t).^7;          %%% NEIL MOD - HIGH KURTOSIS!!! 5 is cusp
t2 = sin(2*pi*3*t).^7;
t3 = sin(2*pi*5*t).^7;  
t4 = sin(2*pi*11*t).^7;         % + 3*ramp; % no ramp for neuroimage 2001 results

t_sources = [t1' t2' t3' t4']; 

mixing = 0.85;

t_sources = t_sources*(mixing*eye(4) + (1-mixing)*randn(4));

%jsize(t_sources,'t_sources')
t_sources = zm_uvar_cols(t_sources);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISPLAY SOURCES.

show4sequences(t_sources,nr,nc,2,'Source');
show4images(s_sources,nr,nc,1,'Source');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INSERT CORRELATIONS BETWEEN SPATIAL SOURCES.

a = 0.0;
sm = (1-a)*eye(4) + a*randn(4,4);
% MIX sources.
sm_sources = s_sources*sm;
p=1; % p=1 => no effect.

% Add Gauss noise to spatial sources.
noise_sources = randn(size(sm_sources));
q=1;
noise_sources = sign(noise_sources).* (abs(noise_sources).^q);
noise_sources = zm_uvar_cols(noise_sources);
noise_fac = 0.0;
sm_sources = (1-noise_fac).*sm_sources + noise_fac.*noise_sources;

sm_sources = sign(sm_sources).* (abs(sm_sources).^p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% INSERT CORRELATIONS BETWEEN TEMPORAL SOURCES.

sm = (1-a)*eye(4) + a*randn(4,4);
tm_sources = t_sources*sm;
p=1/p;
tm_sources = sign(tm_sources).* (abs(tm_sources).^p);

% Add Gauss noise to temporal sources.
% noise_sources = randn(size(tm_sources));
tfac=0.0;
noise_fac = 0;
noise_fac = noise_fac+tfac;
tm_sources = (1-noise_fac).*tm_sources + noise_fac.*noise_sources;

% DISPLAY CORRELATED SOURCES.

show4sequences(tm_sources,nr,nc,3,'Source');
show4images(sm_sources,nr,nc,4,'Source');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKE MIXTURES.

mixt = sm_sources*t_sources';
% Add noise to mixtures ...
mixt = mixt; % + randn(size(mixt))*100;
mixt = zm_uvar_cols(mixt);
% show4images(mixt,nr,nc,4,'Mixture');

show_mixtures=0;
if show_mixtures
jfig(40);
for t=1:nt
	pnshow( reshape(mixt(:,t), nr, nc) );drawnow;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
