% A SIMPLE MATLAB TUTORIAL
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% @ George Kafentzis, SigProcLab, CSD, UOC.
% e-mail : kafentz@csd.uoc.gr

% MATLAB is a convenient software for mathematical (and not only) purposes.
% We will see here its most basic features.
% -------------------------------------------------------------------------




% Notes : 
% -------
% 1)To see the result of each statement, just supress the ; at the end of each line.
% 2)To clear your workspace from all variables, type 'clear all'.





% VARIABLES
% -------------------------------------------------------------------------
% MATLAB handles !every! variable as a matrix, for example:

A = 2; % This is a 1x1 matrix with value 2.
B = 3; % This is a 1x1 matrix with value 3.

% The basic operations of addition, subtraction, multiplication and
% division are the same as in C, as long as the rules of dimensionality are
% applied! The following is quite straightforward, because our variables are
% 1x1 matrices :

C = A + B; % C = 5;
D = A - B; % D = -1;
E = A * B; % E = 6;
F = A / B; % F = 3/2 - careful! It is /, not \ !

% Of course, you can do more than this : C = (A+B)*D - F*C;
% The precedence of operations are the same as in C. Simple, isn't it?





% VECTORS AND MATRICES
% -------------------------------------------------------------------------
% It is very easy in MATLAB to declare and initialize a vector. For
% example:

v = [1 2 3 4];      % This is a 1x4 matrix, which is a row-vector.
v_tr = [1 2 3 4].'; % This is the transpose of row-vector v, a column-vector.
                    % Note the .' notation for transposing a matrix.
x = [4; 1; 0; 2];   % This is a 4x1 matix, which is in fact a column-vector.

% A column vector and a row vector are different in MATLAB (whereas in C
% are the same - x = malloc(3*sizeof(double)); makes an vector of 3 doubles,
% no matter if it is a row-vector or column-vector. By the way, memory
% allocation is automatic in MATLAB, so you don't need 'malloc' or 'free'
% functions.

% Let's make a stop here and say a word about matrix indexing. All
% matrices in MATLAB begin with index 1. For vector x above, this means 
% that x(0) is not defined. x(1) is 4, x(2) is 1, etc. Type 'x(0)', 'x(1)'
% etc. to see what happens.

% Now, check the following operations:

%V1 = v + v_tr;   % This is invalid in MATLAB! You can't add a row-vector 
                 % with a column-vector (or a 1x4 matix with a 4x1 matrix)! 
                 % Write this line to see what happens.

V2 = v.' + v_tr; % This is valid. Do you understand why? (check dimensionality)

% The known rules of matrix addition, multiplication and subtraction in
% matrices must be followed in MATLAB, or else the program will complain...
% Addition and subtraction are done element-by-element, as we all know.
% BUT for multiplication and division, things are different! We all know
% that matrix multiplication follows certain rules. We can have :
%
% -> The 'classic' matrix multiplication, as we all know from Linear
% Algebra :

A = [1 2; 2 1]; % A 2x2 matrix, with first line : 1 2, and second line : 2 1
B = [3 4; 1 2]; % Another 2x2 matrix, first line : 3 4, second line : 1 2

C1 = A*B;        % This is a 2x2 matrix, which is the product of A*B

% -> The element-by-element matrix multiplication :

C2 = A.*B;       % This is different! (Check the dot '.' before the *)
                 % This is an element-by-element multiplication
                 % of the two matrices, which is a 2x2 matrix again, but
                 % with different elements than C1! Check the results.
                 
% The same applies for division :
% -> The 'classic' matrix division, as we know it from Linear Algebra :

C3_1 = A/B;       % This is the same as : A*inv(B);
C3_2 = A\B;       % This is the same as : inv(A)*B;

% -> The element-by-element matrix division, which divides each element of
% A with its respective element in B!

C4 = A./B;      % This is different than C3! Can you understand the difference?

% What is the conclusion here? If we want an element-by-element
% multiplication or division, we should put a '.' before the operand * or /

% Now that we have seen a few things about matrices, let's say a word about
% indexing for non-vector matrices. For matrix A, for example. 
% Let's say we want the first element of the first column : it's 'A(1,1)'.
a = A(1,1);
% Let's say we want the second element of the first column : it's 'A(2,1)'. 
b = A(2,1);
% Let's say we want the third element of the second column : doesn't exist!
% Let's say we want the second element of the second column : it's A(2,2).
c = A(2,2);
% Now let's say we want the first column of the matrix. This should be a
% column vector, right? This is how it works :
v = A(:,1);     % This can be interpreted as "get the elements of all rows
                % of the first column", which is in fact the first column.
% In the same way, we can get the first row of the vector :
vv = A(1,:);    % "get the elements of all columns of the first row" ==
                % == first row.

% You can imagine how it works for an NxM matrix...






% LINEAR SYSTEMS
% -------------------------------------------------------------------------
% As you may have thought, MATLAB requires a good level of understanding of
% Linear Algebra... :-( let's say that we want to solve a linear system, Ax=b.

A = [1 2 1; 3 1 1; 1 1 2];   % This is a 3x3 matrix.
b = [1 1 1].';               % This is the column-vector b. (why is it a column?)
                             % Another way to write it is : b = [1;1;1];

% We want to find the solution of the linear system Ax=b => x = inv(A)*b;
% One way is to use the function inv().

x1 = inv(A)*b;

% Another way is to use MATLAB notation.

x2 = A\b; % \ means left matrix division. Type 'help mldivide' for more details.

% Check x1 and x2 to see that the results are the same. If you notice a
% small difference in the 10th or 12th decimal digit, it is because of the 
% different way x1 and x2 are calculated each time.





% COMPLEX NUMBERS
% -------------------------------------------------------------------------
% This course deals a lot with complex numbers. Let's see how to
% represent them in MATLAB.

% The complex number 1+3j can be declared as :
z = 1 + 3j; % or
z = 1 + 3*j % or
z = 1 + 3i; % or
z = 1 + 3*i;

% All the above declarations are equally valid. j and i are the imaginary 
% unit in MATLAB. Be careful if you are using j or i as an index in your code!

% We can easily extract the real or the imaginary part of a complex number
% z :

r = real(z);    % The real part of z, 1.
im = imag(z);   % The imaginary part of z, 3.

% As we have seen in the lectures, a complex number z can be
% represented in exponential form: z = |z| * exp(jö). Let's see this :

z1 = 4 + 5j;

abs_val = abs(z1);                    % The absolute value of z1, which is |z1|
ang = atan2(imag(z1), real(z1));    % The angle ö of z1, ö = arctan(Im{z1}/Re{z1})
ang = angle(z1);	              % Or using this instead.

% So, we can write :
z2 = abs_val * exp(j * angle);        % Check that z2 = z1

% Addition, subtraction, multiplication and division of complex numbers are
% trivial, using +, -, *, and /.
z3 = z1 + z2;
z4 = z1 - z2;
z5 = z1 * z2;
z6 = z1 / z2;






% (REAL) FUNCTIONS AND PLOTS
% -------------------------------------------------------------------------
% In MATLAB, it is easy to create a function and plot it. The procedure is
% more or less the same for every function :

% We first declare the x-axis (whether it is time or anything else)
t = 0:0.01:4*pi;    % pi is the ð value, ð = 3.1415...

% Let's explain the above statement a little. This notation says : 
% "Create a vector t. Its first value will be zero, and with an
% increasing step of 0.01, stop at 4*pi (including it)". 
% This means that vector t has the values 0, 0.01, 0.02, ... , 4ð. With
% this notation, t is a line-vector. (How can we make it a column-vector?)
% We can choose any incrementing step, positive or negative. So, this is
% our x-axis (time, in this case).

% After that, we define our y-axis, our function. Let's say it is a cosine.
f = 2*cos(t);

% Notice in your workspace that the length of f is the same as the length of t. 
% It makes sense because the cos() function is evaluated at each element (point) 
% of t. Of course, f could be any function, like f = exp(t), f = 2*t + 3,
% f = t.^2 + 2*t - 3, etc.

% Now that we have our function, we want to plot it and see its shape on
% the x-y plane. This is very easy :

plot(t,f); grid;  % Notice that the x-axis goes first, then the y-axis!!
                  % The 'grid' statement just puts the horizontal and
                  % vertical dashed lines on the plane
                  
title('Function : 2*cos(t)'); 
xlabel('time t (sec)'); 
ylabel('2cos(t)');
% If you take a good look at the plot, you will easily understand what
% these three statements do.
                  
% So, theoretically, we expect to see two periods of the cosine 
% (its period is 2ð) and this is exactly what we get.

% Let's see another one, let's say f(t) = 2*(t.^2) - 1;
t = -4:0.01:4;     % We choose this interval
f = 2*(t.^2) - 1;  % Our new function f(t)

% Notice that we can't write 't^2' because this means 't*t', which is
% invalid, because we multiply a 1xN vector with itself, that means a 1xN *
% 1xN multiplication, which is undefined! That's the reason we use the '.'
% before the '^', to tell MATLAB to multiply element-by-element, which is
% what we want : the square of each element of t.

% Let's plot it (it is the same procedure as above):
plot(t,f); grid; title('Function : 2(t^2) - 1');
xlabel('time t (sec)'); ylabel('2(t^2) - 1');


% -----------------------------
% |        PRACTICE: :-)      |
% -----------------------------
% Plot the functions below, each one with step = 0.1 : 
% (just like we did in the examples above)
%
% 1) f(x) = 3*exp(2*x); in [-10,10]
% 2) f(x) = 2*x + x.^2; in [0, 5]
% 3) f(x) = 2*cos(x) - 3*sin(3*x + pi/3) + cos(2*x + pi/8); in [-4ð, 4ð]
% 4) f(x) = (x.^2 - 3)./(x - 1); in [-1,4] (what happens here and why?)
%






% COMPLEX FUNCTIONS AND PLOTS
% -------------------------------------------------------------------------
% In about the same way, we can create complex functions. The plotting
% though is different. We cannot plot the complex function. But what we can
% do is plot its real part and its imaginary part, seperately. Let's see
% how we can do this : 

% Let's define the real part of the function (the same way as we did
% before).
re = -5:0.01:5;

% Let's do the same for the imaginary part.
im = -5:0.01:5;

% As we said, we can't plot the function f(z). So we will plot the real and
% the imaginary part. We need to use 'meshgrid' function in order to plot
% them. Check 'help meshgrid' for further details.
[R,I] = meshgrid(re, im);

% We define our function using R and I, let's say its f(z) = 2*(z.^2) - 1
Z = R + j*I;
f_Z = 2.*(Z.^2) - 1;

% Now we can use 'mesh' to plot the real and the imaginary part. Type 'help
% mesh' for more details.
figure; mesh(R,I,real(f_Z))
figure; mesh(R,I,imag(f_Z))

% You can see now the real part of f(z). Compare it with the function f(x) =
% 2t.^2 - 1 that we plotted earlier. How can we get f(x) from f(z)?
% (theoretically)

% -----------------------------
% |        PRACTICE: :-)      |
% -----------------------------
% Plot the real and the imaginary part of the complex functions below, 
% each one with step = 0.1, for both real and imaginary part:
% (just like we did in the examples above)
%
% 1) f(z) = 3*exp(2*z); in [-10,10] for both axis
% 2) f(z) = 2*z + z.^2; in [0, 5] for both axis
% 3) f(z) = 2*cos(z) - 3*sin(3*z + pi/3) + cos(2*z + pi/8); in [-4ð, 4ð]
% for both axis
% 4) f(z) = (z.^2 - 3)./(z - 1); in [-1,4] for both axis (what happens here and why? Use your imagination! :-) )
