function [STA, t] = stim_sta(T, stim, bounds,dt)
%STIM_STA Spike triggered average of stimulus waveform

% Initialization
t = bounds(1):dt:bounds(2);
STA = zeros(size(t,2),1);

% Get windows
bounds = bounds(ones(size(T)),:);
T = [T, T];
bound_n = ceil((bounds + T)/dt);
K = 0;

% For all spikes occuring in at times T
for i = 1:size(bound_n,1)
    
    % Current window
    n = bound_n(i,1):1:bound_n(i,2);
    
    try
        % Optical power
        STA = STA + stim(n);
        K = K + 1;
        
    catch
        continue
    end
    
end

% Normalize by number of spikes
STA = STA./K;

end

