diff --git a/analysis/lfp/bz_PowerSpectrumSlope.m b/analysis/lfp/bz_PowerSpectrumSlope.m index ddd5a106..5835b5c4 100755 --- a/analysis/lfp/bz_PowerSpectrumSlope.m +++ b/analysis/lfp/bz_PowerSpectrumSlope.m @@ -11,12 +11,14 @@ % % (optional) % 'frange' (default: [4 100]) +% 'channels' subset of channels to calculate PowerSpectrumSlope +% (default: all) % 'showfig' true/false - show a summary figure of the results % (default:false) -% 'saveMat' put your basePath here to save +% 'saveMat' put your basePath here to save/load % baseName.PowerSpectrumSlope.lfp.mat (default: false) -% 'channels' subset of channels to calculate PowerSpectrumSlope -% (default: all) +% 'Redetect' (default: false) to force redetection even if saved +% file exists % % %DLevenstein 2018 @@ -26,13 +28,27 @@ addParameter(p,'saveMat',false) addParameter(p,'channels',[]) addParameter(p,'frange',[4 100]) +addParameter(p,'Redetect',false) parse(p,varargin{:}) SHOWFIG = p.Results.showfig; saveMat = p.Results.saveMat; channels = p.Results.channels; frange = p.Results.frange; +REDETECT = p.Results.Redetect; +%% +if saveMat + basePath = saveMat; + baseName = bz_BasenameFromBasepath(basePath); + savename = fullfile(basePath,[baseName,'.PowerSpectrumSlope.lfp.mat']); + + if exist(savename,'file') && ~REDETECT + load(savename) + return + end +end + %% For multiple lfp channels if ~isempty(channels) usechans = ismember(lfp.channels,channels); @@ -105,10 +121,7 @@ specslope.channels = lfp.channels; if saveMat - basePath = saveMat; - baseName = bz_BasenameFromBasepath(basePath); - savename = fullfile(basePath,[baseName,'.PowerSpectrumSlope.lfp.mat']); - save(savename,'specslope'); + save(savename,'specslope','spec'); end %% Figure @@ -131,7 +144,7 @@ ylabel('f (Hz)') axis xy xlim(bigsamplewin) - set(gca,'XTickLabel',[]) + bz_ScaleBar('s') subplot(8,1,3) plot(lfp.timestamps,lfp.data,'k') axis tight @@ -145,8 +158,8 @@ axis tight box off xlim(bigsamplewin) - ylabel('PSS');xlabel('Time (s)') - + set(gca,'XTickLabel',[]) + subplot(6,2,7) hist(specslope.data,10) box off @@ -168,8 +181,10 @@ box off xlim(exwin(bb,:)');ylim(lfprange) set(gca,'XTickLabel',[]) + set(gca,'YTicks',[]) ylabel('LFP') end + bz_ScaleBar('s') if saveMat figfolder = [basePath,filesep,'DetectionFigures']; diff --git a/detectors/detectBehavior/detect_pupilDilation/GetPupilDilation.m b/detectors/detectBehavior/detect_pupilDilation/GetPupilDilation.m index 38377497..0eb13be9 100755 --- a/detectors/detectBehavior/detect_pupilDilation/GetPupilDilation.m +++ b/detectors/detectBehavior/detect_pupilDilation/GetPupilDilation.m @@ -418,7 +418,7 @@ pupildilation.pupilxy = pupcoords; pupildilation.detectorname = 'GetPupilDilation'; -pupildilation.detectiondate = today('datetime'); +pupildilation.detectiondate = datetime('today'); pupildilation.detectorparms.mask = eyemask; pupildilation.detectorparms.intensitythresh = intensitythresh; pupildilation.detectorparms.unstablewindow = unstablewindow; diff --git a/detectors/detectBehavior/detect_whisking/GetWhiskFromEMG.m b/detectors/detectBehavior/detect_whisking/GetWhiskFromEMG.m index 4470d2b3..d3219099 100644 --- a/detectors/detectBehavior/detect_whisking/GetWhiskFromEMG.m +++ b/detectors/detectBehavior/detect_whisking/GetWhiskFromEMG.m @@ -20,11 +20,11 @@ % 'minNWh' Min. nonwhisking duration (s) (default: 0.1) % 'whiskmerge' Min. interwhisking duration (s) (default: 0.1) % 'NWhmerge' Min. internonwhisking duration (s) (default: 0.02) -% (options - not yet implemented) % 'showfig' show a detection figure (default: true) % 'saveMat' save the output in a baseName.EMGwhisk.behavior.mat file -% 'PulseChannel' -% 'EMGChannel' +% 'PulseChannel' channel with camera pulses (default: 1) +% 'EMGChannel' channel with EMG/Piezo signal (default: 2) +% 'analysis' 'EMG' or 'touch' (default: 'EMG') % % Note: 'auto' threshold detection gradient descents to a local trough in % the smoothed EMG envelope from an initial guess of 0.5 @@ -34,12 +34,9 @@ % 'saveMat' Creates file: % basePath/baseName.EMGwhisk.behavior.mat % -% -% %To Add: option for no camera pulses, which just align time to intan/pupil % -% -%DLevenstein 2017 +%DLevenstein, WMunoz 2017 %% DEV % basePath = '/mnt/proraidDL/Database/WMProbeData/'; % baseName = 'Layers_LFP_Test02_170323_151411'; @@ -54,6 +51,9 @@ addParameter(p,'NWhmerge',0.02,@isnumeric) addParameter(p,'showfig',true,@islogical) addParameter(p,'saveMat',true,@islogical) +addParameter(p,'PulseChannel',1) +addParameter(p,'EMGChannel',2) +addParameter(p,'EMGanalysis','EMG',@islogical) parse(p,varargin{:}) SHOWFIG = p.Results.showfig; @@ -65,7 +65,9 @@ EMGparms.minNWh = p.Results.minNWh; EMGparms.whiskmerge = p.Results.whiskmerge; EMGparms.NWhmerge = p.Results.NWhmerge; - +EMGparms.PulseChannel = p.Results.PulseChannel; +EMGparms.EMGChannel = p.Results.EMGChannel; +EMGanalysis = p.Results.EMGanalysis; if ~exist('basePath','var') basePath = pwd; @@ -74,7 +76,14 @@ %% abfname = fullfile(basePath,[baseName,'.abf']); analogName = fullfile(basePath,['analogin.dat']); -savefile = fullfile(basePath,[baseName,'.EMGwhisk.behavior.mat']); + +switch EMGanalysis + case 'EMG' + savefile = fullfile(basePath,[baseName,'.EMGwhisk.behavior.mat']); + case 'touch' + savefile = fullfile(basePath,[baseName,'.Piezotouch.behavior.mat']); +end + figfolder = fullfile(basePath,'DetectionFigures'); if ~exist(abfname,'file') @@ -86,10 +95,8 @@ %% Clampex File -timechan = 1; -emgchan = 2; -sf_abf = 20000; %Sampling Frequency of the .abf file +sf_abf = 20000; %Sampling Frequency of the .abf file [abffile,si,file_info] = abfload(abfname); @@ -101,8 +108,8 @@ % %Replace this with prompt file_info.recChNames % end -pulse_abf = abffile(:,timechan); -EMG = abffile(:,emgchan); +pulse_abf = abffile(:,EMGparms.PulseChannel); +EMG = abffile(:,EMGparms.EMGChannel); t_abf = [1:length(EMG)]'./sf_abf; @@ -149,12 +156,17 @@ %plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g') %plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r') axis tight - xlabel('EMG Envelope (modZ)'); + xlabel('EMG/Piezo Envelope (modZ)'); xlim([-1.5 max(log10(EMGsm))]) LogScale('x',10) - NiceSave('FAILEDWhiskingDetection',figfolder,baseName) + switch EMGanalysis + case 'EMG' + NiceSave('FAILEDWhiskingDetection',figfolder,baseName) + case 'touch' + NiceSave('FAILEDTouchDetection',figfolder,baseName) + end EMGwhisk = []; return end @@ -180,10 +192,15 @@ axis tight xlim([-1.5 max(log10(EMGsm))]) - xlabel('EMG Envelope (modZ)'); + xlabel('EMG/Piezo Envelope (modZ)'); LogScale('x',10) - NiceSave('FAILEDWhiskingDetection',figfolder,baseName) + switch EMGanalysis + case 'EMG' + NiceSave('FAILEDWhiskingDetection',figfolder,baseName) + case 'touch' + NiceSave('FAILEDTouchDetection',figfolder,baseName) + end EMGwhisk = []; return end @@ -192,7 +209,6 @@ end end %% - % figure % bar(EMGbins,EMGhist) % hold on @@ -205,8 +221,6 @@ % LogScale('x',10) % xlim([-1.5 max(log10(EMGsm))]) - - %% Identify Whisking on/offsets: EMG envelope crosses threshold wh_thresh = EMGsm > EMGparms.Whthreshold; wh_on = find(wh_thresh(2:end)>wh_thresh(1:end-1))+1; %whisking onsets (si) @@ -240,20 +254,15 @@ nwh_on = nwh_on./sf_down; nwh_off = nwh_off./sf_down; - - %Merge brief interruptions [ NWhints ] = MergeSeparatedInts( [nwh_on,nwh_off],EMGparms.NWhmerge ); [ Whints ] = MergeSeparatedInts( [wh_on,wh_off],EMGparms.whiskmerge ); - %Drop nonwhisk epochs smaller than a minimum [nwh_on,nwh_off] = MinEpochLength(NWhints(:,1),NWhints(:,2),EMGparms.minNWh,1); %Drop whisking epochs smaller than a minimum [wh_on,wh_off] = MinEpochLength(Whints(:,1),Whints(:,2),EMGparms.minwhisk,1); - - %% Durations Whdur = wh_off-wh_on; NWhdur = nwh_off-nwh_on; @@ -306,7 +315,6 @@ warning('Frame rate is not constant...') end - firstpulstime_lfp = pulset(1); %% Reset time to align to LFP @@ -317,81 +325,163 @@ %% Figure if SHOWFIG - figure - subplot(4,1,1) - plot(t_align,EMGz,'k') - - hold on - plot(t_align,EMGsm,'b','linewidth',2) - plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) - plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) - axis tight - ylim([-100 100]) - ylabel('EMG (modZ)'); - - subplot(4,1,2) - plot(t_align,EMGz,'k') - - hold on - plot(t_align,EMGsm,'b','linewidth',2) - plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) - plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) - xlim([100 160]) - ylim([-20 40]) - ylabel('EMG (modZ)'); - - subplot(4,2,6) - hist(log10(EMGsm),100) - hold on - plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g') - plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r') - - axis tight - xlabel('EMG Envelope (modZ)'); - LogScale('x',10) - xlim([-1.5 max(log10(EMGsm))]) - - subplot(4,2,8) - plot(durhist.bins,durhist.NWh,'r','linewidth',2) - hold on - plot(durhist.bins,durhist.Wh,'g','linewidth',2) - LogScale('x',10) - xlabel('Duration (s)') - ylabel('# Epochs') - legend('NWh','Wh') - - subplot(4,2,5) - plot(t_abf-firstpulstime_abf+firstpulstime_lfp,pulse_abf,'k') - hold on - plot(pulset,pulsethreshold_abf.*ones(size(pulset)),'r+') - xlim(firstpulstime_lfp+[-0.2 0.5]) - ylabel('Clampex Pulse Onset') - - subplot(4,2,7) - - plot(t_pulse,timepulses,'k') - hold on - plot(pulset,pulsethreshold.*ones(size(pulset)),'r+') - xlim(firstpulstime_lfp+[-0.2 0.5]) - ylabel('Intan Pulse Onset') - - NiceSave('WhiskingDetection',figfolder,baseName) - + switch EMGanalysis + case 'EMG' + + figure + subplot(4,1,1) + plot(t_align,EMGz,'k') + + hold on + plot(t_align,EMGsm,'b','linewidth',2) + plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) + plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) + axis tight + ylim([-100 100]) + ylabel('EMG (modZ)'); + + subplot(4,1,2) + plot(t_align,EMGz,'k') + + hold on + plot(t_align,EMGsm,'b','linewidth',2) + plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) + plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) + xlim([100 160]) + ylim([-20 40]) + ylabel('EMG (modZ)'); + + subplot(4,2,6) + hist(log10(EMGsm),100) + hold on + plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g') + plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r') + + axis tight + xlabel('EMG Envelope (modZ)'); + LogScale('x',10) + xlim([-1.5 max(log10(EMGsm))]) + + subplot(4,2,8) + plot(durhist.bins,durhist.NWh,'r','linewidth',2) + hold on + plot(durhist.bins,durhist.Wh,'g','linewidth',2) + LogScale('x',10) + xlabel('Duration (s)') + ylabel('# Epochs') + legend('NWh','Wh') + + subplot(4,2,5) + plot(t_abf-firstpulstime_abf+firstpulstime_lfp,pulse_abf,'k') + hold on + plot(pulset,pulsethreshold_abf.*ones(size(pulset)),'r+') + xlim(firstpulstime_lfp+[-0.2 0.5]) + ylabel('Clampex Pulse Onset') + + subplot(4,2,7) + + plot(t_pulse,timepulses,'k') + hold on + plot(pulset,pulsethreshold.*ones(size(pulset)),'r+') + xlim(firstpulstime_lfp+[-0.2 0.5]) + ylabel('Intan Pulse Onset') + + NiceSave('WhiskingDetection',figfolder,baseName) + + case 'touch' + figure + subplot(4,1,1) + plot(t_align,EMGz,'k') + + hold on + plot(t_align,EMGsm,'b','linewidth',2) + plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) + plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) + axis tight + ylim([-100 100]) + ylabel('Piezo (modZ)'); + + subplot(4,1,2) + plot(t_align,EMGz,'k') + + hold on + plot(t_align,EMGsm,'b','linewidth',2) + plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2) + plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2) + xlim([100 160]) + ylim([-20 40]) + ylabel('Piezo (modZ)'); + + subplot(4,2,6) + hist(log10(EMGsm),100) + hold on + plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g') + plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r') + + axis tight + xlabel('Piezo Envelope (modZ)'); + LogScale('x',10) + xlim([-1.5 max(log10(EMGsm))]) + + subplot(4,2,8) + plot(durhist.bins,durhist.NWh,'r','linewidth',2) + hold on + plot(durhist.bins,durhist.Wh,'g','linewidth',2) + LogScale('x',10) + xlabel('Duration (s)') + ylabel('# Epochs') + legend('NonTouch','Touch') + + subplot(4,2,5) + plot(t_abf-firstpulstime_abf+firstpulstime_lfp,pulse_abf,'k') + hold on + plot(pulset,pulsethreshold_abf.*ones(size(pulset)),'r+') + xlim(firstpulstime_lfp+[-0.2 0.5]) + ylabel('Clampex Pulse Onset') + + subplot(4,2,7) + + plot(t_pulse,timepulses,'k') + hold on + plot(pulset,pulsethreshold.*ones(size(pulset)),'r+') + xlim(firstpulstime_lfp+[-0.2 0.5]) + ylabel('Intan Pulse Onset') + + NiceSave('TouchDetection',figfolder,baseName) + end end %% -EMGwhisk.timestamps = t_align; -EMGwhisk.EMG = EMGz; -EMGwhisk.EMGsm = EMGsm; -EMGwhisk.ints.Wh = Whints; -EMGwhisk.ints.NWh = NWhints; -EMGwhisk.samplingRate = sf_down; -EMGwhisk.detectorparms = EMGparms; -EMGwhisk.detectorname = 'GetWhiskFromEMG'; -EMGwhisk.detectiondate = today('datetime'); - -if saveMat - save(savefile,'EMGwhisk') +switch EMGanalysis + case 'EMG' + EMGwhisk.timestamps = t_align; + EMGwhisk.EMG = EMGz; + EMGwhisk.EMGsm = EMGsm; + EMGwhisk.ints.Wh = Whints; + EMGwhisk.ints.NWh = NWhints; + EMGwhisk.samplingRate = sf_down; + EMGwhisk.detectorparms = EMGparms; + EMGwhisk.detectorname = 'GetWhiskFromEMG'; + EMGwhisk.detectiondate = today('datetime'); + case 'touch' + Piezotouch.timestamps = t_align; + Piezotouch.Piezo = EMGz; + Piezotouch.Piezosm = EMGsm; + Piezotouch.ints.Touch = Whints; + Piezotouch.ints.NoTouch = NWhints; + Piezotouch.samplingRate = sf_down; + Piezotouch.detectorparms = EMGparms; + Piezotouch.detectorname = 'GetWhiskFromEMG'; + Piezotouch.detectiondate = today('datetime'); +end + +if saveMat + switch EMGanalysis + case 'EMG' + save(savefile,'EMGwhisk') + case 'touch' + save(savefile,'Piezotouch') + end end end diff --git a/externalPackages/FMAToolbox/General/Restrict.m b/externalPackages/FMAToolbox/General/Restrict.m index af8cd838..38bf6704 100755 --- a/externalPackages/FMAToolbox/General/Restrict.m +++ b/externalPackages/FMAToolbox/General/Restrict.m @@ -49,6 +49,8 @@ end % Check parameters +intervals = double(intervals); +samples = double(samples); if ~isdmatrix(intervals) || size(intervals,2) ~= 2, error('Incorrect intervals (type ''help Restrict'' for details).'); end diff --git a/preprocessing/bz_ConcatenateBehavior.m b/preprocessing/bz_ConcatenateBehavior.m new file mode 100644 index 00000000..946c3866 --- /dev/null +++ b/preprocessing/bz_ConcatenateBehavior.m @@ -0,0 +1,130 @@ +function [ behavior ] = bz_ConcatenateBehavior( behaviorName, basePath, varargin) +%[ behavior ] = bz_ConcatenateBehavior( behaviorName, basePath) +%Concatenates buzcode behavior.mat files from subfolders, using the +%MergePoints from bz_ConcatenateDats. Automatically adjusts timestamps from +%each of the subfiles to reflect the timestamps in the merged .lfp file. +%Saves a new file: basePath/baseName.behaviorName.behavior.mat +% +%INPUTS +% behaviorName assumes that one or more subfolders have a file called +% subfolder.behaviorName.behavior.mat, following buzcode +% behavior.mat guidelines. +% Automatically adjusts timestamps stored in +% behaviorName.timestamps and subfields of +% behaviorName.ints. +% basePath upper level basePath into which you have already run +% bz_ConcatenateDats. +% Requires basePath/baseName.MergePoints.events.mat +% (options) +% 'timefields' any fields that contain timestamps that you'd like to +% adjust to match the concatenation MergePoints +% (note: this feature is not yet implemented, let me know +% if you need it and we'll add it to the function...) +% +%note: automatically adjusts behaviorName.timestamps and any fields in +%behaviorName.ints to account for merge points. Any additional fields can +%be indicated with optional 'timefields' input +%DLevenstein 2019 +%% input Parsing + +p = inputParser; +addParameter(p,'timefields',{}); +addParameter(p,'remerge',false); + +parse(p,varargin{:}) + +timefields = p.Results.timefields; +REMERGE = p.Results.remerge; + +%% DEV +% basePath = '/Users/dlevenstein/Desktop/180530_KO_EM1M3'; +% behaviorName = 'pupildiameter'; + +%% +baseName = bz_BasenameFromBasepath(basePath); +savename = fullfile(basePath,[baseName,'.',behaviorName,'.behavior.mat']); + +if exist(savename,'file') && ~REMERGE + RESULT = questdlg([savename,' already exists, overwrite? ',... + '(Use ''remerge'',true to avoid this message in the future)'],... + 'Yes','No'); + switch RESULT + case {'No','Cancel'} + return + end +end + +%% +%Load the mergepoints metadata file +MergePoints = bz_LoadEvents(basePath,'MergePoints'); +if isempty(MergePoints) + error(['No MergePoints.events.mat file found in ',basePath]) +end + + + +%% + +subfolders = cellfun(@(X) fullfile(basePath,X),MergePoints.foldernames,... + 'UniformOutput',false); +behaviorfilenames = cellfun(@(X,Y) fullfile(X,[Y,'.',behaviorName,'.behavior.mat']),... + subfolders,MergePoints.foldernames,'UniformOutput',false); + +%Check that all the subfolders, and which behavior files, exist +existsubfolders = cellfun(@(X) logical(exist(X,'dir')),subfolders); +existbehfiles = cellfun(@(X) logical(exist(X,'file')),behaviorfilenames); +existbehfiles = find(existbehfiles); + + +numsubfolders = length(subfolders); + +if any(~existbehfiles) + display(['Can''t find ',behaviorName,'.behavior.mat in ', MergePoints.foldernames(~existbehfiles),... + '. Merging from ',MergePoints.foldernames(existbehfiles)]) +end + +%% +offsets = MergePoints.timestamps(:,1); +subbeh = struct([]); +for ff = 1:length(existbehfiles) + thisidx = existbehfiles(ff); + newbeh = bz_LoadBehavior(subfolders{thisidx},behaviorName); + + % check that the fields match, if not add the new field + %(add notification here if not alinged) + [subbeh,newbeh] = bz_Matchfields(subbeh,newbeh,'add'); + subbeh(ff) = newbeh; + + % Required fields: timestamps, data. Should prompt user for other + % timestamp fields. + subbeh(ff).timestamps = subbeh(ff).timestamps+offsets(thisidx); + + if isfield(subbeh(ff),'ints') + intfields = fieldnames(subbeh(ff).ints); + for gg = 1:length(intfields) + subbeh(ff).ints.(intfields{gg}) = subbeh(ff).ints.(intfields{gg})+offsets(thisidx); + end + end + +end + + +%Collapse and concatenate the fields from the structures +behavior = bz_CollapseStruct(subbeh,'match','justcat',true); + +tol = 1e-4; %tolerance for difference in sample rates +if range(behavior.samplingRate)>tol + warning('Sampling Rates are not the same...... taking the average') +end +behavior.samplingRate = mean(behavior.samplingRate); + +behavior.MergeInfo.MergedFiles = behaviorfilenames(existbehfiles); +behavior.MergeInfo.MergedDate = datetime('today'); + + +%% Save the new behavior file +eval([behaviorName,'=behavior;']) +save(savename,behaviorName) + +end + diff --git a/preprocessing/bz_ConcatenateDats.m b/preprocessing/bz_ConcatenateDats.m index 55c7a10a..417d8386 100755 --- a/preprocessing/bz_ConcatenateDats.m +++ b/preprocessing/bz_ConcatenateDats.m @@ -351,232 +351,3 @@ function bz_ConcatenateDats(basepath,deleteoriginaldatsbool,sortFiles) end - - - - - - - - - - - - - - - - - - - -%%%% OLD BRENDON VERSION ASSUMING SESSIONMETADATA.MAT %%% -% -% % bz_ConcatenateDats - Concatenate raw .dat files found in a session folder -% % - for either intan type systems or amplipex -% % -% % ALGORITHM OUTLINE: looks for .dat files in a folder (or in subfolders) to -% % concatenate together. The concatenation happens via system commands -% % ("cat" command for linux/mac, "copy" command if windows/pc). Uses -% % different assumptions to find and recognize relevant .dats depending on -% % the acquisition system. -% % -% % REQUIREMENTS: Assumes you are in or pointed to a directory containing -% % subdirectories for various recording files from a single session. *It is -% % assumed that an earlier-acquired data file/folder will have a name that -% % is sorted alphanumerically earlier. Alphanumeric sorting order is -% % assumed to be the recording temporal sequence. -% % Works with acquisition systems: Intan and Amplipex - -% % 1) intan: wherein subfolders are inside the session folder. Each -% % subfolder contains simultaneously-recorded .dat files recorded for a -% % continuous period of time. Start/stop recording commands each create a -% % new folder. *It is assumed that the alphanumeric sorting of these -% % folders corresponds with their sequence in acquisiton time.* -% % These folders contain -% % - info.rhd files with metadata about the recording. -% % - amplifier.dat - int16 file with usually neural data from the -% % headstage -% % - auxiliary.dat (optional) - uint16 file from auxiliary channels on -% % the headstage - often accelerometer -% % - analogin.dat (optional) - uint16 file from analogin channels on -% % main board -% % - digitalin.dat (optional) - uint16 file recording all 16 digital -% % channels on main board -% % - time.dat - int32 file giving recording sample indicies (e.g. -% % 0,1,2,3...) for each sample recorded in other channels -% % - supply.dat - uint16 file showing voltage supplied to (?preamp?) -% % 2) Amplipex: wherein there are a series of sequientially named .dat -% % files, where the first might be called [basename]-01.dat, the next -% % [basename]-02.dat etc. *It is assumed that the alphanumeric sorting of -% % these folders corresponds with their sequence in acquisiton time.* -% % There are not assumed to be any sub-dat files, these .dats are assumed -% % to carry all information. -% % -% % USAGE -% % -% % bz_ConcatenateDats(basepath,deletedats) -% % -% % INPUTS -% % -% % basepath computer path to session folder. Defaults to -% % current folder if no input given -% % deleteoriginaldatsbool - boolean denoting whether to delete (1) or -% % not delete (0) original .dats after -% % concatenation. Default = 0. -% % -% % OUTPUT -% % Operates on files in specified folder. No output variable -% % -% % EXAMPLES -% % Can be called directly or via bz_PreprocessExtracellEphysSession.m -% % -% % Copyright (C) 2017 by Brendon Watson -% % -% % -% % %% Input and directory handling -% % if ~exist('basepath','var') -% % basepath = cd; -% % elseif isempty(basepath) -% % basepath = cd; -% % end -% % -% % basename = bz_BasenameFromBasepath(basepath); -% % -% % if ~exist('deletedats','var') -% % deleteoriginaldatsbool = 0; -% % end -% % -% % %% -% load(fullfile(basepath,[basename,'.SessionMetadata.mat'])); -% newdatpath = fullfile(basepath,[basename,'.dat']); -% switch SessionMetadata.ExtracellEphys.RecordingSystem -% case 'Intan' -% otherdattypes = {'analogin';'digitalin';'auxiliary';'time';'supply'}; -% bad_otherdattypes = []; -% for odidx = 1:length(otherdattypes) -% eval(['new' otherdattypes{odidx} 'path = fullfile(basepath,''' otherdattypes{odidx} '.dat'');']) -% end -% -% for fidx = 1:length(SessionMetadata.ExtracellEphys.Files.Names)% go through each folder in basepath to look for various .dats -% datpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},'amplifier.dat');%int16 -% %datbytes already saved in SessionMetadata -% -% for odidx = 1:length(otherdattypes)%loop through other .dat types found here -% eval([otherdattypes{odidx} 'datpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},''' otherdattypes{odidx} '.dat'');']) -% eval(['d = dir(' otherdattypes{odidx} 'datpaths{fidx});']) -% if isempty(d) -% bad_otherdattypes(odidx) = 1; -% else -% eval([otherdattypes{odidx} 'datsizes(fidx) = d(1).bytes;']) -% end -% end -% % -% % analogindatpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},'analogin.dat'); -% % d = dir(timedatpaths{fidx}); -% % analogindatsizes(fidx) = d(1).bytes; -% % -% % auxiliarydatpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},'auxiliary.dat'); -% % d = dir(timedatpaths{fidx}); -% % auxiliarydatsizes(fidx) = d(1).bytes; -% % -% % timedatpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},'time.dat');%int32 -% % d = dir(timedatpaths{fidx}); -% % timedatsizes(fidx) = d(1).bytes; -% % -% % supplydatpaths{fidx} = fullfile(basepath,SessionMetadata.ExtracellEphys.Files.Names{fidx},'supply.dat');%uint16 -% % d = dir(timedatpaths{fidx}); -% % supplydatsizes(fidx) = d(1).bytes; -% end -% otherdattypes(find(bad_otherdattypes)) = []; -% case 'Amplipex'%As of 4/9/2017 - never tested -% for fidx = 1:length(SessionMetadata.ExtracellEphys.Files.Names) -% datpaths{fidx} = fullfile(basepath,[SessionMetadata.ExtracellEphys.Files.Names{fidx} '.dat']); -% end -% end -% -% %% Concatenate the main data files -% if isunix -% cs = strjoin(datpaths); -% catstring = ['! cat ', cs, ' > ',newdatpath]; -% elseif ispc%As of 4/9/2017 - never tested -% if length(datpaths)>1 -% for didx = 1:length(datpaths)-1; -% datpathsplus{didx} = [datpaths{didx} '+']; -% end -% else -% datpathsplus = datpaths; -% end -% cs = strjoin(datpathsplus); -% catstring = ['! copy /b ', cs, ' ',newdatpath]; -% end -% -% eval(catstring)%execute concatention -% -% % Check that size of resultant .dat is equal to the sum of the components -% t = dir(newdatpath); -% recordingbytes = SessionMetadata.ExtracellEphys.Files.Bytes; -% if t.bytes ~= sum(recordingbytes) -% error('New .dat size not right. Exiting') -% return -% else -% disp(['Primary .dats concatenated successfully']) -% end -% -% %% if intan, also concatenate the other .dats -% if strcmp(SessionMetadata.ExtracellEphys.RecordingSystem,'Intan') -% for odidx = 1:length(otherdattypes) -% eval(['tdatpaths = ' otherdattypes{odidx} 'datpaths;']); -% eval(['tnewdatpath = new' otherdattypes{odidx} 'path;']); -% if isunix -% cs = strjoin(tdatpaths); -% catstring = ['! cat ', cs, ' > ',tnewdatpath]; -% elseif ispc%As of 4/9/2017 - never tested -% if length(tdatpaths)>1 -% for didx = 1:length(tdatpaths)-1; -% datpathsplus{didx} = [tdatpaths{didx} '+']; -% end -% else -% datpathsplus = tdatpaths; -% end -% cs = strjoin(datpathsplus); -% catstring = ['! copy /b ', cs, ' ',tnewdatpath]; -% end -% -% eval(catstring)%execute concatenation -% -% % Check that size of resultant .dat is equal to the sum of the components -% t = dir(tnewdatpath); -% eval(['recordingbytes = ' otherdattypes{odidx} 'datsizes;']) -% if t.bytes ~= sum(recordingbytes) -% error(['New ' otherdattypes{odidx} '.dat size not right. Exiting after .dats converted. Not deleting']) -% deleteoriginaldatsbool = 0; -% else -% disp([otherdattypes{odidx} ' concatenated successfully']) -% end -% end -% end -% -% -% %% Delete original dats, if that option is chosen -% if deleteoriginaldatsbool -% %for other .dats -% for odidx = 1:length(otherdattypes) -% eval(['tdatpaths = ' otherdattypes{odidx} 'datpaths;']); -% for didx = 1:length(tdatpaths) -% if isunix -% eval(['! rm ' tdatpaths{didx}]) -% elseif ispc%As of 4/9/2017 - never tested -% eval(['! del ' tdatpaths{didx}]) -% end -% end -% end -% %for main .dat -% for didx = 1:length(datpaths) -% if isunix -% eval(['! rm ' datpaths{didx}]) -% elseif ispc%As of 4/9/2017 - never tested -% eval(['! del ' datpaths{didx}]) -% end -% end -% end - diff --git a/utilities/bz_FindCatableDims.m b/utilities/bz_FindCatableDims.m index b6729a5b..640b1ed5 100755 --- a/utilities/bz_FindCatableDims.m +++ b/utilities/bz_FindCatableDims.m @@ -10,6 +10,9 @@ % %DLevenstein 2018 %% +%Remove any empty arrays +tryarrays(cellfun(@(X) isempty(X),tryarrays)) = []; + arraysizes = cellfun(@size,tryarrays,'uniformoutput',false); numdims = cellfun(@length,arraysizes); diff --git a/visualization/LogScale.m b/visualization/LogScale.m index 650f7176..09bc9a84 100755 --- a/visualization/LogScale.m +++ b/visualization/LogScale.m @@ -21,14 +21,14 @@ ticks = [range(1):0.5:range(2)]; end - if length(ticks)>=5 + if length(ticks)>5 ticks = ticks(1:2:end); end set(gca,'YTick',ticks) - set(gca,'YTickLabels',round(logbase.^ticks,1,'significant')) + set(gca,'YTickLabels',round(logbase.^ticks,2,'significant')) - if max(abs(ticks))>=3 + if max(abs(ticks))>=3 && logbase==10 tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false); tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false); tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false); @@ -51,9 +51,9 @@ set(gca,'XTick',ticks) - set(gca,'XTickLabels',round(logbase.^ticks,1,'significant')) + set(gca,'XTickLabels',round(logbase.^ticks,2,'significant')) - if max(abs(ticks))>=3 + if max(abs(ticks))>=3 && logbase==10 tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false); tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false); tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false); @@ -69,12 +69,12 @@ ticks = [range(1):0.5:range(2)]; end - if length(ticks)>=5 + if length(ticks)>5 ticks = ticks(1:2:end); end set(gca,'ZTick',ticks) - set(gca,'ZTickLabels',round(logbase.^ticks,1,'significant')) + set(gca,'ZTickLabels',round(logbase.^ticks,2,'significant')) if max(abs(ticks))>=2 tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false); diff --git a/visualization/bz_MultiLFPPlot.m b/visualization/bz_MultiLFPPlot.m index 05625a90..1c344ea7 100755 --- a/visualization/bz_MultiLFPPlot.m +++ b/visualization/bz_MultiLFPPlot.m @@ -20,6 +20,7 @@ % 'axhandle' axes handle in which to put the plot % 'scaleLFP' multiplicative factor to scale the y range of LFP % 'scalespikes' size of spike points (default:5) +% 'spikeside' 'top' (default) or 'bottom' % % %DLevenstein 2017 @@ -39,6 +40,7 @@ addParameter(p,'scaleLFP',1,@isnumeric) addParameter(p,'scalespikes',5,@isnumeric) addParameter(p,'plotcells',nan,@isnumeric) +addParameter(p,'spikeside','top') parse(p,varargin{:}) timewin = p.Results.timewin; channels = p.Results.channels; @@ -49,6 +51,7 @@ ax = p.Results.axhandle; scaleLFP = p.Results.scaleLFP; scalespikes = p.Results.scalespikes; +spikeside = p.Results.spikeside; if isempty(spikes) spikes = spikedefault; @@ -106,18 +109,29 @@ %Space based on median absolute deviation over entire recording - robust to outliers. randtimes = randsample(size(lfp.data,1),1000); -channelrange = 12.*mad(single(lfp.data(randtimes,chindex)),1); +channelrange = 10.*mad(single(lfp.data(randtimes,chindex)),1); lfpmidpoints = -cumsum(channelrange); lfp.plotdata = (bsxfun(@(X,Y) X+Y,single(lfp.data(windex,chindex)).*scaleLFP,lfpmidpoints)); -spikeplotrange = [1 -lfpmidpoints(1)]; +switch spikeside + case 'top' + spikeplotrange = [0 -lfpmidpoints(1)]; + case 'bottom' + spikeplotrange = lfpmidpoints(end)+[1.5 0.5].*lfpmidpoints(1); +end spikes.plotdata = spikes.spindices(winspikes,:); -spikes.plotdata(:,2) = (spikes.plotdata(:,2)./max(spikes.spindices(:,2))).*(diff(spikeplotrange)); - +%spikes.plotdata(:,2) = (spikes.plotdata(:,2)./max(spikes.spindices(:,2))).*(diff(spikeplotrange)); +spikes.plotdata(:,2) = bz_NormToRange(spikes.plotdata(:,2),spikeplotrange); %% Do the plot ywinrange = fliplr(lfpmidpoints([1 end])+1.*[1 -1].*max(channelrange)); if ~isnan(spikes.spindices) - ywinrange(2) = ywinrange(2)+max([spikes.plotdata(:,2);0]); + switch spikeside + case 'top' + ywinrange(2) = ywinrange(2)+max([spikes.plotdata(:,2);0]); + case 'bottom' + ywinrange(1) = min(spikes.plotdata(:,2)); + end + end plot(ax,lfp.timestamps(windex),lfp.plotdata,'k','linewidth',0.5) diff --git a/visualization/bz_NormToRange.m b/visualization/bz_NormToRange.m new file mode 100644 index 00000000..e4e7f724 --- /dev/null +++ b/visualization/bz_NormToRange.m @@ -0,0 +1,30 @@ +function [ normdata ] = bz_NormToRange(data,range,databounds) +%[ normdata ] = NormToRange(data,range,databounds) normalizes some data to +%fit in some range. +% +%INPUTS +% data the data you want to normalize +% range [min max] you would like to normalize it to. +% use 'ylim' to set to min/max of current plot (default) +% databounds [min max] of the data (optional) +% +%OUTPUTS +% normdata normalized data +% +%DLevenstein 2019 +%% +if ~exist('range','var') || strcmp(range,'ylim') + range = get(gca,'ylim'); +end + +if ~exist('databounds','var') + databounds(1) = min(data); databounds(2) = max(data); +end +dataspan = diff(databounds); +rangespan = diff(range); + +normdata = (data-databounds(1))./dataspan; +normdata = normdata.*rangespan+range(1); + +end + diff --git a/visualization/bz_ScaleBar.m b/visualization/bz_ScaleBar.m new file mode 100644 index 00000000..1fdd0758 --- /dev/null +++ b/visualization/bz_ScaleBar.m @@ -0,0 +1,15 @@ +function [ ] = bz_ScaleBar( units ) +%UNTITLED Summary of this function goes here +% Detailed explanation goes here + + % units = 'AU'; + yrange = get(gca,'ylim'); + xrange = get(gca,'xlim'); + xticks = get(gca,'xtick'); + xbar = xticks(1:2);xbarsize = diff(xbar); + set(gca,'xtick',xrange-0.5.*xbarsize) + set(gca,'xticklabels',[num2str(xbarsize),' ',units]) + plot(xrange,[1 1].*yrange(1),'w','linewidth',3) + plot(xrange(2)-[xbarsize 0],[1 1].*yrange(1),'k','linewidth',3) +end +