%% Single on-resonance signal
% The physical quantity actually measured during an MR experiment is the
% current or voltage in the RF coil - the bulk magnetization excited in the
% coil as the sample precesses and, in doing so, induces a small voltage in 
% the receive coil.  
%
% In a uniform sample that is well shimmed, when no gradients are applied, 
% all the protons resonate very close to the resonant frequency.  The
% signal still decays after 50 or 100ms, however, because the protons
% slowly get out of phase with each other and when they are uniformly
% distriubted in phase, they cancel each other out and the signal goes to
% zero. (Phase is measured as the orientation of the bulk magnetization 
% vector projected on the transverse plane, which is the x-y plane, 
% perpendicular to the main magnetic field which is z by convention).  This
% dephasing, even in a perfect sample, is called relaxation and will be the
% topic of Lecture 4.

% For simplicity, we'll model the signal in the coil as an exponential
% decay of a single bulk magnetization vector that is "on-resonance" (at
% the reference frequency of our rotation frame).
dt = .001; % temporal resolution, in s
time = 0:dt:.200; % milliseconds
tau1 = .080; % milliseconds - the decay time constant
signal1 = exp(-time/tau1);
subplot(3,2,1); plot(time,signal1,'k','linewidth',2); axis([0 .2 0 1])
 xlabel('Time (ms)'); ylabel('Signal (uV)'); box off
 title('Time domain')
 set(gcf,'Color','w')
 
% Although this is a very simple experiment, what we want to know is the
% central resonance frequency of the sample, and the width of the peak.  To
% calculate this, we take the Fourier transform.
freq = -.5/dt:1/max(time):0.5/dt;
spectrum1 = fftshift(fft(fftshift(signal1)));
subplot(3,2,2); plot(freq,abs(spectrum1),'k','linewidth',2); axis tight
 xlabel('Frequency (Hz)'); ylabel('Amplitude (a.u.)'); box off
 title('Frequency domain')
 
%% Single poorly-shimmed signal

% If the spins in the sample dephased faster, the signal will decay faster.
% This is what happens when the shim is not perfect, and one side of the
% sample, for example, is at a higher frequency than the other.  For
% comparison, we'll plot the poorly-shimmed signal with the well-shimmed
% signal.
tau2 = .010; % milliseconds - the decay time constant
signal2 = exp(-time/tau2);
subplot(3,2,3); plot(time,signal1,'k','linewidth',1); hold on; 
  plot(time,signal2,'b','linewidth',2); hold off; axis([0 .2 0 1])
 xlabel('Time (ms)'); ylabel('Signal (uV)'); box off
 legend('well-shimmed','poorly-shimmed')
 
spectrum2 = fftshift(fft(fftshift(signal2)));
subplot(3,2,4); plot(freq,abs(spectrum1),'k','linewidth',1); hold on
  plot(freq,abs(3*spectrum2),'b','linewidth',2); hold off; axis tight
 xlabel('Frequency (Hz)'); ylabel('Amplitude (a.u.)'); box off
 
% Note that I cheated there and multiplied spectrum2 by a factor of 3 when
% I plotted it.  The integral of the frequency spectrum is always equal to
% the amplitude of the exponential decay - so since the second spectrum is
% broader, it is actually much shorter.  This makes sense - in the time domain,
% the initial signal is the same because the spins are always in phase right
% after excitation.  However, in the frequency domain, the signal is
% spread out over more frequencies, so there is less of it at the center
% frequency.

%% Off-resonance signals
% Finally, we're going to consider the possibility that my reference
% frequency is not the right frequency.  In this example, we will pick a
% reference frequency that is 50Hz off for my well-shimmed sample, and 1000Hz
% off from my poorly-shimmed sample. Now the signals will be both
% decaying and oscillating
off1 = 50; % Hz
signal1off = exp(-time/tau1).*cos(2*pi*time*off1);
off2 = 100;
signal2off = exp(-time/tau2).*cos(2*pi*time*off2);
subplot(3,2,5); plot(time,signal1off,'k','linewidth',1); hold on; 
  plot(time,signal2off,'b','linewidth',2); hold off; axis tight
 xlabel('Time (ms)'); ylabel('Signal (uV)');
  legend([num2str(off1) 'Hz off'],[num2str(off2) 'Hz off'])
  
% The Fourier transform is separately sensitive to the off-resonance
% shifts and linewidths of the two "samples."
spectrum1off = fftshift(fft(fftshift(signal1off)));
spectrum2off = fftshift(fft(fftshift(signal2off)));
subplot(3,2,6); plot(freq,abs(spectrum1off),'k','linewidth',1); hold on
  plot(freq,abs(3*spectrum2off),'b','linewidth',2); hold off; axis tight
 xlabel('Frequency (Hz)'); ylabel('Amplitude (a.u.)')
 
%% We've been ignoring so far the fact that the frequency axis is negative;
% now this is impossible to ignore because we suddenly have 2 peaks - a
% positive and a negative one - for each off-resonance frequency.  The
% simplest thing to do for now is not to plot the negative portions of the
% freuqency axis.
subplot(3,2,6); plot(0:1/max(time):.5/dt,abs(spectrum1off((1+.5*max(time)/dt):end)),'k','linewidth',1); hold on
  plot(0:1/max(time):.5/dt,abs(spectrum2off((1+.5*max(time)/dt):end)),'b','linewidth',2); hold off; axis tight
 xlabel('Frequency (Hz)'); ylabel('Amplitude (a.u.)')
 
% More information on the Fourier transform is provided as caoFFT.m (with the
% accompanying zebras.mat file). It is possible to have negative
% frequencies because of the way in which phase information is represented
% in the Fourier transform.

%% To turn in
%  If the signal did not decay at all, what would the line in
%  subplot(3,2,2) (the upper right plot) look like -- would it be fatter or
%  even skinnier?  Why?
