%%% HW #2

%%% Due: October 9, 2008


%% 1. Boxcars and Sincs

% The boxcar function (or rect function) is

% rect(t) = 0 if t < -1
%         = 1 if -1 < t < 1
%         = 0 if t > 1

% a. Sample the rect function on an interval -10 to 10 and
%    its Fourier transform with 1024 points. Plot both functions.

% HINT: Be sure to plot the reciprocal domain (say m^-1 or s^-1) correctly
%       with the correct axis.

x1=linspace(-10,10,1024);
S1=zeros(1,1024);
S1( (x1<1)&(x1>-1) ) = 1;

w1=linspace(-1024/10,1024/10,1024);
F1=fftshift(fft(ifftshift(S1)));    
%%% Note the use of ifftshift here - the rect is centered in the window and
%%% so the Fourier transform must be corrected for this.
plot(w1,real(F1),w1,imag(F1));

%%% Sometimes the resolution of the original plot is not sufficient to
%%% determine the FWHM automatically. In these cases, we can interpolate to
%%% generate the missing points.

w1i=linspace(-1024/10,1024/10,30000);
F1i = interp1(w1,F1,w1i);


% b. Sample the rect function on an interval -100 to 100 and its Fourier
%    transform with 1024 points. Plot both functions.

x2=linspace(-100,100,1024);
S2=zeros(1,1024);
S2( (x2<1)&(x2>-1) ) = 1;

w2=linspace(-1024/100,1024/100,1024);
F2=fftshift(fft(ifftshift(S2)));
plot(w2,real(F2),w2,imag(F2));

w2i=linspace(-1024/100,1024/100,30000);
F2i = interp1(w2,F2,w2i);


% c. The full width at half maximum (FWHM) is a measure of linewidth and is
% often defined by the maximum height divided by the width at half of the
% maximum. Measure the FWHM for both plots in (a) and (b). Discuss the
% results. What is Fourier transform of a rect function called?

%%% CORRECTION: The full width at half maximum is the width of the center
%%% lobe at half the maximum height. We could measure this for the real or
%%% imaginary parts or the absolute value.

%%% Locate points that are at half the height.
maxheight1=max(abs(F1i));
minheight1=min(abs(F1i));
points1=find( (F1i<((maxheight1-minheight1)/2+1))&(F1i>((maxheight1-minheight1)/2-1) ));
width1=w1i(points1(size(points1,2)))-w1i(points1(1));
plot(w1i,real(F1i)); hold on;
plot(w1i(points1(size(points1,2))),real(F1i(points1(1))),'ro');
plot(w1i(points1(1)),real(F1i(points1(1))),'ro');

maxheight2=max(F2i);
minheight2=min(F2i);
points2=find( (F2i<( (maxheight2-minheight2)/2+0.01))&(F2i>( (maxheight2-minheight2)/2-0.01) ));
width2=w2i(points2(size(points2,2)))-w2i(points2(1));
plot(w2i,real(F2i)); hold on;
plot(w2i(points2(size(points2,2))),real(F2i(points2(1))),'ro');
plot(w2i(points2(1)),real(F2i(points2(1))),'ro');

%%% The width of the sinc function increases as the width of the rect
%%% function decreases. The two are inversely related in frequency and time
%%% domains.

%% 2. Exponential Decay and Lorentzians

% (a) Plot one second of an exponential decay with T2 = 0.1 s oscillating
%     at a frequency 100 Hz. Take the Fourier transform and plot this also
%     with the correct axis.

T2=0.1;                 %s
v = 100;                %Hz
t=linspace(0,1,1000);   %s
S1=exp(t*(i*2*pi*v-1/T2));

w1=linspace(-500,500,1000);
F1=fftshift(fft(S1));
plot(w1,real(F1));

w1i=linspace(-150,150,30000);
F1i = interp1(w1,F1,w1i);

% (b) Over the same period, plot an exponential decay with T2=0.01 s
%     oscillating with a frequency 100 Hz. Take the Fourier transform and
%     plot this with the correct axis.

T2=0.01;                 %s
v = 100;                %Hz
t=linspace(0,1,1000);   %s
S2=exp(t*(i*2*pi*v-1/T2));

w2=linspace(-500,500,1000);
F2=fftshift(fft(S2));
plot(w2,real(F2));

w2i=linspace(-500,500,30000);
F2i = interp1(w2,F2,w2i);

% (c) Again, measure the FWHM of the lines in (a) and (b) and discuss the
%     results. What is the Fourier transform of an exponential decay
%     called?

maxheight1=max(F1i);
minheight1=min(F1i);
points1=find( (F1i<((maxheight1-minheight1)/2+0.1))&(F1i>((maxheight1-minheight1)/2-0.1) ));
width1=w1i(points1(size(points1,2)))-w1i(points1(1));
plot(w1i,real(F1i)); hold on;
plot(w1i(points1(size(points1,2))),real(F1i(points1(1))),'ro');
plot(w1i(points1(1)),real(F1i(points1(1))),'ro');


maxheight2=max(F2i);
minheight2=min(F2i);
points2=find( (F2i<( (maxheight2-minheight2)/2+0.01))&(F2i>( (maxheight2-minheight2)/2-0.01) ));
width2=w2i(points2(size(points2,2)))-w2i(points2(1));
plot(w2i,real(F2i)); hold on;
plot(w2i(points2(size(points2,2))),real(F2i(points2(1))),'ro');
plot(w2i(points2(1)),real(F2i(points2(1))),'ro');

%%% The shape is a Lorentzian. The relaxation times T2 which
%%% are fast (T2=0.01 s) have much broader linewidths than relaxation times which are small.
%%% T2 relaxation causes homogeneous broadening of the spectral line.


% (d) In this case, how is the spectral line width (FWHM) related to T2?
%     Give an expression for this.

FWHM = 1/pi/T2; %%% although writing this as a statement of proportionality is ok

%% 3. Parts-Per-Million and Alcohols

% Parts-Per-Million (ppm) is a useful independent variable for plotting NMR
% spectra because it is independent of field strength.

% ppm = (frequency - reference frequency) / reference frequency

% For example, with respect to the reference TMS (tetramethysilane), 1H in
% water resonates at 4.7 ppm.

% For example, to calculate the water shift from TMS at 9.4 T...

4.7                 %%% ppm
42.57               %%% gamma MHz/T
9.4                 %%% T
4.7 * 42.57 * 9.4   %%% Hz


% a. Simulate the time domain signal and 1H NMR spectra of 2-phenylethanol
%    C6H5-CH2-CH2-OH.

% %%% HINT:: 2-phenylethanol has 4 groups of protons which resonate at
% different frequencies. For this problem, assume that there is no
% interaction between the protons, so that each group resonates at only one
% frequency (no J-coupling). Also, neglect relaxation.

% %%% Phenyl 1H ~ 7.5 ppm
% %%% OH        ~ 4.5 ppm
% %%% CH2       ~ 4.0 ppm
% %%% CH2       ~ 3.1 ppm

%%% Convert ppm to frequency
factor=42.57*9.4;
v1=-7.5*factor; %%% negative because ppm are downfield (deshielded) from TMS
v2=-4.5*factor;
v3=-4.0*factor;
v4=-3.1*factor;

t=linspace(0,0.2,2048);
S=5*exp(i*2*pi*v1*t)+exp(i*2*pi*v2*t)+2*exp(i*2*pi*v3*t)+2*exp(i*2*pi*v4*t);
plot(t,real(S));

F=fftshift(fft(S));
w=linspace(-1024/0.2/factor,1024/0.2/factor,2048);
plot(w,abs(F));