diff --git a/+nigeLab/+defaults/+SD/ART_HardThresh.m b/+nigeLab/+defaults/+SD/ART_HardThresh.m new file mode 100644 index 00000000..abb9ef6a --- /dev/null +++ b/+nigeLab/+defaults/+SD/ART_HardThresh.m @@ -0,0 +1,9 @@ +function pars = ART_HardThresh() +%% function defining defualt parameters for HARD THRESHOLD spike detection algorithm + +pars.Thresh = 70; % [uV] Fixed voltage threshold for detection; +pars.Samples = 4; % [ms] Window to ignore around artifact (suggest: 4 ms MIN for stim rebound) + +pars.Polarity = 1; + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/FEAT_ica.m b/+nigeLab/+defaults/+SD/FEAT_ica.m new file mode 100644 index 00000000..4c0b417d --- /dev/null +++ b/+nigeLab/+defaults/+SD/FEAT_ica.m @@ -0,0 +1,5 @@ +function pars = FEAT_ica() +%% function defining defualt parameters for ICA based feature extraction +pars = struct(); + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/FEAT_pca.m b/+nigeLab/+defaults/+SD/FEAT_pca.m new file mode 100644 index 00000000..83c650aa --- /dev/null +++ b/+nigeLab/+defaults/+SD/FEAT_pca.m @@ -0,0 +1,10 @@ +function pars = FEAT_pca() +%% function defining defualt parameters for PCA based feature extraction +pars.ExplVar = .95; % Explained Variance to retain durin PCA + % decomposition. + % Takes precedence over NOut. Set to Inf to + % use NOut. + +pars.NOut = 12; % Number of feature inputs for clustering. + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/FEAT_wavelet.m b/+nigeLab/+defaults/+SD/FEAT_wavelet.m new file mode 100644 index 00000000..beb7398c --- /dev/null +++ b/+nigeLab/+defaults/+SD/FEAT_wavelet.m @@ -0,0 +1,10 @@ +function pars = FEAT_wavelet() +%% function defining defualt parameters for wavelet feature extraction + +pars.WaveName = 'bior1.3'; % 'haar' 'bior1.3' 'db4' 'sym8' all examples +[pars.LoD,pars.HiD] = wfilters(pars.WaveName); % get wavelet decomposition parameters +pars.NOut = 12; % Number of feature inputs for clustering +pars.NScales = 3; % Number of scales for wavelet decomposition +% + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/SD_AdaptThresh.m b/+nigeLab/+defaults/+SD/SD_AdaptThresh.m new file mode 100644 index 00000000..42c6572e --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_AdaptThresh.m @@ -0,0 +1,10 @@ +function pars = SD_AdaptThresh() +%% function defining defualt parameters for SNEO spike detection algorithm + +pars.FilterLength = 60; % [ms] Length of the adaptive filter window +pars.Polarity = -1; % polarity of the detection. If positive looks for positive crossings. Negative otherwise. +pars.MinThresh = 15; % [uV] Fixed minimum voltage threshold for detection; +pars.MultCoeff = 4.5; % moltiplicative factor for the adaptive threshold (signal absolute median); +pars.RefrTime = 0.5; % [ms] Refractory time. +pars.PeakDur = 1; % [ms] Peak duration or pulse lifetime period +end diff --git a/+nigeLab/+defaults/+SD/SD_HardThresh.m b/+nigeLab/+defaults/+SD/SD_HardThresh.m new file mode 100644 index 00000000..c6fe1e3a --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_HardThresh.m @@ -0,0 +1,9 @@ +function pars = SD_HardThresh() +%% function defining defualt parameters for HARD THRESHOLD spike detection algorithm + +pars.Polarity = -1; % polarity of the detection. If positive looks for positive crossings. Negative otherwise. +pars.Thresh = 50; % [uV] Fixed voltage threshold for detection; +pars.RefrTime = 0.5; % [ms] Refractory time. +pars.PeakDur = 1; % [ms] Peak duration or pulse lifetime period +pars.NSaround = 7; +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/SD_PTSD.m b/+nigeLab/+defaults/+SD/SD_PTSD.m new file mode 100644 index 00000000..589e3be1 --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_PTSD.m @@ -0,0 +1,9 @@ +function pars = SD_PTSD() +%% function defining defualt parameters for SNEO spike detection algorithm + +% pars.MultCoeff = 4.5; % Multiplication coefficient for noise +pars.Thresh = 50; +pars.RefrTime = 0.5; % [ms] Refractory time. +pars.PeakDur = 2; % [ms] Peak duration or pulse lifetime period +pars.AlignFlag = 0; +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/SD_SNEO.m b/+nigeLab/+defaults/+SD/SD_SNEO.m new file mode 100644 index 00000000..3d0546db --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_SNEO.m @@ -0,0 +1,10 @@ +function pars = SD_SNEO() +%% function defining defualt parameters for SNEO spike detection algorithm + +pars.MultCoeff = 4.5; % Multiplication coefficient for noise +pars.SmoothN = 5; % Number of samples to use for smoothed nonlinear energy operator window +pars.NSaround = 7; % Number of samples around the peak to "look" for negative peak +pars.RefrTime = 0.5; % [ms] Refractory time. +pars.PeakDur = 1; % [ms] Peak duration or pulse lifetime period +pars.Polarity = -1; +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/SD_SWTTEO.m b/+nigeLab/+defaults/+SD/SD_SWTTEO.m new file mode 100644 index 00000000..8573651a --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_SWTTEO.m @@ -0,0 +1,16 @@ +function pars = SD_SWTTEO() +%% function defining defualt parameters for SWTTEO spike detection algorithm +pars.wavLevel = 4; % Wavelet decomposition level +pars.waveName = 'sym5'; % wavelet type + +pars.winType = @hamming; % function handle for the smoothing window type; This is fed to window function +pars.smoothN = 40; % Number of samples for the smoothing operator. Set to 0 to turn off smoothing +pars.winPars = {'symmetric'}; % Optional parameters for the smoothing window + +pars.RefrTime = 1; % [ms] refractory time +pars.MultCoeff = 3.5; % Moltiplication coefficient for SWTTEO thresholding +pars.Polarity = -1; +pars.PeakDur = 2; % [ms] Peak duration or pulse lifetime period + + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/SD_TIFCO.m b/+nigeLab/+defaults/+SD/SD_TIFCO.m new file mode 100644 index 00000000..00d65ac7 --- /dev/null +++ b/+nigeLab/+defaults/+SD/SD_TIFCO.m @@ -0,0 +1,15 @@ +function pars = SD_TIFCO() +%% function defining defualt parameters for TIFCO spike detection algorithm +pars.fMin = 300; % [Hz] lower bound for the gabor based time frequency decomposition +pars.fMax = 3500; % [Hz] higher bound for the gabor based time frequency decomposition + +pars.winType = @hamming; % function handle for the smoothing window type; This is fed to window function +pars.winL = 1; % [s] Length for the time-frequency window decomposition +pars.winPars = {'symmetric'}; % Optional parameters for the smoothing window + +pars.RefrTime = 1; % [ms] refractory time +pars.MultCoeff = 4.5; % Multiplication coefficient for noise +pars.Polarity = -1; % Detection polarity (look for positive or negative peaks) +pars.PeakDur = 1; % [ms] Peak duration or pulse lifetime period + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/WIP_SD_MTEO.m b/+nigeLab/+defaults/+SD/WIP_SD_MTEO.m new file mode 100644 index 00000000..21ccffc6 --- /dev/null +++ b/+nigeLab/+defaults/+SD/WIP_SD_MTEO.m @@ -0,0 +1,8 @@ +function pars = MTEO() +%% function defining defualt parameters for MTEO spike detection algorithm + +pars.multcoeff = 4.5; % Multiplication coefficient for noise +pars.smoothN = 5; % Number of samples to use for smoothed nonlinear energy operator window +pars.nsaround = 7; % Number of samples around the peak to "look" for negative peak + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/+SD/WIP_SD_WTEO.m b/+nigeLab/+defaults/+SD/WIP_SD_WTEO.m new file mode 100644 index 00000000..c6377bda --- /dev/null +++ b/+nigeLab/+defaults/+SD/WIP_SD_WTEO.m @@ -0,0 +1,8 @@ +function pars = SD_WTEO() +%% function defining defualt parameters for WTEO spike detection algorithm + +pars.multcoeff = 4.5; % Multiplication coefficient for noise +pars.smoothN = 5; % Number of samples to use for smoothed nonlinear energy operator window +pars.nsaround = 7; % Number of samples around the peak to "look" for negative peak + +end \ No newline at end of file diff --git a/+nigeLab/+defaults/AutoClustering.m b/+nigeLab/+defaults/AutoClustering.m index 311c1c60..1ba8005d 100644 --- a/+nigeLab/+defaults/AutoClustering.m +++ b/+nigeLab/+defaults/AutoClustering.m @@ -5,8 +5,15 @@ %% PARAMS YOU MIGHT CHANGE pars = struct; -pars.MethodName = 'KMEANS'; % Can be: 'KMEANS' or 'SPC' -pars.NMaxClus = 9; % Maximum # of clusters +pars.MethodName = 'KMEANS'; % Can be: 'KMEANS' or 'SPC' +pars.NMaxClus = 9; % Maximum # of clusters +pars.ID.Clusters = 'kmeans'; % Attached to the cluster folder + + + + + + %% UNLIKELY TO CHANGE % Parameters for each type stored as individual files in ~/+Autoclusters diff --git a/+nigeLab/+defaults/SD.m b/+nigeLab/+defaults/SD.m index 71fa7d41..85cfb62f 100644 --- a/+nigeLab/+defaults/SD.m +++ b/+nigeLab/+defaults/SD.m @@ -1,4 +1,4 @@ -function pars = SD(varargin) +function varargout = SD(varargin) %% defaults.SD Initialize parameters for spike detection % % pars = nigeLab.defaults.SD('NAME',value,...); @@ -11,129 +11,109 @@ % By: MAECI 2018 collaboration (Max Murphy & Federico Barban) % 01/09/2019 :: 4.0.0 -> 4.0.1 = Fix FEAT_NAMES generation -%% DEFAULTS pars = struct; -% General settings -pars.VERSION = 'v4.0.1'; % Version, to be passed with parameters -pars.LIBDIR = 'C:\MyRepos\_SD\APP_Code'; % Location of sub-functions + +% %% OLD SPIKE DETECTION PARAMETERS +% % Parameters +% pars.ARTIFACT_THRESH = 450; % Threshold for artifact -% Folder tags -pars.USE_CAR = true; % By def. use common spatial reference +% pars.PRE_STIM_BLANKING = 0.5; % Window to blank before specifieid stim times (ms) +% pars.POST_STIM_BLANKING = 1.5; % Window to blank after specified stim times (ms) +% pars.ARTIFACT_SPACE = 4; % Window to ignore around artifact (suggest: 4 ms MIN for stim rebound) +% pars.MULTCOEFF = 4.5; % Multiplication coefficient for noise +% pars.PKDURATION = 1.0; % Pulse lifetime period (suggest: 2 ms MAX) +% pars.REFRTIME = 0.5; % Refractory period (suggest: 2 ms MAX). +% pars.PKDETECT = 'sneo';% 'both' or 'pos' or 'neg' or 'adapt' or 'sneo' for peak type +% pars.ALIGNFLAG = 1; % Alignment flag for detection +% % [0 -> highest / 1 -> most negative] +% pars.P2PAMP = 60; % Minimum peak-to-peak amplitude -% File tags -pars.DELETE_OLD_PATH = false; % Remove old files -pars.USE_EXISTING_SPIKES = false; % Use existing spikes on directory -pars.DO_AUTO_CLUSTERING = true; % If false, skips "clustering" portion - % (v2017a in order to use Isilon cluster) +% pars.ART_DIST = 1/35; % Max. time between stimuli (sec) +% pars.NWIN = 120; % Number of windows for automatic thresholding +% pars.WINDUR = 200*1e-3; % Minimum window length (msec) +% pars.INIT_THRESH = 50; % Pre-adaptive spike threshold (micro-volts) +% pars.ART_RATE = 0.0035; % Empirically determined rate for artifacts based on artifact rejection +% pars.M = (-7/3); % See ART_RATE +% pars.B = 1.05; % See ART_RATE +% +% +% % Parameters +% pars.MIN_SPK = 100; % Minimum spikes before sorting +% pars.TEMPSD = 3.5; % Cluster template max radius for template matching +% pars.TSCALE = 3.5; % Scaling for timestamps of spikes as a feature -% % Isilon cluster settings -pars.USE_CLUSTER = false; % Must already have clustering done -%% SPIKE DETECTION PARAMETERS -% Parameters -pars.ARTIFACT_THRESH = 450; % Threshold for artifact +%% User defined parameters for spike detection + pars.STIM_TS = []; % Pre-specified stim times pars.ARTIFACT = []; % Pre-specified artifact times -pars.PRE_STIM_BLANKING = 0.5; % Window to blank before specifieid stim times (ms) -pars.POST_STIM_BLANKING = 1.5; % Window to blank after specified stim times (ms) -pars.ARTIFACT_SPACE = 4; % Window to ignore around artifact (suggest: 4 ms MIN for stim rebound) -pars.MULTCOEFF = 4.5; % Multiplication coefficient for noise -pars.PKDURATION = 1.0; % Pulse lifetime period (suggest: 2 ms MAX) -pars.REFRTIME = 0.5; % Refractory period (suggest: 2 ms MAX). -pars.PKDETECT = 'sneo';% 'both' or 'pos' or 'neg' or 'adapt' or 'sneo' for peak type -pars.ADPT_N = 60; % Number of ms to use for adaptive filter -pars.SNEO_N = 5; % Number of samples to use for smoothed nonlinear energy operator window -pars.NS_AROUND = 7; % Number of samples around the peak to "look" for negative peak -pars.ADPT_MIN = 15; % Minimum for adaptive threshold (fixed) -pars.ALIGNFLAG = 1; % Alignment flag for detection -% [0 -> highest / 1 -> most negative] -pars.P2PAMP = 60; % Minimum peak-to-peak amplitude -pars.W_PRE = 0.4; % Pre-spike window (ms) -pars.W_POST = 0.8; % Post-spike window (ms) -pars.ART_DIST = 1/35; % Max. time between stimuli (sec) -pars.NWIN = 120; % Number of windows for automatic thresholding -pars.WINDUR = 200*1e-3; % Minimum window length (msec) -pars.INIT_THRESH = 50; % Pre-adaptive spike threshold (micro-volts) -pars.PRESCALED = true; % Whether data has been pre-scaled to micro-volts. -pars.FIXED_THRESH = 50; % If alignment is 'neg' or 'pos' this can be set to fix the detection threshold level -pars.ART_RATE = 0.0035; % Empirically determined rate for artifacts based on artifact rejection -pars.M = (-7/3); % See ART_RATE -pars.B = 1.05; % See ART_RATE - -% Spike features and sorting settings (SPC pars in SPIKECLUSTER_SPC) -pars.SC_VER = 'SPC'; % Version of spike clustering - -% Parameters -pars.N_INTERP_SAMPLES = 250; % Number of interpolated samples for spikes -pars.MIN_SPK = 100; % Minimum spikes before sorting -pars.TEMPSD = 3.5; % Cluster template max radius for template matching -pars.TSCALE = 3.5; % Scaling for timestamps of spikes as a feature -pars.USE_TS_FEATURE = false; % Add timestamp as an additional feature for SPC? -pars.FEAT = 'wav'; % 'wav' or 'pca' or 'ica' for spike features -pars.WAVELET = 'bior1.3';% 'haar' 'bior1.3' 'db4' 'sym8' all examples -[pars.LoD,pars.HiD] = wfilters(pars.WAVELET); % get wavelet decomposition parameters -pars.NINPUT = 12; % Number of feature inputs for clustering -pars.NSCALES = 3; % Number of scales for wavelet decomposition - -%% PARSE VARARGIN -if numel(varargin)==1 - varargin = varargin{1}; - if numel(varargin) ==1 - varargin = varargin{1}; - end -end -for iV = 1:2:length(varargin) - pars.(upper(varargin{iV}))=varargin{iV+1}; -end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% Spike Detection +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +pars.WPre = 0.4; % Pre-spike window (ms) +pars.WPost = 0.8; % Post-spike window (ms) -%% PARSE OTHER PARAMETERS BASED ON SELECTED PARAMETERS -switch pars.PKDETECT - case 'neg' - pars.SD_VER = [pars.FEAT '-neg' num2str(pars.FIXED_THRESH)]; - case 'pos' - pars.SD_VER = [pars.FEAT '-pos' num2str(pars.FIXED_THRESH)]; - case 'adapt' - pars.SD_VER = [pars.FEAT '-adapt']; - case 'both' - pars.SD_VER = [pars.FEAT '-PT']; - case 'sneo' - pars.SD_VER = [pars.FEAT '-sneo']; - otherwise - pars.SD_VER = [pars.FEAT '-new']; -end +pars.SDMethodName = 'SWTTEO'; +pars.ID.Spikes = 'SWTTEO'; % implemented to date (2020/06/16): + % SNEO, SWTTEO, WTEO, TIFCO, SWT, + % PTSD, fixed and variable Threshold. + % See documentation for references. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% Artefact Rejection +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +pars.ArtefactRejMethodName = 'HardThresh'; +pars.ID.Artifact = 'HardThresh'; % implemented to date (2020/06/16): + % HardThresh (hard threshold) -if pars.USE_CAR - pars.ID.Artifact = [pars.SD_VER '_CAR']; - pars.ID.Spikes = [pars.SD_VER '_CAR']; - pars.ID.SpikeFeatures = [pars.SD_VER '_CAR']; - pars.ID.Clusters = [pars.SD_VER '_' pars.SC_VER '_CAR']; - pars.ID.Sorted = [pars.SD_VER '_' pars.SC_VER '_CAR']; -else - pars.ID.Artifact = pars.SD_VER; - pars.ID.Spikes = pars.SD_VER; - pars.ID.SpikeFeatures = pars.SD_VER; - pars.ID.Clusters = [pars.SD_VER '_' pars.SC_VER]; - pars.ID.Sorted = [pars.SD_VER '_' pars.SC_VER]; + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% Feature Extraction +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +pars.ID.SpikeFeatures = 'wavelet'; +pars.InterpSamples = 250; % Number of interpolated + % samples per spike. Can help + % with feature extraction. +pars.FeatureExtractionMethodName = 'wavelet'; % implemented to date (2020/06/16): + % wavelet, pca + + + + + + + +%% UNLIKELY TO CHANGE +% Parameters for each type stored as individual files in ~/+SD +SDPath = fullfile(nigeLab.utils.getNigelPath,... + '+nigeLab','+defaults','+SD'); +% Load all the method-specific parameters: +SDConfigFiles = dir(fullfile(SDPath,'*.m')); +for ff = SDConfigFiles(:)' + parName = ff.name(1:end-2); % dropping .m + pars.(parName) = eval(sprintf('nigeLab.defaults.SD.%s',... + parName)); end -switch pars.FEAT - case 'raw' - pars.FEAT_NAMES = cell(1,2*pars.NINPUT + 1); - for ii = 1:pars.NINPUT - pars.FEAT_NAMES{ii} = sprintf('wav-%02d',ii); - end - for ik = 1:pars.NINPUT - pars.FEAT_NAMES{ii+ik} = sprintf('raw-%02d',ik); - end - pars.FEAT_NAMES{ii+ik+1} = 'raw-mean'; - otherwise - pars.FEAT_NAMES = cell(1,pars.NINPUT); - for ii = 1:pars.NINPUT - pars.FEAT_NAMES{ii} = sprintf('%s-%02d',pars.FEAT,ii); +%% Parse output +if nargin < 1 + varargout = {pars}; +else + varargout = cell(1,nargin); + f = fieldnames(pars); + for i = 1:nargin + idx = ismember(lower(f),lower(varargin{i})); + if sum(idx) == 1 + varargout{i} = pars.(f{idx}); end + end end -end \ No newline at end of file +end diff --git a/+nigeLab/+defaults/Sort.m b/+nigeLab/+defaults/Sort.m index 92103937..ac9fa79f 100644 --- a/+nigeLab/+defaults/Sort.m +++ b/+nigeLab/+defaults/Sort.m @@ -7,9 +7,12 @@ pars = struct; % carries all parameter variables % Defaults for UI usability -pars.InFileFilt = {'*_Block.mat';'*_Animal.mat';'*_Tank.mat'}; -pars.InFilePrompt = 'Select BLOCK(S), ANIMAL(S), or TANK'; -pars.InFileDefDir = 'P:\Rat'; +pars.InFileFilt = {'*_Block.mat';'*_Animal.mat';'*_Tank.mat'}; +pars.InFilePrompt = 'Select BLOCK(S), ANIMAL(S), or TANK'; +pars.InFileDefDir = 'P:\Rat'; + +pars.ID.Sorted = 'Manual'; % Label attached to the Sorted folder + pars.ForceNext = false; % Automatically jump to next channel on "confirm" pars.Debug = false; % Set to TRUE to move handles to base workspace diff --git a/+nigeLab/+defaults/nigelColors.m b/+nigeLab/+defaults/nigelColors.m index e2dd0b10..83644fa3 100644 --- a/+nigeLab/+defaults/nigelColors.m +++ b/+nigeLab/+defaults/nigelColors.m @@ -88,6 +88,7 @@ warning('%s is not a recognized nigelColors option.',input); end case 2 + maxColorSteps = 25; switch source case 'cubehelix' if not(isnumeric(input)) @@ -95,8 +96,13 @@ Col = []; return; else - Clrs = cubehelix(25,2.8,1.9,1.45,0.6,[0.05 0.95],[0 .7]); - Col = Clrs(input,:); + idx = 1:maxColorSteps; + while numel(idx) < numel(input) + idx = [idx fliplr(idx(1:end-1))]; + end + idx = idx(1:numel(input)); + Clrs = cubehelix(maxColorSteps,1.45,1.57,2.20,0.78,[0.25 0.57],[0.28 .82]); + Col = Clrs(idx,:); end end end diff --git a/+nigeLab/+evt/dataScrolled.m b/+nigeLab/+evt/dataScrolled.m new file mode 100644 index 00000000..2105f0be --- /dev/null +++ b/+nigeLab/+evt/dataScrolled.m @@ -0,0 +1,11 @@ +classdef (ConstructOnLoad) dataScrolled < event.EventData + properties + ROI + end + + methods + function data = dataScrolled(newRoi) + data.ROI = newRoi; + end + end +end \ No newline at end of file diff --git a/+nigeLab/+evt/saveData.m b/+nigeLab/+evt/saveData.m new file mode 100644 index 00000000..e0961bc7 --- /dev/null +++ b/+nigeLab/+evt/saveData.m @@ -0,0 +1,32 @@ +classdef (ConstructOnLoad) saveData < event.EventData +% SAVEDATA Event for to signal saving of data has begun +% +% Tipycally issued during spike sorting (Sort Interface) +% +% ASSIGNCLUS Properties: +% +% ch - Subscripts of channels to save +% +% ASSIGNCLUS Methods: +% assignClus - Event data constructor. +% evt = nigeLab.evt.assignClus(ch); + + properties + ch % Subscripts of channels to save + end + + methods + function evt = saveData(ch) + % SAVEDATA Cluster assignment event class constructor. + % + % evt = nigeLab.evt.assignClus(subs,class,otherClassToUpdate); + % + % Issued when a subset of spikes from one cluster is assigned to + % a new cluster, typically in the `nigeLab.Sort` spike sorting + % interface + + evt.ch = ch; + end + end + +end \ No newline at end of file diff --git a/+nigeLab/+libs/@ChannelUI/ChannelUI.m b/+nigeLab/+libs/@ChannelUI/ChannelUI.m index 4abd3a9e..6a698667 100644 --- a/+nigeLab/+libs/@ChannelUI/ChannelUI.m +++ b/+nigeLab/+libs/@ChannelUI/ChannelUI.m @@ -108,7 +108,8 @@ function Open(obj) 'MenuBar', 'none', ... 'Units','Normalized',... 'Color',obj.FIG_COL, ... - 'Position',obj.FIG_POS); + 'Position',obj.FIG_POS,... + 'CloseRequestFcn',[]); % it has to be destructed from the parent obj.Menu = uicontrol(obj.Figure,... 'Style','popupmenu',... diff --git a/+nigeLab/+libs/@DataScrollerAxis/DataScrollerAxis.m b/+nigeLab/+libs/@DataScrollerAxis/DataScrollerAxis.m new file mode 100644 index 00000000..56c93be4 --- /dev/null +++ b/+nigeLab/+libs/@DataScrollerAxis/DataScrollerAxis.m @@ -0,0 +1,197 @@ +classdef DataScrollerAxis < handle + + properties + Field + UI + channel + ThisBlock + Channels + MainAxPixelSize + sigLenght + fs + ROIpos = zeros(1,4) + LinePlotExplorer + ReducedPlot + + Listeners + end + + events + roiChanged + end + + methods + + function obj = DataScrollerAxis(nigelObj,Field) + % DATASCROLLERAXIS Axes that allows selecting intresting + % portions of data from a channel. + % this makes use of nigeLab.utils.LinePlotReducer to reduce the + % number of points plotted to the bare minimum to represent all + % the feature of the signal + % + % ax = nigeLab.libs.TimeScrollerAxes(nigelObj,Field); + % nigelObj can be any nigelObj but it will always sample down + % a block to plot the actual data + % Field the actual field to plot. + + if nargin < 2 + Field = 'Filt'; + end + + switch nigelObj.Type + case 'Tank' + anIdx = randi(numel(nigelObj.Children)); + obj = nigeLab.libs.DataScrollerAxis(nigelObj{anIdx},Field); + return + case 'Animal' + blIdx = randi(numel(nigelObj.Children)); + obj = nigeLab.libs.DataScrollerAxis(nigelObj{blIdx},Field); + return + case 'Block' + obj.ThisBlock = nigelObj; + end + + + obj.buildGui(); + changeDataType(obj,Field); + obj.buildROI(); + obj.addListeners(); + end + + function addListeners(obj) + %% Add all needed listeners here + obj.Listeners = [obj.Listeners, ... + addlistener(obj.UI.ChannelSelector,'NewChannel',... + @obj.setChannel)]; + end + + function delete(obj) + delete(obj.UI.ChannelSelector); + delete(obj.UI.Fig); + end + + function buildGui(obj) + obj.UI.Fig = figure('Units','normalized',... + 'Position',[.1 .1 .8 .25],... + 'DeleteFcn',@(~,~)obj.delete,... + 'MenuBar','figure',... + 'NumberTitle','off',... + 'WindowButtonUpFcn',@(~,~)obj.RoiChanged); + + obj.UI.MainAx = axes(obj.UI.Fig,'Units','normalized','Position',[.02 .1 .94 .85]); + obj.MainAxPixelSize = getpixelposition(obj.UI.MainAx); + + obj.Channels.Name = {obj.ThisBlock.Channels.name}; + obj.Channels.Selected = 1; + obj.UI.Parent = obj; + obj.UI.channel = 1; + + obj.UI.ChannelSelector = nigeLab.libs.ChannelUI(obj.UI); + + end + + function buildROI(obj,pos) + if nargin < 2 + yl = ylim(obj.UI.MainAx); + xl = xlim(obj.UI.MainAx); + xmax = min(obj.sigLenght,60*obj.fs)./obj.fs; + obj.ROIpos = [xl(1) yl(1) xmax diff(yl)]; + else + yl = ylim(obj.UI.MainAx); + obj.ROIpos = pos; + end + % builds ROI overlay on top of the datascroller + + obj.UI.ROI = imrect(obj.UI.MainAx,obj.ROIpos,... + 'PositionConstraintFcn',@obj.roiResizeFcn); %#ok + ylim(obj.UI.MainAx,yl); + + end + + function newPos = roiResizeFcn(obj,pos) + newPos = pos; + yl = ylim(obj.UI.MainAx); + xl = xlim(obj.UI.MainAx); + + maxWidth = xl(2) - pos(1); + + % fix ROI height + newPos([2 4]) = [yl(1) diff(yl)]; + + + % fix ROI width + newPos(3) = min(maxWidth,newPos(3)); + + % fix ROI starting x + newPos(1) = max(xl(1),newPos(1)); % minimum + a = newPos(1) - (xl(2)-newPos(3)) >= 0; + newPos = a.*obj.ROIpos + (1-a).*newPos; +% newPos(1) = min(xl(2)-newPos(3),newPos(1)); % maximum + + obj.ROIpos = newPos; + + end + + function plotData(obj) + obj.fs = obj.ThisBlock.SampleRate; + obj.sigLenght = obj.ThisBlock.Samples; + + data = obj.ThisBlock.Channels(obj.Channels.Selected).(obj.Field)(:); + if obj.ThisBlock.getStatus('Time') && ~isempty(obj.ThisBlock.Time) + tt = obj.ThisBlock.Time(:); + else + tt = linspace(0,obj.sigLenght./obj.fs,obj.sigLenght); + end + + if isempty(tt) || length(tt)~=obj.sigLenght + %failsafe + tt = linspace(0,obj.sigLenght./obj.fs,obj.sigLenght); + end + [x_reduced, y_reduced] = nigeLab.utils.reduce_to_width(tt, data, obj.MainAxPixelSize(4) ,[tt(1) tt(end)]); + cla(obj.UI.MainAx); + L = line(obj.UI.MainAx,x_reduced,y_reduced); + obj.ReducedPlot = nigeLab.utils.LinePlotReducer(L,tt,data); + xlim(obj.UI.MainAx,[tt(1) tt(end)]); +% obj.LinePlotExplorer = nigeLab.utils.LinePlotExplorer(obj.UI.Fig); + + end + + function RoiChanged(obj) + evt = nigeLab.evt.dataScrolled(obj.ROIpos); + notify(obj,'roiChanged',evt); + end + + function setChannel(obj,src,~) + obj.Channels.Selected = src.Channel; + plotData(obj); + obj.buildROI(obj.ROIpos); + end + + function flag = set(obj,field,value) + % compatibility function foir the channel selector + + flag = true; + + switch field + case 'channel' + obj.channel = value; + flag = true; + end + end + + function flag = changeDataType(obj,type) + flag = false; + stepDone = obj.ThisBlock.getStatus; + if ~any(strcmp(stepDone,type)) + error('%s field hasn''t been computed yet.',obj.Field); + end + + obj.Field = type; + obj.plotData(); + obj.buildROI(obj.ROIpos); + flag = true; + end + + end + +end \ No newline at end of file diff --git a/+nigeLab/+libs/@ExploreData/ExploreData.m b/+nigeLab/+libs/@ExploreData/ExploreData.m new file mode 100644 index 00000000..576ee32d --- /dev/null +++ b/+nigeLab/+libs/@ExploreData/ExploreData.m @@ -0,0 +1,154 @@ +classdef ExploreData < handle + + properties + ThisBlock + UI + ROI + end + + methods + + function obj = ExploreData(nigelObj) + + switch nigelObj.Type + case 'Block' + obj.ThisBlock = nigelObj; + case 'Animal' + blIdx = randi(numel(nigelObj.Children)); + nigeLab.libs.ExploreData(nigelObj{blIdx}); + case 'Tank' + anIdx = randi(numel(nigelObj.Children)); + nigeLab.libs.ExploreData(nigelObj{anIdx}); + end + + obj.BuildGui(); + obj.addListeners(); + + end + + function addListeners(obj) + addlistener(obj.UI.DataScroller,'roiChanged',@obj.plotSnippets); + end + + function plotSnippets(obj,src,evt) + fs = obj.ThisBlock.SampleRate; + obj.ROI = cumsum(floor(evt.ROI([1 3])*fs))+1; + plotData(obj); + end + + function BuildGui(obj) + + % Parse plottable fields + Fields = obj.ThisBlock.getStatus; + Fields_ = obj.ThisBlock.Fields; + FieldsType_ = obj.ThisBlock.FieldType; + Fields = Fields_(strcmp(FieldsType_,'Channels') & ismember(Fields_,Fields) ); + + % Build gui + obj.UI.DataSelector = obj.buildDataTypeSelector(Fields); + obj.UI.DataScroller = nigeLab.libs.DataScrollerAxis(obj.ThisBlock,'Raw'); + fs = obj.ThisBlock.SampleRate; + + obj.ROI = cumsum(floor(obj.UI.DataScroller.ROIpos([1 3])*fs))+1; % DataScroller.ROIpos is [xpos ypos width height] + + + obj.UI.Fig = figure('Name','Multi-Channel Raw Snippets', ... + 'Units','Normalized', ... + 'Position',[0.05*rand+0.1,0.05*rand+0.1,0.8,0.8],... + 'Color','w','NumberTitle','off',... + 'CloseRequestFcn',@(~,~)obj.delete); + + obj.UI.MainAx = axes(obj.UI.Fig ,'NextPlot','add'); + + obj.plotData; + + end + + function str_box = buildDataTypeSelector(obj,Fields) + + % Create handle to store data and build graphics + obj.UI.DataSelectorFig = figure('Name','Data type selector', ... + 'Units', 'Normalized', ... + 'Position',[0.3 0.5 0.3 0.3],... + 'MenuBar','none',... + 'ToolBar','none',... + 'NumberTitle','off'); + fig = obj.UI.DataSelectorFig; + + p = nigeLab.libs.nigelPanel(fig,... + 'String','Data type selector',... + 'Tag','uidropdownbox',... + 'Units','normalized',... + 'Position',[0 0 1 1],... + 'Scrollable','off',... + 'PanelColor',nigeLab.defaults.nigelColors('surface'),... + 'TitleBarColor',nigeLab.defaults.nigelColors('tertiary'),... + 'TitleColor',nigeLab.defaults.nigelColors('ontertiary')); + + prompt_text = uicontrol('Style','text',... + 'Units','Normalized', ... + 'Position',[0.2 0.775 0.6 0.20],... + 'FontSize', 22, ... + 'FontWeight','bold',... + 'FontName','Droid Sans',... + 'BackgroundColor',nigeLab.defaults.nigelColors('surface'),... + 'ForegroundColor',nigeLab.defaults.nigelColors('tertiary'),... + 'String','Select a field.'); + p.nestObj(prompt_text); + + str_box = uicontrol('Style','popupmenu',... + 'Units', 'Normalized', ... + 'Position',[0.05 0.4 0.9 0.30],... + 'FontSize', 16, ... + 'FontName','Droid Sans',... + 'String',Fields,... + 'Callback',@(~,~)obj.changeDataTyoe); + + p.nestObj(str_box); + end + + function delete(obj) + + delete(obj.UI.DataSelector); + delete(obj.UI.DataSelectorFig); + delete(obj.UI.DataScroller); + delete(obj.UI.MainAx); + delete(obj.UI.Fig); + end + + function changeDataTyoe(obj) + idx = obj.UI.DataSelector.Value; + Field = obj.UI.DataSelector.String{idx}; + + switch Field + case 'Spikes' + if all(obj.ThisBlock.getStatus({'CAR'})) + Field = 'CAR'; + elseif all(obj.ThisBlock.getStatus({'Filt'})) + Field = 'Filt'; + else + error('CAR or Filt are needed to overlay spikes!'); + end + otherwise + ... do nothing for now + end + obj.UI.DataScroller.changeDataType(Field); + + obj.plotData(); + end + + + function plotData(obj) + + idx = obj.UI.DataSelector.Value; + Field = obj.UI.DataSelector.String{idx}; + cla(obj.UI.MainAx,'reset'); + obj.UI.Fig.Name = sprintf('Multi-Channel %s Snippets',Field); + obj.ThisBlock.plotWaves(obj.UI.MainAx,... + Field,obj.ROI(1):obj.ROI(2)); + + end + + end + +end \ No newline at end of file diff --git a/+nigeLab/+libs/@FeaturesUI/FeaturesUI.m b/+nigeLab/+libs/@FeaturesUI/FeaturesUI.m index cd3a139e..6fd86b15 100644 --- a/+nigeLab/+libs/@FeaturesUI/FeaturesUI.m +++ b/+nigeLab/+libs/@FeaturesUI/FeaturesUI.m @@ -41,6 +41,8 @@ HighDimsParams Features2D + PropLinker + ClusterQuality Silhouette QualityIndx @@ -83,6 +85,9 @@ [1,0,1];[0.930000000000000,0.690000000000000,0.130000000000000];... [0.300000000000000,0.950000000000000,0.950000000000000];... [0,0.450000000000000,0.750000000000000]}; + + CurrKeyPress + end % PROTECTED @@ -278,7 +283,10 @@ function Init(obj) 'Position',[0.050,0.075,0.4,0.450],... 'Color',nigeLab.defaults.nigelColors('background'),... 'SizeChangedFcn',@obj.ChangeSize,... - 'CloseRequestFcn',@(~,~)obj.ExitFeatures); + 'WindowScrollWheelFcn',@obj.WindowMouseWheel,... + 'CloseRequestFcn',@(~,~)obj.ExitFeatures,... + 'WindowKeyPressFcn',@obj.WindowKeyPress,... + 'WindowKeyReleaseFcn',@obj.clearKurrKey); featureComboSelectorPanel = uipanel(obj.Figure, ... 'Units', 'Normalized', ... @@ -324,6 +332,8 @@ function Init(obj) 'Units','Normalized',... 'FontSmoothing','off',... 'nextplot','add',... + 'XLim',[-obj.sdMax obj.sdMax],... + 'YLim',[-obj.sdMax obj.sdMax],... 'Position',[0.6 0.125 0.4 0.8],... 'View',[30 30]); @@ -389,9 +399,13 @@ function Init(obj) 'Units','Normalized',... 'nextplot','add',... 'UserData',1,... + 'XLim',[-obj.sdMax obj.sdMax],... + 'YLim',[-obj.sdMax obj.sdMax],... 'Position',[0.1 0.125 0.4 0.8],... 'ButtonDownFcn',@obj.ButtonDownFcnSelect2D); + obj.PropLinker = linkprop([obj.Features2D obj.Features3D],{'XLim','YLim'}); + rgb2Hex = ( @(rgbColour) reshape( dec2hex( rgbColour, 2 )',1, 6)); post = ''; @@ -445,7 +459,7 @@ function Init(obj) end % Initializes mesh map for the 2D featurespace - function InitNewMeshMap(obj,nSD,nMeshEdges) + function InitNewMeshMap(obj,nMeshEdges,xSpan,ySpan) %INITNEWMESHMAP Initialize mesh for 2D featurespace % % obj.InitNewMeshMap(nSD,nMeshEdges); @@ -455,19 +469,24 @@ function InitNewMeshMap(obj,nSD,nMeshEdges) % co-registration with features meeting values % for the selected features. + if nargin < 4 + ySpan = [-1 1]*obj.SD_MAX_DEF; + end + if nargin < 3 - nSD = obj.SD_MAX_DEF; + xSpan = [-1 1]*obj.SD_MAX_DEF; end if nargin < 2 nMeshEdges = obj.SD_MESH_PTS_DEF; end - obj.sdMax = nSD; + obj.sdMax = obj.SD_MAX_DEF; obj.sdMeshPts = nMeshEdges; - obj.sdMeshEdges = linspace(-nSD,nSD,nMeshEdges); - [obj.sdMesh.X,obj.sdMesh.Y] = meshgrid(obj.sdMeshEdges,... - obj.sdMeshEdges); + obj.sdMeshEdges.X = linspace(xSpan(1),xSpan(2),nMeshEdges); + obj.sdMeshEdges.Y = linspace(ySpan(1),ySpan(2),nMeshEdges); + [obj.sdMesh.X,obj.sdMesh.Y] = meshgrid(obj.sdMeshEdges.X,... + obj.sdMeshEdges.Y); end @@ -513,12 +532,58 @@ function ButtonDownFcnSelect2D(obj,src,~) % obj.SetAxesWhereSpikesGo(ax); case 'alt' % Do "cluster cutting" (R-Click) obj.GetSpikesToMove(ax); - otherwise - return; + case 'normal' end end + % CALLBACK: Zoom in or out on featAx by mouse-wheel scroll + function WindowMouseWheel(obj,~,evt) + %WINDOWMOUSEWHEEL Zoom in or out on all plots + % + % fig.WindowScrollWheelFcn = @obj.WindowMouseWheel; + yl = ylim(obj.Features2D); + ySpan = diff(yl); + xl = xlim(obj.Features2D); + xSpan = diff(xl); + if strcmp(obj.CurrKeyPress,'alt') + newXSpan = xl; + newYSpan = yl + ySpan*.05*-evt.VerticalScrollCount; + elseif strcmp(obj.CurrKeyPress,'control') + newXSpan = xl + xSpan*.05*-evt.VerticalScrollCount; + newYSpan = yl; + else + mousepos = get (obj.Features2D, 'CurrentPoint'); + % only alt is pressed + x = mousepos(1); + y = mousepos(2); + newSpan = @(x,xl,xSpan) x - ((((x - xl).*[1 -1])./xSpan)... + .*xSpan*0.9^-evt.VerticalScrollCount).*[1 -1]; + newYSpan = newSpan(y,yl,ySpan); + newXSpan = newSpan(x,xl,xSpan); + end + + + set(obj.Features2D,'YLim',newYSpan,'XLim',newXSpan); + obj.InitNewMeshMap(obj.SD_MESH_PTS_DEF,newXSpan,newYSpan); + end + + function flag = clearKurrKey(obj,src,evt) + if any(strcmp(evt.Key,{'control','alt'})) + obj.CurrKeyPress={}; + else + idx = strcmp(obj.CurrKeyPress,evt.Key); + obj.CurrKeyPress(idx) = []; + end + flag = true; + end + + % CALLBACK: Execute keyboard shortcut on keyboard button press + function WindowKeyPress(obj,src,evt) + %WINDOWKEYPRESS Issue different events on keyboard presses + obj.CurrKeyPress = [obj.CurrKeyPress {evt.Key}]; + end + % CALLBACK: Rotate button click function RotateBtnPress(~,src,~) %ROTATEBTNPRESS Callback when "rotate button" is pressed @@ -585,14 +650,14 @@ function GetSpikesToMove(obj,ax) fi = ismember( obj.Data.class{curCh},find(obj.VisibleClusters)); X = feat*obj.projVecs'; [~,~,~,binX,binY] = histcounts2(X(:,1),X(:,2),... - obj.sdMeshEdges,obj.sdMeshEdges); + obj.sdMeshEdges.X,obj.sdMeshEdges.Y); % Draw polygon set(obj.Figure,'Pointer','circle'); % snipped_region = drawfreehand(ax,'Smoothing',5); axes(ax); - [h,x,y]=nigeLab.utils.freehanddraw(ax); + [h,x,y]=nigeLab.utils.freehanddraw(ax,'color',nigeLab.defaults.nigelColors('onsurface')); % pos = snipped_region.Position; % delete(snipped_region); delete(h) diff --git a/+nigeLab/+libs/@FeaturesUI/PlotFeatures.m b/+nigeLab/+libs/@FeaturesUI/PlotFeatures.m index 09b414d3..6d38eb11 100644 --- a/+nigeLab/+libs/@FeaturesUI/PlotFeatures.m +++ b/+nigeLab/+libs/@FeaturesUI/PlotFeatures.m @@ -1,6 +1,8 @@ function PlotFeatures(obj) %% PLOTFEATURES Plot cluster features from SORT UI in 3D and 2D. - +if ~isvalid(obj.Figure) + return +end % Loop through each subset of 3 features curCh = obj.ChannelSelector.Channel; @@ -9,9 +11,24 @@ function PlotFeatures(obj) % Get TIME binning vector for plotting features vs time (mins) ts = obj.Data.ts{curCh}; +xl = xlim(obj.Features2D); +yl = ylim(obj.Features2D); +cla(obj.Features2D); +set(obj.Features2D,... + 'XLim',xl,... + 'YLim',yl); + +xl = xlim(obj.Features3D); +yl = ylim(obj.Features3D); +zl = ylim(obj.Features3D); cla(obj.Features3D); -cla(obj.Features2D); +set(obj.Features3D,... + 'XLim',xl,... + 'YLim',yl,... + 'ZLim',zl); + +obj.PropLinker = linkprop([obj.Features2D obj.Features3D],{'XLim','YLim'}); % randomly select spikes to plot obj.rsel = false(size(obj.Data.feat{curCh},1),1); @@ -56,7 +73,7 @@ function PlotFeatures(obj) 'Visible',state,... 'UserData',iC); - line(obj.Features2D,... + l=line(obj.Features2D,... X(:,1), ... X(:,2), ... 'LineStyle', 'none',... @@ -92,11 +109,12 @@ function PlotFeatures(obj) 'UserData',iC,... 'ButtonDownFcn',@obj.ButtonDownFcnSelect2D); end + zlim(obj.Features3D,[0 max([zlim(obj.Features3D),ts(fi)'])]); end end drawnow; -obj.ResetFeatureAxes; - +% obj.ResetFeatureAxes; +obj.CountExclusions(obj.ChannelSelector.Channel); end \ No newline at end of file diff --git a/+nigeLab/+libs/@FeaturesUI/PlotQuality.m b/+nigeLab/+libs/@FeaturesUI/PlotQuality.m index 82e63b8d..288874a1 100644 --- a/+nigeLab/+libs/@FeaturesUI/PlotQuality.m +++ b/+nigeLab/+libs/@FeaturesUI/PlotQuality.m @@ -1,6 +1,8 @@ function PlotQuality(obj) %% PLOTFEATURES Plot cluster features from SORT UI in 3D and 2D. - +if ~isvalid(obj.Figure) + return; +end set(obj.Figure,'Pointer','watch'); obj.UpdateSil; % indx = obj.QualityIndx.UserData{2}; diff --git a/+nigeLab/+libs/@FeaturesUI/ResetFeatureAxes.m b/+nigeLab/+libs/@FeaturesUI/ResetFeatureAxes.m index 45fa2ef9..2818e13a 100644 --- a/+nigeLab/+libs/@FeaturesUI/ResetFeatureAxes.m +++ b/+nigeLab/+libs/@FeaturesUI/ResetFeatureAxes.m @@ -32,7 +32,4 @@ function ResetFeatureAxes(obj) 'climmode','manual',... 'clim',[0 1]); -obj.CountExclusions(obj.ChannelSelector.Channel); - - end \ No newline at end of file diff --git a/+nigeLab/+libs/@FeaturesUI/UpdateSil.m b/+nigeLab/+libs/@FeaturesUI/UpdateSil.m index d8eed353..8633ee40 100644 --- a/+nigeLab/+libs/@FeaturesUI/UpdateSil.m +++ b/+nigeLab/+libs/@FeaturesUI/UpdateSil.m @@ -11,6 +11,7 @@ % are also zero % sil = zeros(numel(activeToUpdate),length(feat)); % for ii=1:numel(activeToUpdate) +NNoiseIdx = cl(obj.rsel); scores = silhouette(feat(obj.rsel,:),cl(obj.rsel),measuretype{indx}); % end diff --git a/+nigeLab/+libs/@HighDimsUI/HighDimsUI.m b/+nigeLab/+libs/@HighDimsUI/HighDimsUI.m index 974167fb..7115cf5d 100644 --- a/+nigeLab/+libs/@HighDimsUI/HighDimsUI.m +++ b/+nigeLab/+libs/@HighDimsUI/HighDimsUI.m @@ -10,6 +10,8 @@ Figure matlab.ui.Figure panels + btnPanel + alpha num_dims selected_conds @@ -45,30 +47,47 @@ function PlotFig(obj) 'MenuBar','none',... 'CloseRequestFcn',@(~,~)obj.closeF,'Name','High-dimensional Navigator','numbertitle', 'off'); - p1 = uipanel( obj.Figure,... - 'BackgroundColor',nigeLab.defaults.nigelColors(0.1),... - 'Position',[.01 .01 .48 .98],... - 'BorderType','none'); - - p2 = uipanel( obj.Figure,... - 'BackgroundColor',nigeLab.defaults.nigelColors(0.1),... - 'Position',[.51 .01 .48 .98],... - 'BorderType','none'); + p1 = nigeLab.libs.nigelPanel( obj.Figure,... + 'PanelColor',nigeLab.defaults.nigelColors(0.1),... + 'Position',[.01 .22 .48 .78]); + + p2 = nigeLab.libs.nigelPanel( obj.Figure,... + 'PanelColor',nigeLab.defaults.nigelColors(0.1),... + 'Position',[.51 .22 .48 .78]); + + p3 = nigeLab.libs.nigelPanel( obj.Figure,... + 'PanelColor',nigeLab.defaults.nigelColors(0.1),... + 'Position',[.01 .01 .98 .2],... + 'MinTitleBarHeightPixels',0); + btnAx = axes('Position',[0 0 1 1],'Visible','off','XLim',[0 1],'YLim',[0 1]); + p3.nestObj(btnAx,'btnAx'); + btn1 = nigeLab.libs.nigelButton(btnAx,[0.01 0.1 0.2 0.3],'Best proj.',... + {@(~,~)obj.setViewCallback},'FontSize',.8 ); + dropdown = uicontrol('Style','popupmenu',... + 'Units','normalized',... + 'Position',[0.02 0.5 0.18 0.3],... + 'String',{'LDA','PCA','Random'}); + + p3.nestObj(dropdown,'projDropDown'); for ii = 1:obj.maxDim % max 15 features - ax1 = subplot(5,3,ii,'Parent',p1,'UserData',[ii],'xtick',[],'ytick',[],... - 'ButtonDownFcn', {@obj.Panel_ButtonDownFcn}); - ax2 = subplot(5,3,ii,'Parent',p2,'UserData',[15+ii],'xtick',[],'ytick',[],... + ax1 = subplot(5,3,ii,'Parent',obj.Figure,'UserData',[ii],'xtick',[],'ytick',[],... 'ButtonDownFcn', {@obj.Panel_ButtonDownFcn}); + p1.nestObj(ax1); + ax2 = copy(ax1); + ax2.UserData = 15+ii; + set(ax2,'ButtonDownFcn', {@obj.Panel_ButtonDownFcn}); + p2.nestObj(ax2); ax(ii) = ax1;ax(ii+15) = ax2; end - obj.panels = ax; + obj.panels = ax; InitUI(obj); obj.update_panels; obj.FeaturesUI.FeatX.Enable = 'off'; obj.FeaturesUI.FeatY.Enable = 'off'; - + + obj.btnPanel = p3; end function InitUI(obj) @@ -326,6 +345,208 @@ function closeF(obj) delete(obj.Figure); end + function SetLDAproj(obj) + % The desired vectors are from Fisher's linear discriminant analysis + % + % If only two conditions (should only be a one-d projection), qr gives us a + % second projection vector (that can be considered random) anyway, so it's + % ok + + D = obj.D; + idx = ismember(obj.D.class,find(obj.selected_conds)); + conditions = unique(D.class(idx)); + if numel(conditions)==1 + return; + end + + % TODO trajectories + % if (ismember('cluster', {D.type})) + % D = D(ismember({D.type}, 'cluster')); % get rid of any trajs + % + % if (length(conditions) <= 1) % LDA should do nothing for one condition + % return; + % end + % else + % if (length(conditions) == 1) % goal: only one cond, so make points in the traj far apart + % if (length(D) == 1) % only one trajectory, so LDA will fail + % return; + % end + % newData = []; + % index = 1; + % for itrial = 1:length(D) + % D(itrial).epochStarts = [D(itrial).epochStarts size(D(itrial).data,2)]; %change epochStarts to include the end + % for j = 1:length(D(1).epochStarts) + % newData(index).condition = num2str(j); + % newData(index).data = D(itrial).data(:,D(itrial).epochStarts(j)); + % index = index+1; + % end + % end + % D = newData; + % conditions = unique({D.condition}); + % end + % % else, just use the full trajectory + % end + % + + sigma = zeros(obj.num_dims); + m = []; + + for cond = 1:length(conditions) + thisCond = ismember(D.class, conditions(cond)); + sigma = sigma + cov(D.feat(thisCond,:),1); + m(:,cond) = mean(D.feat(ismember(D.class, conditions(cond)),:),1); + end + + Sigma_within = sigma ./ length(conditions); + Sigma_between = cov(m',1); + + [w e] = eig(Sigma_within \ Sigma_between); + + % LDA does *not* return orthogonal eigenvectors + % perform Gram-Schmidt to find nearest 2nd orthogonal vector + + [q r] = qr(w); + + rotate_to_desired(obj, q(:,1:2)'); + + end + + function SetRandProj(obj) + % rotates to a random set of projection vectors + + + [q r] = qr(randn(obj.num_dims)); + + rotate_to_desired(obj, q(:,1:2)'); + + end + + + function SetPCAProj(obj) + % The desired vectors are the first two PCs of the data + % (trajectories' datapoints are concatenated) + + D = obj.D; + + idx = ismember(obj.D.class,find(obj.selected_conds)); + + [u] = pca([D.feat(idx,:)]); + + rotate_to_desired(obj, u(:,1:2)'); + + end + + + end + + % Methods to precompute projections + methods + + function setViewCallback(obj) + dropdown = obj.btnPanel.Children{strcmp(obj.btnPanel.ChildName,'projDropDown')}; + switch dropdown.String{dropdown.Value} + case 'LDA' + SetLDAproj(obj); + case 'Random' + SetRandProj(obj) + case 'PCA' + SetPCAProj(obj) + end + + end + + + % % Third idea, what GGobi uses + function rotate_to_desired(obj, desired_vecs) +% helper function that modifies the current axes +% rotates the current projection vectors by a small angle until that +% the current vectors equal the desired vectors + + + + + % idea: + % I found a nice idea in (Section 2.2 Cook, Buja, Lee, and Wickham: Grand Tours, + % Projection Pursuit Guided TOurs and Manual Controls) which + % calculates the principal angles between the projection matrices + % + % don't use qr here...we need to keep exact vectors, not just the span + + + projVecs = obj.FeaturesUI.projVecs; + + + + % make sure desired_vecs are normalized + desired_vecs(1,:) = desired_vecs(1,:) ./ norm(desired_vecs(1,:)); + desired_vecs(2,:) = desired_vecs(2,:) ./ norm(desired_vecs(2,:)); + + check_desired = desired_vecs(1,:) * desired_vecs(2,:)'; + if (check_desired > 0) % desired vecs are not orthogonal + desired_vecs(2,:) = desired_vecs(2,:) - desired_vecs(1,:)*desired_vecs(2,:)' * desired_vecs(1,:); % subtract the parallel component + desired_vecs(2,:) = desired_vecs(2,:) ./ norm(desired_vecs(2,:)); + end + + check = projVecs * desired_vecs'; + if (check(1,1) > .99 && check(2,2) > .99) % this is the current view + return; % so no need to move, do nothing + end + + + % Find shortest path between spaces using SVD + [Va lambda Vz] = svd(projVecs * desired_vecs'); + + + % find principal directions in each space + Ba = projVecs' * Va; + Bz = desired_vecs' * Vz; + + + % orthonormalize Bz to get B_star, ensuring projections are orthogonal + % to Ba + % don't use qr...screwed me up, as it negates some vectors + B_star(:,1) = Bz(:,1) - Ba(:,1)'*Bz(:,1)*Ba(:,1) - Ba(:,2)'*Bz(:,1)*Ba(:,2); + B_star(:,1) = B_star(:,1) ./ norm(B_star(:,1)); + B_star(:,2) = Bz(:,2) - B_star(:,1)'*Bz(:,2)*B_star(:,1) - Ba(:,1)'*Bz(:,2)*Ba(:,1) - Ba(:,2)'*Bz(:,2)*Ba(:,2); + B_star(:,2) = B_star(:,2) ./ norm(B_star(:,2)); + + % calculate the principal angles + tau = acos(diag(lambda)); + + % increment angles + for t = linspace(0,1,60) + + % compute the rotation vector comprised of Ba and B_star + % as t-->1, Bt should converge to Bz + % also note that projection vectors are always orthogonal + Bt = []; + + Bt(:,1) = cos(tau(1)*t)*Ba(:,1) + sin(tau(1)*t)*B_star(:,1); + Bt(:,2) = cos(tau(2)*t)*Ba(:,2) + sin(tau(2)*t)*B_star(:,2); + + + + % need to rotate principal vectors back to original basis, + % so that initial projection begins with the old projection vectors + % (if not, you will still get same desired vectors, but the first + % projection will start rotated from the original projection) + projVecs = Va * Bt'; % transform back into the original coordinates + + % normalize projection vectors + projVecs(1,:) = projVecs(1,:) ./ norm(projVecs(1,:)); + projVecs(2,:) = projVecs(2,:) ./ norm(projVecs(2,:)); + + obj.FeaturesUI.projVecs = real(projVecs); + obj.FeaturesUI.PlotFeatures; + end + +% stop main axes from rotating + % Set the projection to exactly the desired vecs + obj.update_panels; +end + + + end end diff --git a/+nigeLab/+libs/@SortUI/SortUI.m b/+nigeLab/+libs/@SortUI/SortUI.m index e22e2ce7..e0f44dda 100644 --- a/+nigeLab/+libs/@SortUI/SortUI.m +++ b/+nigeLab/+libs/@SortUI/SortUI.m @@ -232,7 +232,7 @@ function addSpikeImage(obj) addlistener(obj.SpikeImage,'MainWindowClosed',... @(~,~)sortObj.delete), ... addlistener(obj.SpikeImage,'SaveData',... - @(~,~)sortObj.saveData)]; + @(src,evt)sortObj.saveData(evt.ch))]; end featLabel = parseFeatLabels(obj) % Parse pairs of feature names for dropdown labels list diff --git a/+nigeLab/+libs/@SpikeImage/SpikeImage.m b/+nigeLab/+libs/@SpikeImage/SpikeImage.m index e7b7b644..3cf6acc7 100644 --- a/+nigeLab/+libs/@SpikeImage/SpikeImage.m +++ b/+nigeLab/+libs/@SpikeImage/SpikeImage.m @@ -45,7 +45,8 @@ VisibleToggle % checkbox for selecting visiblity in the feature panel Axes % Axes containers for images Parent % Only set if called by nigeLab.Sort class object - + PropLinker + PlotCB; NumClus_Max = 9; CMap; @@ -61,8 +62,10 @@ 'progPctThresh',2,... 'nImg',11); - UnconfirmedChanges + ConfirmedChanges UnsavedChanges + + CurrKeyPress = {}; end % TRANSIENT,PROTECTED @@ -187,8 +190,8 @@ function delete(obj) obj.UpdateChannel; end - function flag = checkForUnconfirmedChanges(obj,askUser) - % CHECKFORUNCONFIRMEDCHANGES cheacks if there are any unconfirmed changes. + function flag = checkForConfirmedChanges(obj,askUser) + % CHECKFORConfirmedChanges cheacks if there are any unconfirmed changes. % also asks the suer if it's ok to proceed anyway. % returns true if we unconfirmed changes are present, therefore @@ -199,7 +202,8 @@ function delete(obj) flag = false; % Check if it's okay to lose changes if there are any - if obj.UnconfirmedChanges + if ~obj.ConfirmedChanges(obj.Parent.UI.ChannelSelector.Channel) && ... + obj.UnsavedChanges(obj.Parent.UI.ChannelSelector.Channel) if askUser str = questdlg('Unconfirmed changes will be lost. Proceed anyways?',... 'Discard Sorting on this Channel?','Yes','No','Yes'); @@ -223,10 +227,8 @@ function UpdateChannel(obj,~,~) obj.Flatten; % Construct figure - obj.Build; - - % New channel; no changes exist here yet - obj.UnconfirmedChanges = false; + obj.Build; + end function Refresh(obj) @@ -234,14 +236,11 @@ function Refresh(obj) if isa(obj.Parent,'nigeLab.Sort') % Set spike classes - obj.Assign(obj.Parent.spk.class{obj.Parent.UI.ch}); + evt = nigeLab.evt.assignClus(1:numel(obj.Spikes.Class),... + obj.Parent.spk.class{obj.Parent.UI.ch},... + obj.Spikes.Class); + obj.UpdateClusterAssignments(nan,evt); end - - % Flatten spike image - obj.Flatten; - - % Construct figure - obj.Build; end % Add a listener for cluster assignments @@ -345,8 +344,8 @@ function Init(obj,fs) % fs: Sample rate % No changes have been made yet - obj.UnconfirmedChanges = false; - obj.UnsavedChanges = false; + obj.ConfirmedChanges = false(1,numel(obj.Parent.UI.ChannelSelector.Menu.String)); + obj.UnsavedChanges = false(1,numel(obj.Parent.UI.ChannelSelector.Menu.String)); % Add sampling frequency obj.Spikes.fs = fs; @@ -437,11 +436,13 @@ function Build(obj) 'Color',nigeLab.defaults.nigelColors('background'),... 'WindowKeyPressFcn',@obj.WindowKeyPress,... 'WindowScrollWheelFcn',@obj.WindowMouseWheel,... - 'CloseRequestFcn',@obj.CloseSpikeImageFigure); + 'CloseRequestFcn',@obj.CloseSpikeImageFigure,... + 'WindowKeyReleaseFcn',@obj.clearKurrKey); else set(obj.Figure,'CloseRequestFcn',@obj.CloseSpikeImageFigure); set(obj.Figure,'WindowScrollWheelFcn',@obj.WindowMouseWheel); set(obj.Figure,'WindowKeyPressFcn',@obj.WindowKeyPress); + set(obj.Figure,'WindowKeyReleaseFcn',@obj.clearKurrKey); end % Set figure focus figure(obj.Figure); @@ -469,6 +470,7 @@ function Build(obj) end end + % Superimpose the image on everything fprintf(1,'->\tPlotting spikes'); for iC = 1:obj.NumClus_Max @@ -476,8 +478,21 @@ function Build(obj) obj.Draw(iC); end fprintf(1,'complete.\n\n'); + obj.PropLinker = [linkprop([obj.Axes{:}],{'YLim'})... + linkprop([obj.Images{:}],{'YData'});]; end + function flag = clearKurrKey(obj,src,evt) + if any(strcmp(evt.Key,{'control','alt'})) + obj.CurrKeyPress={}; + else + idx = strcmp(obj.CurrKeyPress,evt.Key); + obj.CurrKeyPress(idx) = []; + end + flag = true; + end + + % Initialize checkboxes on axes function InitCheckBoxes(obj,iC) %INITCHECKBOXES Initialize all checkboxes on axes @@ -521,22 +536,17 @@ function Draw(obj,plotNum) if nargin < 2 plotNum = 1:obj.NumClus_Max; else - plotNum = reshape(plotNum,1,numel(plotNum)); - end - - for iPlot = plotNum - set(obj.Images{iPlot},'CData',obj.Spikes.C{iPlot}); - set(obj.Axes{iPlot}.Title,'String',obj.PlotNames{iPlot}); - if obj.Spikes.CurClass == iPlot - obj.SetAxesHighlight(obj.Axes{iPlot},nigeLab.defaults.nigelColors('primary'),20); - else - obj.SetAxesHighlight(obj.Axes{iPlot},nigeLab.defaults.nigelColors('onsurface'),16); - end - set(obj.Axes{iPlot},'YLim',obj.YLim); - set(obj.Images{iPlot},'YData',obj.YLim); - drawnow; + plotNum = plotNum(:)'; end - + isThisAxes = plotNum == obj.Spikes.CurClass; + + cellfun(@(i,c)set(i,'CData',c),obj.Images(plotNum),obj.Spikes.C(plotNum)) + cellfun(@(a,t)set(a.Title,'String',t),obj.Axes(plotNum),obj.PlotNames(plotNum)) + cellfun(@(a)obj.SetAxesHighlight(a,nigeLab.defaults.nigelColors('primary'),20),obj.Axes(plotNum(isThisAxes))); + cellfun(@(a)obj.SetAxesHighlight(a,nigeLab.defaults.nigelColors('onsurface'),16),obj.Axes(plotNum(~isThisAxes))); + set(obj.Axes{1},'YLim',obj.YLim); + set(obj.Images{1},'YData',obj.YLim); + drawnow; end % Initializes a given axes (container for `SpikeImage` image) @@ -607,46 +617,30 @@ function Flatten(obj,plotNum) %FLATTEN Condense spikes into matrix scaled from 0 to 1 if nargin < 2 - + plotNum = 1:obj.NumClus_Max; obj.Spikes.C = cell(obj.NumClus_Max,1); % Colors (spike image) obj.Spikes.A = cell(obj.NumClus_Max,1); % Assignments - for iC = 1:obj.NumClus_Max - % Get bin edges - y_edge = linspace(obj.YLim(1),obj.YLim(2),obj.YPoints); - - % Pre-allocate - clus = obj.Spikes.Waves(obj.Spikes.Class==iC,:); - obj.Spikes.C{iC} = zeros(obj.YPoints-1,obj.XPoints); - obj.Spikes.A{iC} = nan(size(clus)); - for ii = 1:obj.XPoints - [obj.Spikes.C{iC}(:,ii),~,obj.Spikes.A{iC}(:,ii)] = ... - histcounts(clus(:,ii),y_edge); - end - - % Normalize - obj.Spikes.C{iC} = obj.Spikes.C{iC}./... - max(max(obj.Spikes.C{iC})); - end - else - plotNum = reshape(plotNum,1,numel(plotNum)); - for iC = plotNum - % Get bin edges - y_edge = linspace(obj.YLim(1),obj.YLim(2),obj.YPoints); - - % Pre-allocate - clus = obj.Spikes.Waves(obj.Spikes.Class==iC,:); - obj.Spikes.C{iC} = zeros(obj.YPoints-1,obj.XPoints); - obj.Spikes.A{iC} = nan(size(clus)); - for ii = 1:obj.XPoints - [obj.Spikes.C{iC}(:,ii),~,obj.Spikes.A{iC}(:,ii)] = ... - histcounts(clus(:,ii),y_edge); - end + end - % Normalize - obj.Spikes.C{iC} = obj.Spikes.C{iC}./... - max(max(obj.Spikes.C{iC})); - end + plotNum = plotNum(:)'; + % Get bin edges + y_edge = linspace(obj.YLim(1),obj.YLim(2),obj.YPoints); + + for iC = plotNum + % Pre-allocate + clus = obj.Spikes.Waves(obj.Spikes.Class==iC,:); + obj.Spikes.C{iC} = zeros(obj.YPoints-1,obj.XPoints); + obj.Spikes.A{iC} = nan(size(clus)); + for ii = 1:obj.XPoints + [obj.Spikes.C{iC}(:,ii),obj.Spikes.A{iC}(:,ii)] = ... + matlab.internal.math.histcounts(clus(:,ii),y_edge); + end + + % Normalize + obj.Spikes.C{iC} = obj.Spikes.C{iC}./... + max(obj.Spikes.C{iC}(:)); end + end % CALLBACK: Triggered when figure window closes @@ -731,7 +725,7 @@ function GetSpikesToMove(obj,curAxes) set(obj.Figure,'Pointer','circle'); % snipped_region = imfreehand(curAxes); - [h,x,y]=nigeLab.utils.freehanddraw(curAxes); + [h,x,y]=nigeLab.utils.freehanddraw(curAxes,'color',nigeLab.defaults.nigelColors('onsurface')); % pos = getPosition(snipped_region); % delete(snipped_region); delete(h); @@ -770,7 +764,14 @@ function GetSpikesToMove(obj,curAxes) % CALLBACK: Execute keyboard shortcut on keyboard button press function WindowKeyPress(obj,src,evt) %WINDOWKEYPRESS Issue different events on keyboard presses + obj.CurrKeyPress = [obj.CurrKeyPress {evt.Key}]; switch evt.Key + case 'h' + thisClass = obj.Spikes.CurClass; + thisAx = obj.Axes{thisClass}'; + val = obj.VisibleToggle{thisClass}.Value; + obj.VisibleToggle{thisClass}.Value = ~val; + obj.SetVisibleFeatures(thisClass,~val); case {'n','0'} if strcmpi(evt.Modifier,'control') thisClass = obj.Spikes.CurClass; @@ -795,11 +796,18 @@ function WindowKeyPress(obj,src,evt) obj.ConfirmChanges; case 'z' if strcmpi(evt.Modifier,'control') - obj.UndoChanges; +% obj.UndoChanges; + uiundo(obj.Figure,'execUndo'); + end + case 'y' + if strcmpi(evt.Modifier,'control') + uiundo(obj.Figure,'execRedo'); + end case 's' if strcmpi(evt.Modifier,'control') - if obj.UnconfirmedChanges + if obj.UnsavedChanges(obj.Parent.UI.ChannelSelector.Channel) &&... + ~obj.ConfirmedChanges(obj.Parent.UI.ChannelSelector.Channel) str = questdlg('Confirm current changes before save?',... 'Use Most Recent Scoring?','Yes','No','Yes'); else @@ -812,8 +820,6 @@ function WindowKeyPress(obj,src,evt) obj.SaveChanges; end - case 'escape' - notify(obj,'MainWindowClosed'); case {'x','c'} if strcmpi(evt.Modifier,'control') notify(obj,'MainWindowClosed'); @@ -822,7 +828,12 @@ function WindowKeyPress(obj,src,evt) if strcmpi(evt.Modifier,'control') obj.Recluster; else - obj.Refresh; + ButtonName = questdlg(sprintf('Reload data from disk?\nAll unsaved changes will be lost.'), ... + 'Reload', ... + 'Yes', 'No', 'No'); + if strcmp(ButtonName,'Yes') + obj.Refresh; + end end case 'q' obj.Parent.UI.FeaturesUI.ReopenWindow; @@ -861,8 +872,15 @@ function SaveChanges(obj) % % obj.SaveChanges(); - notify(obj,'SaveData'); - obj.UnsavedChanges = false; + confirmedChans = find(obj.ConfirmedChanges); + if isempty(confirmedChans) + fprintf(1,'No changes were done to the scoring.\nNothing was saved!\n'); + return; + end + evt = nigeLab.evt.saveData(confirmedChans); + notify(obj,'SaveData',evt); + obj.UnsavedChanges(confirmedChans) = false; + obj.ConfirmedChanges(confirmedChans) = false; disp('Scoring saved.'); end @@ -877,7 +895,8 @@ function UndoChanges(obj) else obj.Spikes.Class = obj.Parent.spk.class{1}; end - obj.UnconfirmedChanges = false; + obj.ConfirmedChanges(obj.Parent.UI.ChannelSelector.Channel) = false; + obj.UnsavedChanges(obj.Parent.UI.ChannelSelector.Channel) = false; obj.Flatten; obj.SetPlotNames; obj.Draw; @@ -900,7 +919,7 @@ function ConfirmChanges(obj) else obj.Parent.spk.class{1} = obj.Spikes.Class; end - obj.UnconfirmedChanges = false; + obj.ConfirmedChanges(obj.Parent.UI.ChannelSelector.Channel) = true; notify(obj,'ChannelConfirmed'); fprintf(1,'Scoring for channel %d confirmed.\n',... obj.Parent.UI.ChannelSelector.Channel); @@ -911,39 +930,79 @@ function WindowMouseWheel(obj,~,evt) %WINDOWMOUSEWHEEL Zoom in or out on all plots % % fig.WindowScrollWheelFcn = @obj.WindowMouseWheel; - - obj.YLim(1) = min(obj.YLim(1) + 20*evt.VerticalScrollCount,-20); - obj.YLim(2) = max(obj.YLim(2) - 10*evt.VerticalScrollCount,10); - obj.Spikes.Y = linspace(obj.YLim(1),obj.YLim(2),obj.YPoints-1); + modifier = obj.CurrKeyPress; + if ~isempty(modifier) && all(strcmp(modifier,'alt')) + % only alt is pressed + offset = 20*evt.VerticalScrollCount; + obj.YLim = obj.YLim + offset; + obj.Spikes.Y = obj.Spikes.Y + offset; + + else + obj.YLim(1) = min(obj.YLim(1) + 20*evt.VerticalScrollCount,-20); + obj.YLim(2) = max(obj.YLim(2) - 10*evt.VerticalScrollCount,10); + obj.Spikes.Y = linspace(obj.YLim(1),obj.YLim(2),obj.YPoints-1); + end obj.Flatten; obj.Draw; end % CALLBACK: Updates cluster assignments - function UpdateClusterAssignments(obj,~,evt) + function UpdateClusterAssignments(obj,~,evt,logundo) %UPDATECLUSTERASSIGNMENTS Update cluster assigns and notify - + if nargin<4 + logundo = true; + end if ~isprop(evt,'subs') return; + elseif isempty(evt.subs) + return; end + % identify new classes to assign + newClasses = unique(evt.class); + newClasses = newClasses(:)'; + % log status quo for undo/redo + oldEvt = nigeLab.evt.assignClus(evt.subs,... + obj.Spikes.Class(evt.subs),... + unique(evt.class)); + % Identify plots to update plotsToUpdate = unique(obj.Spikes.Class(evt.subs)); plotsToUpdate = reshape(plotsToUpdate,1,numel(plotsToUpdate)); - plotsToUpdate = unique([plotsToUpdate, obj.Spikes.CurClass]); % in case + plotsToUpdate = unique([plotsToUpdate, newClasses]); % in case if ~isnan(evt.otherClassToUpdate) plotsToUpdate = unique([plotsToUpdate, evt.otherClassToUpdate]); end - % Assign and redo graphics - obj.Assign(obj.Spikes.CurClass,evt.subs); + for class = newClasses + % Assign and redo graphics + idx = evt.class == class; + obj.Assign(class,evt.subs(idx)); + end obj.SetPlotNames(plotsToUpdate); obj.Flatten(plotsToUpdate); obj.Draw(plotsToUpdate); % Indicate that there have been some changes in class ID - obj.UnconfirmedChanges = true; - obj.UnsavedChanges = true; + obj.ConfirmedChanges(obj.Parent.UI.ChannelSelector.Channel) = false; + obj.UnsavedChanges(obj.Parent.UI.ChannelSelector.Channel) = true; notify(obj,'ClassAssigned',evt); + + % Add undo redo functionality + % Prepare an undo/redo action + source = sprintf('%d,',plotsToUpdate); + destination = sprintf('%d,',newClasses); + cmd.Name = sprintf('Cluster assignment (%g to %g)',source(1:end-1),destination(1:end-1)); + + cmd.Function = @obj.UpdateClusterAssignments; % Redo action + cmd.Varargin = {nan,evt,false}; + cmd.InverseFunction = @obj.UpdateClusterAssignments; % Undo action + + cmd.InverseVarargin = {nan,oldEvt,false}; +% % Register the undo/redo action with the figure + if logundo + uiundo(obj.Figure,'function',cmd); + end + end % Assign spikes to a given class diff --git a/+nigeLab/+libs/@configSD/configSD.m b/+nigeLab/+libs/@configSD/configSD.m new file mode 100644 index 00000000..0a5c399c --- /dev/null +++ b/+nigeLab/+libs/@configSD/configSD.m @@ -0,0 +1,504 @@ +classdef configSD < handle + %% CONFIGSD class. It summons a UI where you can configure the SD parameters prior to SD execution + properties + UI + Listeners + Channels + + Fig + DataAx + artPlot + spkPlot + + SpikeAx + SDParsPanel + ArtRejParsPanel + BtnPanel + SliderLbl + durSlider + SDBtn + ARTBtn + ZoomBtn + ResampleBtn + ExportBtn + + AllTextBoxes + + ExBlock + SDMethods + ARTMethods + Pars + + startIdx + endIdx + + data + artRejData + end + + + methods + + function obj = configSD(nigelObj) + switch nigelObj.Type + case 'Tank' + idxA = randi(numel(nigelObj.Children)); + idxB = randi(numel(nigelObj.Children(idxA).Children)); + obj.ExBlock = nigelObj{idxA,idxB}; + case 'Animal' + idx = randi(numel(nigelObj.Children)); + obj.ExBlock = nigelObj{idx}; + case 'Block' + obj.ExBlock = nigelObj; + end + + ff = fieldnames(obj.ExBlock.Pars.SD); + obj.SDMethods = ff(cellfun(@(fName)numel(fName)>2 && strcmp(fName(1:3),'SD_') ,ff)); + obj.SDMethods = cellfun(@(MName) MName(4:end),obj.SDMethods, 'UniformOutput',false); + + obj.ARTMethods = ff(cellfun(@(fName)numel(fName)>2 && strcmp(fName(1:4),'ART_') ,ff)); + obj.ARTMethods = cellfun(@(MName) MName(5:end),obj.ARTMethods, 'UniformOutput',false); + + obj.Pars = obj.ExBlock.Pars.SD; + + % Build UI + obj.buildGUI(); + + obj.Channels.Name = {obj.ExBlock.Channels.name}; + obj.Channels.Selected = 1; + obj.UI.Parent = obj; + obj.UI.channel = 1; + + obj.UI.ChannelSelector = nigeLab.libs.ChannelUI(obj.UI); + obj.Listeners = [obj.Listeners, ... + addlistener(obj.UI.ChannelSelector,'NewChannel',... + @obj.setChannel)]; + + + obj.sampleData() + + + end + + function flag = set(~,field,~) + % This does not serve any purpose. Just to keep the channel + % selector happy. For whatever reason was requiring this to be + % in place + if strcmp(field,'channel') + flag = true; + else + flag = false; + end + end + + function setChannel(obj,src,evt) + obj.Channels.Selected = src.Channel; + sampleData(obj); + renderData(obj) + end + + function renderData(obj) + cla(obj.DataAx,'reset'); + cla(obj.SpikeAx,'reset'); + + thisBlock = obj.ExBlock; + fs = thisBlock.SampleRate; + t = (obj.startIdx:obj.endIdx)./fs; + + + plot(obj.DataAx,t,obj.data); + xlim(obj.DataAx,[t(1) min(t(1)+60,t(end))]); + obj.DataAx.YAxis.Color = [1,1,1];obj.DataAx.XAxis.Color = [1,1,1]; + title(obj.DataAx,sprintf('Filtered Data, %.2f minutes',numel(t)./fs./60),'Color',[1 1 1]); + + end + + function sampleData(obj) + + waitFig = plotWaitFigure(obj,'Sampling data...'); + lockedObjs = obj.lockUnlockGui([],'off'); + + thisBlock = obj.ExBlock; + obj.data = []; + obj.artRejData = []; + + L = thisBlock.Samples; + fs = thisBlock.SampleRate; + Samples = floor(obj.durSlider.Value*60*fs); % length of ther selected recording + obj.startIdx = randi(max(1,L-Samples)); + obj.endIdx = min(L,obj.startIdx+Samples); + chanIdx = obj.Channels.Selected; + if L<60*fs + % removing durSlider from the list of objs to reenable + lockedObjs(lockedObjs == obj.durSlider) = []; + lockedObjs(lockedObjs == obj.SliderLbl) = []; + obj.durSlider.Enable = 'off'; + obj.SliderLbl.Enable = 'off'; + obj.SliderLbl.String = 'less then 1min data.'; + end + + isCar = thisBlock.getStatus('CAR'); + isFilt = thisBlock.getStatus('Filt'); + if isCar + obj.data = thisBlock.Channels(chanIdx).CAR(obj.startIdx:obj.endIdx); + elseif isFilt + obj.data = thisBlock.Channels(chanIdx).Filt(obj.startIdx:obj.endIdx); + else + error('Niether Filt nor Car detected.%cYou need to have ate least your signal filtered for spike detection!',newline); + end + renderData(obj); + + obj.lockUnlockGui(lockedObjs,'on'); + delete(waitFig); + end + + function delete(obj) + delete(obj.UI.ChannelSelector); + delete(obj.UI.Fig); + end + + function buildGUI(obj,fig) + % BUILDGUI Build the graphical interface + % + % fig = obj.buildGUI(); Constructs new figure and adds panels + % + % fig = obj.buildGUI(fig); Adds the panels only + + if nargin < 2 + obj.UI.Fig = figure(... + 'Toolbar','none',... + 'MenuBar','none',... + 'NumberTitle','off',... + 'Units','pixels',... + 'Position',[1500 200 800 700],... + 'Color',nigeLab.defaults.nigelColors('bg'),... + 'Visible','on',... + 'CloseRequestFcn',@(~,~) obj.delete); + fig = obj.UI.Fig; + else + obj.UI.Fig = fig; + end + + + obj.DataAx = axes(fig,'Units','normalized','Position',[.05 .7 .9 .25],'Box','off'); + obj.ZoomBtn = uicontrol(fig,'Style','togglebutton',... + 'String','Z',... + 'Units','normalized','Position',[.05 .96 .03 .03],'Callback',@(~,~)zoom(fig)); + jh = nigeLab.utils.findjobj(obj.ZoomBtn); + jh.setBorderPainted(false); + jh.setContentAreaFilled(false); + zbtnImgPath = fullfile(matlabroot,'toolbox\matlab\appdesigner\web\release\images\figurefloatingpalette','zoomin_cursor3D.png'); + if exist(zbtnImgPath,'file') + img = imread(zbtnImgPath); + set(obj.ZoomBtn,'CData',img); + end + obj.SpikeAx = axes(fig,'Units','normalized','Position',[.75 .35 .2 .28],'Box','off'); + title(obj.SpikeAx,'Spikes','Color',[1 1 1]); + + obj.SpikeAx.YAxis.Color = [1,1,1];obj.SpikeAx.XAxis.Color = [1,1,1]; + obj.BtnPanel =uipanel(fig,'Position',[.75 .05 .2 .25]); + + + obj.SliderLbl = uicontrol(obj.BtnPanel,'Style','text',... + 'Units','normalized',... + 'Position',[.1 .7 .8 .1],... + 'String','1 min',... + 'HorizontalAlignment','left'); + obj.durSlider = uicontrol(obj.BtnPanel,'Style','slider',... + 'Units','normalized',... + 'Position',[.1 .8 .8 .15],... + 'Value',1,'Min',1,'Max',10,... + 'Callback',@(hObj,~,~)set(obj.SliderLbl,'String',sprintf('%.2f Min',hObj.Value))); + obj.ResampleBtn = uicontrol(obj.BtnPanel,'Style','pushbutton','String','Resample',... + 'Units','normalized',... + 'Position',[.5 .5 .35 .2],... + 'Callback',@(~,~)obj.sampleData()); + + + obj.SDBtn = uicontrol(obj.BtnPanel,'Style','pushbutton','String','SD',... + 'Units','normalized',... + 'Position',[.1 .3 .35 .2],... + 'Callback',@(~,~)obj.StartSD); + + obj.ARTBtn = uicontrol(obj.BtnPanel,'Style','pushbutton','String','ArtRej',... + 'Units','normalized',... + 'Position',[.5 .3 .35 .2],... + 'Callback',@(~,~)obj.StartArtRej); + + obj.ExportBtn = uicontrol(obj.BtnPanel,'Style','pushbutton','String','Export',... + 'Units','normalized',... + 'Position',[.1 .05 .75 .2],... + 'Callback',@(~,~)obj.ExportPars); + + obj.SDParsPanel = uitabgroup(fig,... + 'Units','normalized',... + 'Position',[.05 .25 .65 .4]); + + obj.ArtRejParsPanel = uitabgroup(fig,... + 'Units','normalized',... + 'Position',[.05 .05 .65 .18]); + + % fill SD panel + for ii = 1:numel(obj.SDMethods) + % Add all uitabs + thisTab = uitab(obj.SDParsPanel,'Title',obj.SDMethods{ii}); + thisPars = (obj.ExBlock.Pars.SD.(['SD_' obj.SDMethods{ii}])); + thisParsNames = fieldnames(thisPars); + thisTab.Units = 'pixels'; + Width = thisTab.InnerPosition(3); + thisTab.Units = 'normalized'; + for jj=1:numel(thisParsNames) + Hoff = floor(Width/3)*mod(jj-1,3); + Woff = -40*idivide(int16(jj-1),3); + thisParText = uicontrol('Style','edit',... + 'Parent',thisTab,... + 'Units','pixels',... + 'Position',[38 + Hoff 215 + Woff 37 18],... + 'String',nigeLab.utils.ToString(thisPars.(thisParsNames{jj})),... + 'Callback',@(hObj,~,~)textBoxCallback(obj,hObj,['SD_' obj.SDMethods{ii}],thisParsNames{jj})); + + thisParLbl = uicontrol('Style','text',... + 'Parent',thisTab,... + 'Units','pixels',... + 'Position',[38 + Hoff 235 + Woff 70 18],... + 'String',thisParsNames{jj},... + 'HorizontalAlignment','left'); + + end + end + obj.SDParsPanel.SelectedTab = findobj(obj.SDParsPanel,... + 'Title',obj.Pars.SDMethodName); + + % fill artefact rejection panel + for ii = 1:numel(obj.ARTMethods) + % Add all uitabs + thisTab = uitab(obj.ArtRejParsPanel,'Title',obj.ARTMethods{ii}); + thisPars = (obj.ExBlock.Pars.SD.(['ART_' obj.ARTMethods{ii}])); + thisParsNames = fieldnames(thisPars); + thisTab.Units = 'pixels'; + Width = thisTab.InnerPosition(3); + thisTab.Units = 'normalized'; + for jj=1:numel(thisParsNames) + Hoff = floor(Width/3)*mod(jj-1,3); + Woff = -40*idivide(int16(jj-1),3); + thisParText = uicontrol('Style','edit',... + 'Parent',thisTab,... + 'Units','pixels',... + 'Position',[38 + Hoff 60 + Woff 37 18],... + 'String',nigeLab.utils.ToString(thisPars.(thisParsNames{jj})),... + 'Callback',@(hObj,~,~)textBoxCallback(obj,hObj,['ART_' obj.ARTMethods{ii}],thisParsNames{jj})); + + thisParLbl = uicontrol('Style','text',... + 'Parent',thisTab,... + 'Units','pixels',... + 'Position',[38 + Hoff 78 + Woff 70 18],... + 'String',thisParsNames{jj},... + 'HorizontalAlignment','left'); + + end + end + obj.ArtRejParsPanel.SelectedTab = findobj(obj.ArtRejParsPanel,... + 'Title',obj.Pars.ArtefactRejMethodName); + + end + + function ExportPars(obj) + SDAlgName = obj.SDParsPanel.SelectedTab.Title; + ArtRejAlgName = obj.ArtRejParsPanel.SelectedTab.Title; + + obj.Pars.ArtefactRejMethodName = ArtRejAlgName; + obj.Pars.SDMethodName = SDAlgName; + + obj.ExBlock.Pars.SD = obj.Pars; + obj.ExBlock.saveParams; + end + + function textBoxCallback(obj,hObj,MethodName,ParName) + val = get(hObj,'String'); + thiPar = obj.ExBlock.Pars.SD.(MethodName).(ParName); + + if isnumeric(thiPar) + convertedVal = str2double(val); + + elseif isa(thiPar,'function_handle') + if ismember(exist(thiPar),[2 3 5 6]) + convertedVal = str2func(val); + else + error('The defined function does not exist!/n'); + end + + elseif ischar(thiPar) + convertedVal = val; + + end + + obj.Pars.(MethodName).(ParName) = convertedVal; + end + + function AllObjs = lockUnlockGui(obj,objs,status) + if nargin<2 + % no additional inputs. Looks for all objs and toggles the + % status + AllObjs = findobj(obj.UI.Fig,'type','uicontrol'); + lockUnlockGui(obj,AllObjs); + + elseif nargin <3 + % All inputs are provided. Sets the enabled property of the + % objs provided to status + onObjs = findobj(objs,'type','uicontrol','enable','on'); + offObjs = findobj(objs,'type','uicontrol','enable','off'); + lockUnlockGui(obj,onObjs,'off'); + lockUnlockGui(obj,offObjs,'on'); + AllObjs = [onObjs(:);offObjs(:)]; + elseif nargin <4 && isempty(objs) + % Second input is []. Looks for all objs with enable = + % ~status and sets the status to status. + % Ex obj.lockUnlockGui([],'on') enables all disabled controls + switch status + case {true,'on'} + Currentstatus = 'off'; + case {false,'off'} + Currentstatus = 'on'; + end + AllObjs = findobj(obj.UI.Fig,'type','uicontrol','enable',Currentstatus); + lockUnlockGui(obj,AllObjs,status); + else + switch status + case {true,'on'} + status = 'on'; + case {false,'off'} + status = 'off'; + end + set(objs,'Enable',status); + end + end + + function waitFig = plotWaitFigure(obj,message) + waitFig = figure('Units','pixels','MenuBar','none','ToolBar','none','NumberTitle','off','CloseRequestFcn',[]); + waitFig.Position = [obj.UI.Fig.Position(1:2) + obj.UI.Fig.Position(3:4)./2 272 40]; + waitFig.Name = 'Wait! I''m thinking...'; + uicontrol(waitFig,'Style','text','Units','normalized','String',message,'Position',[.05 .05 .9 .9]); + drawnow; + figAlwaysOnTop(obj,waitFig,true) + end + + function StartSD(obj) + if isempty(obj.artRejData) + obj.StartArtRej; + end + + waitFig = plotWaitFigure(obj,'Detecting spikes...'); + lockedObjs = obj.lockUnlockGui([],'off'); + + AlgName = obj.SDParsPanel.SelectedTab.Title; + SDFun = ['SD_' AlgName]; + SDPars = obj.Pars.(SDFun); + SDPars.fs = obj.ExBlock.SampleRate; + SDargsout = obj.ExBlock.testSD(SDFun,obj.artRejData,SDPars); + + tIdx = SDargsout{1}; + peak2peak = SDargsout{2}; + peakAmpl = SDargsout{3}; + peakWidth = SDargsout{4}; + + plotSpikes(obj,tIdx,peakAmpl,peakWidth) + + obj.lockUnlockGui(lockedObjs,'on'); + delete(waitFig); + end + + function StartArtRej(obj) + waitFig = plotWaitFigure(obj,'Detecting artefacts...'); + lockedObjs = obj.lockUnlockGui([],'off'); + + AlgName = obj.ArtRejParsPanel.SelectedTab.Title; + + ArtFun = ['ART_' AlgName]; + ArtPars = obj.Pars.(ArtFun); + ArtPars.fs = obj.ExBlock.SampleRate; + Artargsout = obj.ExBlock.testSD(ArtFun,obj.data,ArtPars); + obj.artRejData = Artargsout{1}; + artifact = Artargsout{2}; + + plotArt(obj,artifact) + + + obj.lockUnlockGui(lockedObjs,'on'); + delete(waitFig); + end + + function plotArt(obj,art) + hold(obj.DataAx,'on'); + fs = obj.ExBlock.SampleRate; + t = (obj.startIdx:obj.endIdx)./fs; + delete([obj.artPlot,obj.spkPlot]); + cla(obj.SpikeAx); + obj.artPlot = plot(obj.DataAx,t(art),obj.data(art),'og'); + legend(obj.DataAx,{'Data','Artifacts'}); + end + + function plotSpikes(obj,tIdx,peakAmpl,peakWidth) + fs = obj.ExBlock.SampleRate; + t = (obj.startIdx:obj.endIdx)./fs; + + WindowPreSamples = obj.Pars.WPre * 1e-3 * fs; + WindowPostSamples = obj.Pars.WPost * 1e-3 * fs; + out_of_record = tIdx <= WindowPreSamples+1 | tIdx >= length(obj.data) - WindowPostSamples - 2; + + peakAmpl(out_of_record) = []; + peakWidth(out_of_record) = []; + tIdx(out_of_record) = []; + + + hold(obj.DataAx,'on'); + delete(obj.spkPlot); + obj.spkPlot = plot(obj.DataAx,t(tIdx),peakAmpl,'*r'); + legend(obj.DataAx,{'Data','Artifacts','Spikes'}); + % BUILD SPIKE SNIPPET ARRAY AND PEAK_TRAIN + tIdx = tIdx(:); % make sure it's vertical + cla(obj.SpikeAx); + if (any(tIdx)) % If there are spikes in the current signal + snippetIdx = (-WindowPreSamples : WindowPostSamples) + tIdx; + spikes = obj.data(snippetIdx); + plot(obj.SpikeAx,(-WindowPreSamples : WindowPostSamples)./fs,spikes); + obj.SpikeAx.YAxis.Color = [1,1,1];obj.SpikeAx.XAxis.Color = [1,1,1]; + else + spikes = []; + end + + title(obj.SpikeAx,sprintf('%d spikes',size(spikes,1)),'Color',[1 1 1]); + + end + + function figAlwaysOnTop(obj,fig,mode) + if nargin < 2 + mode = false; + end + % toggle figure alway on top. Has to be changed before + % the visibility property + % drawnow; + warningTag = 'MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame'; + warning('off',warningTag); + + jFrame = get(handle(fig),'JavaFrame'); + jFrame_fHGxClient = jFrame.fHG2Client; + jFrame_fHGxClientW=jFrame_fHGxClient.getWindow; + tt = tic; + while(isempty(jFrame_fHGxClientW)) + if toc(tt) > 5 + % this is not critical, no need to lock everythong here if + % it's not working + warning('on',warningTag); + return; + end + jFrame_fHGxClientW=jFrame_fHGxClient.getWindow; + end + + jFrame_fHGxClientW.setAlwaysOnTop(mode); + warning('on',warningTag); + end + + end + + +end \ No newline at end of file diff --git a/+nigeLab/+libs/@nigelButton/nigelButton.m b/+nigeLab/+libs/@nigelButton/nigelButton.m index 4cdff5ec..d1a6611d 100644 --- a/+nigeLab/+libs/@nigelButton/nigelButton.m +++ b/+nigeLab/+libs/@nigelButton/nigelButton.m @@ -430,7 +430,15 @@ function delete(b) b.Label.Color = b.SelectedColor; b.Position = b.SelectedPosition; uistack(b.Group,'top'); - set(b.Figure,'WindowButtonUpFcn',@(obj,~,~)ButtonUpFcn(b)); + if iscell(b.Figure.WindowButtonUpFcn) + oldfun = b.Figure.WindowButtonUpFcn; + else + oldfun = {b.Figure.WindowButtonUpFcn}; + end + fcnList = {oldfun,... + {@(obj,~,~)ButtonUpFcn(b)}}; + set(b.Figure,'WindowButtonUpFcn',... + {@nigeLab.utils.multiCallbackWrap, fcnList}); case 'off' if strcmp(b.Hovered,'on') b.Hovered = 'off'; @@ -727,11 +735,14 @@ function ButtonUpFcn(b) return; end + fcnList = b.Figure.WindowButtonUpFcn{2}; + originalFcn = fcnList{1}; if strcmpi(b.Selected,'off') - set(b.Figure,'WindowButtonUpFcn',[]); + % revert the modification to windowbuttonupfcn + set(b.Figure,'WindowButtonUpFcn',originalFcn); return; elseif strcmpi(b.Enable,'off') - set(b.Figure,'WindowButtonUpFcn',[]); + set(b.Figure,'WindowButtonUpFcn',originalFcn); return; end % If button is released, turn off "highlight border" diff --git a/+nigeLab/+libs/@nigelPanel/nigelPanel.m b/+nigeLab/+libs/@nigelPanel/nigelPanel.m index d8051f3d..473fc261 100644 --- a/+nigeLab/+libs/@nigelPanel/nigelPanel.m +++ b/+nigeLab/+libs/@nigelPanel/nigelPanel.m @@ -434,7 +434,6 @@ function buildTitleBar(obj) % BUILDTITLEBAR Build graphics for "nice header box" % % obj.buildTitleBar; - if isempty(obj.TitleBar) obj.TitleBar = struct; obj.TitleBar.axes = axes(obj.OutPanel,... % formerly "a" diff --git a/+nigeLab/+utils/@LinePlotReducer/reduce_to_width.m b/+nigeLab/+utils/@LinePlotReducer/reduce_to_width.m new file mode 100644 index 00000000..88f5ae42 --- /dev/null +++ b/+nigeLab/+utils/@LinePlotReducer/reduce_to_width.m @@ -0,0 +1,113 @@ +function [x_reduced, y_reduced] = reduce_to_width(x, y, width, lims) + +% [x_reduced, y_reduced] = reduce_to_width(x, y, width, lims) +% +% (This function is primarily used by LinePlotReducer, but has been +% provided as a stand-alone function outside of that class so that it can +% be used potentially in other projects. See help LinePlotReducer for +% more.) +% +% Reduce the data contained in x and y for display on width pixels, +% stretching from lims(1) to lims(2). +% +% For example, if x and y are 20m points long, but we only need to plot +% them over 500 pixels from x=5 to x=10, this function will return 1000 +% points representing the first point of the window (x=5), the last point +% of the window (x=10), and the maxima and minima for each pixel between +% those two points. The result is that x_reduced and y_reduced can be +% plotted and will look exactly like x and y, without all of the waste of +% plotting too many points. +% +% x can be n-by-1 or n-by-m with n samples of m columns (that is, there can +% be 1 x for all y or 1 x for each y. +% +% y must be n-by-m with n samples of m columns +% +% [xr, yr] = reduce_to_width(x, y, 500, [5 10]); % Reduce the data. +% +% plot(xr, yr); % This contains many fewer points than plot(x, y) but looks +% the same. +% +% Tucker McClure +% Copyright 2013, The MathWorks, Inc. + + % We'll need the first point to the left of the limits, the first point + % to the right to the right of the limits, and the min and max at every + % pixel inbetween. That's 1 + 1 + 2*(width - 2) = 2*width total points. + n_points = 2*width; + + % If the data is already small, there's no need to reduce. + if size(x, 1) <= n_points + x_reduced = x; + y_reduced = y; + return; + end + + % Reduce the data to the new axis size. + x_reduced = nan(n_points, size(y, 2)); + y_reduced = nan(n_points, size(y, 2)); + for k = 1:size(y, 2) + + % Find the starting and stopping indices for the current limits. + if k <= size(x, 2) + + % Rename the column. This actually makes stuff substantially + % faster than referencing x(:, k) all over the place. On my + % timing trials, this was 20x faster than referencing x(:, k) + % in the loop below. + xt = x(:, k); + + % Map the lower and upper limits to indices. + nx = size(x, 1); + lower_limit = binary_search(xt, lims(1), 1, nx); + [~, upper_limit] = binary_search(xt, lims(2), lower_limit, nx); + + % Make the windows mapping to each pixel. + x_divisions = linspace(x(lower_limit, k), ... + x(upper_limit, k), ... + width + 1); + + end + + % Create a place to store the indices we'll need. + indices = [lower_limit, zeros(1, n_points-2), upper_limit]; + + % For each pixel... + right = lower_limit; + for z = 1:width-1 + + % Find the window bounds. + left = right; + [~, right] = binary_search(xt, ... + x_divisions(z+1), ... + left, upper_limit); + + % Get the indices of the max and min. + yt = y(left:right, k); + [~, max_index] = max(yt); + [~, min_index] = min(yt); + + % Record those indices. + indices(2*z:2*z+1) = sort([min_index max_index]) + left - 1; + + end + + % Sample the original x and y at the indices we found. + x_reduced(:, k) = xt(indices); + y_reduced(:, k) = y(indices, k); + + end + +end + +% Binary search to find boundaries of the ordered x data. +function [L, U] = binary_search(x, v, L, U) + while L < U - 1 % While there's space between them... + C = floor((L+U)/2); % Find the midpoint + if x(C) < v % Move the lower or upper bound in. + L = C; + else + U = C; + end + end +end diff --git a/+nigeLab/+utils/ToString.m b/+nigeLab/+utils/ToString.m new file mode 100644 index 00000000..ec1b4d3d --- /dev/null +++ b/+nigeLab/+utils/ToString.m @@ -0,0 +1,28 @@ +function str = ToString(in1) +%TOSTRING Universal to string convertor + +if iscell(in1) + if isscalar(in1) + str = {nigeLab.utils.ToString(in1{1})}; + elseif isempty(in1) + str = ''; + else + str = [{nigeLab.utils.ToString(in1{1})},... + nigeLab.utils.ToString(in1(2:end))]; + end + return; +end + +if isnumeric(in1) || islogical(in1) + str = num2str(in1); + +elseif isa(in1,'function_handle') + str = func2str(in1); + +elseif ischar(in1) + str = in1; + +end + +end + diff --git a/+nigeLab/+utils/buildWorkerConfigScript.m b/+nigeLab/+utils/buildWorkerConfigScript.m index b22d0a62..bcb0f911 100644 --- a/+nigeLab/+utils/buildWorkerConfigScript.m +++ b/+nigeLab/+utils/buildWorkerConfigScript.m @@ -36,7 +36,7 @@ case {'remote','fromremote'} % Just makes sure nigeLab is added to path (didn't work for MM) % This works if the method is queued FROM repo on remote location: - p = varargin{1}; + p = getNigelPath('UNC'); if ~iscell(p) p = {p}; end @@ -45,7 +45,7 @@ case {'local','fromlocal'} % "Wrap" everything in a script that is executed by the worker, % instead of the `doAction` - p = getNigelPath('UNC'); + p = varargin{1}; if ~iscell(p) p = {p}; end diff --git a/+nigeLab/@Block/private/SNEOThreshold.m b/+nigeLab/+utils/fastsmooth.m similarity index 70% rename from +nigeLab/@Block/private/SNEOThreshold.m rename to +nigeLab/+utils/fastsmooth.m index dc6c5d34..449b2e41 100644 --- a/+nigeLab/@Block/private/SNEOThreshold.m +++ b/+nigeLab/+utils/fastsmooth.m @@ -1,118 +1,3 @@ -function [p2pamp,ts,pmin,dt,E] = SNEOThreshold(data,pars,art_idx) -%% SNEOTHRESHOLD Smoothed nonlinear energy operator thresholding detect -% -% [p2pamp,ts,pmin,dt,E] = SNEOTHRESHOLD(data,pars,art_idx) -% -% -------- -% INPUTS -% -------- -% data : 1 x N double of bandpass filtered data, preferrably -% with artifact excluded already, on which to perform -% monopolar spike detection. -% -% pars : Parameters structure from SPIKEDETECTCLUSTER with -% the following fields: -% -% -> SNEO_N \\ number of samples for smoothing window -% -> MULTCOEFF \\ factor to multiply NEO noise threshold by -% -% art_idx : Indexing vector for artifact rejection periods, -% which are temporarily removed so that thresholds -% are not underestimated. -% -% -------- -% OUTPUT -% -------- -% p2pamp : Peak-to-peak amplitude of spikes. -% -% ts : Timestamps (sample indices) of spike peaks. -% -% pmin : Value at peak minimum. (pw) in SPIKEDETECTIONARRAY -% -% dt : Time difference between spikes. (pp) in -% SPIKEDETECTIONARRAY -% -% E : Smoothed nonlinear energy operator value at peaks. -% -% By: Max Murphy 1.0 01/04/2018 Original version (R2017a) - -%% GET NONLINEAR ENERGY OPERATOR SIGNAL AND SMOOTH IT -Y = data - mean(data); -Yb = Y(1:(end-2)); -Yf = Y(3:end); -Z = [0, Y(2:(end-1)).^2 - Yb .* Yf, 0]; % Discrete nonlinear energy operator -% Zs = fastsmooth(Z,pars.SNEO_N); -kern = ones(1,pars.SNEO_N)./pars.SNEO_N; -Zs = fliplr( conv( fliplr(conv(Z,kern,'same')) ,kern,'same')); % the same as the above tri option,(default one here used) but 10x faster -clear('Z','Y','Yb','Yf'); -%% CREATE THRESHOLD FILTER -tmpdata = data; -tmpdata(art_idx) = []; -tmpZ = Zs; -tmpZ(art_idx) = []; - -th = pars.MULTCOEFF * median(abs(tmpZ)); -data_th = pars.MULTCOEFF * median(abs(tmpdata)); -clear('tmpZ','tmpdata'); -%% PERFORM THRESHOLDING -pk = Zs > th; - -if sum(pk) <= 1 - p2pamp = []; - ts = []; - pmin = []; - dt = []; - E = []; - return -end - -%% REDUCE CONSECUTIVE CROSSINGS TO SINGLE POINTS -z = zeros(size(data)); -% pkloc = repmat(find(pk),pars.NS_AROUND*2+1,1) + (-pars.NS_AROUND:pars.NS_AROUND).'; -% pkloc(pkloc < 1) = 1; -% pkloc(pkloc > numel(data)) = numel(data); -% pkloc = unique(pkloc(:)); -% z(pkloc) = data(pkloc); - -%%%%%%%%%%%%%%%% FB, 5/28/2019 optimized for speed. The above process took 2.9s -%%%%%%%%%%%%%%%% now it takes more or less 0.85s. Same result -pkloc = conv(pk,ones(1,pars.NS_AROUND*2+1),'same')>0; -z(pkloc) = data(pkloc); - - -minTime = 1e-3*pars.REFRTIME; % parameter in milliseconds -[ts,pmin] = nigeLab.libs.peakseek(-z,minTime*pars.FS,data_th); -E = Zs(ts); - - -%% GET PEAK-TO-PEAK VALUES -tloc = repmat(ts,2*pars.PLP+1,1) + (-pars.PLP:pars.PLP).'; -tloc(tloc < 1) = 1; -tloc(tloc > numel(data)) = numel(data); -pmax = max(data(tloc)); - -p2pamp = pmax + pmin; - -%% EXCLUDE VALUES OF PMAX <= 0 -pm_ex = pmax<=0; -ts(pm_ex) = []; -p2pamp(pm_ex) = []; -pmax(pm_ex) = []; -pmin(pm_ex) = []; -E(pm_ex) = []; - -%% GET TIME DIFFERENCES -if numel(ts)>1 - dt = [diff(ts), round(median(diff(ts)))]; -elseif numel(ts)==1 - dt = 0; -else - dt = []; -end - - -end - function Y=fastsmooth(X,N,varargin) %% FASTSMOOTH Smooths vector X % @@ -352,4 +237,3 @@ end - diff --git a/+nigeLab/+utils/ginput.m b/+nigeLab/+utils/ginput.m index 5c4faa10..9c56027e 100644 --- a/+nigeLab/+utils/ginput.m +++ b/+nigeLab/+utils/ginput.m @@ -297,7 +297,7 @@ function updateCrossHair(fig, crossHair) % Create thin uicontrols with black backgrounds to simulate fullcrosshair pointer. % 1: horizontal left, 2: horizontal right, 3: vertical bottom, 4: vertical top for k = 1:4 - crossHair(k) = uicontrol(fig, 'Style', 'text', 'Visible', 'off', 'Units', 'pixels', 'BackgroundColor', [0 0 0], 'HandleVisibility', 'off', 'HitTest', 'off'); %#ok + crossHair(k) = uicontrol(fig, 'Style', 'text', 'Visible', 'off', 'Units', 'pixels', 'BackgroundColor', nigeLab.defaults.nigelColors('onsecondary'), 'HandleVisibility', 'off', 'HitTest', 'off'); %#ok end end diff --git a/+nigeLab/+utils/multiCallbackWrap.m b/+nigeLab/+utils/multiCallbackWrap.m index 79eae325..33383fe1 100644 --- a/+nigeLab/+utils/multiCallbackWrap.m +++ b/+nigeLab/+utils/multiCallbackWrap.m @@ -3,7 +3,7 @@ function multiCallbackWrap(ObjH, EventData,fcnList) % usage : % fcnList = {@myfunc0, @myfunc1, @myfunc2}; % obj = uicontrol(... -% 'Callback', {@nigeLab.Utils.multiCallbackWrap, fcnList}); +% 'Callback', ); % % or if input arguments are needed % @@ -14,23 +14,20 @@ function multiCallbackWrap(ObjH, EventData,fcnList) % obj = uicontrol(... % 'Callback', {@nigeLab.Utils.multiCallbackWrap, fcnList}); -if ~isscalar(fcnList) - for iFcn = 1:numel(fcnList) - nigeLab.utils.multiCallbackWrap(ObjH,EventData,fcnList{iFcn}); - end - return; +if ~iscell(fcnList) + error('nigeLab:BadInput','Third input argument must be a cell array!') +elseif isempty(fcnList) + return; elseif iscell(fcnList{1}) - nigeLab.utils.multiCallbackWrap(ObjH,EventData,fcnList{1}); + nigeLab.utils.multiCallbackWrap(ObjH,EventData,fcnList{1}); + nigeLab.utils.multiCallbackWrap(ObjH,EventData,fcnList(2:end)); + return; end -if iscell(fcnList) - thisFunction = fcnList{1}; - argsIn = fcnList; - argsIn(1) = []; -else - thisFunction = fcnList; - argsIn = {}; -end + +thisFunction = fcnList{1}; +argsIn = fcnList; +argsIn(1) = []; % If function handle is empty, continue if isempty(thisFunction) diff --git a/+nigeLab/@Block/Block.m b/+nigeLab/@Block/Block.m index ea1a4227..e74fc471 100644 --- a/+nigeLab/@Block/Block.m +++ b/+nigeLab/@Block/Block.m @@ -1021,7 +1021,7 @@ function updateVideosFolder(blockObj,newFolderPath) rms_out = analyzeRMS(blockObj,type,sampleIndices) % Compute RMS for channels % Methods for visualizing data: - flag = plotWaves(blockObj) % Plot stream snippets + flag = plotWaves(blockObj,ax,field,idx,computeRMS) % Plot stream snippets flag = plotSpikes(blockObj,ch) % Show spike clusters for a single channel flag = plotOverlay(blockObj) % Plot overlay of values on skull @@ -1214,6 +1214,11 @@ function unlockData(blockObj,fieldType) % SEALED,HIDDEN,PUBLIC methods (Sealed,Hidden,Access=public) varargout = testbench(blockObj,varargin); % Testbench with access to protected methods + + function SDargsout = testSD(blockobj,SDFun,data,SDPars) + SDargsout = cell(1,nargout(SDFun)); + [SDargsout{:}] = feval(SDFun,data,SDPars); + end end % % % % % % % % % % END METHODS% % % end \ No newline at end of file diff --git a/+nigeLab/@Block/analyzeRMS.m b/+nigeLab/@Block/analyzeRMS.m index b3ddd436..a8864cd3 100644 --- a/+nigeLab/@Block/analyzeRMS.m +++ b/+nigeLab/@Block/analyzeRMS.m @@ -30,9 +30,13 @@ for iT = 1:numel(type) data = blockObj.Channels(iCh).(type{iT})(:); x.(type{iT}) = rms(data); - rms_out(iCh).(type{iT}) = rms(data(sampleIndices)); + rms_out(iCh).(type{iT}) = rms(data); + end + if isempty(blockObj.RMS) + blockObj.RMS = [struct2table(x)]; + else + blockObj.RMS = [blockObj.RMS; struct2table(x)]; end - blockObj.RMS = [blockObj.RMS; struct2table(x)]; pct = 100 * (iCh / blockObj.NumChannels); fprintf(1,'\b\b\b\b\b%.3d%%\n',floor(pct)) end diff --git a/+nigeLab/@Block/doReReference.m b/+nigeLab/@Block/doReReference.m index 4663ae88..9036dd1e 100644 --- a/+nigeLab/@Block/doReReference.m +++ b/+nigeLab/@Block/doReReference.m @@ -57,8 +57,8 @@ curCh = curCh + 1; if ~doSuppression % Filter and and save amplifier_data by probe/channel - iProbe = blockObj.Channels(iCh).probe; - nChanPb = sum(iProbe == [blockObj.Channels.probe]); + iProbe = probe==blockObj.Channels(iCh).probe; + nChanPb = sum(probe(iProbe) == [blockObj.Channels.probe]); data = blockObj.Channels(iCh).Filt(:); refMean(iProbe,:)=refMean(iProbe,:) + data ./ nChanPb; else @@ -93,10 +93,11 @@ end curCh = 0; for iCh = blockObj.Mask + iProbe = probe==blockObj.Channels(iCh).probe; curCh = curCh + 1; % Do re-reference data = doCAR(blockObj.Channels(iCh),... - refMean(blockObj.Channels(iCh).probe)); + refMean(iProbe,:)); % Get filename pNum = num2str(blockObj.Channels(iCh).probe); @@ -108,7 +109,7 @@ blockObj.Channels(iCh).CAR = nigeLab.libs.DiskData(... 'MatFile',fName,data,'access','w','overwrite',true); lockData(blockObj.Channels(iCh).CAR); - blockObj.Channels(iCh).refMean = refMeanFile{blockObj.Channels(iCh).probe}; + blockObj.Channels(iCh).refMean = refMeanFile{iProbe}; % Update user diff --git a/+nigeLab/@Block/doSD.m b/+nigeLab/@Block/doSD.m index d9237989..a0465f60 100644 --- a/+nigeLab/@Block/doSD.m +++ b/+nigeLab/@Block/doSD.m @@ -32,8 +32,8 @@ return; end -[~,pars] = blockObj.updateParams('SD'); -pars.FS = blockObj.SampleRate; +[~,pars] = blockObj.updateParams('SD','KeepPars'); +pars.fs = blockObj.SampleRate; % UPDATE STATUS FOR THESE STAGES blockObj.updateStatus('Spikes',false,blockObj.Mask); @@ -51,9 +51,6 @@ curCh = 0; for iCh = blockObj.Mask curCh = curCh + 1; - % Parse file-name information - pNum = num2str(blockObj.Channels(iCh).probe); - chNum = blockObj.Channels(iCh).chStr; % No longer need check for CAR since checkActionIsValid does this data = blockObj.Channels(iCh).CAR(:); @@ -140,21 +137,22 @@ % Adapted by: MAECI 2018 Collaboration (Federico Barban & Max Murphy) % CONVERT PARAMETERS - pars.w_pre = double(round(pars.W_PRE / 1000 * pars.FS)); % Samples before spike - pars.w_post = double(round(pars.W_POST / 1000 * pars.FS)); % Samples after spike - pars.ls = double(pars.w_pre+pars.w_post); % Length of spike - pars.art_dist = double(pars.ART_DIST*pars.FS); % Maximum Stimulation frequency - pars.PLP = double(floor(pars.PKDURATION*1e-3*pars.FS)); % Pulse lifetime period [samples] - pars.RP = double(floor(pars.REFRTIME*1e-3*pars.FS)); % Refractory period [samples] - pars.nc_artifact = double(floor(pars.ARTIFACT_SPACE*1e-3*pars.FS)); % PLP [samples] - pars.npoints = double(length(data)); % Sample length of record - if pars.PRESCALED - pars.th_artifact = pars.ARTIFACT_THRESH; - else - pars.th_artifact = pars.ARTIFACT_THRESH * 1e-6; - end +% pars.w_pre = double(round(pars.W_PRE / 1000 * pars.FS)); % Samples before spike +% pars.w_post = double(round(pars.W_POST / 1000 * pars.FS)); % Samples after spike +% pars.ls = double(pars.w_pre+pars.w_post); % Length of spike +% pars.art_dist = double(pars.ART_DIST*pars.FS); % Maximum Stimulation frequency +% pars.PLP = double(floor(pars.PKDURATION*1e-3*pars.FS)); % Pulse lifetime period [samples] +% pars.RP = double(floor(pars.REFRTIME*1e-3*pars.FS)); % Refractory period [samples] +% pars.npoints = double(length(data)); % Sample length of record +% if pars.PRESCALED +% pars.th_artifact = pars.ARTIFACT_THRESH; +% else +% pars.th_artifact = pars.ARTIFACT_THRESH * 1e-6; +% end - % REMOVE ARTIFACT + + +% REMOVE ARTIFACT if ~isempty(pars.STIM_TS) data_ART = RemoveStimPeriods(data,pars); else @@ -166,174 +164,189 @@ else art_idx = []; end - [data_ART,artifact] = HardArtifactRejection(data_ART,pars); - % COMPUTE SPIKE THRESHOLD AND DO DETECTION - % SpikeDetection_PTSD_core.cpp; - if mod(pars.PLP,2)>0 - pars.PLP = pars.PLP + 1; % PLP must be even or doesn't work... - end - data_ART = double(data_ART); + ArtFun = ['ART_' pars.ArtefactRejMethodName]; + ArtPars = pars.(ArtFun); + ArtPars.fs = pars.fs; + Artargsout = cell(1,nargout(ArtFun)); + [Artargsout{:}] = feval(ArtFun,data_ART,ArtPars); + data_ART = Artargsout{1}; + artifact = Artargsout{2}; - switch pars.PKDETECT - case 'both' % (old, probably not using any more -MM 8/3/2017) - tmpdata = data_ART; - tmpdata(art_idx) = []; - pars.thresh = PreciseTimingThreshold(tmpdata,pars); - [spkValues, spkTimeStamps] = SpikeDetection_PTSD_core(data_ART, ... - pars.thresh, ... - pars.PLP, ... - pars.RP, ... - pars.ALIGNFLAG); - - % +1 added to accomodate for zero- (c) or one-based (matlab) array indexing - ts = 1 + spkTimeStamps( spkTimeStamps > 0); - p2pamp = spkValues( spkTimeStamps > 0); - pw = nan(size(p2pamp)); - pp = nan(size(p2pamp)); - - clear spkValues spkTimeStamps; - - case 'neg' % (probably use this in future -MM 8/3/2017) - - pars.thresh = pars.FIXED_THRESH; - - [p2pamp,ts,pw,pp] = ThresholdDetection(data_ART,pars,-1); - - case 'pos' - - pars.thresh = pars.FIXED_THRESH; - - [p2pamp,ts,pw,pp] = ThresholdDetection(data_ART,pars,1); - - case 'adapt' % Use findpeaks in conjunction w/ adaptive thresh -12/13/17 - pars.thresh = pars.MULTCOEFF; - [p2pamp,ts,pw,pp] = AdaptiveThreshold(data_ART,pars); - - case 'sneo' % Use findpeaks in conjunction w/ SNEO - 1/4/17 - pars.thresh = pars.MULTCOEFF; - [p2pamp,ts,pw,pp,E] = SNEOThreshold(data_ART,pars,art_idx); - otherwise - error('Invalid PKDETECT specification.'); + + + + SDFun = ['SD_' pars.SDMethodName]; + SDPars = pars.(SDFun); + SDPars.fs = pars.fs; + SDargsout = cell(1,nargout(SDFun)); + [SDargsout{:}] = feval(SDFun,data_ART,SDPars); + + tIdx = SDargsout{1}; + peak2peak = SDargsout{2}; + peakAmpl = SDargsout{3}; + peakWidth = SDargsout{4}; + + if numel(SDargsout)>4 + pTransformed = SDargsout{5}; end + dt = [diff(tIdx), round(median(diff(tIdx)))]; + + % COMPUTE SPIKE THRESHOLD AND DO DETECTION + % SpikeDetection_PTSD_core.cpp; +% if mod(pars.PLP,2)>0 +% pars.PLP = pars.PLP + 1; % PLP must be even or doesn't work... +% end +% data_ART = double(data_ART); +% +% switch pars.PKDETECT +% case 'both' % (old, probably not using any more -MM 8/3/2017) +% tmpdata = data_ART; +% tmpdata(art_idx) = []; +% pars.thresh = PreciseTimingThreshold(tmpdata,pars); +% [spkValues, spkTimeStamps] = SpikeDetection_PTSD_core(data_ART, ... +% pars.thresh, ... +% pars.PLP, ... +% pars.RP, ... +% pars.ALIGNFLAG); +% +% % +1 added to accomodate for zero- (c) or one-based (matlab) array indexing +% ts = 1 + spkTimeStamps( spkTimeStamps > 0); +% p2pamp = spkValues( spkTimeStamps > 0); +% pw = nan(size(p2pamp)); +% pp = nan(size(p2pamp)); +% +% clear spkValues spkTimeStamps; +% +% case 'neg' % (probably use this in future -MM 8/3/2017) +% +% pars.thresh = pars.FIXED_THRESH; +% +% [p2pamp,ts,pw,pp] = ThresholdDetection(data_ART,pars,-1); +% +% case 'pos' +% +% pars.thresh = pars.FIXED_THRESH; +% +% [p2pamp,ts,pw,pp] = ThresholdDetection(data_ART,pars,1); +% +% case 'adapt' % Use findpeaks in conjunction w/ adaptive thresh -12/13/17 +% pars.thresh = pars.MULTCOEFF; +% [p2pamp,ts,pw,pp] = AdaptiveThreshold(data_ART,pars); +% +% case 'SNEO' % Use findpeaks in conjunction w/ SNEO - 1/4/17 +% +% % [p2pamp,ts,pw,pp,E] = SNEOThreshold(data_ART,pars,art_idx); +% otherwise +% error('Invalid PKDETECT specification.'); +% end % ENSURE NO SPIKES REMAIN FROM ARTIFACT PERIODS if any(artifact) - [ts,ia]=setdiff(ts,(artifact./pars.FS)); - p2pamp=p2pamp(ia); - pw = pw(ia); - pp = pp(ia); - if exist('E','var')~=0 - E = E(ia); + [tIdx,ia]=setdiff(tIdx,(artifact./pars.fs)); + peak2peak=peak2peak(ia); + peakAmpl = peakAmpl(ia); + peakWidth = peakWidth(ia); + if exist('pTransformed','var')~=0 + pTransformed = pTransformed(ia); end end % EXCLUDE SPIKES THAT WOULD GO OUTSIDE THE RECORD - out_of_record = ts <= pars.w_pre+1 | ts >= pars.npoints-pars.w_post-2; - p2pamp(out_of_record) = []; - pw(out_of_record) = []; - pp(out_of_record) = []; - ts(out_of_record) = []; - if exist('E','var')~=0 - E(out_of_record) = []; + WindowPreSamples = pars.WPre * 1e-3 * pars.fs; + WindowPostSamples = pars.WPost * 1e-3 * pars.fs; + out_of_record = tIdx <= WindowPreSamples+1 | tIdx >= length(data) - WindowPostSamples - 2; + peak2peak(out_of_record) = []; + peakAmpl(out_of_record) = []; + peakWidth(out_of_record) = []; + tIdx(out_of_record) = []; + if exist('pTransformed','var')~=0 + pTransformed(out_of_record) = []; end % BUILD SPIKE SNIPPET ARRAY AND PEAK_TRAIN - if (any(ts)) % If there are spikes in the current signal + tIdx = tIdx(:); % make sure it's vertical + if (any(tIdx)) % If there are spikes in the current signal + snippetIdx = (-WindowPreSamples : WindowPostSamples) + tIdx; + spikes = data(snippetIdx); +% [peak_train,spikes] = BuildSpikeArray(data,ts,peak2peak,pars); - [peak_train,spikes] = BuildSpikeArray(data,ts,p2pamp,pars); +% %No interpolation in this case +% if length(spikes) > 1 +% %eliminates borders that were introduced for interpolation +% spikes(:,end-1:end)=[]; +% spikes(:,1:2)=[]; +% end + + %% Extract spike features + features = []; + pars.FEAT_NAMES = {}; + if ~strcmp(pars.FeatureExtractionMethodName,{'none',''}) + + %Interpolate spikes + pars.InterpSamples = max(size(spikes,2),pars.InterpSamples); + Ispikes = interp1(1:size(spikes,2),... + spikes.',... + linspace(1,size(spikes,2),pars.InterpSamples),... + 'spline').'; + + FeatFun = ['FEAT_' pars.FeatureExtractionMethodName]; + FeatPars = pars.(FeatFun); + FeatPars.fs = pars.fs; + Featargsout = cell(1,nargout(FeatFun)); + [Featargsout{:}] = feval(FeatFun,Ispikes,FeatPars); + features = Featargsout{1}; + if nargout(FeatFun) == 1 + % if feature labels are not provided just use the method + % name and numbers + pars.FEAT_NAMES = arrayfun(@(ii)... + sprintf('%s-%.2d',... + pars.FeatureExtractionMethodName, ii),... + 1:size(features,2),... + 'UniformOutput',false); + else + pars.FEAT_NAMES = Featargsout{2}; + end + end + + + normalize = @(x)(x - mean(x))./std(x,[],1); + features = [features, normalize(peakAmpl(:))]; + if ~ismember({'pk-amp'},pars.FEAT_NAMES) + pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-amp'}]; + end - %No interpolation in this case - if length(spikes) > 1 - %eliminates borders that were introduced for interpolation - spikes(:,end-1:end)=[]; - spikes(:,1:2)=[]; + features = [features, normalize(peak2peak(:))]; + if ~ismember({'pk-prom'},pars.FEAT_NAMES) + pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-prom'}]; end - % Extract spike features - if size(spikes,1) > pars.MIN_SPK % Need minimum number of spikes - features = WaveFeatures(spikes,pars); - features = features./std(features); - if ~any(isnan(p2pamp)) - tmp = (reshape(p2pamp,size(features,1),1)./max(p2pamp)-0.5)*3.0; - features = [features, tmp]; - if ~ismember({'pk-amp'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-amp'}]; - end - else - tmp = rand(size(features,1),1); - features = [features, tmp]; - end - if ~any(isnan(pp)) - tmp = (reshape(pp,size(features,1),1)./max(pp)-0.5)*3.0; - features = [features, tmp]; - if ~ismember({'pk-prom'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-prom'}]; - end - else - tmp = rand(size(features,1),1); - features = [features, tmp]; - end - if ~any(isnan(pw)) - tmp = (reshape(pw,size(features,1),1)./max(pw)-0.5)*3.0; - features = [features, tmp]; - if ~ismember({'pk-width'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-width'}]; - end - else - tmp = rand(size(features,1),1); - features = [features, tmp]; - end - if exist('E','var')~=0 - tmp = (reshape(E,size(features,1),1)./max(E)-0.5)*3.0; - features = [features, tmp]; - if ~ismember({'pk-energy'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-energy'}]; - end - else - tmp = rand(size(features,1),1); - features = [features, tmp]; - end - else - - if ~ismember({'pk-amp'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-amp'}]; - end - - - if ~ismember({'pk-prom'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-prom'}]; - end - - - if ~ismember({'pk-width'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-width'}]; - end - - - if ~ismember({'pk-energy'},pars.FEAT_NAMES) - pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-energy'}]; - end - - % Just make features reflect poor quality of (small) cluster - features = randn(size(spikes,1),numel(pars.FEAT_NAMES)) * 10; + features = [features, normalize(peakWidth(:))]; + if ~ismember({'pk-width'},pars.FEAT_NAMES) + pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-width'}]; end + + features = [features, normalize(pTransformed(:))]; + if ~ismember({'pk-energy'},pars.FEAT_NAMES) + pars.FEAT_NAMES = [pars.FEAT_NAMES, {'pk-energy'}]; + end + + else % If there are no spikes in the current signal - peak_train = sparse(double(pars.npoints) + double(pars.w_post), double(1)); - spk = ones(0,pars.ls + 4); + spk = ones(0,-WindowPreSamples + WindowPostSamples + 4); feat = ones(0,numel(pars.FEAT_NAMES)+4); art = ones(0,5); return; end % GENERATE OUTPUT VECTORS - tIdx = find(peak_train); - tIdx = reshape(tIdx,numel(tIdx),1); +% tIdx = find(peak_train); +% tIdx = reshape(tIdx,numel(tIdx),1); type = zeros(size(tIdx)); value = tIdx; % Store in value for now tag = zeros(size(tIdx)); - ts = tIdx./pars.FS; % Get values in seconds + ts = tIdx./pars.fs; % Get values in seconds % CONCATENATE OUTPUT MATRICES if isempty(spikes) @@ -352,17 +365,15 @@ if isempty(artifact) art = ones(0,5); else - artifact = artifact(:); % make sure is column oriented + artifact = artifact(:); % make sure is column oriented iStart = artifact([true, diff(artifact)' ~= 1]); iStop = artifact(fliplr([true, diff(fliplr(artifact))' ~= 1])); - artifact = reshape(artifact,numel(artifact),1); - value = reshape(iStart,numel(iStart),1); tag = reshape(iStop,numel(iStop),1); type = zeros(size(value)); - ts = value./pars.FS; - snippet = tag./pars.FS; + ts = value./pars.fs; + snippet = tag./pars.fs; art = [type,value,tag,ts,snippet]; end diff --git a/+nigeLab/@Block/getClus.m b/+nigeLab/@Block/getClus.m index 5c1f2032..ff2ad26f 100644 --- a/+nigeLab/@Block/getClus.m +++ b/+nigeLab/@Block/getClus.m @@ -51,14 +51,15 @@ end % CHECK TO BE SURE THAT THIS BLOCK/CHANNEL HAS BEEN SORTED +fName = fullfile(sprintf(strrep(blockObj.Paths.Clusters.file,'\','/'),... + num2str(blockObj.Channels(ch).probe),... + blockObj.Channels(ch).chStr)); + if getStatus(blockObj,'Clusters',ch) - clusterIndex = blockObj.Channels(ch).Clusters.value; + clusterIndex = getCIFromExistingFile(blockObj,ch); ts = getSpikeTimes(blockObj,ch); n = numel(ts); - if isempty(clusterIndex) % If the file does not exist - clusterIndex = zeros(n,1); % Assignment for output - return; - elseif numel(clusterIndex)~=n + if numel(clusterIndex)~=n warning(['[GETCLUS]::[%s] Mismatch between number of Cluster '... 'assignments (%g) and number of spikes (%g) for ' ... '%s::P%g-%s\n\t->\tAssigning all spikes to cluster zero.\n'],... @@ -67,35 +68,105 @@ clusterIndex = zeros(n,1); return; end +elseif exist(fName,'file')~=0 + load(fName,'data'); + clusterIndex = data(2,:); + ts = getSpikeTimes(blockObj,ch); + n = numel(ts); + if numel(clusterIndex)~=n + warning(['[GETCLUS]::[%s] Mismatch between number of Cluster '... + 'assignments (%g) and number of spikes (%g) for ' ... + '%s::P%g-%s\n\t->\tAssigning all spikes to cluster zero.\n'],... + blockObj.Name,numel(clusterIndex),n,... + blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); + clusterIndex = zeros(n,1); + return; + end + blockObj.Channels(ch).Clusters = nigeLab.libs.DiskData(fType,... + fName); + updateStatus(blockObj,'Clusters',true,ch); else % If it doesn't exist + nigeLab.utils.cprintf('*SystemCommands',true,'[%s Channel %d] *Cluster* data not found.Initializing empty Clusters file.\n',blockObj.Name,ch); if getStatus(blockObj,'Spikes',ch) % but spikes do - ts = getSpikeTimes(blockObj,ch); - n = numel(ts); - clusterIndex = zeros(n,1); % Assignment for output - data = [zeros(n,3) clusterIndex zeros(n,1)]; - - % initialize the 'Sorted' DiskData file - fType = blockObj.getFileType('Clusters'); - fName = fullfile(sprintf(strrep(blockObj.Paths.Clusters.file,'\','/'),... - num2str(blockObj.Channels(ch).probe),... - blockObj.Channels(ch).chStr)); - if exist(blockObj.Paths.Clusters.dir,'dir')==0 - mkdir(blockObj.Paths.Clusters.dir); - end - if exist(fName,'file')==0 - blockObj.Channels(ch).Clusters = nigeLab.libs.DiskData(fType,... - fName,data,'access','w'); - - if ~suppressText - fprintf(1,'Initialized Clusters file for P%d: Ch-%s\n',... - blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); - end - end + initClusters(blockObj,ch,suppressText); + clusterIndex = getCIFromExistingFile(blockObj,ch); +% updateStatus(blockObj,'Clusters',true,ch); return; else - clusterIndex = zeros(0,1); - return; + error(['[GETCLUS]::[%s Channel %d] Trying to access *Clusters* without *Spikes*.\n'... + 'To access spike sorting funcionalities run doSD first.'],... + blockObj.Name,ch); end end +if ~any(clusterIndex) % if all indexes are zero + nigeLab.utils.cprintf('*SystemCommands',true,'[%s Channel %d] *Cluster* data contains only zeros.\n',blockObj.Name,ch); +end +end + +function clusterIndex = getCIFromExistingFile(blockObj,ch) +%%GETCIFROMEXISTINGFILE Get cluster index from existing file +% For backwards compatibility, make sure "tags" is not a cell + +ftype = getFileType(blockObj,'Clusters'); +info = getInfo(blockObj.Channels(ch).Clusters); +names = {info.Name}; +tagIdx = find(strcmpi(names,'tag'),1,'first'); + +if isempty(tagIdx) % Everything is fine, return it the normal way + clusterIndex = blockObj.Channels(ch).Clusters.value; +elseif strcmp(info(tagIdx).class,'cell') + tag = blockObj.Channels(ch).Clusters.tag(:); + tag = parseSpikeTagIdx(blockObj,tag); + fname = getPath(blockObj.Channels(ch).Clusters); + + clusters = zeros(numel(tag),5); + clusters(:,2) = blockObj.Channels(ch).Clusters.class(:); + clusters(:,3) = tag; + clusters(:,4) = getSpikeTimes(blockObj,ch); + + blockObj.Channels(ch).Clusters = ... + nigeLab.libs.DiskData(ftype,fullfile(fname),... + clusters,'access','w'); + clusterIndex = clusters(:,2); + +elseif numel(info) > 1 + fname = getPath(blockObj.Channels(ch).Clusters); + + clusters = zeros(numel(info(1).size(1)),5); + clusters(:,2) = blockObj.Channels(ch).Clusters.class(:); + clusters(:,3) = blockObj.Channels(ch).Clusters.tag(:); + clusters(:,4) = getSpikeTimes(blockObj,ch); + + blockObj.Channels(ch).Clusters = ... + nigeLab.libs.DiskData(ftype,fullfile(fname),... + clusters,'access','w'); + clusterIndex = clusters(:,2); + +else % Everything is fine, return it the normal way + clusterIndex = blockObj.Channels(ch).Sorted.value; +end +end -end \ No newline at end of file +function initClusters(blockObj,ch,suppressText) + if getStatus(blockObj,'Spikes',ch) % if spikes are present + ts = getSpikeTimes(blockObj,ch); + n = numel(ts); + clusterIndex = zeros(n,1); % Assignment for output + data = [zeros(n,3) clusterIndex zeros(n,1)]; + + % initialize the 'Sorted' DiskData file + fType = blockObj.getFileType('Clusters'); + fName = fullfile(sprintf(strrep(blockObj.Paths.Clusters.file,'\','/'),... + num2str(blockObj.Channels(ch).probe),... + blockObj.Channels(ch).chStr)); + if exist(blockObj.Paths.Clusters.dir,'dir')==0 + mkdir(blockObj.Paths.Clusters.dir); + end + blockObj.Channels(ch).Clusters = nigeLab.libs.DiskData(fType,... + fName,data,'access','w'); + if ~suppressText + fprintf(1,'Initialized Clusters file for P%d: Ch-%s\n',... + blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); + end + end +end diff --git a/+nigeLab/@Block/getSort.m b/+nigeLab/@Block/getSort.m index 6f16c54d..16be8b8c 100644 --- a/+nigeLab/@Block/getSort.m +++ b/+nigeLab/@Block/getSort.m @@ -29,7 +29,7 @@ end if nargin < 3 - suppressText = false; + suppressText = ~blockObj.Verbose; end %% USE RECURSION TO ITERATE ON MULTIPLE CHANNELS @@ -62,8 +62,8 @@ clusterIndex = blockObj.getClus(ch); end elseif exist(fName,'file')~=0 % Technically, files could exist but Status not updated... - - clusterIndex = getCIFromExistingFile(blockObj,ch); + load(fName,'data'); + clusterIndex = data(:,2); if ~any(clusterIndex) % if all indexes are zero clusterIndex = blockObj.getClus(ch); ts = getSpikeTimes(blockObj,ch); @@ -73,53 +73,57 @@ mkdir(blockObj.Paths.Sorted.dir); end blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... - fName,data,'access','w'); + fName,data,'overwrite', true); if ~suppressText fprintf(1,'Initialized Sorted file for P%d: Ch-%s\n',... blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); end + else + blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... + fName); end updateStatus(blockObj,'Sorted',true,ch); -elseif getStatus(blockObj,'Clusters',ch) +else + nigeLab.utils.cprintf('*SystemCommands',true,'*Sorted* data were not found. Using *Clusters* instead.\n'); clusterIndex = blockObj.getClus(ch); - ts = getSpikeTimes(blockObj,ch); - n = numel(ts); - data = [zeros(n,2) clusterIndex ts zeros(n,1)]; - if exist(blockObj.Paths.Sorted.dir,'dir')==0 - mkdir(blockObj.Paths.Sorted.dir); - end - blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... - fName,data,'access','w'); - if ~suppressText - fprintf(1,'Initialized Sorted file for P%d: Ch-%s\n',... - blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); - - end - - updateStatus(blockObj,'Sorted',true,ch); -elseif getStatus(blockObj,'Spikes',ch) % but spikes do - % initialize the 'Sorted' DiskData file - - - - ts = getSpikeTimes(blockObj,ch); - n = numel(ts); - clusterIndex = zeros(n,1); - data = [zeros(n,2) clusterIndex ts zeros(n,1)]; - if exist(blockObj.Paths.Sorted.dir,'dir')==0 - mkdir(blockObj.Paths.Sorted.dir); - end - blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... - fName,data,'access','w'); - if ~suppressText - fprintf(1,'Initialized Sorted file for P%d: Ch-%s\n',... - blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); - + if ~isempty(clusterIndex) + ts = getSpikeTimes(blockObj,ch); + n = numel(ts); + data = [zeros(n,2) clusterIndex ts zeros(n,1)]; + if exist(blockObj.Paths.Sorted.dir,'dir')==0 + mkdir(blockObj.Paths.Sorted.dir); + end + blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... + fName,data,'access','w'); + if ~suppressText + fprintf(1,'Initialized Sorted file for P%d: Ch-%s\n',... + blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); + + end end -else - clusterIndex = []; +% elseif getStatus(blockObj,'Spikes',ch) % but spikes do +% % initialize the 'Sorted' DiskData file +% +% +% +% ts = getSpikeTimes(blockObj,ch); +% n = numel(ts); +% clusterIndex = zeros(n,1); +% data = [zeros(n,2) clusterIndex ts zeros(n,1)]; +% if exist(blockObj.Paths.Sorted.dir,'dir')==0 +% mkdir(blockObj.Paths.Sorted.dir); +% end +% blockObj.Channels(ch).Sorted = nigeLab.libs.DiskData(fType,... +% fName,data,'access','w'); +% if ~suppressText +% fprintf(1,'Initialized Sorted file for P%d: Ch-%s\n',... +% blockObj.Channels(ch).probe,blockObj.Channels(ch).chStr); +% +% end +% else +% clusterIndex = []; end end diff --git a/+nigeLab/@Block/getSpikeTimes.m b/+nigeLab/@Block/getSpikeTimes.m index 17528512..1c394771 100644 --- a/+nigeLab/@Block/getSpikeTimes.m +++ b/+nigeLab/@Block/getSpikeTimes.m @@ -79,7 +79,7 @@ elseif strcmpi(clusterIndex,'all') ts = getEventData(blockObj,'Spikes','ts',ch); else - ts = getEventData(blockObj,'Spikes','ts',ch,'tag',clusterIndex); + ts = getEventData(blockObj,'Sorted','ts',ch,'value',clusterIndex); end if isempty(ts) ts = zeros(0,1); diff --git a/+nigeLab/@Block/plotSpikes.m b/+nigeLab/@Block/plotSpikes.m index 5f893f5e..180d9667 100644 --- a/+nigeLab/@Block/plotSpikes.m +++ b/+nigeLab/@Block/plotSpikes.m @@ -34,14 +34,14 @@ %% PARSE INPUT ARGUMENTS flag = false; -if ~ismember('Spikes',blockObj.Fields(blockObj.Status)) +if ~all(blockObj.getStatus('Spikes')) warning('No spikes detected yet.'); return; end if nargin < 3 class = nan; -elseif ~ismember('Sorted',blockObj.Fields(blockObj.Status)) +elseif ~all(blockObj.getStatus('Sorted')) || ~all(blockObj.getStatus('Clusters')) class = nan; end diff --git a/+nigeLab/@Block/plotWaves.m b/+nigeLab/@Block/plotWaves.m index d4c83fa9..60c34ed3 100644 --- a/+nigeLab/@Block/plotWaves.m +++ b/+nigeLab/@Block/plotWaves.m @@ -1,4 +1,4 @@ -function flag = plotWaves(blockObj) +function flag = plotWaves(blockObj,ax,field,idx,computeRMS) %% PLOTWAVES Plot multi-channel waveform snippets for BLOCK % % flag = PLOTWAVES(blockObj); @@ -19,70 +19,142 @@ %% DEFAULTS flag = false; -blockObj.PlotPars = nigeLab.defaults.Plot(); - +% blockObj.Pars.Pars.Plot = nigeLab.defaults.Plot(); %% FIGURE OUT WHAT TO PLOT str_in = blockObj.getStatus; -[~,idx] = nigeLab.utils.uidropdownbox('Choose Wave Type',... - 'Select type of waveform to plot:',... - str_in); - -str = str_in{idx}; -if strcmp(str,'Spikes') - warning('Spikes overlay not yet supported.'); - return; -end -%% GET INDEXING VECTOR -if strcmp(str,'LFP') - fs = blockObj.LFPPars.DownSampledRate; +%% Parse first input for axes reference +if isa(ax,'matlab.graphics.axis.Axes') + fig = ax.Parent; + nargs = nargin-1; else - fs = blockObj.SampleRate; + field = ax; + nargs = nargin; + + fig = figure('Name',sprintf('Multi-Channel %s Snippets',field), ... + 'Units','Normalized', ... + 'Position',[0.05*rand+0.1,0.05*rand+0.1,0.8,0.8],... + 'Color','w','NumberTitle','off'); + + ax = axes(fig,'NextPlot','add'); +end + +%% Parse field input and change options accordingly +if nargs < 2 + [~,idx] = nigeLab.utils.uidropdownbox('Choose Wave Type',... + 'Select type of waveform to plot:',... + str_in); + field = str_in{idx}; end -tStart = (blockObj.PlotPars.DefTime - blockObj.PlotPars.PreAlign)/1000; -iStart = max(1,round(tStart * fs)); -tStop = (blockObj.PlotPars.DefTime + blockObj.PlotPars.PostAlign)/1000; -iStop = min(blockObj.Samples,round(tStop * fs)); +plotSpikesOverlay = false; +switch field + case 'LFP' + fs = blockObj.LFPPars.DownSampledRate; + case 'Spikes' + fs = blockObj.SampleRate; + plotSpikesOverlay = true; + if all(blockObj.getStatus({'CAR'})) + field = 'CAR'; + elseif all(blockObj.getStatus({'Filt'})) + field = 'Filt'; + else + error('CAR or Filt are needed to overlay spikes!'); + end + otherwise + fs = blockObj.SampleRate; +end -vec = iStart:iStop; -t = linspace(tStart,tStop,numel(vec)); +%% Get correct index +if nargs < 3 + tStart = (blockObj.Pars.Plot.DefTime - blockObj.Pars.Plot.PreAlign)/1000; + iStart = max(1,round(tStart * fs)); + + tStop = (blockObj.Pars.Plot.DefTime + blockObj.Pars.Plot.PostAlign)/1000; + iStop = min(blockObj.Samples,round(tStop * fs)); + + idx = iStart:iStop; +end +t = linspace(idx(1)./fs,idx(end)./fs,numel(idx)); dt = mode(diff(t)); +%% RMS input +if nargs < 4 + computeRMS = false; +end + + + %% ASSIGN CHANNEL COLORS BASED ON RMS -load(blockObj.PlotPars.ColorMapFile,'cm'); -if isempty(blockObj.RMS) - analyzeRMS(blockObj); -elseif ~ismember(str,blockObj.RMS.Properties.VariableNames) - analyzeRMS(blockObj); +load(blockObj.Pars.Plot.ColorMapFile,'cm'); +if isempty(blockObj.RMS) && computeRMS + analyzeRMS(blockObj); + r = blockObj.RMS.(field); + ic = assignColors(r); +elseif ~ismember(field,blockObj.RMS.Properties.VariableNames) && computeRMS + analyzeRMS(blockObj); + r = blockObj.RMS.(field); + ic = assignColors(r); +else + %Multicolored version +% cm = nigeLab.defaults.nigelColors(1:blockObj.NumChannels,'cubehelix'); +% ic = 1:blockObj.NumChannels; + + %Monochrome version + cm = nigeLab.defaults.nigelColors('primary'); + ic = ones(1,blockObj.NumChannels); + + r = nan(1,blockObj.NumChannels); end -r = blockObj.RMS.(str); -ic = assignColors(r); -%% MAKE FIGURE AND PLOT -fig = figure('Name',sprintf('Multi-Channel %s Snippets',str), ... - 'Units','Normalized', ... - 'Position',[0.05*rand+0.1,0.05*rand+0.1,0.8,0.8],... - 'Color','w'); -ax = axes(fig,'NextPlot','add'); +%% MAKE FIGURE AND PLOT tickLabs = cell(blockObj.NumChannels,1); -tickLocs = 0:blockObj.PlotPars.VertOffset:(blockObj.PlotPars.VertOffset*(... - blockObj.NumChannels-1)); +tickLocs = 1:blockObj.Pars.Plot.VertOffset:(blockObj.Pars.Plot.VertOffset*(... + blockObj.NumChannels)); +pixelPos = getpixelposition(ax); for iCh = 1:blockObj.NumChannels tickLabs{iCh} = blockObj.Channels(iCh).custom_channel_name; - y = blockObj.Channels(iCh).(str)(vec)+tickLocs(iCh); - plot(ax,t,y, ... + y = blockObj.Channels(iCh).(field)(idx)+tickLocs(iCh); + [t_reduced, y_reduced] = nigeLab.utils.reduce_to_width(t, y, pixelPos(end) ,[t(1) t(end)]); + line(ax,t_reduced,y_reduced, ... 'Color',cm(ic(iCh),:), ... 'LineWidth',1.75,... 'UserData',iCh); %#ok + + if plotSpikesOverlay + fs = blockObj.SampleRate; + Wpre = blockObj.Pars.SD.WPre * 1e-3; + Wpost = blockObj.Pars.SD.WPost * 1e-3; + + tSpk = blockObj.getSpikeTimes(iCh); + spkToKeep = tSpk < t(end) & tSpk > t(1); + tSpk = tSpk(spkToKeep,:); + tSpk = (-Wpre : 1/fs : Wpost) + tSpk; + spkIndex = logical(sum( ... + tSpk(:,1) <= t_reduced & tSpk(:,end) >= t_reduced, 1)); + + spkBound = conv(spkIndex,[-1 1],'same')>0; + spkBound(end) = false; + + tSpk = t_reduced; + tSpk(~spkIndex) = nan; + + spk = y_reduced; + spk(~spkIndex) = nan; + + line(ax,tSpk,spk, ... + 'Color','k', ... + 'LineWidth',1.75,... + 'UserData',iCh); %#ok + end - text(ax,max(t)+dt,tickLocs(iCh),... + text(ax,max(t)*1.01,tickLocs(iCh),... sprintf('RMS: %.3g',r(iCh)),... 'FontName','Arial',... 'FontWeight','bold',... 'Color','k',... - 'FontSize',14,... + 'FontSize',10,... 'UserData',iCh); end @@ -90,15 +162,17 @@ ax.YTickLabel = tickLabs; xlabel('Time (sec)','FontName','Arial','FontSize',14,'Color','k'); ylabel('Channel','FontName','Arial','FontSize',14,'Color','k'); -title(sprintf('Multi-Channel %s Snippets',str),'FontName','Arial',... +title(sprintf('Multi-Channel %s Snippets',field),'FontName','Arial',... 'Color','k','FontSize',20); -fname = fullfile(blockObj.paths.TW,... - [blockObj.Name sprintf(blockObj.PlotPars.SnippetString,str)]); - -savefig(fig,[fname '.fig']); -saveas(fig,[fname '.png']); -blockObj.Graphics.Waves = fig; +xlim(ax,t([1,end])) +ylim(ax,tickLocs([1,end])+ diff(tickLocs([1,end]))*[-.05 .05]) +% fname = fullfile(blockObj.paths.TW,... +% [blockObj.Name sprintf(blockObj.Pars.Plot.SnippetString,str)]); +% +% savefig(fig,[fname '.fig']); +% saveas(fig,[fname '.png']); +% blockObj.Graphics.Waves = fig; flag = true; diff --git a/+nigeLab/@Block/private/HardArtifactRejection.m b/+nigeLab/@Block/private/ART_HardThresh.m similarity index 78% rename from +nigeLab/@Block/private/HardArtifactRejection.m rename to +nigeLab/@Block/private/ART_HardThresh.m index 8cf628c1..feb34cb7 100644 --- a/+nigeLab/@Block/private/HardArtifactRejection.m +++ b/+nigeLab/@Block/private/ART_HardThresh.m @@ -1,4 +1,4 @@ -function [data_ART,art_idx] = HardArtifactRejection(data,pars) +function [data_ART,art_idx] = ART_HardThresh(data,pars) %% HARDARTIFACTREJECTION Automatically changes data to zero on samples and samples within a pre-specified window around it upon crossing a pre-specified threshold. % % -------- @@ -10,9 +10,9 @@ % pars : Parameter structure from SPIKEDETECTCLUSTER with % the following fields: % -% -> th_artifact \\ Artifact threshold (micro-volts). +% -> Thresh \\ Artifact threshold (micro-volts). % -% -> nc_artifact \\ Number of samples around threshold to replace +% -> Samples \\ Number of samples around threshold to replace % with zeros to cancel out any ripple effects of % amplifier saturation and rebound. % @@ -37,14 +37,15 @@ % Kelly RM v1.0 11/04/2015 Original version. %% FIND THRESHOLD CROSSINGS -segm = find(abs(data) >= pars.th_artifact); +segm = find(pars.Polarity*data >= pars.Thresh); art_idx = []; +Nsamples = double(floor(pars.Samples*1e-3*pars.fs)); % from ms to samples %% REMOVE (ZERO) SECTIONS AROUND ARTIFACT if any(segm) - if (pars.nc_artifact) - rm_pre = max(segm - pars.nc_artifact,1); % Start index >= 1 - rm_post = min(length(data),segm + pars.nc_artifact); % End index <= total samples in data + if (Nsamples) + rm_pre = max(segm - Nsamples,1); % Start index >= 1 + rm_post = min(length(data),segm + Nsamples); % End index <= total samples in data for in = 1:length(segm) % art_idx = [art_idx ,rm_pre(in):segm(in)]; @@ -59,7 +60,7 @@ else data_ART = data; - art_idx = 0; + art_idx = []; end end \ No newline at end of file diff --git a/+nigeLab/@Block/private/FEAT_ica.m b/+nigeLab/@Block/private/FEAT_ica.m new file mode 100644 index 00000000..3f6cd868 --- /dev/null +++ b/+nigeLab/@Block/private/FEAT_ica.m @@ -0,0 +1,58 @@ +function feat = FEAT_ica(spikes,pars) +%% WAVEFEATURES Calculates the spike features +% +% [inspk,K] = WAVEFEATURES(spikes,pars) +% +% -------- +% INPUTS +% -------- +% spikes : N x M matrix of putative spike waveforms. N is the +% number of spikes, M is the number of samples in +% each spike "snippet." +% +% pars : Parameters structure for feature decomposition. +% +% -------- +% OUTPUT +% -------- +% feat : N x K matrix of spike features to be input to the +% SPC algorithm. N is the number of spikes, K is the +% user-defined number of features (which depends on +% the feature-extraction algorithm selected). +% +% +% Reformatted for use in nigeLab by FB 2020/06/19 (what a great year btw) +% Based on previous design by Max Murphy + +%% GET VARIABLES +error('nigeLab:WIP','ICA based decomposition is not yet implemented.\nSorry, the cat ate my homework...'); + +% N =size(spikes,1); % # spikes +% +% if N <= 10 +% K = 3; +% feat = zeros(N,K); % Just put everything into same cluster +% return; +% end +% +% +% M =size(spikes,2); % # samples per spike +% +% %% CALCULATES FEATURES +% K = 3; +% Z = fastICA(spikes.',K); +% cc = Z.'; +% coeff = 1:K; +% +% +% +% %% CREATES INPUT MATRIX FOR SPC +% feat=zeros(N,K); +% for iN=1:N +% for iK=1:K +% feat(iN,iK)=cc(iN,coeff(iK)); +% end +% end + +end + diff --git a/+nigeLab/@Block/private/FEAT_pca.m b/+nigeLab/@Block/private/FEAT_pca.m new file mode 100644 index 00000000..d4364b14 --- /dev/null +++ b/+nigeLab/@Block/private/FEAT_pca.m @@ -0,0 +1,48 @@ +function feat = FEAT_pca(spikes,pars) +%% WAVEFEATURES Calculates the spike features +% +% [inspk,K] = WAVEFEATURES(spikes,pars) +% +% -------- +% INPUTS +% -------- +% spikes : N x M matrix of putative spike waveforms. N is the +% number of spikes, M is the number of samples in +% each spike "snippet." +% +% pars : Parameters structure from feature decomposition. +% +% -------- +% OUTPUT +% -------- +% feat : N x K matrix of spike features to be input to the +% SPC algorithm. N is the number of spikes, K is the +% user-defined number of features (which depends on +% the feature-extraction algorithm selected). +% +% +% Reformatted for use in nigeLab by FB 2020/06/19 (what a great year btw) +% Based on previous design by Max Murphy + + +%% GET VARIABLES +N =size(spikes,1); % # spikes + +if N <= 10 + K = pars.NOut; + feat = rand(N,K); % Just put everything into same cluster + return; +end + +[~, SCORE, LATENT] = pca(spikes); +if isinf(coeff.ExplVar) + K = pars.NOut; +else + K = find( cumsum(LATENT)./sum(LATENT) > pars.ExplVar,1); +end +%% CREATES INPUT MATRIX FOR SPC +feat = SCORE(:,1:K); + + +end + diff --git a/+nigeLab/@Block/private/FEAT_wavelet.m b/+nigeLab/@Block/private/FEAT_wavelet.m new file mode 100644 index 00000000..fae41363 --- /dev/null +++ b/+nigeLab/@Block/private/FEAT_wavelet.m @@ -0,0 +1,147 @@ +function [feat,featName] = FEAT_wavelet(spikes,pars) +%% WAVEFEATURES Calculates the spike features +% +% [inspk,K] = WAVEFEATURES(spikes,pars) +% +% -------- +% INPUTS +% -------- +% spikes : N x M matrix of putative spike waveforms. N is the +% number of spikes, M is the number of samples in +% each spike "snippet." +% +% pars : Parameters structure from feature decomposition. +% +% -------- +% OUTPUT +% -------- +% feat : N x K matrix of spike features to be input to the +% SPC algorithm. N is the number of spikes, K is the +% user-defined number of features (which depends on +% the feature-extraction algorithm selected). +% +% +% Reformatted for use in nigeLab by FB 2020/06/19 (what a great year btw) +% Based on previous design by Max Murphy + +%% GET VARIABLES +N =size(spikes,1); % # spikes + +if N <= 10 + + K = pars.NOut; + feat = rand(N,K); % Just put everything into same cluster + featName = arrayfun(@(ii) sprintf('rand-%.2d',ii),1:K,'UniformOutput',false); + return; +end + +M =size(spikes,2); % # samples per spike + +%% CALCULATES Wavelt decomposition features + K = pars.NOut; + cc = zeros(N,floor(M/2)); + for iN=1:N % Wavelet decomposition (been using 3 scales, 'bior1.3') + [C,L] = wavedec(spikes(iN,:), ... + pars.NScales, ... + pars.LoD,... + pars.HiD); + cc(iN,:) = C((sum(L(1:pars.NScales))+2):(end-1)); + end + + % Remove columns (features) that are mostly 0 + aux = []; + for iC = 1:size(cc,2) + if sum(abs(cc(:,iC)) + end + end + + % Normalize features + cc = cc - mean(cc,1); + cc = cc./std(cc,[],1); + + % Find kurtosis peaks in the time-series distribution + y = kurtosis(cc); + [k_pk,k_loc] = findpeaks(y); + try + [~,ind] = sort(k_pk,'descend'); + loc = k_loc(ind); + catch + loc = []; + end + + if (numel(loc) >= K) + coeff = loc(1:K); + else % Not enough coefficients, look for skewness + y = skewness(cc); + try + [p_pk, p_loc] = findpeaks(y); + [p_pk,ind] = sort(p_pk,'descend'); + p_loc = p_loc(ind); + catch + p_loc = []; + end + + try + [n_pk, n_loc] = findpeaks(-y); + [n_pk,ind] = sort(n_pk,'descend'); + n_loc = n_loc(ind); + catch + n_loc = []; + end + + % Decide which one to include first + if ~isempty(p_loc) && ~isempty(n_loc) + flag = p_loc(1) > n_loc(1); + elseif isempty(p_loc) && isempty(n_loc) + flag = true; + else + if isempty(p_loc) + flag = false; + elseif isempty(n_loc) + flag = true; + end + end + + % Add coeffs based on skew until enough coefficients + pk_count = 0; + while ((pk_count <= max(numel(p_pk),numel(n_pk))) && ... + (numel(loc) < K)) + pk_count = pk_count + 0.5; + + if (flag && (pk_count <= numel(p_pk))) + loc = [loc, p_loc(ceil(pk_count))]; %#ok + elseif (~flag && (pk_count <= numel(n_pk))) + loc = [loc, n_loc(ceil(pk_count))]; %#ok + end + loc = unique(loc); + end + + % If still need coefficients, add random ones + if numel(loc) < K + vec = 1:size(cc,2); + vec = setdiff(vec,loc); + n_remain = K - numel(loc); + + loc = [loc, vec(randperm(numel(vec),n_remain))]; + end + + coeff = loc(1:K); + end + + + +%% Format Output +feat=zeros(N,K); +for iN=1:N + for iK=1:K + feat(iN,iK)=cc(iN,coeff(iK)); + end +end + +featName = arrayfun(@(ii) sprintf('wave-%.2d',ii),1:K,'UniformOutput',false); + + + +end + diff --git a/+nigeLab/@Block/private/AdaptiveThreshold.m b/+nigeLab/@Block/private/SD_AdaptThresh.m similarity index 62% rename from +nigeLab/@Block/private/AdaptiveThreshold.m rename to +nigeLab/@Block/private/SD_AdaptThresh.m index 7179f4df..658ffcfb 100644 --- a/+nigeLab/@Block/private/AdaptiveThreshold.m +++ b/+nigeLab/@Block/private/SD_AdaptThresh.m @@ -1,4 +1,4 @@ -function [p2pamp,ts,pmin,dt] = AdaptiveThreshold(data,pars) +function [ts,p2pamp,pmin,pW] = SD_AdaptThresh(data,pars) %% ADAPTIVE_THRESHOLD Set adaptive threshold based on local noise % % [pwpamp,ts,pmin,dt] = ADAPTIVE_THRESHOLD(data,pars,art_idx) @@ -13,7 +13,7 @@ % pars : Parameters structure from SPIKEDETECTCLUSTER with % the following fields: % -% -> ADPT_N \\ Number of samples in local noise detection filter +% -> FilterLength \\ Number of samples in local noise detection filter % % -------- % OUTPUT @@ -29,30 +29,36 @@ % By: Max Murphy 1.0 12/13/2017 Original version (R2017b) %% CREATE THRESHOLD FILTER -n = round((pars.ADPT_N/1e3)*pars.FS); -th = fastsmooth(abs(data),pars.ADPT_N,'abs_med',1); -th(th (th * pars.MULTCOEFF); +pk = (pars.Polarity * data) > th; if sum(pk) <= 1 p2pamp = []; ts = []; pmin = []; dt = []; + pW = []; return end %% REDUCE CONSECUTIVE CROSSINGS TO SINGLE POINTS z = zeros(size(data)); -z(pk) = data(pk); -[pmin,ts] = findpeaks(-z); +z(pk) = pars.Polarity * data(pk); +[pmin,ts] = findpeaks(z); +pmin = pmin * pars.Polarity; %% GET PEAK-TO-PEAK VALUES -tloc = repmat(ts,pars.PLP,1) + (1:pars.PLP).'; -pmax = max(data(tloc)); - +PLP = pars.PeakDur*1e-3*pars.fs; % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); p2pamp = pmax + pmin; %% EXCLUDE VALUES OF PMAX <= 0 @@ -61,22 +67,7 @@ p2pamp(pm_ex) = []; pmax(pm_ex) = []; pmin(pm_ex) = []; - -%% GET TIME DIFFERENCES -dt = [diff(ts), round(median(diff(ts)))]; - -%% PERFORM SECONDARY EXCLUSIONS -a = p2pamp/(2*pars.ARTIFACT_THRESH); -b = pmin/pars.ARTIFACT_THRESH; -c = dt/pars.FS; - -dt_ex = c <= pars.ART_RATE & b >= (pars.M.*a + pars.B); - -ts(dt_ex) = []; -p2pamp(dt_ex) = []; -pmax(dt_ex) = []; -pmin(dt_ex) = []; -dt(dt_ex) = []; +pW(pm_ex) = []; diff --git a/+nigeLab/@Block/private/ThresholdDetection.m b/+nigeLab/@Block/private/SD_HardThresh.m similarity index 66% rename from +nigeLab/@Block/private/ThresholdDetection.m rename to +nigeLab/@Block/private/SD_HardThresh.m index bddccfa0..41eb5642 100644 --- a/+nigeLab/@Block/private/ThresholdDetection.m +++ b/+nigeLab/@Block/private/SD_HardThresh.m @@ -1,4 +1,4 @@ -function [v,ts,w,p] = ThresholdDetection(data,pars,polarity) +function [ts,p2pamp,pmin,pW] = SD_HardThresh(data,pars) %% THRESHOLDDETECTION Use monopolar threshold-crossing to get spike times % % [v,ts,w,p] = THRESHOLDDETECTION(data,thresh,PLP,RP,polarity) @@ -45,19 +45,37 @@ % v1.0 08/02/2017 Original version (R2017a) %% SET THRESHOLD -% Min of flat threshold or adaptive from Quiroga 2004 -pars.thresh = min(pars.FIXED_THRESH,... - pars.MULTCOEFF*median(abs(data))/0.6745); - -%% JUST USE FINDPEAKS... -[v,ts,w,p] = findpeaks(polarity*data, ... - 'MinPeakHeight',pars.thresh, ... - 'MinPeakProminence',pars.thresh,... - 'MaxPeakWidth',pars.PLP/3, ... - 'WidthReference','halfheight', ... - 'MinPeakDistance',pars.RP); - -%% MAKE SURE SIGN IS CORRECT -v = v * polarity; + +pk = pars.Polarity * data > pars.Thresh; + +%% REDUCE CONSECUTIVE CROSSINGS TO SINGLE POINTS +z = zeros(size(data)); +pkloc = conv(pk,ones(1,pars.NSaround*2+1),'same')>0; +z(pkloc) = pars.Polarity .* data(pkloc); + + +minTime = 1e-3*pars.RefrTime; % parameter in milliseconds +[ts,pmin] = nigeLab.libs.peakseek(z,minTime*pars.fs,pars.Thresh); +pmin = pmin .* pars.Polarity; + + +%% GET PEAK-TO-PEAK VALUES +PLP = pars.PeakDur*1e-3*pars.fs; % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); +p2pamp = pmax + pmin; + +%% EXCLUDE VALUES OF PMAX <= 0 +pm_ex = pmax<=0; +ts(pm_ex) = []; +p2pamp(pm_ex) = []; +pmax(pm_ex) = []; +pmin(pm_ex) = []; +pW(pm_ex) = []; + + end \ No newline at end of file diff --git a/+nigeLab/@Block/private/SD_MTEO.m b/+nigeLab/@Block/private/SD_MTEO.m new file mode 100644 index 00000000..f1a13af2 --- /dev/null +++ b/+nigeLab/@Block/private/SD_MTEO.m @@ -0,0 +1,64 @@ +function [spikepos,y] = SD_MTEO(f,pars) +%MTEO computes the timestamps of detected spikes in timedomain using a +%multiresolution teager energy operator. +% +% Input parameters: +% in_struc: Input structure which contains +% M: Matrix with data, stored columnwise +% SaRa: Sampling frequency +% k: Resolutation parameter +% optional input parameters: +% none +% Output parameters: +% spikepos: Timestamps of the detected spikes stored columnwise +% +% Description: + % This method is based on the work of J.H.Choi, H.K.Jung, T.Kim "A New Action Potential Detector Using the MTEO + % and Its Effects on Spike Sorting Systems at Low Signal-to-Noise Ratios". It describes a modified TEO. Therefor + % the parameter k is used to vary resolution and is related to frequency of actionpotentials. The input signal is split + % into three TEO outputs followed by a smoothing and a maximum filter. For every time instance the maximum form the + % three smoothed TEO outputs is picked and is indicated in spikepos via variable max. + % + % % Dependencies: +% +% +% Author: F. Lieb, September 2016 + + +fs = pars.fs; +k = pars.k; +L = length(s); + +if size(s, 1) > size(s, 2) + s = s'; +end + + + +%make it zero mean +s2 = f;% - mean(f); + +%this implements fig.4 in the paper +out = zeros(length(s2),length(k)); +for ik = 1:length(k) + %k-teo + out(:,ik) = s2.^2 - circshift(s2,[0, -k(ik)]).*circshift(s2,[0, k(ik)]); + %smoothing window (hamming of length 4*k+1) normalized by noise power + wind = hamming(4*k(ik)+1); + wind = wind./(sqrt(3*sum(wind.^2) + sum(wind)^2)); + + out(:,ik) = conv(out(:,ik),wind,'same'); +end + +%max filter +out = max(out, [], 2); +y = out; + +% extract spike positions from mteo output +switch params.method + case 'numspikes' + spikepos = getSpikePositions(out,fs,s,params); + otherwise + error('unknown method specified'); +end + diff --git a/+nigeLab/@Block/private/SD_PTSD.m b/+nigeLab/@Block/private/SD_PTSD.m new file mode 100644 index 00000000..a5e80390 --- /dev/null +++ b/+nigeLab/@Block/private/SD_PTSD.m @@ -0,0 +1,39 @@ +function [ts,p2pamp,pmin,pW] = SD_PTSD(data,pars) + + + +% PRECISION TIMINIG SPIKE DETECTION +[spkValues, spkTimeStamps] = SpikeDetection_PTSD_core( double(data(:)), pars.Thresh , pars.PeakDur, pars.RefrTime, pars.AlignFlag); + +% %%%%%%%%%%%%%%% Valentina - end +ts = 1 + spkTimeStamps( spkTimeStamps > 0 )'; % +1 added to accomodate for zero- (c) or one-based (matlab) array indexing +pmin = data( ts ); +% % %%%%%%%%%%%% Valentina - begin +% spikesValue(spikesTime<=w_pre+1 | spikesTime>=length(data)-w_post-2)=[]; +% spikesTime(spikesTime<=w_pre+1 | spikesTime>=length(data)-w_post-2)=[]; +% nspk = length(spikesTime); +% +% minTime = 1e-3*pars.RefrTime; % parameter in milliseconds +% [ts,pmin] = nigeLab.libs.peakseek(spkValues,minTime*pars.fs,pars.Thresh ); +% pmin = pmin .* pars.Polarity; + +%% GET PEAK-TO-PEAK VALUES +PLP = floor(pars.PeakDur*1e-3*pars.fs); % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); +p2pamp = pmax + pmin; + +%% EXCLUDE VALUES OF PMAX <= 0 +pm_ex = pmax<= 0 | pmin >= 0; +ts(pm_ex) = []; +p2pamp(pm_ex) = []; +pmax(pm_ex) = []; +pmin(pm_ex) = []; +pW(pm_ex) = []; + + +end + diff --git a/+nigeLab/@Block/private/SD_SNEO.m b/+nigeLab/@Block/private/SD_SNEO.m new file mode 100644 index 00000000..5b87fe51 --- /dev/null +++ b/+nigeLab/@Block/private/SD_SNEO.m @@ -0,0 +1,109 @@ +function [ts,p2pamp,pmin,pW,E] = SD_SNEO(data,pars) +%% SNEOTHRESHOLD Smoothed nonlinear energy operator thresholding detect +% +% [ts,p2pamp,pmin,pW,E] = SNEOTHRESHOLD(data,pars,art_idx) +% +% -------- +% INPUTS +% -------- +% data : 1 x N double of bandpass filtered data, preferrably +% with artifact excluded already, on which to perform +% monopolar spike detection. +% +% pars : Parameters structure from SPIKEDETECTCLUSTER with +% the following fields: +% +% -> SNEO_N \\ number of samples for smoothing window +% -> MULTCOEFF \\ factor to multiply NEO noise threshold by +% +% art_idx : Indexing vector for artifact rejection periods, +% which are temporarily removed so that thresholds +% are not underestimated. +% +% -------- +% OUTPUT +% -------- +% p2pamp : Peak-to-peak amplitude of spikes. +% +% ts : Timestamps (sample indices) of spike peaks. +% +% pmin : Value at peak minimum. (pw) in SPIKEDETECTIONARRAY +% +% dt : Time difference between spikes. (pp) in +% SPIKEDETECTIONARRAY +% +% E : Smoothed nonlinear energy operator value at peaks. +% +% pW : peakWidth, distance between pmax and pmin +% +% By: Max Murphy 1.0 01/04/2018 Original version (R2017a) + +%% GET NONLINEAR ENERGY OPERATOR SIGNAL AND SMOOTH IT +Y = data - mean(data); +Yb = Y(1:(end-2)); +Yf = Y(3:end); +Z = [0, Y(2:(end-1)).^2 - Yb .* Yf, 0]; % Discrete nonlinear energy operator +% Zs = fastsmooth(Z,pars.SNEO_N); +kern = ones(1,pars.SmoothN)./pars.SmoothN; +Zs = fliplr( conv( fliplr(conv(Z,kern,'same')) ,kern,'same')); % the same as the above tri option,(default one here used) but 10x faster +clear('Z','Y','Yb','Yf'); +%% CREATE THRESHOLD FILTER +tmpdata = data; +% tmpdata(art_idx) = []; +tmpZ = Zs; +% tmpZ(art_idx) = []; + +th = pars.MultCoeff * median(abs(tmpZ)); +data_th = pars.MultCoeff * median(abs(tmpdata)); +clear('tmpZ','tmpdata'); +%% PERFORM THRESHOLDING +pk = Zs > th; + +if sum(pk) <= 1 + p2pamp = []; + ts = []; + pmin = []; + dt = []; + E = []; + return +end + +%% REDUCE CONSECUTIVE CROSSINGS TO SINGLE POINTS +z = zeros(size(data)); +% pkloc = repmat(find(pk),pars.NS_AROUND*2+1,1) + (-pars.NS_AROUND:pars.NS_AROUND).'; +% pkloc(pkloc < 1) = 1; +% pkloc(pkloc > numel(data)) = numel(data); +% pkloc = unique(pkloc(:)); +% z(pkloc) = data(pkloc); + +%%%%%%%%%%%%%%%% FB, 5/28/2019 optimized for speed. The above process took 2.9s +%%%%%%%%%%%%%%%% now it takes more or less 0.85s. Same result +pkloc = conv(pk,ones(1,pars.NSaround*2+1),'same')>0; +z(pkloc) = pars.Polarity .* data(pkloc); + + +minTime = 1e-3*pars.RefrTime; % parameter in milliseconds +[ts,pmin] = nigeLab.libs.peakseek(z,minTime*pars.fs,data_th); +pmin = pmin .* pars.Polarity; +E = Zs(ts); + + +%% GET PEAK-TO-PEAK VALUES +PLP = pars.PeakDur*1e-3*pars.fs; % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); +p2pamp = pmax + pmin; + +%% EXCLUDE VALUES OF PMAX <= 0 +pm_ex = pmax<=0; +ts(pm_ex) = []; +p2pamp(pm_ex) = []; +pmax(pm_ex) = []; +pmin(pm_ex) = []; +E(pm_ex) = []; +pW(pm_ex) = []; + +end diff --git a/+nigeLab/@Block/private/SD_SWT2012.m b/+nigeLab/@Block/private/SD_SWT2012.m new file mode 100644 index 00000000..290a8b38 --- /dev/null +++ b/+nigeLab/@Block/private/SD_SWT2012.m @@ -0,0 +1,133 @@ +function [spikepos,y] = SD_SWT2012(in,params) +%SWTEO computes the timestamps of detected spikes in timedomain using a +%stationary wavelet transform algorithm. +% +% Input parameters: +% in_struc: Input structure which contains +% M: Matrix with data, stored columnwise +% SaRa: Sampling frequency +% optional input parameters: +% none +% Output parameters: +% spikepos: Timestamps of the detected spikes stored columnwise +% +% Description: + % This method is based on the work of V.Shalchyan, W.Jensen, and D.Farina "Spike detection and clustering + % with unsupervised wavelet optimization in extracellular neural recordings". + % The algorithem splits signal into 5 levels of wavelet transforms. The coefficients are each thresholded by + % a level depended threshold. The values of the three highest Coefficients are sumed up and smoothed by a Bartlett + % window. The indicator signal is indicated to spikepos where the spikes can be detected. + + % + % % Dependencies: +% +% +% Author: F. Lieb, September 2016 + +s = in.M(:); +L = length(s); +fs = in.SaRa; + + +%prefilter signal +if params.filter + if ~isfield(params,'F1') + params.Fstop = 100; + params.Fpass = 200; + Apass = 0.2; + Astop = 80; + params.F1 = designfilt( 'highpassiir',... + 'StopbandFrequency',params.Fstop ,... + 'PassbandFrequency',params.Fpass,... + 'StopbandAttenuation',Astop, ... + 'PassbandRipple',Apass,... + 'SampleRate',fs,... + 'DesignMethod','butter'); + end + f = filtfilt(params.F1,s); +else + f = s; +end + + +spikelength = 2e-3*fs; + +%thfactor +thfactor = 0.8; + +%do zero padding if the L is not divisible by a power of two +level = 5; + +pow = 2^level; +if rem(L,pow) > 0 + Lok = ceil(L/pow)*pow; + Ldiff = Lok - L; + f = [s; zeros(Ldiff,1)]; +else + Lok = L; +end + +%stationary wavelet transform with 5 levels and symlet +%swdec = swt(s,level,'sym4'); +[swta,swtd] = swt(f,level,'sym5'); + +spd = swtd; + +%estimate noise for each level +%tmp = swdec(1:level,:)'; +tmp = spd'; +thr = thfactor*sqrt(2*log(Lok)).*mad(tmp,1)/0.6745; + +% Denoise. +%swdec2 = zeros(size(swdec)); +swdec2 = zeros(size(spd)); +for k = 1:size(spd,1) + %swdec2(k,:) = wthresh(swdec(k,:),'h',thr(k)); + swdec2(k,:) = wthresh(spd(k,:),'h',thr(k)); +end + +%reconstruction not needed +%sigDEN = iswt(swdec2,'sym4'); + + + + +%compute sum of level +%temp = swdec2(1:level,:)'; +temp = swdec2'; +EW = sum((temp-repmat(mean(temp),Lok,1)).^2); +[~,idx] = sort(EW,'descend'); + +%manifestation variable +Sn = abs(temp(:,idx(1))) + abs(temp(:,idx(2))) + abs(temp(:,idx(3))); + +%smoothing window +w = bartlett(ceil(spikelength/2)); +Tn = conv(Sn,w,'same'); + +y = Tn; + + + + +switch params.method + case 'numspikes' + out = y; + np = 0; + spikepos = zeros(1,params.numspikes); + while (np < params.numspikes) + [~, idxmax] = max(out); + indx = idxmax-ceil(spikelength/2) : idxmax + floor(spikelength/2); + indx = max(1,indx); + indx = min(L,indx); + out(indx) = 0; + + idxx = find(abs(s(indx)) == max(abs(s(indx))),1,'first'); + + %spikepos(np+1) = idxmax; + spikepos(np+1) = indx(1)+idxx-1; + np = np + 1; + end + otherwise + error('unknown method specified'); +end \ No newline at end of file diff --git a/+nigeLab/@Block/private/SD_SWTTEO.m b/+nigeLab/@Block/private/SD_SWTTEO.m new file mode 100644 index 00000000..ea2df87f --- /dev/null +++ b/+nigeLab/@Block/private/SD_SWTTEO.m @@ -0,0 +1,155 @@ +function [ts,p2pamp,pmin,pW,E] = SD_SWTTEO(data,pars) +%SWTTEO Detects Spikes Location using a modified WTEO approach +% Usage: spikepos = swtteo(in); +% spikepos = swtteo(in,params); +% +% Input parameters: +% in_struc: Input structure which contains +% M: Matrix with data, stored columnwise +% SaRa: Sampling frequency +% optional input parameters: +% none +% Output parameters: +% spikepos: Timestamps of the detected spikes stored columnwise +% +% Description: +% swtteo(in,params) computes the location of action potential in +% noisy MEA measurements. This method is based on the work of N. +% Nabar and K. Rajgopal "A Wavelet based Teager Engergy Operator for +% Spike Detection in Microelectrode Array Recordings". The algorithm +% therein was further improved by using a stationary wavelet +% transform and a different thresholding concept. +% For an unsupervised usage the sensitivity of the algorithm can be +% adapted by changing the value of the variable global_fac in line +% 108. A larger value results in fewer detected spikes but also the +% number of false positives decrease. Decreasing this factor makes it +% more sensitive to detect spikes. +% +% References: +% tbd. +% +% +% Author: F. Lieb, February 2016 +% + + +%parse inputs +fs = pars.fs; +TEO = @(x,k) (x.^2 - myTEOcircshift(x,[-k, 0]).*myTEOcircshift(x,[k, 0])); +L = length(data); +data = data(:); % ensure in is column + + +%do zero padding if the L is not divisible by a power of two +% TODO use nextpow2 here +pow = 2^pars.wavLevel; +if rem(L,pow) > 0 + Lok = ceil(L/pow)*pow; + Ldiff = Lok - L; + data = [data; zeros(Ldiff,1)]; +end + + + +%vectorized version: +lo_D = wfilters(pars.waveName); +out_ = zeros(size(data)); +ss = data; +for k=1:pars.wavLevel + %Extension + lf = length(lo_D); + ss = extendswt(ss,lf); + %convolution + swa = conv2(ss,lo_D','valid'); + swa = swa(2:end,:); %even number of filter coeffcients + %apply teo to swt output + + + temp = abs(TEO(swa,1)); + + + + if pars.smoothN + wind = window(pars.winType,pars.smoothN,pars.winPars{:}); + temp2 = conv2(temp,wind','same'); + else + temp2 = temp; + end + + out_ = out_ + temp2; + + + %dyadic upscaling of filter coefficients + lo_D = dyadup(lo_D,0,1); + %updates + ss = swa; +end +clear('ss'); + +% Standard detection + +lambda_swtteo = prctile(out_,99); +lambda_data = pars.MultCoeff*median(abs(data)); +data_th = zeros(size(data)); +data_th(out_>lambda_swtteo) = pars.Polarity .* data(out_>lambda_swtteo); + +minTime = 1e-3*pars.RefrTime; % parameter in milliseconds +[ts,pmin] = nigeLab.libs.peakseek(data_th,minTime*pars.fs,lambda_data); +pmin = pmin .* pars.Polarity; +E = out_(ts); + +%% GET PEAK-TO-PEAK VALUES +PLP = pars.PeakDur*1e-3*pars.fs; % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); +p2pamp = pmax + pmin; + +%% EXCLUDE VALUES OF PMAX <= 0 +pm_ex = pmax<=0; +ts(pm_ex) = []; +p2pamp(pm_ex) = []; +pmax(pm_ex) = []; +pmin(pm_ex) = []; +E(pm_ex) = []; +pW(pm_ex) = []; + + + + +function y = extendswt(x,lf) +%EXTENDSWT extends the signal periodically at the boundaries +[r,c] = size(x); +y = zeros(r+lf,c); +y(1:lf/2,:) = x(end-lf/2+1:end,:); +y(lf/2+1:lf/2+r,:) = x; +y(end-lf/2+1:end,:) = x(1:lf/2,:); + + +function X = myTEOcircshift(Y,k) +%circshift without the boundary behaviour... + +colshift = k(1); +rowshift = k(2); + +temp = circshift(Y,k); + +if colshift < 0 + temp(end+colshift+1:end,:) = flipud(Y(end+colshift+1:end,:)); +elseif colshift > 0 + temp(1:1+colshift-1,:) = flipud(Y(1:1+colshift-1,:)); +else + +end + +if rowshift<0 + temp(:,end+rowshift+1:end) = fliplr(Y(:,end+rowshift+1:end)); +elseif rowshift>0 + temp(:,1:1+rowshift-1) = fliplr(Y(:,1:1+rowshift-1)); +else +end + +X = temp; + diff --git a/+nigeLab/@Block/private/SD_TIFCO.m b/+nigeLab/@Block/private/SD_TIFCO.m new file mode 100644 index 00000000..047ed389 --- /dev/null +++ b/+nigeLab/@Block/private/SD_TIFCO.m @@ -0,0 +1,304 @@ +function [ts,p2pamp,pmin,pW,E] = SD_TIFCO(data,pars) +%GABOR computes the timestamps of detected spikes in timedomain using a +%Gabor Transform based spike detection. +% +% Input parameters: +% in_struc: Input structure which contains +% M: Matrix with data, stored columnwise +% SaRa: Sampling frequency +% optional input parameters: +% sx: +% Output parameters: +% spikepos: Timestamps of the detected spikes stored columnwise +% +% Description: + % This method is based on the work F.Lieb "...". This algorithm uses + % a special case of Short-Time-Fourier Transform using a filter. As spikes + % can be specified by a certain time-frequency behavior, frequencies below + % and above a certain level will be ignored by specifying them within the + % GABOR transform. By applying a moving average to the time-frequency + % coefficients the spike form is enforced. Then a STF is used on the singal + % generating an indicator signal indicated in spikepos. + % % Dependencies: +% +% +% Author: f. Lieb, September 2016 + +fs = pars.fs; +L = length(data); +W = window(pars.winType, pars.winL*fs ,pars.winPars{:}); + +a = 1; +M = 100; + +[c,freq] = dgtsf(data,W,a,pars.fMin,pars.fMax,M,fs); + +numt = 1; +numf = size(c,2); +%numf = 100; + +W = convFreqWeights(c,numf,numt); + + +sx1 = sum(W,2); + + +% Standard detection + +lambda_tifco = prctile(sx1,99); +lambda_data = pars.MultCoeff*median(abs(data)); +data_th = zeros(size(data)); +data_th(sx1>lambda_tifco) = pars.Polarity .* data(sx1>lambda_tifco); + +minTime = 1e-3*pars.RefrTime; % parameter in milliseconds +[ts,pmin] = nigeLab.libs.peakseek(data_th,minTime*pars.fs,lambda_data); +pmin = pmin .* pars.Polarity; +E = sx1(ts); + +%% GET PEAK-TO-PEAK VALUES +PLP = pars.PeakDur*1e-3*pars.fs; % from ms to samples +tloc = repmat(ts,2*PLP+1,1) + (-PLP:PLP).'; +tloc(tloc < 1) = 1; +tloc(tloc > numel(data)) = numel(data); +[pmax,Imax] = max(data(tloc)); +pW = abs(Imax-PLP); +p2pamp = pmax + pmin; + +%% EXCLUDE VALUES OF PMAX <= 0 +pm_ex = pmax<=0; +ts(pm_ex) = []; +p2pamp(pm_ex) = []; +pmax(pm_ex) = []; +pmin(pm_ex) = []; +E(pm_ex) = []; +pW(pm_ex) = []; + +function [c, freq] = dgtsf(f,g,a,fmin,fmax,M,fs,tinv) +%DGTSF Discrete Gabor transform with specified frequency range +% Usage: c = dgtsf(f,g,a,fmin,fmax,M,fs); +% c = dgtsf(f,g,a,fmin,fmax,M,fs,'timeinv'); +% [c,freq] = dgtsf(f,g,a,fmin,fmax,M,fs); +% +% Input parameters: +% f : Signal +% a : Length of time shifts +% M : Number of frequency channels +% g : Window function (currently supported: hanning window only) +% fmin : Minimum frequency +% fmax : Maximum frequency +% tinv : Option to compute the time invariant version +% +% Output parameters: +% c : Time-Frequency respresentation +% freq : Frequency axis ranging from [fmin fmax] +% +% Description: +% dgtsf(f,g,a,fmin,fmax,M,fs) computes the time-frequency +% represenation of f very efficiently if M is relatively small. The +% option tinv computes the time-invariant phase (shift and modulation +% operator are reversed). This phase option is requiered in the +% spikeDetection Algorithm. +% +% The window g can be a vector containing the window, or a cell array +% where the first entry is the window type as a string. If the ltfat +% toolbox is installed a variety of different windows are available, +% if not only the hanning window is supported so far. The second +% parameter in the cell array is the size of the window. For example +% {'hann', 100} gives a hanning window with 100 datapoints. Please be +% reminded that the windowing is done in the frequency domain, so a +% narrow window will give a good frequency resolution. +% +% +% Dependencies: +% findFreqIndx +% ~gabwin if ltfat exists on path +% +% Author: F. Lieb, Janurary 2016 +% +if nargin<8 + tinv = ''; +end + +%get signal length +L=size(f,1); +if L == 1 + L = size(f,2); + f = f'; +end + +numelectrodes = size(f,2); + +%get window function +if (iscell(g)) + if exist('gabwin','file') == 2 + if strcmp(g{1},'gauss') + winsize = L; + glh = L/2; + g2 = fftshift(gabwin(g,a,L,L)); + else + winsize = ceil(g{2}); + glh=floor(winsize/2); + g2 = fftshift(gabwin(g,a,winsize)); + end + + else + winsize = ceil(g{2}); + glh = floor(winsize/2); + g2 = hann(g{2},'periodic'); + g2 = g2./norm(g2); + end + offset = 0; +else + g2 = g; + winsize = length(g2); + glh=floor(winsize/2); + [~,ymaxidx] = max(g2); + offset = ymaxidx; +end + +%number of time channels +N=L/a; + +%find the most approriate M +[M, freq,dfindx] = findFreqIndx(L,fs,fmin,fmax,M); + +%container for result +if isa(f,'single') + c=zeros(N,M,numelectrodes,'single'); +else + c=zeros(N,M,numelectrodes,'double'); +end + +%map to frq axis: +frqdst = fs/L; +fminidx = floor(fmin/frqdst); +if (iscell(g)) + win_range = fminidx-glh+1:fminidx+glh; +else + win_range = fminidx-offset+1:fminidx+winsize-offset ; +end +win_range(win_range<1) = win_range(win_range<1)+L; + +df = diff(dfindx); +df = df(1); + +%fft of input signal +fhat = fft(f); + +%loop over all frequency bins +for k = 1:M + %pointwise multiplication + c([end-floor(winsize/2)+1:end,1:ceil(winsize/2)],k,:) = bsxfun(@times,fhat(win_range,:),g2); + %update window range + win_range = win_range + df; + win_range(win_range>L) = win_range(win_range>L)-L; +end +%ifft of output matrix +c = ifft(c); + +%todo: do this computation inside the loop to save some comp time... +% (couldnt get this to work...) +if (strcmp(tinv,'timeinv')) + %apply time invariant factor: + TimeInd = (0:(L/a-1))*a; + phase = exp(2*pi*1i*(TimeInd'*dfindx)/L); + c = bsxfun(@times,c,phase); +end + +function [M_out, freq, freqindx] = findFreqIndx(L,fs,fmin,fmax,M) +%FINDFREQINDX finds the corresponding frequency indices for specific set of +% parameters +% +% If fmin and fmax as well as L (signal length) and fs are specified, not +% all M are possible due to the discretization and the limits due to fs +% and L. +% This function looks for a possible number of frequency bins M_out which +% is close the specified number M +% +% Author: F. Lieb, February 2016 +% + +%get minimum freq resolution: +dfreq = fs/L; + +%get indx of minimum frequency: +fminidx = floor(fmin/dfreq); +%get indx of maximum frequency: +fmaxidx = ceil(fmax/dfreq); + +M_out = findFreqIndx_helper(fmaxidx,fminidx,M); +freq = linspace(fminidx*dfreq,fmaxidx*dfreq,M_out); +freqindx = linspace(fminidx,fmaxidx, M_out); + +function W2 = convFreqWeights(coeff, numfn, numtn) +%CONVFREQWEIGHTS sliding window neighborhood +% Usage: W = convFreqWeights(coeff,numfn, numtn); +% +% Input parameters: +% coeff : Time-Frequency matrix with columnwise frequency content +% numfn : Number of neighbors in frequency direction +% numtn : Number of neighbors in time direction +% +% Output parameters: +% W2 : Convolved Time Frequency respresentation +% +% Description: +% convFreqWeights(coeff,numfn,numtn) computes the convolution of a +% mean-window with neighboring tf-coefficients +% +% Author: F. Lieb, February 2016 +% + +%get inputsize +[cm, cn] = size(coeff); + +if cn > cm + coeff = coeff.'; +end + +factor = 2; +%get kernel +neigh = ones(numtn,numfn); +neigh = neigh./(norm(neigh(:),1)); +windowcenter = ceil(numtn/2); +windowcenter2= ceil(numfn/2); + +%container +if isa(coeff,'single') + W = zeros(cm+numtn-1,cn + numfn-1,'single'); +else + W = zeros(cm+numtn-1,cn + numfn-1,'double'); +end + +%extend the borders +W(windowcenter:cm+windowcenter-1,windowcenter2:cn+windowcenter2-1)=abs(coeff).^factor;%abs(coeff).^2; +W(1:windowcenter-1,:) = flipud( W(windowcenter:2*(windowcenter-1),:) ); +W(cm+windowcenter:end,:) = flipud( W(cm-numtn+2*windowcenter :cm+windowcenter-1,:) ); +W(:,1:windowcenter2-1)= fliplr( W(:,windowcenter2:2*(windowcenter2-1))); +W(:,cn+windowcenter2:end)= fliplr( W(:,cn-numfn+2*windowcenter2:cn+windowcenter2-1)); + +%do convolution +W2 = (conv2(W,neigh,'valid')).^(1/factor); + +function out = findFreqIndx_helper(io,iu,M) +%FINDFREQINDX_HELPER helper function for findFreqIndx +% +% looks for an equidistant discretization of the intervall [iu io] with M +% elements or some number close to M. +% +% Author: F. Lieb, February 2016 +% + +isize = io-iu; + +%alldiv(isize) +d=isize./(isize:-1:2); +d=d(d==round(d)); + +tmp = isize./d + 1; + +%find element closest to M +[~,idx] = min(abs(tmp-M)); +out = tmp(idx(1)); + + diff --git a/+nigeLab/@Block/private/SD_WTEO.m b/+nigeLab/@Block/private/SD_WTEO.m new file mode 100644 index 00000000..610efa8b --- /dev/null +++ b/+nigeLab/@Block/private/SD_WTEO.m @@ -0,0 +1,93 @@ +function [spikepos,y] = SD_WTEO(s,pars) +%MTEO computes the timestamps of detected spikes in timedomain using a +%wavelet based teager energy operator. +% +% Input parameters: +% in_struc: Input structure which contains +% M: Matrix with data, stored columnwise +% SaRa: Sampling frequency +% optional input parameters: +% none +% Output parameters: +% spikepos: Timestamps of the detected spikes stored columnwise +% +% Description: + % This method is based on the work of N.Naber and F.Franke "a wavelet + % based Teager Energy operator for Spike Detection + % in Microelectrode Array Recordings". It descriebes a wavelet + % based TEO. Hence a low-pass filter using first and second level + % approximation coefficients of the discrete wavelet transform is + % apllied to each sub-band, followed by TEO. Each output is + % thresholded independently and then up sampled to the orignial + % signal lenght indicated in spikepos. + + % + % % Dependencies: +% +% +% Author: F. Lieb, September 2016 + + +%parse inputs +fs = pars.fs; +L = length(s); +wavLevel = 2; +wavelet = 'sym7'; + +TEO = @(x) (x.^2 - circshift(x,[0, -1]).*circshift(x,[0, 1])); + + + +if size(s, 1) > size(s, 2) + s = s'; +end + +%do zero padding if the L is not divisible by a power of two +pow = 2^wavLevel; +if rem(L,pow) > 0 + Lok = ceil(L/pow)*pow; + Ldiff = Lok - L; + s = [s; zeros(Ldiff,1)]; + L = Lok; +end + + + +f = s; + + + + +%My interpretation +%[SWTa,SWTd] = swt(s,wavLevel,wavelet); +%out = TEO(SWTa); +lo = wfilters(wavelet); +zz = wconv1(f,lo,'same'); +cA1 = zz(1:2:end); +out1 = TEO(cA1); + +zz = wconv1(cA1,lo,'same'); +cA2 = zz(1:2:end); +out2 = TEO(cA2); + +%o1 = interp1(1:length(out1),out1,linspace(1,length(out1),L),'nearest'); +o1 = dyadup(out1,1); +o1 = o1(1:end-1); +%o2 = interp1(1:length(out2),out2,linspace(1,length(out2),L),'nearest'); +o2 = dyadup(dyadup(out2,1),1); +o2 = o2(1:end-3); + +out = [o1; o2]; + + %sum over approximation coefficients +out = sum(out,1); +y = out; +switch pars.method + case 'numspikes' + spikepos = getSpikePositions(out,fs,in.M,pars); + case 'lambda' + thout = wthresh(out,'h',pars.lambda); + spikepos = getSpikePositions(thout,fs,s,pars); + otherwise + error('unknown method specified'); +end \ No newline at end of file diff --git a/+nigeLab/@Sort/Sort.m b/+nigeLab/@Sort/Sort.m index a359348d..68ce6bd1 100644 --- a/+nigeLab/@Sort/Sort.m +++ b/+nigeLab/@Sort/Sort.m @@ -161,7 +161,7 @@ function delete(sortObj) return; end - if checkForUnconfirmedChanges(sortObj.UI.SpikeImage,true) + if checkForConfirmedChanges(sortObj.UI.SpikeImage,true) return; end @@ -209,7 +209,7 @@ function delete(sortObj) methods (Access=public) setChannel(sortObj,~,~) % Set the current channel in the UI setClass(sortObj,class) % Set the current sort class - saveData(sortObj) % Save the sorting + saveData(sortObj,ch) % Save the sorting end % PROTECTED diff --git a/+nigeLab/@Sort/saveData.m b/+nigeLab/@Sort/saveData.m index 0707c274..8b63c813 100644 --- a/+nigeLab/@Sort/saveData.m +++ b/+nigeLab/@Sort/saveData.m @@ -1,4 +1,4 @@ -function saveData(sortObj) +function saveData(sortObj,ch) %% SAVEDATA Save data on the disk % % sortObj.saveData; @@ -14,15 +14,21 @@ function saveData(sortObj) % Writes the sorting for all channels to the disk. %% -fprintf(1,'Saving...%03g%%\n',0); -iTotal = numel(sortObj.Blocks) * numel(sortObj.Channels.Mask); +if nargin < 2 + ch = sortObj.Channels.Mask; +else + ch = ch(:)'; +end +chans = sprintf('%d,',ch); +fprintf(1,'Saving channels %s...%03g%%\n',chans(1:end-1),0); +iTotal = numel(sortObj.Blocks) * numel(ch); iCount = 0; for iBlock = 1:numel(sortObj.Blocks) blockObj = sortObj.Blocks(iBlock); - for iCh = sortObj.Channels.Mask + for iCh = ch idx = sortObj.spk.block{iCh} == iBlock; - blockObj.Channels(iCh).Sorted.unlockData; - blockObj.Channels(iCh).Sorted.value = sortObj.spk.class{iCh}(idx); + blockObj.Channels(sortObj.Channels.Idx(iCh,iBlock)).Sorted.unlockData; + blockObj.Channels(sortObj.Channels.Idx(iCh,iBlock)).Sorted.value = sortObj.spk.class{iCh}(idx); iCount = iCount + 1; fprintf(1,'\b\b\b\b\b%03g%%\n',round((iCount/iTotal)*100)); end diff --git a/+nigeLab/@nigelObj/nigelObj.m b/+nigeLab/@nigelObj/nigelObj.m index 3b6186d5..1b30c432 100644 --- a/+nigeLab/@nigelObj/nigelObj.m +++ b/+nigeLab/@nigelObj/nigelObj.m @@ -148,7 +148,7 @@ ChildContainer % Container for .Children Dependent property ChildListener event.listener % Listens for changes in object Children GUIContainer nigeLab.libs.DashBoard % Container for handle to GUI (nigeLab.libs.DashBoard) - TreeNodeContainer uiw.widget.TreeNode % Container for handle to Tree nodes in the nigelTree obj + TreeNodeContainer % Container for handle to Tree nodes in the nigelTree obj ParentListener event.listener % Listens for changes in object Parent PropListener event.listener % Listens for changes in properties of this object SortGUIContainer nigeLab.Sort % Container for handle to Spike Sorting GUI (nigeLab.Sort) @@ -547,21 +547,29 @@ function delete(obj) n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); end case '.' - switch class(nigelObj) - case 'nigeLab.Block' - if ~isempty(nigelObj(1).Pars) - Ffields = nigelObj(1).Pars.Block.Fields; - idx = strcmpi(Ffields,s(1).subs); - if any(idx) - n = numel(nigelObj); + meta = metaclass(nigelObj); + methods = meta.MethodList; + methodIdx = strcmp({methods.Name},s(1).subs); + if any(methodIdx) + thisMethod = methods(methodIdx); + n = numel(thisMethod.OutputNames); + else + switch class(nigelObj) + case 'nigeLab.Block' + if ~isempty(nigelObj(1).Pars) + Ffields = nigelObj(1).Pars.Block.Fields; + idx = strcmpi(Ffields,s(1).subs); + if any(idx) + n = numel(nigelObj); + else + n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); + end else n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); end - else + otherwise n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); - end - otherwise - n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); + end end otherwise n = builtin('numArgumentsFromSubscript',nigelObj,s,indexingContext); @@ -3159,13 +3167,28 @@ function addChild(obj,childData,idx) % Get the current "Spike Detection method," which gets % added onto the front part of the _Spikes and related % folders + + % default. No attribute in teh name + p.Folder = sprintf(strrep(p.Folder,'\','/'),... + ''); if ~isfield(obj.HasParsInit,'SD') obj.updateParams('SD'); elseif ~obj.HasParsInit.SD obj.updateParams('SD'); end - p.Folder = sprintf(strrep(p.Folder,'\','/'),... - obj.Pars.SD.ID.(F{iF})); + + % Look for ID in Pars + ParsFields = fieldnames(obj.Pars); + ParsWithID = find(cellfun(@(F) isfield(obj.Pars.(F),'ID'),ParsFields))'; + for ii=ParsWithID + if isfield(obj.Pars.(ParsFields{ii}).ID,(F{iF})) + % if found put the ID in the foldername path + p.Folder = sprintf(strrep(p.Folder,'\','/'),... + obj.Pars.(ParsFields{ii}).ID.(F{iF})); + break; + end + end + end % Set folder name for this particular Field