%% Examples of vectors and matrices

% sound file with vector
% load Matlab example file
load('handel')
% play sound file
sound(y, Fs)

% plot time course
plot(y)

% plot with correct x-axis (in milliseconds)
t = 0:1000/Fs:(length(y)-1)*1000/Fs;
plot(t,y);

% manipulate the vector
y2 = y;
y2(1:2:end) = 0;
plot(t,y2);
sound(y2, Fs);

% how to silence the section between 1-2s?
% how to delete the section between 1-2s? (take care when plotting...)


% image with matrix
img = imread('greens.jpg');
image(img)
% changing colormap won't work for colour image:
colormap(gray)

% convert to greyscale and try again
img2 = rgb2gray(img);
image(img2)
colormap(winter)
colormap(summer)
colormap(hot)
colormap(cool)

size(img)  % colour - 3 dimensions
size(img2) % greyscale - 2 dimensions

% changing rows
img2(1:2:end,:) = 0;
image(img2)
img2(1:2:end,:) = 255;
image(img2)

% how to change columns?

%% Vector basics

%% a vector is an ordered sequence of numbers
x = [2 4 6 8]

y = [2:2:8]
z = [8:-2:2]

% getting an element out of a vector
y(3)

% basic operations
x/2
2*x
whos

% Element-wise operations
x.^2 % x^2 won't work here

% the number of elements of a vector
length(y)
% note: this is the number of elements, not the "Euclidean length"

% loop over a vector
for i=1:length(y)
    2*y(i)
end

% plotting vectors
plot(y)
hold on   % don't overwrite this
plot(z, 'r')  % plot additional line in red colour

whos  % look at your workspace
clear all  % clear your workspace
whos  % look at your (now empty) workspace
close all  % close your figures

% modulate time course of sound wave
load('handel')
t = 0:1000/Fs:(length(y)-1)*1000/Fs;
plot(t, y);

y2 = y;
x = exp((1:length(y))/(Fs*3));
x = x/max(x);
y2 = y.*x';
plot(t, y2);

sound(y2, Fs)

%% plotting functions of vectors

x = [-5:0.1:5]

% square every element of x
y = x.^2
plot(y)

% look at the x-axis of the plot - what's wrong?
% this is better
plot(x, y);
plot(x, y, 'x');

% a few more things
hold on
plot(x, (x-1).^2, 'r')
hold on
plot(x, x.^2+5, 'k')
plot(x, 2*(x-1).^2, 'g')

close all

% plot sine and cosine
sinx = sin(x);
cosx = cos(x);

plot(x, sinx, 'k', x, cosx, 'r')

% correlation
corr(sinx, cosx)
% produces lots of NaNs because dimension not appropriate, better:
corr(sinx', cosx')
% sine and cosine are orthogonal - why is this correlation not exactly 0?


%% angle between vectors
x = [0 1];
y = [1 1];

% try multiplying x and y in different ways (use transposes)

% what is the angle between these two vectors?

% cosine of angle between vectors
cosxy = x*y' / (norm(x)*norm(y));

% get the angle from cosine
acosxy = acos( cosxy )

% Why does is this not what we expected?
% tranform radians to degrees:
acosxy * (180/pi)

%% Functions for two vectors

% correlation between two vectors
% try
corr(x, y)
% and you'll get a strange result
% the input should be column vectors, so do it like this

corr(x', y')

% ttest
ttest(x', y')
% produces a suspiciously significant result...
% check how the ttest works:
help ttest

% this produces the p-value
[h, p] = ttest(x', y')

% if you want the t-value:
[h, p, ci, stats] = ttest(x', y')
stats.tstat  % "stats" is a structure, with multiple elements

% for simple measures
sum(x)
mean(x)
std(x)

% Händel again, look at effect of smoothing
load('handel')
y2 = smooth(y, 5);
sound(y2, Fs)
corr(y, y2);
y2 = smooth(y, 5);
corr(y, y2);

%% Matrices

% Define a matrix
x = [1 1 1; 2 2 2]

% Getting elements out of matrices
x(1,2)
x(:,1)
x(2,:)

% Getting matrix dimensions
size(x) % returns number of rows and columns, respectively
[m,n] = size(x)

% Matrix transpose
x
y = x'

% visualise a matrix
imagesc(x)
imagesc(y)

% adding/multiplying values etc.
x + 2
2 * x
x.^2

% Define another matrix
z = [1 2 3; 4 5 6]

% Subtract z from x
x - z

% Add z to x
x + z

% !!! Elementwise multiplication (NOT "matrix multiplication": "*")
x .* z

% Try... (matrix multiplication)
x * z
% doesn't work because dimensions don't match (we'll get to this later)

% Elementwise division (NOT "matrix division/inversion": "/")
x ./ z

% Try... (matrix inversion)
x / z

% if you cannot get enough...
y3 = repmat(y, 3, 1)
% try also:
y3 = repmat(y, 1, 3)

%% Multiplying matrix and vector
% Define a vector
v = [1 1 1]'

% Mutliply matrix x with vector v (and obtain sum within each row of x)
x
v
x * v

% matrix and vector dimensions must agree
v * x % does not work

% Modulating time course of sound (again)
load('handel')
y2 = y(1:2*Fs); % take fewer samples (or might be out of memory)
m = 1:2*Fs;     % modulate with linear function
m = m/max(m);   % normalise the maximum
M = diag(m);    % make it diagonal of matrix
y3 = M * y2;    % multiply matrix with sound vector
sound(y3, Fs);

%% interpretation as sum of column vectors
% x * v results in the sum of column vectors in x
x * v
% produces the same result as
1*x(:,1) + 1*x(:,2) + 1*x(:,3)

%% Matrix multiplication
% Multiply two matrices
y = x'
x * y

% Multiply two matrices
y * x

% !!! Note: x*y is not necessarily y*x!!!

% Define new matrix
z = [1 1; 2 2]

% Works:
y*z
z*x

% Doesn't work:
z*y
x*z



%% Make your own matrices

% Identity matrix
id = eye(3)

% Example:
x*id

id*x'

% diagonal matrix
di = eye(3)
di(1,1) = 1
di(2,2) = 2
di(3,3) = 3

% Multiplies columns by a factor:
x*di

% Multiplies rows by a factor
di*y

%% Some useful functions for matrices

% sum/average across columns of a matrix
% note: most Matlab function operate column-by-column
sum( x )
mean( x )

% sum across rows of a matrix
sum( x' )
mean( x' )

% get the diagonal of a (square) matrix
mydiag = diag(di)

% dimensions of a matrix, size()
[m,n] = size(x)

% dimension of a vector
diaglen = length(mydiag)

% create a matrix/vector with ones
ones(10)   % creates 10x10 matrix
ones(10,1) % creates 10x1 column vector

% create a matrix/vector with zeros
zeros(10)   % creates 10x10 matrix
zeros(1,10) % creates 1x10 row vector

% create matrix/vector with standard uniform random numbers in interval 0-1
rand(3)

% create matrix/vector with standard normal random numbers in interval 0-1
randn(3)

% correlation among columns of a matrix
R = rand(5,3)
% correlation among columns
corr(R)
% correlation among rows
corr(R')
% visualise a matrix
imagesc( corr(R') )
colorbar


% "reshaping" a matrix or a vector
x = 1:10
reshape(x,5,2)
reshape(x,2,5)

% Kronecker product ("concatenates" matrices/vectors)
kron(id,y)

%% example for "image processing"
img2 = double(img2); % convert to double format (so that matrix operations work)
X = diff( eye(300) );
img3 = X * img2;
image(img3);


%% EMEG example
load EMEGdata.mat
whos
data
data.evoked
data.evoked.epochs

% get only EEG data
mydat = data.evoked.epochs(307:376,:)
plot( mydat' )

% find bad channel
max(max(mydat))
max(max(mydat'))
[a,b] =max(max(mydat'))

% remove bad channel
mydat(63,:) = []
whos

plot( mydat' )

[m,n] = size( mydat );

% create "average reference operator"
avg_op = (eye(m)-ones(m)/m);

% this subtracts the mean across electrodes at each sample (column)
mydatref = (eye(69)-ones(69)/69)*mydat;

% this should look different:
plot( mydatref' )

% a different way of plotting 2D data
imagesc(mydatref)
colorbar

% plot correlation matrix among electrodes
imagesc(corr(mydatref'))
colorbar