% Module 4:  SNR vs. resolution
%
% Working with new SNR flag in (fixed) EPI code.
%
% Download the most recent version of EPI.m and make sure it's in your
% path.  
% The pineapple's inhomogeneous enough that it's a pain to do SNR stuff.
% Same's true for human brains, so this week it's a boring sphere ...
load sphere

% Problem 1.  SNR vs. bandwidth
% We'll simulate a "normal" 3T EPI experiment
TR = 1500; TE = 30;
FOV = .064; 
BW = 100000;

% We'll look at the effect of resolution and field of view, starting with:
FOV = .064;
[hiBW kData] = EPI(FOV,16,pd,T1,T2,TR,TE,[],[],[],BW,1);

% Three ways to estimate SNR. 
% 1) Take the average intensity out on the side, away from ghosts in the PE
% direction
noiseEst = 1.65*std(abs(hiBW(:,1)));
% 2) Take the standard deviation of a "uniform" image region (but this is
% bogus because the variation of intensity in our image is due to Gibbs
% ringing, a truncation artifact that shows up at strong intensity
% boundaries.
noiseEst = std(reshape(abs(hiBW(7:11,7:11)),[1 25]));
% 3) If you're going after thermal SNR, you can turn off the RF pulse and
% run the experiment the same way:
alpha = 0;
[noise kData] = EPI(FOV,16,pd,T1,T2,TR,TE,alpha,[],[],BW,1);
noiseEst = 1.65*std(abs(noise(:)));

% I'm going with the last ... (since, secretly, nothing worked out right
% using the first 2, which I suspect is a monstrous bug in how I
% conceptualized the code ...)
signal = mean(mean(abs(hiBW(7:11,7:11))));
SNR_hiBW = signal/noiseEst

% Repeat for low bandwidth
BW = 50000;
[loBW kData] = EPI(FOV,16,pd,T1,T2,TR,TE,[],[],[],BW,1);
[noise kData] = EPI(FOV,16,pd,T1,T2,TR,TE,0,[],[],BW,1);
noiseEst = 1.65*std(abs(noise(:)));
signal = mean(mean(abs(loBW(7:11,7:11))));
SNR_loBW = signal/noiseEst

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Explain why the SNR in the high BW case is approximately 1.4X worse
% than in the lowBW case, for fixed resolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Problem 2.  Effect of FOV (keeping resolution the same).
FOV = .064;
[lowRes kData] = EPI(FOV,16,pd,T1,T2,TR,TE,[],[],[],[],1);
noise = EPI(FOV,16,pd,T1,T2,TR,TE,0,[],[],[],1);
noiseEst = 1.65*std(abs(noise(:)));
signal = mean(mean(abs(lowRes(7:11,7:11))));
SNR_smallFOVloRes = signal/noiseEst

FOV = .128;
[lowRes kData] = EPI(FOV,32,pd,T1,T2,TR,TE,[],[],[],[],1);
noise = EPI(FOV,32,pd,T1,T2,TR,TE,0,[],[],[],1);
noiseEst = 1.65*std(abs(noise(:)));
signal = mean(mean(abs(lowRes(15:19,15:19))));
SNR_lgFOVloRes = signal/noiseEst

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Explain why the SNR in the large FOV case is approximately 2X better
% than in the small FOV case, for fixed resolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Problem 3.  Resolution
FOV = .064;
[lowRes kData] = EPI(FOV,16,pd,T1,T2,TR,TE,[],[],[],[],1);
noise = EPI(FOV,16,pd,T1,T2,TR,TE,0,[],[],[],1);
noiseEst_lowRes = mean(abs(noise(:)));
signal_lowRes = mean(mean(abs(lowRes(7:11,7:11))));
SNR_lowRes = signal_lowRes/noiseEst_lowRes

[hiRes kData] = EPI(FOV,32,pd,T1,T2,TR,TE,[],[],[],[],1);
noise = EPI(FOV,32,pd,T1,T2,TR,TE,0,[],[],[],1);
noiseEst_hiRes = mean(abs(noise(:)));
signal_hiRes = mean(mean(abs(hiRes(15:19,15:19))));
SNR_hiRes = signal_hiRes/noiseEst_hiRes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Again, explain the relative SNRs, low res and high res (this time, the
% high resolution SNR is about half the low-res).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%