echo off;
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology
%
% hw6p1.m - Homework #6, Problem #1
%
% Stochastic resonance using Integrate & Fire Neuron
% 
% UNITS: R in MOhm, C in nF, T in msec, V in mV, I in nA
%
% M. Nelson, 2004 
%----------------------------------------------------------------------
clear;                          % clear MATLAB variables
figure(1); clf;                 % clear MATLAB graphcis
orient tall;                    % use the full page

%%%%%%%%%%%%%%%%%%
% INITIALIZATION
%%%%%%%%%%%%%%%%%%

C = 1;   % nF
R = 10; % MOhm

VTHR = 5;  % mV
VSPK = 70; % mV

TMAX = 5000; %msec   
DT = 1;                         % integration time step
t=(0:DT:TMAX)';                 % time vector
I = 0.3*sin(2*pi*t/1000);      % current injection vector 
Inoise = randn(size(t));

AnoiseArray = [0 .1 .2 .5 1 2 5 10];

for j = 1:length(AnoiseArray)
    Anoise = AnoiseArray(j);
    
    v = zeros(length(t),1);         % pre-allocate voltage array
    
    %%%%%%%%%%%%%%%
    % MAIN LOOP
    %%%%%%%%%%%%%%%
    for i=2:length(t);                 
        
        if(v(i-1) == VSPK)                  % if we just spiked, 
            v(i) = 0;                       %   reset the voltage to 0
        else                                % otherwise,
            dvdt = (-v(i-1)/R + I(i) + Anoise*Inoise(i))/C;    %   calculate dv/dt using
            v(i) = v(i-1) + dvdt*DT;        %   Euler integration
        end
        
        if(v(i) >= VTHR)                    % check for threshold crossing
            v(i) = VSPK;                    % generate a spike
        end
        
    end   
    
    plot(t,100*(j-1)+v,'b-');                     % voltage vs. time
    hold on;
    str = sprintf('\\sigma_{noise} = %3.1f', Anoise);
    text(-1200,100*(j-1) + 20, str, 'color','r');
    
end

xlabel('time (msec)');
ylabel('voltage (mV)');
axis([-1500 TMAX -20 100*length(AnoiseArray)])    
title('Stochastic Resonance Example');
grid on;