% Demonstrate an example of stochastic resonance

% % Sampling rate
sr = 64; %Hz NB: if you don't choose a power of 2 here it gets messy

% Time basis
t0 = 0;
t1 = 1;
t = linspace(t0,t1,sr*(t1-t0));  % 1 second

% Frequency basis
nu = (-length(t)/2):(length(t)/2-1);

% Signal parameters
f = 4;
A = 2;
sig = A*sin(2*pi*f*t);

% Noise
nA = 3; %Noise amplitude
nsig = nA*2*(rand(size(t)));
nsig  = nsig - mean(nsig); % remove DC


% Threshold 
thresh = 2.2;

figure(1)
clf
subplot(3,1,1)
plot(t,sig)
title('Signal');
xlabel('t (s)');
ylabel('sig')

subplot(3,1,2)
plot(t,nsig)
title('Noise');
xlabel('t (s)');
ylabel('noise')

subplot(3,1,3)
sigtn = sig+nsig;
sigtn(find(sigtn<thresh))=thresh;
plot(t,sigtn)
title('Sig + Noise after threshold');
xlabel('t (s)');
ylabel('noise')

figure(2)
clf

% First DFT with no noise
SIGnn = fftshift(fft(sig));
subplot(3,1,1);
plot(nu,abs(SIGnn)/max(abs(SIGnn)))
title('No noise, no threshold');
xlabel('\nu (Hz)');
ylabel('FFT(sig)')

% Next DFT with no noise with threshold
sigt = sig;
sigt(find(sigt<thresh))=thresh;
sigt=sigt-thresh; % subtract DC
SIGnnt = fftshift(fft(sigt));
subplot(3,1,2);
plot(nu,abs(SIGnnt)/max(abs(SIGnnt)))
title('No noise, threshold');
xlabel('\nu (Hz)');
ylabel('FFT(sig)')

% Finally DFT with  noise with threshold
sigtn = sig+nsig;
sigtn(find(sigtn<thresh))=thresh;
sigtn=sigtn-mean(sigtn); % subtract DC
SIGnt = fftshift(fft(sigtn));
subplot(3,1,3);
plot(nu,abs(SIGnt)/max(abs(SIGnt)))
title('Noise, threshold');
xlabel('\nu (Hz)');
ylabel('FFT(sig)')

% Calculate SNR
[discard,Isig] = max(abs(SIGnn((length(t)/2):end))) % Just check for positive frequencies

hold on
line([nu(length(t)/2+Isig-1), nu(length(t)/2+Isig-1)],[0,1]);
hold off

SNR = abs(SIGnt((length(t)/2)+Isig-1))/mean(abs(SIGnt((length(t)/2):end)))

