% Show an SNR curve for the stochastic resonance example


% % 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 = logspace(-1,2,400); %Noise amplitude
Navg = 20; % Number of averages over noise

% Threshold 
thresh = 2.05;

for jj = 1:length(nA)
    
    for kk=1:Navg

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

        % DFT of signal + noise with threshold
        sigtn = sig+nsig;
        sigtn(find(sigtn<thresh))=thresh;
        sigtn=sigtn-mean(sigtn); % subtract DC
        SIGnt = fftshift(fft(sigtn));

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

        SNRtemp(kk) = abs(SIGnt((length(t)/2)+Isig-1))/mean(abs(SIGnt((length(t)/2):end)));
    end
    SNR(jj) = mean(SNRtemp);
    disp(jj)
end

figure(1)
clf
semilogx(nA,10*log10(SNR))
xlabel('Noise amplitude')
ylabel('SNR (dB)')
title('SNR as a function of noise')
axes('position',[0.35,0.25,0.3,0.2])
plot(t,sig)
hold on
plot(t,thresh*ones(size(t)),'g-')
hold off
xlabel('t(s)')
ylabel('sig');