% ICE6 - Signal-to-noise ratio and contrast-to-noise ratio


%% This is unreasonably high temporal resolution ... really fake ... but I
%% wanted 1000 data points so I could get clean numbers
exptLength = 1000; stimStart = 330; stimEnd = 660;
% Make a constant signal with unit itensity
S = ones([1 exptLength]); 
S(stimStart:stimEnd) = S(stimStart:stimEnd) + C;
plot(S); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])

% Add white noise, SNR = 100;
SNR = 100;
% Generate random noise
noise =  rand([1 length(S)]);
% Scale it so the SNR will be 100
noise = (1/SNR)/std(noise(:))*noise;
% verify that it's uniformly distributed
hist(noise,20)
% noise spectrum should be flat
plot(abs(fftshift(fft(fftshift(noise))))); xlabel('freq'); ylabel('amp')
plot(abs(fftshift(fft(fftshift(noise - mean(noise(:))))))); xlabel('freq'); ylabel('amp')
% add noise to signal and plot
SplusN = S + noise - mean(noise(:));
plot(SplusN); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])
% calculate SNR
signal = mean(SplusN)
noise = std(SplusN)
SNR = signal/noise


%%%% Now do CNR
exptLength = 1000; stimStart = 330; stimEnd = 660;
S = ones([1 exptLength]); 
% Add a 1% modulation on top of signal
C = 0.01;
S(stimStart:stimEnd) = S(stimStart:stimEnd) + C;
plot(S); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])

% Add white noise, SNR = 100 again
SNR = 100;
noise =  rand([1 length(S)]);
noise = (1/SNR)/std(noise(:))*noise;
SplusN = S + noise - mean(noise(:));
plot(SplusN); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])

% calculate CNR, or contrast-to-noise ratio.  It should be about 1.
contrast = (mean(SplusN(stimStart:stimEnd)) - mean(SplusN([1:32 67:100])))/ ...
    mean(SplusN([1:(stimStart-1) (stimEnd+1):exptLength]))
noise = std(SplusN([1:(stimStart-1) (stimEnd+1):exptLength]))
CNR = contrast/noise

% repeat with Gaussian noise (same story), just for constant signal
exptLength = 1000; stimStart = 330; stimEnd = 660;
S = ones([1 exptLength]); 
% Add normally distributed white noise, with SNR = 100;
SNR = 100;
noise =  randn([1 length(S)]);
noise = (1/SNR)/std(noise(:))*noise;
hist(noise,20)
% noise spectrum should be flat
plot(abs(fftshift(fft(fftshift(noise))))); xlabel('freq'); ylabel('amp')
plot(abs(fftshift(fft(fftshift(noise - mean(noise(:))))))); xlabel('freq'); ylabel('amp')

SplusN = S + noise - mean(noise(:));
plot(SplusN); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])
signal = mean(SplusN)
noise = std(SplusN)
SNR = signal/noise

% repeat with correlated noise, again just added to constant signal
exptLength = 1000; stimStart = 330; stimEnd = 660;
S = ones([1 exptLength]); 
% Add white noise, SNR = 10;
SNR = 100;
noise =  randn([1 length(S)]);
noise = conv(noise,[0 1 1 1 0]); % convolution with a top-hat is FT multiplication by sinc
noise = noise(1:exptLength);
noise = (1/SNR)/std(noise(:))*noise;
hist(noise,20)
% noise spectrum should be flat
plot(abs(fftshift(fft(fftshift(noise - mean(noise(:))))))); xlabel('freq'); ylabel('amp')

SplusN = S + noise - mean(noise(:));
plot(SplusN); xlabel('time'); ylabel('signal intensity')
axis([0 exptLength 0 1.1])
signal = mean(SplusN)
noise = std(SplusN)
SNR = signal/noise
% SNR is the same, noise distribution looks the same, but averaging will
% have less of an effect on beating down the SNR.  Prove this?
