/MatlabCode/branches/Greg's Branch/Analysis/WhiteNoise/WNAnalysis.m
MATLAB | 1351 lines | 1227 code | 47 blank | 77 comment | 64 complexity | 88de9d76cb6fc4a79c34fe2e4283db7c MD5 | raw file
Possible License(s): GPL-2.0, AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- % Scripts to analyze white noise data
- %
- % Loading the data
-
- WN=nex2stro;
- framerate = WN.sum.exptParams.framerate;
- nstixperside = WN.sum.exptParams.nstixperside;
- ntrials = length(WN.sum.absTrialNum);
- stimonidx = find(strcmp(WN.sum.trialFields(1,:),'stim_on'));
- stimoffidx = find(strcmp(WN.sum.trialFields(1,:),'all_off'));
- nframesidx = find(strcmp(WN.sum.trialFields(1,:),'num_frames'));
- noisetypeidx = find(strcmp(WN.sum.trialFields(1,:),'noise_type'));
- sigmaidxs = strmatch('sigma',WN.sum.trialFields(1,:));
-
- hepidx = find(strcmp(WN.sum.rasterCells(1,:),'AD11'));
- vepidx = find(strcmp(WN.sum.rasterCells(1,:),'AD12'));
- anlgStartTimeidx = find(strcmp(WN.sum.rasterCells(1,:),'anlgStartTime'));
- eyestart_t = [WN.ras{:,anlgStartTimeidx}]';
- eyesampperiod = 1/WN.sum.analog.storeRates{1};
- gammaTable = WN.sum.exptParams.gamma_table;
- gammaTable = reshape(gammaTable, length(gammaTable)/3, 3);
- gammaTable1 = interp1(linspace(0,255,256),gammaTable,linspace(0,255,65536), 'spline');
- invgamma = InvertGamma(gammaTable, 0);
-
- % Reconstructing the M matrix
- fundamentals = WN.sum.exptParams.fundamentals;
- fundamentals = reshape(fundamentals,[length(fundamentals)/3,3]);
- mon_spd = WN.sum.exptParams.mon_spd;
- mon_spd = reshape(mon_spd,[length(mon_spd)/3, 3]);
- mon_spd = SplineRaw([380:4:780]', mon_spd, [380:5:780]');
- M = fundamentals'*mon_spd;
-
- % Getting the background rgb/lms
- ridx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_r'));
- gidx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_g'));
- bidx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_b'));
- bkgndRGB = [mode(WN.trial(:,ridx)), mode(WN.trial(:,gidx)), mode(WN.trial(:,bidx))];
- bkgndrgb = [gammaTable(bkgndRGB(1)+1,1); gammaTable(bkgndRGB(2)+1,2); gammaTable(bkgndRGB(3)+1,3)];
- bkgndlms = M*bkgndrgb;
-
- %%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Stuff for looking at spike-triggered averages, etc.
- spikename = getSpikenum(WN);
- spikeidx = find(strcmp(WN.sum.rasterCells(1,:),spikename));
- maxT = 9;
- Lgunnoise = WN.trial(:,noisetypeidx) == 1;
- Lconenoise = WN.trial(:,noisetypeidx) == 2;
- if (sum(Lconenoise) == 0)
- whichnoisetype = 1;
- disp('Gun noise only');
- elseif (sum(Lgunnoise) == 0)
- whichnoisetype = 2;
- disp('Cone noise only');
- else
- querystr = ['Which noisetype? 1=gun (n=',num2str(sum(Lgunnoise)),') 2=cone (n=',num2str(sum(Lconenoise)),'): '];
- whichnoisetype = input(querystr);
- end
- tmpstro = WN;
- if (whichnoisetype == 2) && (any(Lgunnoise)) % Removing gun noise
- disp(['Getting rid of ',num2str(sum(Lgunnoise)),' gun noise trials!']);
- tmpstro.ras(Lgunnoise,:) = [];
- tmpstro.trial(Lgunnoise,:) = [];
- elseif (whichnoisetype == 1) && (any(Lconenoise)) % Removing cone noise
- disp(['Getting rid of ',num2str(sum(Lconenoise)),' cone noise trials!']);
- tmpstro.ras(Lconenoise,:) = [];
- tmpstro.trial(Lconenoise,:) = [];
- end
- out = getWhtnsStats(tmpstro,maxT,'STCOVmex', {nstixperside^2, 3, maxT}, spikename);
- tmpstro = [];
- STAs = out{1};
- STCs = out{2};
- nspikes = out{3};
- %%
- % Plotting the STA and PCs
- if (~exist('mask'))
- mask = zeros(nstixperside, nstixperside); % Pixel mask. Someday make a tool to make this non-zero.
- end
-
- Lmask = logical(repmat(~mask(:),[3 1]));
- PCs = [];
- for i = 1:size(STCs,2)
- STC = reshape(STCs(:,i), 3*nstixperside^2, 3*nstixperside^2);
- subSTC = STC(Lmask, Lmask);
- subSTA = STAs(Lmask,i);
- P = eye(size(subSTC))-subSTA*inv(subSTA'*subSTA)*subSTA';
- subSTC = P*subSTC*P';
- [tmp,d] = eig(subSTC);
- v = repmat(double(Lmask),[1 size(tmp,2)]);
- v(Lmask,:) = tmp;
- [newd, idxs] = sort(diag(d));
- v = v(:,idxs);
- v = v(:,[end end-1 end-2 2]); % Not taking last (trivially zero) vector
- PCs = cat(3, PCs, v);
- end
- % Normalizing images
- nstixperside = WN.sum.exptParams.nstixperside;
- template = reshape([1:nstixperside^2],nstixperside,nstixperside);
- edgepixels = [template(:,1); template(1,[2:end-1])'; template(end,[2:end-1])'; template(:,end)];
- if (all(mask(edgepixels)))
- PCs = PCs.*max(STAs(:));
- else
- edgepixelidxs = [edgepixels; edgepixels+nstixperside^2; edgepixels+2*(nstixperside^2)];
- PCelements = PCs(edgepixelidxs,:,:);
- PCsds = std(PCelements); % One std calculated per PC
- PCs = PCs.*repmat(std(STAs(:,1))./PCsds,[size(PCs,1),1,1]);
- end
-
- rowidxs = reshape([1:3*nstixperside^2],[nstixperside^2 3]);
- maxes = []; mins = [];
- imagevectors = [STAs, reshape(PCs,[300 size(PCs,2)*size(PCs,3)])];
- for i = 1:3
- maxes = [maxes; max(max(imagevectors(rowidxs(:,i),:)))];
- mins = [mins; min(min(imagevectors(rowidxs(:,i),:)))];
- end
- potentialnormfactors = [(1-[.5; .5; .5]-eps)./(maxes-[.5; .5; .5]); (-[.5; .5; .5]+eps)./(mins-[.5; .5; .5])];
- % 'eps' in above line is a kludge that is required for avoiding
- % out of bounds errors.
- potentialnormfactors(potentialnormfactors < 0) = []; % if min > mu or max < mu
- normfactor = min(potentialnormfactors);
-
- %muvect = reshape(repmat(bkgndrgb',nstixperside^2,1),nstixperside^2*3,1);
- muvect = reshape(repmat([.5 .5 .5],nstixperside^2,1),nstixperside^2*3,1);
-
- % Plotting
- figure;
- for i = 1:size(STAs,2)
- STA = normfactor*(STAs(:,i)-muvect)+muvect;
- STA = reshape(STA,[nstixperside nstixperside 3]);
- subplot(6,size(STAs,2),i);
- image(STA);
- set(gca,'XTick',[],'YTick',[]); axis square;
- if (i == 1)
- ylabel('STA');
- end
- for j = 1:size(v,2)
- PC = normfactor*(PCs(:,j,i)-muvect)+muvect;
- PC = reshape(PC,[nstixperside nstixperside 3]);
- subplot(6,size(STAs,2), j*size(STAs,2)+i);
- image(PC);
- set(gca,'XTick',[],'YTick',[]); axis square;
- end
- STC = reshape(STCs(:,i),[sqrt(length(STCs(:,i))),sqrt(length(STCs(:,i)))]);
- STV = reshape(diag(STC),[nstixperside nstixperside 3]);
- STV = (STV-mean(STV(:)))./(2*range(STV(:)))+.5;
- subplot(6,size(STAs,2),5*size(STAs,2)+i);
- image(STV);
- set(gca,'XTick',[],'YTick',[]); axis square;
- if (i == 1)
- ylabel('STV');
- end
- end
-
- %%
- % Significance testing on the STAs, STVs, broken down by color
- L = WN.trial(:,noisetypeidx)==1;
- mu1idx = find(strcmp(WN.sum.trialFields(1,:),'mu1'));
- mu2idx = find(strcmp(WN.sum.trialFields(1,:),'mu2'));
- mu3idx = find(strcmp(WN.sum.trialFields(1,:),'mu3'));
- sigma1idx = find(strcmp(WN.sum.trialFields(1,:),'sigma1'));
- sigma2idx = find(strcmp(WN.sum.trialFields(1,:),'sigma2'));
- sigma3idx = find(strcmp(WN.sum.trialFields(1,:),'sigma3'));
- muvect = unique(WN.trial(L,[mu1idx mu2idx mu3idx]),'rows')/1000;
- sigmavect = unique(WN.trial(L,[sigma1idx sigma2idx sigma3idx]),'rows')/1000;
- sigmavect(all(any(sigmavect == 0),2),:) = [];
- gausslims = [WN.sum.exptParams.gauss_locut WN.sum.exptParams.gauss_hicut]/1000;
- mumat = repmat(reshape(repmat(muvect,nstixperside^2,1),[nstixperside^2*3, 1]),[1,size(STAs,2)]);
- sigmamat = repmat(reshape(repmat(sigmavect,nstixperside^2,1),[nstixperside^2* 3, 1]),[1,size(STAs,2)]);
- zscoremeans = (STAs-mumat)./(sigmamat*sqrt(nspikes));
-
- % Doing the calculations for the variances
- % Only considering one correction factor per dimension
- % (assuming variances on green and blue guns are same as red gun)
-
- NPOINTS = 65536;
- x = linspace(gausslims(1),gausslims(2),NPOINTS);
- Fx = norminv(x)*sigmavect(1);
- sigmacorrectionfactor = std(Fx)./sigmavect(1);
- for i = 1:size(STCs,2)
- STC = reshape(STCs(:,i),[sqrt(length(STCs(:,i))),sqrt(length(STCs(:,i)))]);
- STVs(:,i) = diag(STC)./nspikes;
- end
- muvar = (sigmavect(1)*sigmacorrectionfactor)^2;
- sigmavar = muvar*sqrt(2/nspikes);
- zscorevars = (STVs-muvar)./sigmavar;
- maxzscore = max(abs([zscoremeans(:);zscorevars(:)]));
- alpha = 0.01;
- crit = norminv(1-alpha,0,1);
- % Plotting
- figure;
- for i = 1:size(STAs,2)
- % First STAs
- zmat = reshape(zscoremeans(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- pmat = logical(abs(zmat(:,:,j))>crit);
- im = zmat(:,:,j)./(2*maxzscore)+.5;
- im = repmat(im,[1 1 3]);
- sigidxs = find(pmat);
- im(sigidxs) = .5; % red to .5 where sig. Looks red on dark and cyan on bright.
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- % Then STVs
- zmat = reshape(zscorevars(:,i),[nstixperside nstixperside 3]);
- for j = 4:6
- pmat = logical(abs(zmat(:,:,j-3))>crit);
- im = zmat(:,:,j-3)./(2*maxzscore)+.5;
- im = repmat(im,[1 1 3]);
- sigidxs = find(pmat);
- im(sigidxs) = .5;
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- end
-
- %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Looking at the responses to synthetic images
- % Finding which trials were synthetic image trials
- imageidx = find(strcmp(WN.sum.rasterCells(1,:),'synth_image'));
-
- allimages = WN.ras(:,imageidx);
- L = [];
- for i = 1:length(allimages)
- L = [L; ~isnan(allimages{i}(1))];
- end
- figure;
- plot(L);
- [x,y] = ginput(2);
- close;
- imageidxs = find(L);
- imageidxs = imageidxs(imageidxs<max(x) & imageidxs>min(x));
- % Pulling out the images and putting them through the gamma functions
- % to get them in normalized intensity units.
- nelements = length(allimages{imageidxs(1)});
- npix = nelements/3;
-
-
- npixperside = sqrt(npix);
- images = nan(nelements,length(imageidxs));
- for i = 1:length(imageidxs)
- im = allimages{imageidxs(i)};
- for j = 1:3
- idxs = [1:npix]+npix*(j-1);
- images(idxs,i) = gammaTable1(im(idxs),j);
- end
- end
-
- % Subtracting the background
- bkgndrgbmat = repmat(bkgndrgb', npix, 1);
- bkgndrgbmat = repmat(bkgndrgbmat(:), 1, length(imageidxs));
- [uniqueimages, i, whichimage] = unique((images-bkgndrgbmat)','rows');
- nuniqueims = max(whichimage);
- spikerates = {};
- meanfr = [];
- stdfr = [];
- ns = [];
- for i = 1:max(whichimage)
- Lcont = whichimage == i;
- stimon_t = WN.trial(imageidxs(Lcont),stimonidx);
- numframes = WN.trial(imageidxs(Lcont),nframesidx)
- stimoff_t = stimon_t+numframes/framerate;
- spikes = WN.ras(imageidxs(Lcont),spikeidx);
- stimdurations = stimoff_t-stimon_t;
- for j = 1:sum(Lcont)
- spikerates{i}(j) = sum((spikes{j}>stimon_t(j)) & (spikes{j} < stimoff_t(j)))/stimdurations(j);
- end
- meanfr = [meanfr; mean(spikerates{i})];
- stdfr = [stdfr; std(spikerates{i})];
- ns = [ns; sum(Lcont)];
- end
- sem = stdfr./sqrt(ns);
- sem(ns < 2) = nan;
-
- [u,s,v] = svd(uniqueimages');
- if ((size(s,2) == 1) || s(1,1)./s(2,2) > 1000) % 1-D images
- [y,contrastsort] = sort(v(:,1));
- figure;
- subplot(2,1,1); hold on;
- plot(v(contrastsort,1)*s(1),meanfr(contrastsort),'ko-', 'Linewidth', 1,'MarkerSize',4)
- plot([v(contrastsort,1) v(contrastsort,1)]*s(1),[meanfr(contrastsort)-sem(contrastsort) meanfr(contrastsort)+sem(contrastsort)],'k-.')
- ylabel('sp/sec');
- xlabel('Projection');
- subplot(2,2,3);
- avgim = v(contrastsort(1))*u(:,1)+bkgndrgbmat(:,1);
- image(reshape(avgim,npixperside, npixperside, 3));
- axis image; set(gca,'Xtick',[],'Ytick',[]);
- subplot(2,2,4);
- avgim = v(contrastsort(end))*u(:,1)+bkgndrgbmat(:,1);
- image(reshape(avgim,npixperside, npixperside, 3));
- axis image; set(gca,'Xtick',[],'Ytick',[]);
-
- figure; % Subplots may not be arranged in a meaningful order
- for i = 1:nuniqueims
- subplot(1,nuniqueims,i);
- scalefactor = .5/max(abs(uniqueimages(:)));
- image(reshape(scalefactor*uniqueimages(i,:)',[nstixperside, nstixperside, 3])+.5);
- title(num2str(meanfr(i),2));
- set(gca,'Xtick',[],'YTick',[]);
- axis image;
- end
-
- else % 2-D images
- % Rotating coordinate axes
- dists = sum(v(:,1).^2+v(:,2).^2,2);
- corneridx = find(dists == max(dists),1);
- thetas = [atan(v(corneridx,2)./v(corneridx,1))+pi/4;...
- atan(v(corneridx,2)./v(corneridx,1))-pi/4];
- theta = thetas(abs(thetas) == min(abs(thetas)));
- rotmat = [cos(theta) sin(theta); -sin(theta) cos(theta)];
- basis = u(:,[1 2])*rotmat' * rotmat*s([1 2], [1 2])*rotmat';
- basis = pinv(basis)'; % im = b'*w so w = im*pinv(b)
- norms = sqrt(sum(basis.^2));
- basis = basis./repmat(norms, size(basis,1),1);
- % This is flawed because the noise in the STA contributes to
- % its norm but we're zeroing out nuisance elements in the PC1 so the
- % when the are constrained to have the same norm, the "signal" in the
- % PC1 vector is huge, and so the weights onto the PC1 vector are
- % artificially small (relative to projections onto the STA).
- %weights = (rotmat*s([1 2], [1 2])*rotmat') * rotmat*v(:,[1 2])';
-
- weights = uniqueimages*basis;
- normweights = [(weights(:,1)-min(weights(:,1)))/(max(weights(:,1))-min(weights(:,1))),...
- (weights(:,2)-min(weights(:,2)))/(max(weights(:,2))-min(weights(:,2)))];
- rankweights = round(normweights*(sqrt(nuniqueims)-1))+1;
- frim = zeros(sqrt(nuniqueims));
- semim = zeros(sqrt(nuniqueims));
- figure;
- for i = 1:nuniqueims;
- frim(rankweights(i,1),rankweights(i,2)) = meanfr(i);
- semim(rankweights(i,1),rankweights(i,2)) = sem(i);
- ind = sub2ind(max(rankweights), rankweights(i,1), rankweights(i,2))
- subplot(max(rankweights(:,1)), max(rankweights(:,2)), ind);
- scalefactor = .5/max(abs(uniqueimages(:)));
- image(reshape(scalefactor*uniqueimages(i,:)',[nstixperside, nstixperside, 3])+.5);
- title(num2str(meanfr(i),2));
- set(gca,'Xtick',[],'YTick',[]);
- axis image;
- end
- frim = frim';
- semim = semim';
- figure; imagesc(frim); axis xy; axis image;
- colormap(gray);
- end
-
- %%
- % For L vs M replays
- % L is on the x axis, M is on the y axis
- LMScontrast = [];
- for i = 1:nuniqueims
- im = reshape(uniqueimages(i,:),[nstixperside^2 3]);
- [u,s,v] = svd(im);
- if (max(abs(u(:,1))) ~= max(u(:,1)))
- s(1) = -s(1);
- end
- LMScontrast(:,i) = (M*im(find(abs(u(:,1))==max(abs(u(:,1)))),:)')./bkgndlms;
- end
- rmscontrast = unique(sqrt(mean(LMScontrast.^2))); % in case we're interested in this
- totalcontrast = sqrt(sum(LMScontrast.^2));
- uniquecontrasts = unique(round(totalcontrast*100))/100;
- data = [];
- for i = 1:length(uniquecontrasts)
- L = abs(totalcontrast-uniquecontrasts(i)) < 1/100;
- if (sum(L) == 8)
- frs = meanfr(L);
- theta = atan2(LMScontrast(2,L), LMScontrast(1,L))';
- [theta,j] = sort(theta);
- data = cat(3,data,[theta, frs(j)]);
- disp('theta contrast');
- disp([theta, totalcontrast(L)']);
- end
- end
- figure; subplot(3,3,5);
- posbColors = {'k', 'b', 'm', 'c', 'g', 'r'};
- for i = size(data,3):-1:1
- theta = data(:,1,i);
- resp = data(:,2,i);
- polar([theta; theta(1)],[resp; resp(1)],[posbColors{i}, '-']);
- hold on;
- end
- ims = uniqueimages(L,:);
- theta = squeeze(data(:,1,end));
- resp = squeeze(data(:,2,end));
- for i = 1:8
- binnedtheta = round(mod(theta(i),2*pi)./(pi/4));
- if (binnedtheta == 8 | binnedtheta == 0)
- plotnum = 6;
- elseif (binnedtheta > 4)
- plotnum = binnedtheta +2;
- elseif (binnedtheta == 4)
- plotnum = 4;
- else
- plotnum = 4-binnedtheta;
- end
- subplot(3,3,plotnum);
- image(reshape(ims(j(i),:)',[nstixperside, nstixperside, 3])+.5);
- title(num2str(resp(i),2));
- set(gca,'Xtick',[],'YTick',[]);
- axis image;
- end
- %%
- % For L vs M replays at multiple contrasts
- if (size(data,3) > 1)
- posbColors = {'k', 'b', 'm', 'c', 'g', 'r'};
- figure; axes; hold on;
- for i = 1:size(data,3);
- plot(data(:,1,i),data(:,2,i),[posbColors{i}, 'o-'])
- end
- end
- set(gca,'YScale','log');
- %%
- % Rasters for the replay images
-
- figure; axes; hold on;
- counter = 0;
- for i = 1:size(uniqueimages,2)
- Lcont = whichimage == i;
- stimon_t = WN.trial(imageidxs(Lcont),stimonidx);
- numframes = WN.trial(imageidxs(Lcont),nframesidx);
- stimoff_t = stimon_t+numframes/framerate;
- spikes = WN.ras(imageidxs(Lcont),spikeidx);
-
- for j = 1:sum(Lcont)
- nsp = length(spikes{j})
- if (nsp > 0)
- plot(0,counter,'g*');
- plot(stimoff_t(j)-stimon_t(j),counter,'r*');
- else
- plot(0,counter,'k*');
- plot(stimoff_t(j)-stimon_t(j),counter,'k*');
- end
- plot([spikes{j} spikes{j}]'-stimon_t(j),[zeros(nsp,1) .5*ones(nsp,1)]'+counter,'k-')
- counter = counter + 1;
- end
- plot([-.2 1],[counter counter],'k-');
- end
- set(gca,'Xlim',[-.2 1]); xlabel('time(s)');
- title(WN.sum.fileName);
-
-
- %%
- % Getting the projections onto the PC1 and STA
- % Just experimenting with the code
-
- whichframe = 4; % in the middle
- nframestoconsider = 2;
- msperframe = 1000/WN.sum.exptParams.framerate;
- %%%%%%%%%%%%
- % Skip this part if using a basis that came from somewhere else
- %%%%%%%%%%
- if (~exist('basis'))
- STC = reshape(STCs(:,whichframe), 3*nstixperside^2, 3*nstixperside^2);
- [v,d] = eig(STC);
- basis = [STAs(:,whichframe) v(:,end)];
- basis = basis./repmat(sqrt(sum(basis.^2)), size(basis,1),1);
- end
-
- nframestotal = sum(WN.trial(:,nframesidx));
- initargs = {mkbasis(repmat(basis,nframestoconsider,1)) 3, nframestotal, [nstixperside^2 3 nframestoconsider]};
- out = getWhtnsStats(WN,whichframe-ceil(nframestoconsider/2),'STPROJmod',initargs);
- projs = out{1};
- Lspike = out{2};
-
- % 1-D firing rate plot
- whichvect = 2;
- nbins = 5;
- x = linspace(prctile(projs(:,whichvect), 5), prctile(projs(:,whichvect), 95),nbins+2);
- Na = hist(projs(:,whichvect), x);
- Ns = zeros(size(Na));
- for i = 1:max(Lspike)
- Ns = Ns + hist(projs(Lspike >=i,whichvect), x);
- end
- figure;
- plot(x(2:end-1),Ns(2:end-1)./Na(2:end-1)/msperframe*1000,'ko');
-
- % 2-D firing rate plot
- if (round(sqrt(nuniqueims)) == sqrt(nuniqueims));
- nbins = sqrt(nuniqueims);
- else
- nbins = 4;
- end
- x = [nbins nbins; prctile(projs, 10); prctile(projs, 90)];
- %x = [nbins nbins; min(weights); max(weights)];
-
- Na = hist2(projs, x);
- Ns = zeros(size(Na));
- for i = 1:max(Lspike)
- Ns = Ns + hist2(projs(Lspike >=i,:), x);
- end
-
- figure;
- imageprops = mkbasis([-1 1; 0 1; 1 1; -1 0; 0 0; 1 0; -1 -1; 0 -1; 1 -1]')';
- for i = [1 2 3 4 6 7 8 9]
- subplot(3,3,i);
- tmpim = basis/2*imageprops(i,:)'+.5; % just for a quick rendering
- image(reshape(tmpim,[nstixperside nstixperside 3]));
- axis image; set(gca,'XTick',[]); set(gca,'Ytick',[]);
- end
- subplot(3,3,5)
- imagesc(Ns./Na);
- axis image; set(gca,'XTick',[]); set(gca,'Ytick',[]);
- colormap gray;
- figure;
- plot(frim,(Ns./Na)/msperframe*1000,'k.');
-
- %%
- % Rotating the basis vectors onto which we project the spike-triggered
- % stimuli. Ideally, we'll see the best correlation between the projected
- % firing rate plot and the empirical firing rate plot when we don't rotate
- % at all.
- nbins = sqrt(length(frim(:)));
- x = [nbins nbins; min(weights); max(weights)];
- data = [];
- for theta = linspace(-pi/4,pi/4,5)
- theta
- rotmat = [cos(theta) sin(theta); -sin(theta) cos(theta)];
- initargs = {basis*rotmat', 3, nframestotal, [nstixperside^2 3 1]};
- out = getWhtnsStats(WN,whichframe,'STPROJmod',initargs);
- projs = out{1};
- Lspike = out{2};
- Na = hist2(projs, x);
- Ns = zeros(size(Na));
- for i = 1:max(Lspike)
- Ns = Ns + hist2(projs(Lspike >=i,:), x);
- end
- [rho,p] = corr(frim(:),Ns(:)./Na(:),'type','Kendall')
- data = [data; theta rho p];
- end
- figure;
- plot(data(:,1)*180/pi,data(:,2),'k-.');
- xlabel('Deg axes rotation')
- ylabel('r')
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Everything below this point is more special-purpose code...
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- %%
- % Looking for missed frames
- % (Late peak is from the synth image trials)
- stimon_t = WN.trial(:,stimonidx);
- stimoff_t = WN.trial(:,stimoffidx);
- framerate = WN.sum.exptParams.framerate;
- predstimoff_t = stimon_t+WN.trial(:,nframesidx)/framerate;
- x = stimoff_t-predstimoff_t;
-
- hist(x(x<0),100);
- plot(WN.trial(:,nframesidx), x,'k.')
-
-
- %%
- % Analyzing fixational saccades
- sacstats = getSacData(WN);
- L = logical(WN.trial(:,noisetypeidx)~= 0);
- amplitudes = [];
- directions = [];
- peakv = [];
- pathlengths = [];
- durations = [];
- BINWIDTH = .01;
- timebins = [-.25:BINWIDTH:.5];
- PSTHmat = zeros(sum(L),length(timebins));
- for i = find(L')
- stimon_t = WN.trial(i,stimonidx);
- numframes = WN.trial(i,nframesidx);
- stimoff_t = stimon_t+numframes/WN.sum.exptParams.framerate;
- st = sacstats.starttimes{i};
- Lsac = (st > stimon_t) & (st < stimoff_t-.1); %.1 to avoid a saccade that leaves the window
- if any(Lsac)
- if any(sacstats.amplitudes{i}(Lsac) > 2)
- keyboard
- end
- amplitudes = [amplitudes; sacstats.amplitudes{i}(Lsac)];
- directions = [directions; sacstats.directions{i}(Lsac)];
- peakv = [peakv; sacstats.peakv{i}(Lsac)];
- pathlengths = [pathlengths; sacstats.pathlengths{i}(Lsac)];
- durations = [durations; sacstats.durations{i}(Lsac)];
- spiketimes = [];
- for j = find(Lsac')
- tmp = WN.ras{i,1}-sacstats.starttimes{i}(j);
- spiketimes = [spiketimes; tmp((tmp > timebins(1)-BINWIDTH/2) & (tmp < timebins(end)+BINWIDTH/2))];
- end
- [n,x] = hist(spiketimes, timebins);
- PSTHmat(i,:) = n./(BINWIDTH*sum(Lsac));
- end
- end
- figure;
- subplot(3,2,1);
- hist(amplitudes,100);
- xlabel('amplitude'); ylabel('count');
- subplot(3,2,3);
- plot(amplitudes,peakv,'k.');
- xlabel('amplitude'); ylabel('peak vel.');
- subplot(3,2,4);
- plot(amplitudes,pathlengths,'k.');
- xlabel('amplitude'); ylabel('traj. length');
- subplot(3,2,5);
- [rho, theta]= hist(directions,20);
- polar([theta theta(1)],[rho rho(1)],'k-')
- subplot(3,2,6); hold on;
- plot(timebins, mean(PSTHmat),'k-');
- plot([0 0],[0 max(mean(PSTHmat))],'b:');
- set(gca,'Xlim',[min(timebins) max(timebins)]);
- xlabel('time wrt saccade onset (s)'); ylabel('sp/sec');
- %%
- % Cross validation to find static nonlinearities
- % 1-D version
-
- niters = 2;
- CVprop = .7;
- whichframe = 4;
- nbins = 5;
- Lgaussnoise = WN.trial(:,noisetypeidx) == 1;
- binlim = 2*mode(WN.trial(Lgaussnoise,9))/1000
- ntrials = size(WN.trial,1);
- ntrialsSTC = round(ntrials*CVprop);
- frs = [];
- vs = [];
- for iter = 1:niters
- trialidxs = randperm(ntrials);
- tmpstro = WN;
- tmpstro.trial = WN.trial(trialidxs(1:ntrialsSTC),:);
- tmpstro.ras = WN.ras(trialidxs(1:ntrialsSTC),:);
- out = getWhtnsStats(tmpstro, maxT, 'STCOVmex', {nstixperside^2, 3, whichframe});
- STAs = out{1};
- STCs = out{2};
- STC = reshape(STCs(:,whichframe), 3*nstixperside^2, 3*nstixperside^2);
- STC = STC.*~Lmask;
- [v,d] = eig(STC);
- d = diag(d);
-
- basis = [STAs(:,whichframe) v(:,end)];
- basis = basis./sqrt(sum(basis.^2));
- if (basis'*synthimvect < 0)% synthimvect is created by the chunk of code above
- basis = -basis;
- end
- tmpstro = WN;
- tmpstro.trial = WN.trial(trialidxs(ntrialsSTC+1:end),:);
- tmpstro.ras = WN.ras(trialidxs(ntrialsSTC+1:end),:);
- nframestotal = sum(WN.trial(:,nframesidx));
- initargs = {basis, whichframe-1, nframestotal, [nstixperside^2 3 1]};
- out = getWhtnsStats(tmpstro,whichframe,'STPROJmod',initargs);
- projs = out{1};
- Lspike = out{2};
-
- x = linspace(-binlim, binlim, nbins+2);
- Na = hist(projs(:,1), x);
- Ns = zeros(size(Na));
- for i = 1:max(Lspike)
- Ns = Ns + hist(projs(Lspike >=i,1), x);
- end
- frs = [frs; Ns./Na]
- vs = [vs,basis];
- end
-
- figure; subplot(2,1,1);
- plot(x,mean(frs));
- subplot(2,2,4);
- image(reshape(mean(vs,2),[nstixperside, nstixperside, 3])+.5); % Not rendered accurately
- axis image; set(gca,'XTick',[],'YTick',[]);
- subplot(2,2,3);
- image(reshape(-mean(vs,2),[nstixperside, nstixperside, 3])+.5); % Not rendered accurately
- axis image; set(gca,'XTick',[],'YTick',[]);
-
-
- %%
- % Cross validation to find static nonlinearities
- % 2-D version
-
- niters = 10;
- CVprop = .7;
- whichframe = 4;
- nbins = 10;
- Lgaussnoise = WN.trial(:,noisetypeidx) == 1;
- binlim = 1*mode(WN.trial(Lgaussnoise,9))/1000;
- ntrials = size(WN.trial,1);
- ntrialsSTC = round(ntrials*CVprop);
- Nas = [];
- Nss = [];
- for iter = 1:niters
- trialidxs = randperm(ntrials);
- tmpstro = WN;
- tmpstro.trial = WN.trial(trialidxs(1:ntrialsSTC),:);
- tmpstro.ras = WN.ras(trialidxs(1:ntrialsSTC),:);
- out = getWhtnsStats(tmpstro, maxT, 'STCOVmex', {nstixperside^2, 3, whichframe});
- STAs = out{1};
- STCs = out{2};
- STC = reshape(STCs(:,whichframe), 3*nstixperside^2, 3*nstixperside^2);
- [v,d] = eig(STC);
-
- basis = mkbasis([STAs(:,whichframe) v(:,end)]);
- basis(:,2) = sign(basis(:,1)'*basis(:,2))*basis(:,2); % PC constrained to have a + proj onto STA
- tmpstro = WN;
- tmpstro.trial = WN.trial(trialidxs(ntrialsSTC+1:end),:);
- tmpstro.ras = WN.ras(trialidxs(ntrialsSTC+1:end),:);
- nframestotal = sum(WN.trial(:,nframesidx));
- initargs = {basis, whichframe-1, nframestotal, [nstixperside^2 3 1]};
- out = getWhtnsStats(tmpstro,whichframe,'STPROJmod',initargs);
- projs = out{1};
- Lspike = out{2};
-
- x = [nbins nbins; -binlim -binlim; binlim binlim];
- Na = hist2(projs, x);
- Ns = zeros(size(Na));
- for i = 1:max(Lspike)
- Ns = Ns + hist2(projs(Lspike >=i,:), x);
- end
- Nas = cat(3,Nas,Na);
- Nss = cat(3,Nss,Ns);
- end
- fr = Nss./Nas;
- figure;
- imagesc(mean(fr,3));
-
- %% Looking at STA from the multielectrode array (no STC analysis because stimulus is too big).
- %
-
- WN=nex2stro;
- LFPsamprate = WN.sum.analog.storeRates{1}; % Assuming all channels are sampled at the same rate
- framerate = WN.sum.exptParams.framerate;
- nstixperside = WN.sum.exptParams.nstixperside;
- ntrials = length(WN.sum.absTrialNum);
- stimonidx = find(strcmp(WN.sum.trialFields(1,:),'stim_on'));
- stimoffidx = find(strcmp(WN.sum.trialFields(1,:),'all_off'));
- nframesidx = find(strcmp(WN.sum.trialFields(1,:),'num_frames'));
- noisetypeidx = find(strcmp(WN.sum.trialFields(1,:),'noise_type'));
- sigmaidxs = strmatch('sigma',WN.sum.trialFields(1,:));
-
- hepidx = find(strcmp(WN.sum.rasterCells(1,:),'AD11'));
- vepidx = find(strcmp(WN.sum.rasterCells(1,:),'AD12'));
- anlgStartTimeidx = find(strcmp(WN.sum.rasterCells(1,:),'anlgStartTime'));
- eyestart_t = [WN.ras{:,anlgStartTimeidx}]';
- eyesampperiod = 1/WN.sum.analog.storeRates{1};
- gammaTable = WN.sum.exptParams.gamma_table;
- gammaTable = reshape(gammaTable, length(gammaTable)/3, 3);
- gammaTable1 = interp1(linspace(0,255,256),gammaTable,linspace(0,255,65536), 'spline');
- invgamma = InvertGamma(gammaTable, 0);
-
- % Reconstructing the M matrix
- fundamentals = WN.sum.exptParams.fundamentals;
- fundamentals = reshape(fundamentals,[length(fundamentals)/3,3]);
- mon_spd = WN.sum.exptParams.mon_spd;
- mon_spd = reshape(mon_spd,[length(mon_spd)/3, 3]);
- mon_spd = SplineRaw([380:4:780]', mon_spd, [380:5:780]');
- M = fundamentals'*mon_spd;
-
- % Getting the background rgb/lms
- ridx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_r'));
- gidx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_g'));
- bidx = find(strcmp(WN.sum.trialFields(1,:),'bkgnd_b'));
- bkgndRGB = [mode(WN.trial(:,ridx)), mode(WN.trial(:,gidx)), mode(WN.trial(:,bidx))];
- bkgndrgb = [gammaTable(bkgndRGB(1)+1,1); gammaTable(bkgndRGB(2)+1,2); gammaTable(bkgndRGB(3)+1,3)];
- bkgndlms = M*bkgndrgb;
- maxT = 9;
-
- potentialspikeidxs = find(strncmp(WN.sum.rasterCells,'sig',3))
- for spikeidx = potentialspikeidxs
- out = getWhtnsStats(WN,maxT,'STAmex', {nstixperside^2, 3, maxT}, WN.sum.rasterCells{spikeidx});
- nspikes = out{3};
- STAs = out{1};
- %STVs = (out{2}-out{1}.^2)./nspikes;
- STVs = out{2}./nspikes;
-
- % Significance testing on the STAs, STVs, broken down by color
- L = WN.trial(:,noisetypeidx)==1; % 1 = gun noise
- mu1idx = find(strcmp(WN.sum.trialFields(1,:),'mu1'));
- mu2idx = find(strcmp(WN.sum.trialFields(1,:),'mu2'));
- mu3idx = find(strcmp(WN.sum.trialFields(1,:),'mu3'));
- sigma1idx = find(strcmp(WN.sum.trialFields(1,:),'sigma1'));
- sigma2idx = find(strcmp(WN.sum.trialFields(1,:),'sigma2'));
- sigma3idx = find(strcmp(WN.sum.trialFields(1,:),'sigma3'));
- muvect = unique(WN.trial(L,[mu1idx mu2idx mu3idx]),'rows')/1000;
- sigmavect = unique(WN.trial(L,[sigma1idx sigma2idx sigma3idx]),'rows')/1000;
- sigmavect(all(any(sigmavect == 0),2),:) = [];
- gausslims = [WN.sum.exptParams.gauss_locut WN.sum.exptParams.gauss_hicut]/1000;
- if (sum(L) > sum(~L)) % If more gun noise than cone noise
- mumat = repmat(reshape(repmat(muvect,nstixperside^2,1),[nstixperside^2*3, 1]),[1,size(STAs,2)]);
- sigmamat = repmat(reshape(repmat(sigmavect,nstixperside^2,1),[nstixperside^2* 3, 1]),[1,size(STAs,2)]);
- zscoremeans = (STAs-mumat)./(sigmamat*sqrt(nspikes));
-
- % Doing the calculations for the variances
- % Only considering one correction factor per dimension
- % (assuming variances on green and blue guns are same as red gun)
-
- NPOINTS = 65536;
- x = linspace(gausslims(1),gausslims(2),NPOINTS);
- Fx = norminv(x)*sigmavect(1);
- sigmacorrectionfactor = std(Fx)./sigmavect(1);
- muvar = (sigmavect(1)*sigmacorrectionfactor)^2;
- sigmavar = muvar*sqrt(2/nspikes);
- zscorevars = (STVs-muvar)./sigmavar;
- else
- zscoremeans = (STAs-mumat)./sqrt(nspikes);
- zscorevars = zeros(size(zscoremeans));
- end
- maxzscore = max(abs([zscoremeans(:);zscorevars(:)]));
- alpha = 0.01;
- crit = norminv(1-alpha,0,1);
-
- % Plotting
- figure;
- for i = 1:size(STAs,2)
- % First STAs
- zmat = reshape(zscoremeans(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- pmat = logical(abs(zmat(:,:,j))>crit);
- im = zmat(:,:,j)./(2*maxzscore)+.5;
- im = repmat(im,[1 1 3]);
- sigidxs = find(pmat);
- im(sigidxs) = .5; % red to .5 where sig. Looks red on dark and cyan on bright.
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- % Then STVs
- zmat = reshape(zscorevars(:,i),[nstixperside nstixperside 3]);
- for j = 4:6
- pmat = logical(abs(zmat(:,:,j-3))>crit);
- im = zmat(:,:,j-3)./(2*maxzscore)+.5;
- im = repmat(im,[1 1 3]);
- sigidxs = find(pmat);
- im(sigidxs) = .5;
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- end
- set(gcf,'Name',WN.sum.rasterCells{spikeidx});
- drawnow;
- end
-
- %%
- % Looking at LFP-triggered averages (filtered, interpolated LFPs)
- centerfreq = 100; % Hz
- cyclespersample = centerfreq./LFPsamprate;
- ncyclesper6sigma = 3; % Controls the bandwidth
- nsamplesper6sigma = ceil(ncyclesper6sigma/cyclespersample);
- ncycles = nsamplesper6sigma*cyclespersample;
- filtkernel1 = normpdf(linspace(-3,3,nsamplesper6sigma),0,1).*cos(linspace(0,2*pi*ncycles,nsamplesper6sigma));
- filtkernel2 = normpdf(linspace(-3,3,nsamplesper6sigma),0,1).*sin(linspace(0,2*pi*ncycles,nsamplesper6sigma));
- filtkernel1 = filtkernel1./norm(filtkernel1);
- filtkernel2 = filtkernel2./norm(filtkernel2);
- for whichADchan = 1:32;
- ntrials = size(WN.trial,1);
- nstixperside = WN.sum.exptParams.nstixperside; %get number of stixels per side
- msperframe = 1000/WN.sum.exptParams.framerate; %calculate msec per frame
- secsperframe = 1/WN.sum.exptParams.framerate;
- gammaTable = WN.sum.exptParams.gamma_table; %get gamma_table
- gammaTable = reshape(gammaTable,length(gammaTable)/3,3); %reshapse gamma_table into three columns
- invgammaTable = InvertGamma(gammaTable,1); %invert gamma_table
- ngammasteps = size(invgammaTable,1); %get number of rows of gamma_table (65536)
- seedidx = strcmp(WN.sum.trialFields(1,:),'seed'); %get seed index from trialFields
- nframesidx = strcmp(WN.sum.trialFields(1,:),'num_frames'); %get nframes index from trialFields
- stimonidx = strcmp(WN.sum.trialFields(1,:),'stim_on'); %get stimon index from trialFields
- muidxs = [find(strcmp(WN.sum.trialFields(1,:),'mu1')),... %get mu indices into vector from trialFields
- find(strcmp(WN.sum.trialFields(1,:),'mu2')),...
- find(strcmp(WN.sum.trialFields(1,:),'mu3'))];
- sigmaidxs = [find(strcmp(WN.sum.trialFields(1,:),'sigma1')),... %get sigma indices into vector from trialFields
- find(strcmp(WN.sum.trialFields(1,:),'sigma2')),...
- find(strcmp(WN.sum.trialFields(1,:),'sigma3'))];
- LFPstarttimes = [WN.ras{:,strcmp(WN.sum.rasterCells,'anlgStartTime')}];
-
- feval('STAmex','init', {nstixperside^2, 3, maxT});
- % This code was largely lifted from getWhtnsStats
- hwait = waitbar(0,'Finding stimuli...'); %display progress bar
- L = WN.trial(:,noisetypeidx)==1; % 1 = gun noise
- for i = 1:ntrials %get values and insert into given column
- if (WN.trial(i,noisetypeidx)==1)
- if (sum(L) > sum(~L))
- noisetype = 1;
- else
- continue
- end
- end
-
- if (WN.trial(i,noisetypeidx)==2)
- if (sum(~L) > sum(L))
- noisetype = 2;
- else
- continue
- end
- end
-
- seed = WN.trial(i,seedidx);
- nframes = WN.trial(i,nframesidx);
- mu = WN.trial(i,muidxs)/1000;
- sigma = WN.trial(i,sigmaidxs)/1000;
-
- if (noisetype == 1) % Gaussian gun
- x = linspace(WN.sum.exptParams.gauss_locut/1000, WN.sum.exptParams.gauss_hicut/1000, ngammasteps);%dividing gauss by 1000 and making equal intervals so that there are 65536 values
- for gun = 1:3
- invnormcdf(:,gun) = norminv(x)*sigma(gun)+mu(gun);
- end
- randnums = getEJrandnums(3*nstixperside^2*nframes, seed);
- randnums = reshape(randnums, [nstixperside^2*3, nframes]);
- for j = 1:3
- idxs = [1:nstixperside^2]+nstixperside^2*(j-1);
- randnums(idxs,:) = reshape(invnormcdf(randnums(idxs,:)+1,j),[length(idxs),nframes]);
- end
- elseif (noisetype == 2) % Binary cone
- colordirlms = sign(fullfact([2,2,2])-1.5);
- randnums = getEJrandnums(nstixperside^2*nframes, seed);
- randnums = reshape(randnums, [nstixperside^2, nframes]);
- randnums = mod(randnums, size(colordirlms,1))+1;
- tmp = colordirlms(randnums,:);
- tmp = reshape(tmp, [nstixperside^2 nframes 3]);
- tmp = permute(tmp,[1 3 2]);
- randnums = reshape(tmp,[nstixperside^2*3 nframes]);
- end
-
- t_stimon = WN.trial(i,stimonidx);
- LFP = WN.ras{i,strcmp(WN.sum.rasterCells,['AD',num2str(16+whichADchan)])};
- filteredLFP = conv(LFP, filtkernel1,'same').^2+conv(LFP, filtkernel2,'same').^2;
- LFPtimes = LFPstarttimes(i)+[0:length(filteredLFP)-1]/LFPsamprate; % sec
- frametimes = t_stimon+(linspace(0, nframes*secsperframe, nframes)+(secsperframe/2)'); % sec
- rectLFP = interp1(LFPtimes, filteredLFP, frametimes);
- % plot(LFPtimes,abs(LFP)); hold on; plot(frametimes,rectLFP,'m-');
- feval('STAmex',randnums(:),rectLFP);
- if (ishandle(hwait))
- waitbar(i/ntrials,hwait);
- else
- break;
- end
- end
- close(hwait);
- out = feval('STAmex','return');
- eval('clear STAmex');
-
- % Now doing some plotting
- STAs = out{1};
- STVs = out{2};
- STAlims = [min(STAs(:)) max(STAs(:))];
- STVlims = [min(STVs(:)) max(STVs(:))];
- STAs = (STAs-STAlims(2))/(STAlims(1)-STAlims(2));
- STVs = (STVs-STVlims(2))/(STVlims(1)-STVlims(2));
-
- figure; colormap(gray(255));
- set(gcf,'Name',['LFP ',num2str(whichADchan)]);
- for i = 1:size(STAs,2)
- % First STAs
- mat = reshape(STAs(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- im = mat(:,:,j);
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im*255);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- % Then STVs
- mat = reshape(STVs(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- im = mat(:,:,j);
- subplot(6,size(STAs,2),(j-1+3)*size(STAs,2)+i);
- image(im*255);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- end
- end
-
- %%
- % Looking at rectified LFPs at the time of stimus onset
- ADchans = find(strncmp('AD',WN.sum.rasterCells,2));
- LFPchans = ADchans(3:end); % Getting rid of the eye position channels
- hwait = waitbar(0,'Computing stimulus-locked LFPs...'); %display progress bar
- poststim_t = 0.5; % secs
- prestim_t = 0.3; % secs
- nprestimsamples = prestim_t*LFPsamprate;
- npoststimsamples = poststim_t*LFPsamprate;
- nsamples = nprestimsamples+npoststimsamples+1;
- PSTHs = zeros(length(LFPchans),nsamples);
- for i = 1:ntrials %get values and insert into given column
- t_stimon = WN.trial(i, stimonidx);
- for ADchan1 = LFPchans
- LFP = WN.ras{i,ADchan1};
- LFPtimes = LFPstarttimes(i)+[0:length(LFP)-1]/LFPsamprate; % sec
- err = (LFPtimes-t_stimon).^2;
- startidx = find(err == min(err),1);
- if (startidx > nprestimsamples & length(LFP)-startidx > npoststimsamples)
- LFPclip = LFP(startidx+[-nprestimsamples:npoststimsamples]);
- PSTHs(ADchan1 == LFPchans,:) = PSTHs(ADchan1 == LFPchans,:) + abs(LFPclip)';
- end
- end
- waitbar(i/ntrials,hwait);
- end
- close(hwait);
- figure; subplot(2,1,1);
- imagesc(PSTHs);
- subplot(2,1,2);
- tvect = linspace(-prestim_t, poststim_t, nsamples);
- plot(tvect,mean(PSTHs));
- % Cool. I think we're seeing the stimulus refresh in the LFP.
- L = tvect>.3;
- freqs = linspace(0, LFPsamprate/2, sum(L)/2);
- [Pxx,F] = periodogram(mean(PSTHs(:,L)),[],[],LFPsamprate)
- figure; plot(F,Pxx); set(gca,'YScale','log');
-
- %%
- % Looking at the auto- and cross-correlation functions of the LFPs
-
- hwait = waitbar(0,'Finding stimuli...'); %display progress bar
- maxlag = 0.2; % sec
- delayfromstimon = 0.25;
- maxlagsamples = maxlag*LFPsamprate;
- nsamples = 2*maxlagsamples+1;
- data = zeros(length(LFPchans),length(LFPchans), nsamples);
- for i = 1:ntrials % get values and insert into given column
- t_stimon = WN.trial(i, stimonidx);
- t_stimoff = WN.trial(i, stimoffidx);
- for idx1 = 1:length(LFPchans)
- for idx2 = idx1:length(LFPchans)
- LFP1 = WN.ras{i,LFPchans(idx1)};
- LFP2 = WN.ras{i,LFPchans(idx2)};
- LFPtimes = LFPstarttimes(i)+[0:length(LFP1)-1]/LFPsamprate; % sec
- err = (LFPtimes-t_stimon-delayfromstimon).^2;
- startidx = find(err == min(err),1);
- err = (LFPtimes-t_stimoff).^2;
- endidx = find(err == min(err),1);
- LFP1 = LFP1(startidx:endidx);
- LFP2 = LFP2(startidx:endidx);
- if (length(LFP1) >= nsamples)
- tmp = xcorr(LFP1,LFP2,maxlagsamples,'coeff');
- data(idx1,idx2,:) = tmp;
- end
- end
- end
- waitbar(i/ntrials,hwait);
- end
- close(hwait);
-
- % First looking at all the autocorrelation functions
- tvect = [-maxlagsamples:maxlagsamples]/LFPsamprate;
- tmp = zeros(length(LFPchans),nsamples);
- for i = 1:length(LFPchans)
- tmp(i,:) = data(i,i,:);
- end
- figure;
- subplot(2,1,1);
- title('Autocorrelations');
- imagesc(tmp); set(gca,'Xtick',[]);
- subplot(2,1,2);
- plot(tvect,tmp');
- xlabel('lag secs');
-
- % Now looking at cross-correlations
- tmp = [];
- for sig1 = 1:length(LFPchans)
- for sig2 = sig1+1:length(LFPchans)
- tmp = [tmp;squeeze(data(sig1,sig2,:))'];
- end
- end
- figure;
- imagesc(tmp);
-
- % % Looking at cross correlograms of adjacent and non-adjacent pairs of electrodes
- % % Only appropriate for most medial bank of 32 (port C).
- tmp1 = []; tmp2 = [];
- for i = 1:length(LFPchans)-1
- if (mod(i,2)) % if i is odd
- tmp1 = [tmp1; squeeze(data(i,i+1,:))'];
- else % if i is even
- tmp2 = [tmp2; squeeze(data(i,i+1,:))'];
- end
- end
- figure;
- subplot(1,2,1); % Nearby pairs
- plot(tvect, tmp1);
- set(gca,'Ylim',[min([tmp1(:);tmp2(:)])*.9 max([tmp1(:);tmp2(:)])*1.1]);
- subplot(1,2,2); % Distant pairs
- plot(tvect, tmp2);
- set(gca,'Ylim',[min([tmp1(:);tmp2(:)])*.9 max([tmp1(:);tmp2(:)])*1.1]);
-
- %%
- % Cranking through a bunch of frequencies to find where we get the best
- % signal in the LFP-TA and LFP-TV
- LFPsamprate = WN.sum.analog.storeRates{1}; % Assuming all channels are sampled at the same rate
- whichADchan = 4; % Pick something good
- nfreqs = 10;
- lowfreq = 30;
- highfreq = 300;
- %highfreq = LFPsamprate/2; % Nyquist
- ncyclesper6sigma = 3; % Controls the bandwidth
- ntrials = size(WN.trial,1);
- nstixperside = WN.sum.exptParams.nstixperside; %get number of stixels per side
- msperframe = 1000/WN.sum.exptParams.framerate; %calculate msec per frame
- secsperframe = 1/WN.sum.exptParams.framerate;
- gammaTable = WN.sum.exptParams.gamma_table; %get gamma_table
- gammaTable = reshape(gammaTable,length(gammaTable)/3,3); %reshapse gamma_table into three columns
- invgammaTable = InvertGamma(gammaTable,1); %invert gamma_table
- ngammasteps = size(invgammaTable,1); %get number of rows of gamma_table (65536)
- seedidx = strcmp(WN.sum.trialFields(1,:),'seed'); %get seed index from trialFields
- nframesidx = strcmp(WN.sum.trialFields(1,:),'num_frames'); %get nframes index from trialFields
- stimonidx = strcmp(WN.sum.trialFields(1,:),'stim_on'); %get stimon index from trialFields
- muidxs = [find(strcmp(WN.sum.trialFields(1,:),'mu1')),... %get mu indices into vector from trialFields
- find(strcmp(WN.sum.trialFields(1,:),'mu2')),...
- find(strcmp(WN.sum.trialFields(1,:),'mu3'))];
- sigmaidxs = [find(strcmp(WN.sum.trialFields(1,:),'sigma1')),... %get sigma indices into vector from trialFields
- find(strcmp(WN.sum.trialFields(1,:),'sigma2')),...
- find(strcmp(WN.sum.trialFields(1,:),'sigma3'))];
- LFPstarttimes = [WN.ras{:,strcmp(WN.sum.rasterCells,'anlgStartTime')}];
-
- for centerfreq = logspace(log10(lowfreq),log10(highfreq),nfreqs); % Hz
- cyclespersample = centerfreq./LFPsamprate;
- nsamplesper6sigma = ceil(ncyclesper6sigma/cyclespersample);
- ncycles = nsamplesper6sigma*cyclespersample;
- filtkernel1 = normpdf(linspace(-3,3,nsamplesper6sigma),0,1).*cos(linspace(0,2*pi*ncycles,nsamplesper6sigma));
- filtkernel2 = normpdf(linspace(-3,3,nsamplesper6sigma),0,1).*sin(linspace(0,2*pi*ncycles,nsamplesper6sigma));
- filtkernel1 = filtkernel1./norm(filtkernel1);
- filtkernel2 = filtkernel2./norm(filtkernel2);
-
- feval('STAmex','init', {nstixperside^2, 3, maxT});
- % This code was largely lifted from getWhtnsStats
- hwait = waitbar(0,'Finding stimuli...'); %display progress bar
- L = WN.trial(:,noisetypeidx)==1; % 1 = gun noise
- for i = 1:ntrials %get values and insert into given column
- if (WN.trial(i,noisetypeidx)==1)
- if (sum(L) > sum(~L))
- noisetype = 1;
- else
- continue
- end
- end
-
- if (WN.trial(i,noisetypeidx)==2)
- if (sum(~L) > sum(L))
- noisetype = 2;
- else
- continue
- end
- end
-
- seed = WN.trial(i,seedidx);
- nframes = WN.trial(i,nframesidx);
- mu = WN.trial(i,muidxs)/1000;
- sigma = WN.trial(i,sigmaidxs)/1000;
-
- if (noisetype == 1) % Gaussian gun
- x = linspace(WN.sum.exptParams.gauss_locut/1000, WN.sum.exptParams.gauss_hicut/1000, ngammasteps);%dividing gauss by 1000 and making equal intervals so that there are 65536 values
- for gun = 1:3
- invnormcdf(:,gun) = norminv(x)*sigma(gun)+mu(gun);
- end
- randnums = getEJrandnums(3*nstixperside^2*nframes, seed);
- randnums = reshape(randnums, [nstixperside^2*3, nframes]);
- for j = 1:3
- idxs = [1:nstixperside^2]+nstixperside^2*(j-1);
- randnums(idxs,:) = reshape(invnormcdf(randnums(idxs,:)+1,j),[length(idxs),nframes]);
- end
- elseif (noisetype == 2) % Binary cone
- colordirlms = sign(fullfact([2,2,2])-1.5);
- randnums = getEJrandnums(nstixperside^2*nframes, seed);
- randnums = reshape(randnums, [nstixperside^2, nframes]);
- randnums = mod(randnums, size(colordirlms,1))+1;
- tmp = colordirlms(randnums,:);
- tmp = reshape(tmp, [nstixperside^2 nframes 3]);
- tmp = permute(tmp,[1 3 2]);
- randnums = reshape(tmp,[nstixperside^2*3 nframes]);
- end
-
- t_stimon = WN.trial(i, stimonidx);
-
- LFP = WN.ras{i,strcmp(WN.sum.rasterCells,['AD',num2str(16+whichADchan)])};
- filteredLFP = conv(LFP, filtkernel1,'same').^2+conv(LFP, filtkernel2,'same').^2;
- LFPtimes = LFPstarttimes(i)+[0:length(filteredLFP)-1]/LFPsamprate; % sec
- frametimes = t_stimon+(linspace(0, nframes*secsperframe, nframes)+(secsperframe/2)'); % sec
- rectLFP = interp1(LFPtimes, abs(filteredLFP), frametimes); % This is probably not the best
- % plot(LFPtimes,abs(LFP)); hold on; plot(frametimes,rectLFP,'m-');
- feval('STAmex',randnums(:),rectLFP);
- if (ishandle(hwait))
- waitbar(i/ntrials,hwait);
- else
- break;
- end
- end
- close(hwait);
- out = feval('STAmex','return');
- eval('clear STAmex');
-
- % Now doing some plotting
- STAs = out{1};
- STVs = out{2};
- STAlims = [-max(abs(STAs(:))) max(abs(STAs(:)))];
- STVlims = [min(STVs(:)) max(STVs(:))];
- STAs = (STAs-STAlims(2))/(STAlims(1)-STAlims(2));
- STVs = (STVs-STVlims(2))/(STVlims(1)-STVlims(2));
-
- figure; colormap(gray(255));
- set(gcf,'Name',['LFP ',num2str(whichADchan),' ',num2str(centerfreq),' Hz']);
- for i = 1:size(STAs,2)
- % First STAs
- mat = reshape(STAs(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- im = mat(:,:,j);
- subplot(6,size(STAs,2),(j-1)*size(STAs,2)+i);
- image(im*255);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- % Then STVs
- mat = reshape(STVs(:,i),[nstixperside nstixperside 3]);
- for j = 1:3
- im = mat(:,:,j);
- subplot(6,size(STAs,2),(j-1+3)*size(STAs,2)+i);
- image(im*255);
- axis image;
- set(gca,'XTick',[],'YTick',[]);
- end
- end
- end
-
- %%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % The stuff below is basically obsolete. It's stuff I was working on just
- % to make sure that could do offline analysis at all.
- %
- % First, looking at eye movements
- figure;
- for i = 1:ntrials
- hep = WN.ras{i,hepidx}*4096/400;
- vep = WN.ras{i,vepidx}*4096/400;
- nsamps = length(hep);
- t = [0:1:nsamps-1]*eyesampperiod+eyestart_t(i);
- subplot(3,1,1); hold on;
- plot(hep,vep);
- % subplot(3,1,2); hold on;
- % plot(t-stimon_t(i), hep,'m-');
- % plot(t-stimon_t(i), vep,'g-');
- end
-
- subplot(3,1,1)
- axis square
- % Theory on the eye movements:
- % in nex2stro.m we multiply the raw AD values by 10/2^12 to get to
- % "millivolts". But we know that on REX 40 AD values is 1 degree and we
- % know that REX and Plexon have similar A/D cards. Both Rex and Plexon
- % appear to want A/D signals from -5V to +5V (they both saturate outside of
- % this range).
- %
- % 10 V/102.4 deg * 4096 steps/10 V = 40 steps/deg
-
- %%
- % Select a pixel for analysis
- [x,y] = ginput;
- whichpix = [round(x) round(y)];
- idx = (whichpix(1)-1)*nstixperside+whichpix(2);
- STA_t = STAs([idx, idx+nstixperside^2, idx+2*nstixperside^2],:);
- figure; axes; hold on;
- plot(STA_t(1,:),'r-','LineWidth',2);
- plot(STA_t(2,:),'g-','LineWidth',2);
- plot(STA_t(3,:),'b-','LineWidth',2);
-
- %%
- % Looking at the gun intensities of STA and PCs
- % at a single frame
- colors = {'red','green','blue'};
- rowidxs = reshape([1:3*nstixperside^2],[nstixperside^2, 3]);
- whichframe = input('Which frame would you like to see?');
- figure;
- subplot(6,1,1); hold on;
- maxval = max(abs([STAs(:,whichframe);PCs(:,1,whichframe)]));
- for gun = 1:3
- plot(STAs(rowidxs(:,gun),whichframe),'Color',colors{gun});
- end
- set(gca,'Ylim', [-1.5 1.5]*maxval);
- for i = 1:4
- subplot(6,1,1+i); hold on;
- for gun = 1:3
- plot(PCs(rowidxs(:,gun),i,whichframe),'Color',colors{gun});
- end
- set(gca,'Ylim', [-1.5 1.5]*maxval);
- end
- STC = reshape(STCs(:,whichframe),[sqrt(length(STCs(:,whichframe))),sqrt(length(STCs(:,whichframe)))]);
- STV = reshape(diag(STC),[nstixperside nstixperside 3]);
- STV = (STV-mean(STV(:)));
- subplot(6,1,6); hold on;
- for gun = 1:3
- plot(reshape(STV(:,:,gun),[nstixperside^2 1]),'Color',colors{gun});
- end
-
- %% Looking trial by trial at eye movements, spikes, events, etc.
- figure; axes; hold on;
- framerate = WN.sum.exptParams.framerate;
- predstimoff_t = stimon_t+WN.trial(:,nframesidx)/framerate;
- for i = 1:ntrials
- timeanchor = stimon_t(i);
- hep = WN.ras{i,hepidx}*4096/400;
- vep = WN.ras{i,vepidx}*4096/400;
- spiketimes = WN.ras{i,spikeidx}-timeanchor;
- nspikes = length(spiketimes);
- nsamps = length(hep);
- t = [0:1:nsamps-1]*eyesampperiod+eyestart_t(i)-timeanchor;
- minep = min([hep; vep]);
-
- plot(t, hep,'m-');
- plot(t, vep…
Large files files are truncated, but you can click here to view the full file