%% CLASS NOTES

%%%% Introduction to Matlab Techniques for NMR
%%%% Date: September 16, 2008

%% User Interface

%% Command Window

% The smallest, but most significant part of Matlab. Matlab executes
% commands like a calculator.

% EXAMPLES:::

% (1) Add two numbers together and print them.

% (2) Add two numbers together and store them as a variable 'x'

% (3) Repeat a command with the up arrow, but add a semicolon to the end of
% the command. What happens?

% (4) Repeat a command with the up arrow, but add a percent sign to the
% beginning of the command.

% (5) Calculate the Larmor frequency for 1H in a 3T magnetic field.

gamma = 42.57     % MHz/T
B0    = 3         % T
v0    = gamma*B0  % MHz

% (6) Calculate the Larmor frequency for 23Na in a 7T magnetic
% field.

gamma = 11.2      % MHz/T
B0    = 7         % T
v0    = gamma*B0  % MHz

%% Editor

% The editor is where you create programs or scripts.

% To create a new document to edit, you can use either the shortcut or
% choose File->New->M-File

% EXAMPLES:::
% Create an M-File script which adds two numbers together and then divides
% the result. Save the M-File script when you are finished.

%% Current Directory

% The current directory is a directory browser, where you can browse the
% hard disc.

%% Workspace

% In the workspace window, a list of variables in memory is displayed.

% EXAMPLE:::
% (1) Create a variable and determine how many bytes it consumes in memory.
% (2) Save your entire workspace, clear memory, and then recall your saved
% workspace from the hard disc.

%% Command History

% This contains a brief and interactive list of previous commands that were
% entered.

%% Other features

% (1) Current directory browser in top window
% (2) Docking and undocking windows
% (3) Data browser, much like a spreadsheet
% (4) Creating a variable of another type
% (5) Matlab Toolkits (statistics, symbolic, etc)


%% Help

% At the top of the GUI, select Help->Product Help

% EXAMPLE::: Permute
% Search for the function 'permute' in the help manual.

%% Matrices

%%% Creating matrices

x = [1 2; 3 4]; %% Create a 2x2 matrix
y = [1,2; 3,4]; %% Create a 2x2 matrix

z = [1; 2; 3; 4]; %% Create a column vector (1x2 matrix)
q = [1 2 3 4]; %% Create a row vector 2x1 matrix

%%% Indexing
%%% indexing makes it possible to change individual elements of a matrix

x(1,1)=3; %% Change element (1,1) to 3.
x(2,1)=6; %% Change element (2,1) to 6.

%%% Preallocate vectors

t1=0:1:100 %%% create a vector going from 0 to 100 in steps of 1.
t2=1:1:100 %%% create a vector going from 1 to 100 in steps of 1.
t3=1:5:100 %%% Create a vector going from 1 to 100 in steps of 5.

%%% Fencepost problems
%%% If you make a fence with 10 panels, how many posts do you need? 11!

size(t1)   %%% t1 from the example above has 101 elements, but I want 100.
t1=linspace(0,1,100);   %%% linspace solves the fencepost problem

%%% Preallocate matrices

w=zeros(2,2);
x=ones(2,2);
y=rand(10000,1);

%%% Create matrices for 2D functions

[X,Y]=meshgrid(1:1:100,1:1:100);


%% Loops

%%% Loops enable you to repeat a single computation multiple times.

for j=1:10
    k(j)=j*10;
end

%% Vector Manipulation

%%% If possible, don't use loops in Matlab. Find another way to code.

%%% Example:::

%%% hint: tic, toc is a stopwatch
%%% tic <- start watch
%%% do something...
%%% toc <- stop watch and display time

%%% instead of...
tic
for j=1:100
    t(j)=j*2;
    S(j)=exp(-t(j)/100);
end
z1=toc
%%% try this...
tic
t=2:2:200;
S=exp(-t/100);
z2=toc

z1/z2 %%% how much faster?

%%% But what if the vector was really big?

%%% instead of...
tic
for j=1:10000
    t(j)=j*2;
    S(j)=exp(-t(j)/100);
end
z1=toc
%%% try this...
tic
t=2:2:20000;
S=exp(-t/100);
z2=toc

z1/z2 %%% how much faster?

%%% Some other useful functions

%%% Circshift
A = [ 1 2 3;4 5 6; 7 8 9]
A =
     1     2     3
     4     5     6
     7     8     9

B = circshift(A,1)
B =
     7     8     9
     1     2     3
     4     5     6
     
%%% Flip and Rotate
A = [ 1 2 3;4 5 6; 7 8 9]
B = flipud(A);
C = fliplr(A);
D = flipdim(A,2);
E = rot90(A);

%%% Reshape
x=1:10:1e5;
size(x,2)/10
x=reshape(x,10,1000);


%%% Sort
y=rand(100,1);
plot(y); hold on;   %%% hint::: hold on keeps a plot up
plot(sort(y),'r');  %%% hint::: 'r' specifies plot color 'red'

%%% Vector operations
x=rand(100,1);
y=rand(100,1);
c=10;
c*x   %%% scalar operations are the same
x+y   %%% so are vector addition...
x-y   %%% and subtraction
x./y  %%% but, you must use period operator for some others
x.*y
x.^y


%% Plotting a Line and Figure Manipulation

m = 0.1       % slope
b = 2        % intercept
x =-100:100  % x
y = m*x+b    % y
plot(x,y)

% To edit your figure, within the figure menu, select View->Property
% Editor

%% Random Numbers

%%% TIP::: This is very useful for modeling noise...

z=randn(10000,1);
hist(z)
hist(z,100)

%%% Shift the mean and standard deviation. 
 z =   0.6  +  2 *   randn(10000,1);
%%%      mean  + std * random numbers


%% Complex Numbers

%%% Complex numbers are represented by two values.
%%% a               + i*b          % i = sqrt(-1)
%%% real part       + imaginary part

%%% real numbers
a = 2
%%% imaginary numbers
b = 2i    %%% i = sqrt(-1)

%%% Euler's Equation
phi = 45/360*2*pi   %%% radians
exp(i*phi)          %%% is the same as...
cos(phi)+i*sin(phi) %%% this...

%%% Representing signals with complex numbers
t = 0:1:100  %%% s
v = 100      %%% Hz
S_real = cos(2*pi*v*t) %%%  real part of complex signal
plot(t,S_real);
S_imag = sin(2*pi*v*t) %%% imaginary part of complex signal
plot(t, S_imag);
S = cos(v*t)+i*sin(2*pi*v*t)
plot(t,S)
plot(t,real(S))
plot(t,imag(S))
plot(S)
S_euler = exp(i*2*pi*v*t)
plot(t,S);       hold on;   %%% is the same as....
plot(t,S_euler);            %%% this...


%% Fourier Transform

%%% A special operation which takes a function from one space to another

%%% In NMR, the FT is used ALL THE TIME!
%%% S(t) -> S(v)   Forward Transform (time to frequency space)
%%% S(v) -> S(t)   Inverse Transform (frequency to time space)

%%% The FT is calculated using a special computational algorithm called the
%%% fast Fourier Transform (or 'FFT'). FFT and FT are sometimes used
%%% interchangeably although FFT refers to the computer algorithm and FT
%%% to the operation.

%%% Later, we will have a whole class on the FT.

%%% Example:::

%%% Suppose I have a signal with multiple frequencies and amplitudes
tmin=0;
tinc=.0005;      % tinc is very important to prevent aliasing
tmax=0.256;
t = tmin:tinc:tmax;         % s
v1 = 750; v2 = -400; v3 = 1; %%% frequencies (Hz) %%% Hint:: multiple lines
a1 = 1;   a2 = 2;   a3 = 3 ; %%% amplitudes
S = a1*exp(2*pi*i*v1*t)+a2*exp(2*pi*i*v2*t)+a3*exp(2*pi*i*v3*t);
plot(t,S);

F = fftshift(fft(S)); %%% Take the FFT and shift (FFTSHIFT) the data N/2.
v = linspace(-1/tinc/2,1/tinc/2,size(t,2));
plot(v,abs(F));
xlabel('Hz');

% try this
subplot(2,1,1); plot(t,S);
subplot(2,1,2); plot(v,abs(F));

%% Find

%%% find() returns the elements of a vector satisfying some condition

%%% repeat the histogram of a normal random distribution
z=randn(100000,1);
x = -4:0.01:4;
hist(z,x);

%%% now find the points lying outside some range
p = find(z < 0.2 | z > 0.4);
whos p;

%%% see that it worked
hist(z(p),x);

%%% FIND example #2:
%%% find() can also work with arrays
x = -100.:1.:100.;
y = -100.:1.:100.;
[x,y] = ndgrid(x,y); 

z = sqrt(x.*x+y.*y);

colormap hot;
imagesc(z);
colorbar;

p = find(z > 100.)
z(p) = 0.
imagesc(z);
colorbar;
%%% also mesh
mesh(z)


%% Rotation Matrices

%%% a special class of 3x3 matrices which rotate vectors in x,y,z

theta=45;              % degrees
theta=45/360*2*pi;     % radians
zvect=[0; 0; 1];

%%% rotation around the x-axis
Rx = [1 0 0; 0 cos(theta) sin(theta); 0 -sin(theta) cos(theta)];
%%% rotation around the y-axis
Ry = [cos(theta) 0 -sin(theta); 0 1 0; sin(theta) 0 cos(theta)];
%%% rotation around the z-axis
Rz = [cos(theta) sin(theta) 0; -sin(theta) cos(theta) 0; 0 0 1];

plot3(zvect(1),zvect(2),zvect(3),'ro'); hold on;
axis([-1 1 -1 1 -1 1]);
ax=linspace(-1,1,30);
line=zeros(1,30);
plot3(ax,line,line,'-.');
text(1,0,0,...
     'x',...
     'FontSize',16);
 text(0,1,0,...
     'y',...
     'FontSize',16);
 text(0,0,1,...
     'z',...
     'FontSize',16);
plot3(line,ax,line,'-.');
plot3(line,line,ax,'-.');

y=Rx*zvect;
plot3(y(1),y(2),y(3),'ro'); hold on;