Photocurrent Code

%% This script sets up a gaussian grating that is constant over the simulation, % then loops over an array of spatial frequencies to find various results % such as photocurrent spectra and amplification. Parameters can be % changed to test different factors for eye movements.

ieInit;

% Initialize parameters

% Gaussian onset and offset of the grating stimWeights = ieScale(fspecial('gaussian',[1,50],15),0,1); % Padded by zeroes weights = [zeros(1, 30), stimWeights, zeros(1, 30)];

% This is the field of view of the scene. sparams.fov = 0.5;

% Initialize the harmonic parameters structure with default % Change entries that are common to uniform and harmonic freq = linspace(1,50,5); ratios = zeros(3,1); fft_bin2_vals = zeros(10,1);

figure; hold on; for a = 1:length(freq)

   clear params
   for ii=2:-1:1
       params(ii) = harmonicP;
       params(ii).GaborFlag = 0.2;
       params(ii).freq      = freq(a);
       params(ii).row = 256;
       params(ii).col = 256;
   end
   
   % params(1) is for the uniform field
   params(1).contrast  = 1.0;  % contrast of the two frequencies
   
   % params(2) is matched and describes the grating
   params(2).contrast = 1.0;
   
   % The call to create the retinal image sequence
   oisH = oisCreate('harmonic','blend',weights,'testParameters',params,'sceneParameters',sparams);
   fov = oiGet(oisH.oiFixed,'fov');
   tSamples = 128;
   
   cMosaic = coneMosaic;
   cMosaic.integrationTime = 0.001;
   cMosaic.setSizeToFOV(fov);
   cMosaic.os.noiseFlag = 'none';
   
   % create em object without movement
   em_noMovement = emCreate;     % Create an eye movement object
   em_noMovement.emFlag = [0 0 0];  % Make sure tremor, draft and saccade are all off
   cMosaic.emGenSequence(tSamples,'em',em_noMovement);  % Generate the sequence
   
   cMosaic.compute(oisH);
   cMosaic.computeCurrent;
   
   deMeanedMosaic = bsxfun(@minus,cMosaic.current,mean(cMosaic.current,3));

   Spectra_noMovements = abs(fft(deMeanedMosaic,128,3));
   avgSpectrum_noMovements = squeeze(sum(sum(Spectra_noMovements)))./93^2;
   
   % create em object with movement
   emF = 1; emA = 93/30;
   x = round(emA*sin(2*pi*emF*(1:tSamples)/tSamples));
   y = zeros(size(x(:)));
   cMosaic.emPositions = [x(:),y(:)];
   cMosaic.name = 'Horizontal em';
   cMosaic.compute(oisH);
   cMosaic.computeCurrent;
   
   deMeanedMosaic = bsxfun(@minus,cMosaic.current,mean(cMosaic.current,3));
   Spectra_withMovements = abs(fft(deMeanedMosaic,128,3));
   avgSpectrum_withMovements = squeeze(sum(sum(Spectra_withMovements)))./93^2;
   
   
   avgSpectraAmp = sum(sum(Spectra_withMovements./Spectra_noMovements))./93^2;
    redCones = find(cMosaic.pattern == 3);
    avgredConesAmp = zeros(128,1);
    for i = 1:length(redCones)
        col = ceil(redCones(i)/93);
        row = mod(redCones(i),93);
        if row == 0
            row = 93;
        end
        coneSpectrum_noMovements = Spectra_noMovements(row,col,:);
        coneSpectrum_withMovements = Spectra_withMovements(row,col,:);
        avgredConesAmp = avgredConesAmp + squeeze(coneSpectrum_withMovements./coneSpectrum_noMovements);
    end

avgredConesAmp = avgredConesAmp./length(redCones); plot((0:tSamples/2-1)./0.128,squeeze(avgredConesAmp(1:64))); end; hold off; xlabel('Temporal Frequency (Hz)'); ylabel('Ratio (EM to no EM)'); title('Ratios for Red Cones'); legend('1 cpd','13.25 cpd','25.5 cpd','37.75 cpd','50 cpd');