Skip to content

Commit

Permalink
added preprocessing functions for eye, accel, and photodiode; clarifi…
Browse files Browse the repository at this point in the history
…ed runAnalysis calling syntax
  • Loading branch information
sserene committed Jul 31, 2017
1 parent 2456f01 commit ea1513c
Show file tree
Hide file tree
Showing 6 changed files with 336 additions and 42 deletions.
57 changes: 44 additions & 13 deletions buildAnalysisParamFile.m
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,8 @@
%note: use Blackrock indexing for unitsToUnsort and unitsToDiscard, so unsorted is 0, first defined unit is 1, etc.
ephysParams.unitsToUnsort = {[],[],[]}; %these units will be re-grouped with u0
ephysParams.unitsToDiscard = {[],[],[]}; %these units will be considered noise and discarded
ephysParams.spikeWaveformPca = 0;
ephysParams.plotSpikeWaveforms = 0; %0, 1 to build then close, 2 to build and leave open
ephysParams.spikeWaveformPca = 1;
ephysParams.plotSpikeWaveforms = 2; %0, 1 to build then close, 2 to build and leave open
ephysParams.shiftSpikeWaveforms = 0;
% see http://www.mathworks.com/help/signal/examples/filter-design-gallery.html
hp1Hz = designfilt('highpassiir', 'FilterOrder',8,'PassbandFrequency',1, ...
Expand All @@ -75,7 +75,7 @@
analogInParams.needAnalogIn = 1;
analogInParams.analogInChannels = [138,139,140];
analogInParams.channelNames = {'eyeX','eyeY','eyeD'};
analogInParams.analogInChannelScaleBy = [1 1 1]; %converts raw values to microvolts
analogInParams.analogInChannelScaleBy = [5/32764 5/32764 5/32764]; %converts raw values to volts
analogInParams.decimateFactorPass1 = 1; %note: product of the two decimate factors should be 30, if 1 khz samples desired
analogInParams.decimateFactorPass2 = 1;
analogInParams.samPerMS = 1; %THIS IS AFTER DECIMATION, and applies to analogIn (should be raw rate/productOfDecimateFactors)
Expand All @@ -85,25 +85,42 @@
analogInParams.filters = {0,0,0};%{butter200Hz_v1;butter200Hz_v1;butter200Hz_v1}; %filter channel i if filters{i} is digital filter or 1x2 numeric array
analogInParams.plotFilteredSignal = 1; %#ok

photodiodeParams.needPhotodiode = 0;
photodiodeParams.channels = [1;2]; %#ok

% parameters preprocessLogFile, see function for details
stimSyncParams.usePhotodiode = 0; %#ok
%
eyeCalParams.method = 'auto';
eyeCalParams.needEyeCal = 1;
eyeCalParams.method = 'zeroEachFixation';
eyeCalParams.makePlots = 0;
eyeCalParams.eyeXChannelInd = 1;
eyeCalParams.eyeYChannelInd = 2;
eyeCalParams.eyeDChannelInd = 3;
eyeCalParams.gainX = 112;
eyeCalParams.gainY = 107;
eyeCalParams.flipX = 1;
eyeCalParams.flipY = 1; %#ok
eyeCalParams.flipY = 1;
eyeCalParams.offsetX = -6.4;
eyeCalParams.offsetY = -5.6;
eyeCalParams.minFixZeroTime = 1000; %#ok

accelParams.needAccelCal = 0;
accelParams.accelChannels = {[4;5;6]};
accelParams.channelGains = {[1/.666 1/.666 1/.666]};
accelParams.calMethods = {'hardcode'}; %other option is 'calFile'; calibration method
accelParams.calFiles = {''}; %if method is 'calFile', an ns2 filename

% parameters for excludeStimuli, see function for details
excludeStimParams.fixPre = 0; %ms
excludeStimParams.fixPost = 0; %ms
excludeStimParams.fixPre = 100; %ms
excludeStimParams.fixPost = 100; %ms
excludeStimParams.flashPre = 0; %ms
excludeStimParams.flashPost = 0; %ms
excludeStimParams.juicePre = 0; % optional, ms
excludeStimParams.juicePost = 0; % optional, ms
excludeStimParams.DEBUG = 0; % makes exclusion criterion plots if true
% additional optional excludeStimParams: accel1, accel2, minStimDur (ms)


psthParams.psthPre = 100; % if e.g. +200, then start psth 200ms before trial onset;
psthParams.psthImDur = 0; % only need to set this for variable length stim runs; else, comes from log file
psthParams.psthPost = 300;
Expand Down Expand Up @@ -188,7 +205,9 @@
plotSwitch.calcLatencyRF = 0;
plotSwitch.calcEvokedPowerRF = 0;
plotSwitch.evokedPsthMuaMultiCh = 0;
plotSwitch.evokedByCategory = 0;
plotSwitch.evokedByCategory = 1;
plotSwitch.analogInByItem = 1;
plotSwitch.analogInDerivativesByItem = 1;
plotSwitch.colorPsthEvoked = 0;
plotSwitch.linePsthEvoked = 0;
plotSwitch.runSummary = 0;
Expand All @@ -205,8 +224,8 @@
plotSwitch.spikeSpectraByCategory = 1;
plotSwitch.SpikeSpectraTfByImage = 1;
plotSwitch.lfpSpectraTfByImage = 1;
plotSwitch.tfSpectraByCategory = 1;
plotSwitch.tfErrs = 1; %#ok
plotSwitch.tfSpectraByCategory = 0;
plotSwitch.tfErrs = 0; %#ok

%%%% note: all analysisGroups cell arrays are nx1, NOT 1xn
analysisGroups.selectivityIndex.groups = {{'face';'nonface'},{'face';'object'},{'face';'body'}};
Expand All @@ -223,6 +242,18 @@
analysisGroups.evokedPotentials.names = {'fobPlus'};
analysisGroups.evokedPotentials.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
%
analysisGroups.analogInPotentials.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
analysisGroups.analogInPotentials.channels = {[1; 2]};
analysisGroups.analogInPotentials.names = {'eyePositions,fobPlus'};
analysisGroups.analogInPotentials.units = {'degrees visual angle'};
analysisGroups.analogInPotentials.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
%
analysisGroups.analogInDerivatives.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
analysisGroups.analogInDerivatives.channels = {[1; 2]};
analysisGroups.analogInDerivatives.names = {'eyeVelocity,fobPlus'};
analysisGroups.analogInDerivatives.units = {'degrees visual angle/sec'};
analysisGroups.analogInDerivatives.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
%
analysisGroups.colorPsthEvoked.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
analysisGroups.colorPsthEvoked.names = {'fobPlus'};
analysisGroups.colorPsthEvoked.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
Expand Down Expand Up @@ -261,11 +292,11 @@
calcSwitch.faceSelectIndexLate = 0;
calcSwitch.inducedTrialMagnitudeCorrection = 0;
calcSwitch.evokedSpectra = 1;
calcSwitch.inducedSpectra = 0;
calcSwitch.inducedSpectra = 1;
calcSwitch.evokedImageTF = 0;
calcSwitch.inducedImageTF = 0;
calcSwitch.evokedCatTF = 1;
calcSwitch.inducedCatTF = 0;
calcSwitch.inducedCatTF = 1;
calcSwitch.meanEvokedTF = 1;
calcSwitch.trialMeanSpectra = 0;
calcSwitch.coherenceByCategory = 1;
Expand Down
17 changes: 17 additions & 0 deletions preprocessAccelSignals.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
function [ analogInData ] = preprocessAccelSignals( analogInData,params )
%Applies volts to m/s^2 conversion for accelerometer signals.
% - Optionally: applies coordinate system conversion based on calibration runs (not yet implemented)
% Note: supports multiple accelerometers
if ~params.needAccelCal
return
end
assert(~any(strcmp(params.calMethods,'calFile')),'Requested file-based accelerometer calibration; not yet implemented');


for accel_i = 1:length(params.accelChannels)
analogInData(params.accelChannels{accel_i}) = analogInData(params.accelChannels{accel_i}) - mean(analogInData(params.accelChannels{accel_i}),2);
analogInData(params.accelChannels{accel_i}) = params.channelGains{accel_i} .* analogInData(params.accelChannels{accel_i});
end

end

91 changes: 91 additions & 0 deletions preprocessEyeSignals.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
function [ analogInData ] = preprocessEyeSignals( analogInData,taskData,params )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
if ~params.needEyeCal
return
end
calibMethod = params.method;
gainX = params.gainX;
gainY = params.gainY;
flipY = params.flipY;
flipX = params.flipX;
if strcmp(calibMethod,'hardcodeZero')
offsetX = params.offsetX;
offsetY = params.offsetY;
end
if strcmp(calibMethod, 'zeroEachFixation')
minFixZeroTime = params.minFixZeroTime;
end

eyeX = analogInData(params.eyeXChannelInd,:);
eyeY = analogInData(params.eyeYChannelInd,:);
eyeD = analogInData(params.eyeDChannelInd,:);

fixInInds = zeros(size(eyeX));
transitionInInds = zeros(size(eyeX));
transitionOutInds = zeros(size(eyeX));
for fix_i = 1:length(taskData.fixationInTimes)
assert(all(taskData.fixationOutTimes > taskData.fixationInTimes));
fixInInds(ceil(taskData.fixationInTimes(fix_i)):(floor(taskData.fixationOutTimes(fix_i))-10)) = 1;
transitionInInds(ceil(taskData.fixationInTimes(fix_i))) = 1;
transitionOutInds(floor(taskData.fixationOutTimes(fix_i))-10) = 1;
end
eyeX = gainX*eyeX*(-1*flipX + ~flipX);
eyeY = gainY*eyeY*(-1*flipY + ~flipY);
if strcmp(calibMethod,'autoZeroSingle')
eyeX = eyeX - mean(eyeX(fixInInds == 1));
eyeY = eyeY - mean(eyeY(fixInInds == 1));
end
if strcmp(calibMethod,'zeroEachFixation')
for fix_i = 1:length(taskData.fixationInTimes)
if fix_i == 1 || (taskData.fixationOutTimes(fix_i) - taskData.fixationInTimes(fix_i)) > minFixZeroTime
offsetX = mean(eyeX(ceil(taskData.fixationInTimes(fix_i)):(floor(taskData.fixationOutTimes(fix_i))-10)));
offsetY = mean(eyeY(ceil(taskData.fixationInTimes(fix_i)):(floor(taskData.fixationOutTimes(fix_i))-10)));
end
if fix_i < length(taskData.fixationInTimes)
eyeX(ceil(taskData.fixationInTimes(fix_i)):floor(taskData.fixationInTimes(fix_i+1))) = ...
eyeX(ceil(taskData.fixationInTimes(fix_i)):floor(taskData.fixationInTimes(fix_i+1))) - offsetX;
eyeY(ceil(taskData.fixationInTimes(fix_i)):floor(taskData.fixationInTimes(fix_i+1))) = ...
eyeY(ceil(taskData.fixationInTimes(fix_i)):floor(taskData.fixationInTimes(fix_i+1))) - offsetY;
else
eyeX(ceil(taskData.fixationInTimes(fix_i)):end) = eyeX(ceil(taskData.fixationInTimes(fix_i)):end) - offsetX;
eyeY(ceil(taskData.fixationInTimes(fix_i)):end) = eyeY(ceil(taskData.fixationInTimes(fix_i)):end) - offsetY;
end
end
end
if strcmp(calibMethod,'hardcodeZero')
eyeX = eyeX - offsetX;
eyeY = eyeY - offsetY;
end
if strcmp(calibMethod, 'ninePoint')
error('ninePoint eye calibration not yet implemented');
end
eyeD = eyeD - mean(eyeD(fixInInds == 1));
fixInX = eyeX(fixInInds == 1);
fixInY = eyeY(fixInInds == 1);
fixInD = eyeD(fixInInds == 1);
if params.makePlots
figure()
scatter(fixInX,fixInY,1,'b');
hold on
scatter(eyeX(~fixInInds),eyeY(~fixInInds),1,'r')
figure()
scatter(eyeX(transitionInInds == 1),eyeY(transitionInInds == 1),20,'b');
hold on
scatter(eyeX(transitionOutInds == 1),eyeY(transitionOutInds == 1),20,'r');
figure()
histogram(fixInD,100,'Normalization','pdf');
hold on
histogram(eyeD(~fixInInds),100,'Normalization','pdf');
title('eyeD')
legend({'In','Out'});
figure()
a1 = subplot(1,2,1);
histogram2(fixInX',fixInY',[50 50],'Normalization','probability','DisplayStyle','tile','BinMethod','integer');
a2 = subplot(1,2,2);
histogram2(eyeX(~fixInInds)',eyeY(~fixInInds)',[50 50],'Normalization','probability','DisplayStyle','tile','BinMethod','integer');
end
analogInData(params.eyeXChannelInd,:) = eyeX;
analogInData(params.eyeYChannelInd,:) = eyeY;
end

10 changes: 10 additions & 0 deletions preprocessPhotodiodeStrobe.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
function [ analogInData ] = preprocessPhotodiodeStrobe( analogInData, params )
%preprocessPhotodiodeStrobe cleans a photodiode signal
% Detailed explanation goes here
if ~params.needPhotodiode
return
end
error('photodiode strobe preprocessing not yet implemented');

end

72 changes: 52 additions & 20 deletions processRun.m
Original file line number Diff line number Diff line change
@@ -1,11 +1,16 @@
function [ spikesByImage, spikesByCategory, lfpByImage, lfpByCategory, categoryList, picFiles ] = processRun( varargin )
%processRun is a top level function to analyze units, MUA, lfps, coupling, and RFs.
%processRun loads, preprocesses, and aligns task events, lfps, muas, units, and analog signals (eg eye, accelerometer,photodiodes)
% - handles single-channel and multi-channel sessions
% - relies on raw visiko (.log) and blackrock (.ns5) files
% - aligns visiko events to blackrock clock (preprocessLogFile.m)
% - excludes trials with broken fixation, fix spot flash, or (optionally)
% juice delivery (via exludeStimuli)
% - relies only on raw task log and physiology fiels: currently visiko (.log) and blackrock (.ns5,.ns2)
% - aligns task events to ephys system clock (preprocessLogFile.m)
% - excludes trials with broken fixation, fix spot flash, and (optionally)
% juice delivery and high head acceleration(via exludeStimuli)
% - decimates and (optionally) filters raw lfp data (default, to 1 kHz)
% - calibrates eye signals (via preprocessEyeSignals)
% - calibrates accelerometer signals (via processAccelSignals)
% - calculates spike isolation measures (via preprocessSpikes)
% - shows task event summary (via excludeTrials)
% - optionally calls an analysis routine; default is runAnalyses.m
% Inputs:
% - varargin can have the following forms:
% - empty (default assignments: buildAnalysisParamFile, runAnalyses
Expand Down Expand Up @@ -40,6 +45,7 @@
channelNames = ephysParams.channelNames;
spikeChannels = ephysParams.spikeChannels;
lfpChannels = ephysParams.lfpChannels;
analogInChannels = analogInParams.analogInChannels;
psthPre = psthParams.psthPre;
psthImDur = psthParams.psthImDur;
psthPost = psthParams.psthPost;
Expand All @@ -62,7 +68,13 @@
analogInData = preprocessAnalogIn(analogInFilename, analogInParams);
[spikesByChannel, taskTriggers, channelUnitNames] = preprocessSpikes(spikeFilename, ephysParams);
lfpData = preprocessLFP(lfpFilename, ephysParams);
analogInData = preprocessPhotodiodeStrobe(analogInData, photodiodeParams);
[ taskData, stimTiming ] = preprocessLogFile(taskFilename, taskTriggers, stimSyncParams); % load visual stimulus data and transform its timestamps to ephys clock reference
analogInData = preprocessEyeSignals(analogInData,taskData,eyeCalParams);
analogInData = preprocessAccelSignals(analogInData, accelParams);
taskDataAll = taskData;
taskData = excludeTrials( taskData, excludeStimParams); %exclude trials for lost fixation etc.

if stimTiming.shortest == stimTiming.longest
psthParams.psthImDur = stimTiming.shortest;
psthImDur = psthParams.psthImDur;
Expand All @@ -79,11 +91,6 @@
fprintf('Variable stimulus length run. Excluding trials shorter than %d\n',psthImDur);
end

[eyeX,eyeY,eyeD] = calibrateEyeSignals(analogInData(1,:),analogInData(2,:),analogInData(3,:),taskData,eyeCalParams);
return
taskDataAll = taskData;
taskData = excludeTrials( taskData, excludeStimParams); %exclude trials for lost fixation etc.

%sort trials by image and image category
tmp = load(stimParamsFilename); %loads variables paramArray, categoryLabels,pictureLabels
picCategories = tmp.paramArray;
Expand Down Expand Up @@ -153,8 +160,8 @@
% align LFP data by trial, sort by image and category, and possibly remove DC and linear components
lfpByImage = alignLFP(lfpData, onsetsByImage, lfpChannels, lfpAlignParams);
lfpByCategory = alignLFP(lfpData, onsetsByCategory, lfpChannels, lfpAlignParams);
analogInByImage = alignAnalogIn(analogInData, onsetsByImage, analogInChannels, analogInParams);
analogInByCategory = alignAnalogIn(analogInData, onsetsByCategory, analogInChannels, analogInParams);
analogInByImage = alignAnalogIn(analogInData, onsetsByImage, analogInChannels, lfpAlignParams);
analogInByCategory = alignAnalogIn(analogInData, onsetsByCategory, analogInChannels, lfpAlignParams);

for cat_i = 1:length(categoryList)
Output.VERBOSE(categoryList{cat_i});
Expand All @@ -168,16 +175,41 @@
'stimTiming', 'picCategories', 'onsetsByImage', 'onsetsByCategory')
end
end

runAnalysisInputs.analysisParamFilename = analysisParamFilename;
runAnalysisInputs.spikesByChannel = spikesByChannel;
runAnalysisInputs.lfpData = lfpData;
runAnalysisInputs.analogInData = analogInData;
runAnalysisInputs.taskData = taskData;
runAnalysisInputs.taskDataAll = taskDataAll;
runAnalysisInputs.psthImDur = psthImDur;
runAnalysisInputs.preAlign = preAlign;
runAnalysisInputs.postAlign = postAlign;
runAnalysisInputs.categoryList = categoryList;
runAnalysisInputs.pictureLabels = pictureLabels;
runAnalysisInputs.jumpsByImage = jumpsByImage;
runAnalysisInputs.spikesByImage = spikesByImage;
runAnalysisInputs.psthEmptyByImage = psthEmptyByImage;
runAnalysisInputs.spikesByCategory = spikesByCategory;
runAnalysisInputs.psthEmptyByCategory = psthEmptyByCategory;
runAnalysisInputs.spikesByImageForTF = spikesByImageForTF;
runAnalysisInputs.spikesByCategoryForTF = spikesByCategoryForTF;
runAnalysisInputs.lfpByImage = lfpByImage;
runAnalysisInputs.lfpByCategory = lfpByCategory;
runAnalysisInputs.analogInByImage = analogInByImage;
runAnalysisInputs.analogInByCategory = analogInByCategory;
runAnalysisInputs.channelUnitNames = channelUnitNames;
runAnalysisInputs.stimTiming = stimTiming;
runAnalysisInputs.picCategories = picCategories;
runAnalysisInputs.onsetsByImage = onsetsByImage;
runAnalysisInputs.onsetsByCategory = onsetsByCategory;



if nargin == 0 || (nargin == 2 && strcmp(varargin{1},'paramBuilder')) || (nargin == 2 && strcmp(varargin{1},'preprocessed'))
runAnalyses( analysisParamFilename, spikesByChannel, lfpData, analogInData, taskData, taskDataAll, psthImDur, preAlign, postAlign,...
categoryList, pictureLabels, jumpsByImage, spikesByImage, psthEmptyByImage, spikesByCategory, psthEmptyByCategory,...
spikesByImageForTF, spikesByCategoryForTF, lfpByImage, lfpByCategory, analogInByImage, analogInByCategory, channelUnitNames, ...
stimTiming, picCategories, onsetsByImage, onsetsByCategory);
runAnalyses(runAnalysisInputs);
else
feval(varargin{end},analysisParamFilename, spikesByChannel, lfpData, analogInData, taskData, taskDataAll, psthImDur, preAlign, postAlign,...
categoryList, pictureLabels, jumpsByImage, spikesByImage, psthEmptyByImage, spikesByCategory, psthEmptyByCategory,...
spikesByImageForTF, spikesByCategoryForTF, lfpByImage, lfpByCategory, analogInByImage, analogInByCategory, channelUnitNames, ...
stimTiming, picCategories, onsetsByImage, onsetsByCategory);
feval(varargin{end},runAnalysisInputs);
end
end

Loading

0 comments on commit ea1513c

Please sign in to comment.