function [sboth] = makeTernaryStim(T,X,pixelRate,pixelWidth,dtinstance,parity,spatialFilterType)
% This function simulates the input to two photoreceptors that are
% observing a correlated noise stimulus.
%
% Inputs:
% T: The amount of time in seconds to simulate
% X: The amount of space in degrees to simulate
% pixelRate: The number of frames per second
% pixelWidth: The size of the pixels in degrees
% dtinstance: The time interval in frames on which correlations will exist
% dtinstance=-1 simulates no correlations
% parity: 1 for positive correlations, -1 for negative correlations
% spatialFilterType: 0 for delta function spatial filters, 1 for 5 degree
% fwhm gaussian filters
%
% Output: the simulated inputs in a matrix with dimensions
% time (in milliseconds) by cell (left and right)

    % A correlated noise stimulus is generated by averaging a random binary
    % stimulus with a version of itself which has been shifted in time and
    % space
    c1 = (rand(ceil(T*pixelRate/1000)+2,ceil(X/pixelWidth)+2)>0.5)*2 - 1;

    % this simulates no correlations in the stimulus
    if dtinstance == -1
        c2 = (rand(ceil(T*pixelRate/1000)+2,ceil(X/pixelWidth)+2)>0.5)*2 - 1;
        xttemp = c1/2 + c2/2;
    else
        xttemp = c1/2 + parity*circshift(c1,[dtinstance,1])/2;
    end

    if X/pixelWidth == 2 && spatialFilterType == 0
        sboth = xttemp(:,2:3);
    else
        % interpolate in space to simulate pixel width
        xt = interp1([0:size(xttemp,2)-1]*pixelWidth,xttemp',[0:X-1]','previous')';
        % Get filter weights and apply them
        filts = make_spatial_filters(spatialFilterType,X);
        sboth = xt*filts;
    end

    % interpolate in time to simulate the duration of a frame
    sboth = interp1([0:size(sboth,1)-1]/pixelRate*1000,sboth,[0:T-1]','previous');
end

function F = make_spatial_filters(spatialFilterType,X)

    dX = 1;
    centers = [12.5 17.5]; %Location of the centers of the receptive fields in degrees

    fwhm = 5;
    sigma = fwhm/sqrt(8*log(2)); % convert to sigma
    xvals = [0.5:dX:X]'; % vector of pixel center locations in degrees

    switch spatialFilterType
        case 0 % delta functions
            F = zeros(length(xvals),2);
            F(round(centers(1)),1)=1;
            F(round(centers(2)),2)=1;
        case 1 % Gaussian filters
            % note: weights are computed for the center of the pixel, not 
            % integrated over the pixel. This approximation is less correct
            % when fwhm is small.
            f1 = exp(-(xvals-centers(1)).^2/2/sigma^2);
            f2 = exp(-(xvals-centers(2)).^2/2/sigma^2);
            f1 = f1/sum(f1);
            f2 = f2/sum(f2);
            F = [f1(:),f2(:)];
    end
end