diff --git a/linearalgebra/simulations/covareigen.m b/linearalgebra/simulations/covareigen.m
new file mode 100644
index 0000000..1acf330
--- /dev/null
+++ b/linearalgebra/simulations/covareigen.m
@@ -0,0 +1,82 @@
+function covareigen( x, y, pm, w, h )
+% computes covariance matrix from the pairs x, y
+% diagonalizes covariance matrix
+% x and y are column vectors
+% pm: 0 - scatter of data, 1 histogram, 2 multivariate gauus from
+% covariance matrix, 1 and 2 require w and h
+% w and h are the width and the height (in data units) for plotting
+% histograms instead of scatter
+
+ % covariance matrix:
+ cv = cov( [ x y ] );
+
+ % eigen values:
+ [ v , d] = eig( cv )
+ s = sign( v );
+ s(1,:) = s(2,:);
+ v = v .* s;
+
+ % plots:
+ subplot( 2, 2, 1 );
+ hold on;
+
+ if (pm > 0) & (nargin == 5)
+ if pm == 1
+ % histogram of data:
+ xp = -w:0.5:w;
+ yp = -h:0.5:h;
+ [n,c] = hist3([x y], { xp, yp } );
+ contourf( c{1}, c{2}, n' );
+ else
+ xp = -w:0.1:w;
+ yp = -h:0.1:h;
+ [xg,yg] = meshgrid( xp, yp );
+ xy = [ xg(:) yg(:) ];
+ gauss = reshape( exp(-0.5*diag(xy*inv(cv)*xy'))/sqrt((2.0*pi)^2.0*det(cv)), size( xg ) );
+ contourf( xp, yp, gauss )
+ end
+ colormap( 'gray' );
+ else
+ % scatter plot:
+ scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
+ end
+
+ % plot eigenvectors:
+ quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
+
+ axis( 'equal' );
+ hold off;
+
+ subplot( 2, 2, 3 );
+ hist( x, 50, 'b' );
+
+ % sort the eigenvalues:
+ [d,inx] = sort( diag(d), 'descend' );
+
+ subplot( 2, 2, 2 );
+ hold on;
+ % subtract means:
+ x = x - mean( x );
+ y = y - mean( y );
+ % project onto eigenvectors:
+ nx = [ x y ] * v(:,inx(1));
+ ny = [ x y ] * v(:,inx(2));
+ cv = cov( [nx, ny] )
+ [ v , d] = eig( cv )
+ if (pm > 0) & (nargin == 5)
+ if pm == 1
+ [n,c] = hist3([nx ny], { xp, yp } );
+ contourf( c{1}, c{2}, n' );
+ else
+ gauss = reshape( exp(-0.5*diag(xy*inv(cv)*xy'))/sqrt((2.0*pi)^2.0*det(cv)), size( xg ) );
+ contourf( xp, yp, gauss )
+ end
+ else
+ scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
+ end
+ axis( 'equal' );
+ hold off;
+
+ subplot( 2, 2, 4 );
+ hist( nx, 50, 'b' );
+end
diff --git a/linearalgebra/simulations/covareigen3.m b/linearalgebra/simulations/covareigen3.m
new file mode 100644
index 0000000..c11f2de
--- /dev/null
+++ b/linearalgebra/simulations/covareigen3.m
@@ -0,0 +1,45 @@
+function covareigen3( x, y, z )
+% computes covariance matrix from the triples x, y, z
+% diagonalizes covariance matrix
+% x, y and z are column vectors
+
+ % covariance matrix:
+ cv = cov( [ x y z ] );
+
+ % eigen values:
+ [ v , d] = eig( cv )
+
+ % plots:
+ subplot( 1, 2, 1 );
+ hold on;
+
+ % scatter plot:
+ view( 3 );
+ scatter3( x, y, z, 0.1, 'b', 'filled', 'MarkerEdgeColor', 'blue' );
+ xlabel( 'x' );
+ ylabel( 'y' );
+ zlabel( 'z' );
+ grid on;
+
+ % plot eigenvectors:
+ quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
+
+ %axis( 'equal' );
+ hold off;
+
+ % sort the eigenvalues:
+ [d,inx] = sort( diag(d), 'descend' );
+
+ subplot( 2, 2, 2 );
+ hold on;
+ % subtract means:
+ x = x - mean( x );
+ y = y - mean( y );
+ % project onto eigenvectors:
+ nx = [ x y z ] * v(:,inx(1));
+ ny = [ x y z ] * v(:,inx(2));
+ scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
+ axis( 'equal' );
+ hold off;
+
+end
diff --git a/linearalgebra/simulations/covareigen3examples.m b/linearalgebra/simulations/covareigen3examples.m
new file mode 100644
index 0000000..3d16628
--- /dev/null
+++ b/linearalgebra/simulations/covareigen3examples.m
@@ -0,0 +1,14 @@
+n = 10000;
+
+% three distributions:
+x = randn( n, 1 );
+y = randn( n, 1 );
+z = randn( n, 1 );
+f = figure( 'Position', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ]);
+d = [ 0 4 -4 8 ];
+for k = 1:4
+ x((k-1)*n/4+1:k*n/4) = x((k-1)*n/4+1:k*n/4) + d(k);
+ y((k-1)*n/4+1:k*n/4) = y((k-1)*n/4+1:k*n/4) - d(k);
+ z((k-1)*n/4+1:k*n/4) = z((k-1)*n/4+1:k*n/4) + d(k);
+end
+covareigen3( x, y, z );
diff --git a/linearalgebra/simulations/covareigenexamples.m b/linearalgebra/simulations/covareigenexamples.m
new file mode 100644
index 0000000..f360e8c
--- /dev/null
+++ b/linearalgebra/simulations/covareigenexamples.m
@@ -0,0 +1,35 @@
+scrsz = get( 0, 'ScreenSize' );
+set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
+
+% correlation coefficients:
+n = 10000;
+x = randn( n, 1 );
+f = figure( 1 );
+for r = 0.01:0.199:1
+ fprintf( 'Correlation = %g\n', r );
+ clf( f );
+ y = r*x + sqrt(1-r^2)*randn( n, 1 );
+ covareigen( x, y, 2, 5.0, 3.0 );
+ key = waitforbuttonpress;
+end
+return
+
+% two distributions:
+n = 10000;
+x1 = randn( n/2, 1 );
+y1 = randn( n/2, 1 );
+x2 = randn( n/2, 1 );
+y2 = randn( n/2, 1 );
+f = figure( 1 );
+pause( 'on' );
+for d = 0:1:5
+ fprintf( 'Distance = %g\n', d );
+ clf( f );
+ d2 = d / sqrt( 2.0 );
+ x = [ x1; x2 ];
+ y = [ y1+d2; y2-d2 ];
+ scrsz = get(0,'ScreenSize');
+ covareigen( x, y, 1, 10.0, 7.0 );
+ %key = waitforbuttonpress;
+ pause( 1.0 );
+end
diff --git a/linearalgebra/simulations/matrix2dtrafos.m b/linearalgebra/simulations/matrix2dtrafos.m
new file mode 100644
index 0000000..da6840e
--- /dev/null
+++ b/linearalgebra/simulations/matrix2dtrafos.m
@@ -0,0 +1,74 @@
+m = [ 1 0; 0 1];
+matrixbox( m );
+title( 'Identity' );
+waitforbuttonpress;
+
+m = [ 2 0; 0 1];
+matrixbox( m );
+title( 'Scale x' );
+waitforbuttonpress;
+
+m = [ 1 0; 0 2];
+matrixbox( m );
+title( 'Scale y' );
+waitforbuttonpress;
+
+m = [ 2 0; 0 2];
+matrixbox( m );
+title( 'Scale both' );
+waitforbuttonpress;
+
+m = [ -1 0; 0 1];
+matrixbox( m );
+title( 'Flip x' );
+waitforbuttonpress;
+
+m = [ 1 0; 0 -1];
+matrixbox( m );
+title( 'Flip y' );
+waitforbuttonpress;
+
+m = [ -1 0; 0 -1];
+matrixbox( m );
+title( 'Flip both' );
+waitforbuttonpress;
+
+m = [ 1 0; 1 1];
+matrixbox( m );
+title( 'Shear x' );
+waitforbuttonpress;
+
+m = [ 1 1; 0 1];
+matrixbox( m );
+title( 'Shear y' );
+waitforbuttonpress;
+
+phi = 0.1667*pi;
+m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
+matrixbox( m );
+title( 'Rotate 30' );
+waitforbuttonpress;
+
+phi = 0.333*pi;
+m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
+matrixbox( m );
+title( 'Rotate 60' );
+waitforbuttonpress;
+
+phi = 0.5*pi;
+m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
+matrixbox( m );
+title( 'Rotate 90' );
+waitforbuttonpress;
+
+phi = 0.75*pi;
+m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
+matrixbox( m );
+title( 'Rotate 135' );
+waitforbuttonpress;
+
+phi = 1.0*pi;
+m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
+matrixbox( m );
+title( 'Rotate 180' );
+waitforbuttonpress;
diff --git a/linearalgebra/simulations/matrixbox.m b/linearalgebra/simulations/matrixbox.m
new file mode 100644
index 0000000..5abbfa0
--- /dev/null
+++ b/linearalgebra/simulations/matrixbox.m
@@ -0,0 +1,24 @@
+function matrixbox( m )
+% visualizes the effect of a matrix m on a set of vectors forming a box
+% m: a 2x2 matrix
+
+ v = [ 0 1 1 0 0; 0 0 1 1 0];
+ w = m*v;
+ clf;
+ hold on;
+ %set(gca, 'Xtick', [], 'Ytick', [], 'box', 'off')
+ % axis:
+ plot( [-2 2], [0 0], 'k', 'LineWidth', 1 );
+ plot( [0 0], [-2 2], 'k', 'LineWidth', 1 );
+ % old box:
+ plot( v(1,:), v(2,:), 'k', 'LineWidth', 1.0 )
+ quiver( [v(1,1)], [v(2,1)], [v(1,3)], [v(2,3)], 1.0, 'k', 'LineWidth', 2.0 );
+ scatter( [v(1,2)], [v(2,2)], 60.0, 'filled', 'k' );
+ % transfomred box:
+ plot( w(1,:), w(2,:), 'b', 'LineWidth', 2.0 )
+ quiver( [w(1,1)], [w(2,1)], [w(1,3)], [w(2,3)], 1.0, 'b', 'LineWidth', 3.0 );
+ scatter( [w(1,2)], [w(2,2)], 100.0, 'filled', 'b' );
+ hold off;
+ xlim( [ -2 2 ] );
+ ylim( [ -2 2 ] );
+end
diff --git a/linearalgebra/simulations/plotvector.m b/linearalgebra/simulations/plotvector.m
new file mode 100644
index 0000000..668818f
--- /dev/null
+++ b/linearalgebra/simulations/plotvector.m
@@ -0,0 +1,22 @@
+function plotvector( v, s, sp )
+% plot an arrow for the 2D-vector v
+% v: the 2D vector
+% s: a string passed to quiver for the color
+% sp: if this string is set it defines the color of some helper lines
+
+ if nargin < 2
+ s = 'b';
+ end
+ quiver( [0.0], [0.0], [v(1)], [v(2)], 1.0, s );
+ grid on;
+ if nargin == 3
+ hh = ishold;
+ if ~hh
+ hold on;
+ end
+ plot( [0.0 v(1) v(1)], [0.0 0.0 v(2)], sp );
+ if ~hh
+ hold off;
+ end
+ end
+end
diff --git a/linearalgebra/simulations/plotvector3.m b/linearalgebra/simulations/plotvector3.m
new file mode 100644
index 0000000..d93e873
--- /dev/null
+++ b/linearalgebra/simulations/plotvector3.m
@@ -0,0 +1,24 @@
+function plotvector3( v, s, sp )
+% plot an arrow for the 3D-vector v
+% v: the 3D vector
+% s: a string passed to quiver for the color
+% sp: if this string is set it defines the color of some helper lines
+
+ if nargin < 2
+ s = 'b';
+ end
+ quiver3( [0.0], [0.0], [0.0], [v(1)], [v(2)], [v(3)], 1.0, s );
+ view( 42, 12 );
+ grid on
+ if nargin == 3
+ hh = ishold;
+ if ~hh
+ hold on;
+ end
+ plot3( [v(1) v(1)], [v(2) v(2)], [0.0 v(3)], sp );
+ plot3( [0.0 v(1) v(1)], [v(2) v(2) 0.0], [0.0 0.0 0.0], sp );
+ if ~hh
+ hold off;
+ end
+ end
+end
diff --git a/linearalgebra/simulations/spikesortingwave.m b/linearalgebra/simulations/spikesortingwave.m
new file mode 100644
index 0000000..369915e
--- /dev/null
+++ b/linearalgebra/simulations/spikesortingwave.m
@@ -0,0 +1,127 @@
+% generate spikes:
+n = 1000;
+misi = 0.01;
+isis = exprnd( misi, n, 1 );
+isis = isis + 0.01;
+spikes = cumsum( isis );
+p = rand( size( spikes ) );
+
+% spike waveforms:
+dt = 0.0001;
+x = -0.01:dt:0.01;
+y1 = 2.0*exp( -(x-0.0003).^2/2.0/0.0005^2 ) - 1.4*exp( -(x-0.0005).^2/2.0/0.0005^2 );
+y2 = exp( -x.^2/2.0/0.002^2 ).*cos(2.0*pi*x/0.005);
+p12 = 0.5;
+
+% voltage trace:
+noise = 0.2;
+t = 0:dt:spikes(end)+3.0*x(end);
+v = noise*randn( 1, length( t ) );
+for k = 1:length( spikes )
+ inx = ceil( spikes(k)/dt );
+ if p(k) < p12
+ v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y1;
+ else
+ v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y2;
+ end
+end
+
+figure( 1 );
+clf;
+plot( t(t<1.0), v(t<1.0) );
+hold on;
+
+% find peaks:
+thresh = 0.7;
+inx = find( v(2:end-1) > thresh & v(1:end-2)<v(2:end-1) & v(2:end-1) > v(3:end) ) + 1;
+spiket = t(inx);
+spikev = v(inx);
+tinx = inx;
+for k=1:2
+ inx = find( ( spiket(2:end-1)-spiket(1:end-2)>0.005 | spikev(2:end-1) > spikev(1:end-2) ) & ( spiket(3:end)-spiket(2:end-1)>0.005 | spikev(2:end-1) > spikev(3:end) ) )+1;
+ spiket = spiket(inx);
+ spikev = spikev(inx);
+ tinx = tinx(inx);
+end
+scatter( spiket(spiket<1.0), spikev(spiket<1.0), 100.0, 'r', 'filled' );
+%scatter( t(tinx), v(tinx), 100.0, 'r', 'filled' );
+hold off;
+
+% spike waveform snippets:
+w = ceil( 0.005/dt );
+vs = [];
+for k=1:length(tinx)
+ vs = [ vs; v(tinx(k)-w:tinx(k)+w-1) ];
+end
+ts = t(1:size(vs,2));
+ts = ts - ts(floor(length(ts)/2));
+figure( 2 );
+clf;
+hold on
+for k=1:size(vs,1)
+ plot( ts, vs(k,:), '-b' );
+end
+hold off
+
+% pca:
+cv = cov( vs );
+[ ev , ed ] = eig( cv );
+[d,dinx] = sort( diag(ed), 'descend' );
+% spectrum of eigenvalues:
+figure( 3 );
+clf;
+subplot( 4, 1, 1 );
+scatter( 1:length(d), d, 'b', 'filled' );
+% project onto eigenvectors:
+nx = vs * ev(:,dinx(1));
+ny = vs * ev(:,dinx(2));
+%scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
+%hold on;
+
+% features:
+subplot( 4, 2, 5 );
+plot( 1000.0*ts, ev(:,dinx(1)), 'g', 'LineWidth', 2 );
+subplot( 4, 2, 6 );
+plot( 1000.0*ts, ev(:,dinx(2)), 'r', 'LineWidth', 2 );
+
+% clustering:
+%kx = kmeans( [ nx, ny ], 2 );
+% nx smaller or greater a threshold:
+kthresh = 1.6;
+kx = ones( size( vs, 1 ), 1 );
+kx(nx>kthresh) = 2;
+
+subplot( 4, 1, 2 );
+scatter( nx(kx==1), ny(kx==1), 'g', 'filled', 'MarkerEdgeColor', 'white' );
+hold on;
+scatter( nx(kx==2), ny(kx==2), 'r', 'filled', 'MarkerEdgeColor', 'white' );
+hold off;
+
+% show sorted spike waveforms:
+subplot( 4, 2, 7 );
+hold on
+kinx1 = find(kx==1);
+for k=1:length(kinx1)
+ plot( 1000.0*ts, vs(kinx1(k),:), '-g' );
+end
+plot( 1000.0*x, y1, '-k', 'LineWidth', 2 );
+xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
+hold off
+subplot( 4, 2, 8 );
+hold on
+kinx2 = find(kx==2);
+for k=1:length(kinx2)
+ plot( 1000.0*ts, vs(kinx2(k),:), '-r' );
+end
+plot( 1000.0*x, y2, '-k', 'LineWidth', 2 );
+xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
+hold off
+
+% spike trains:
+figure( 1 );
+hold on;
+six1 = find( spiket(kinx1)<1.0 );
+scatter( spiket(kinx1(six1)), spikev(kinx1(six1)), 100.0, 'g', 'filled' );
+six2 = find( spiket(kinx2)<1.0 );
+scatter( spiket(kinx2(six2)), spikev(kinx2(six2)), 100.0, 'r', 'filled' );
+hold off;
diff --git a/linearalgebra/simulations/spiketime.m b/linearalgebra/simulations/spiketime.m
new file mode 100644
index 0000000..fb22807
--- /dev/null
+++ b/linearalgebra/simulations/spiketime.m
@@ -0,0 +1,11 @@
+p = 0.0;
+n=0;
+for k = 1:length( spikes )
+ %a= 0.5*(spikes(k+1)+spikes(k));
+ %fprintf( 'set label %d "$T_{%d}$" at %g, -0.5 center\n', k+2, k, a );
+ %fprintf( 'set arrow %d from %g, 0.3 to %g, 0.3 heads\n', k+2, spikes(k), spikes(k+1) );
+ %fprintf( '%g %d\n%g %d\n\n', p, n, spikes(k), n );
+ fprintf( '%g %d\n', spikes(k), n+1 );
+ n = n+1;
+ p = spikes(k);
+end