From 1942e1a20713a103f4b3ef9a5e5e66af6a1480c5 Mon Sep 17 00:00:00 2001 From: Nabarb Date: Mon, 2 Nov 2020 16:24:58 +0100 Subject: [PATCH] Feature/visualization tools1 (#95) * Proposed fix for 2017 nargout bug This tries to fix the bug #86 emerging in matlab 2017 editions (and maybe others) where the builtin numArgumentsFromSubscript is not behaving as expected, Instead of returning the correct nargout it basically returns numel(nigelObj) creating all sorts of problems, such as the too many output args when calling tankObj.nigelDash from base. This solution relies on the number of declared args in the method declaration to determine nargout at this stage. This does not affect the effective number of output arguments which are still controlled by the number of variables on the left of an assignment. Bear in mind that this will effect the number of output args when called without an explicit assignment. * fixes Fixed small bug in buildWorkerCoonfigScript. Now works as intended, it was getting the two modalities mixed up. Fixed ismethod issue in numArgsFromSub. in nigelObj. references #86. * Critical fix in rereferencing Found a critical in rereferencing that was actually preventing the rereferencing from happening. * bugfix getspiketime (when sorted) getspiketime was returning an empty array ehen called with optional "class" input. (Was looking at the wrong field of the event file) fixed. * Feature/new sort (#88) * Sort new features Added a pletora of features to the Sort interface. Fixed and improved the way the SortFiles are saved and loaded, (getClus getSort). Fixed small bugs here and there, HighDimsUi now has the same interface language as dashboard (with nigelpanels etc). A new button was added to find the most convenient projection according to LDA. When pressed moves the plane visualized in featuresUI in the position that best separates the selected cluster. Thanks DataHigh! ginput Changed the cross-hair color to improve visibility. multicallbackwrap wasn't really working. fixed. NigelButton fixed a small bug that was overwriting WindowButtonUpFnc of the figure it was put on. Now it uses multiCallbackWrap to avoid overwriting. FaturesUI The 2D graph now cna be zoomed and panned. Use the MouseWheel to zoom. MouseWheel + alt or control to pan horizontally and vertically. Saving performance greatly improved by saving only the clusters that actually need saving. Each channels starts with a UnsavedChanges target to false and a ConfirmedChanges to false. When a change is done the UnsavedChanges goes to true. When confirmed, ConfirmedChanges goes to true. SpikeImage; SortUi; * Fixed zooming in FeaturesUI Now the two graphs have linked axes that change and shift together. Also, navigation is now homogeneous between the two figures * Fixed bug in Sort/savedata fixed a bug in savedata. Program was crashing when a channel was turned off in mask * performance update Small performance update in plotting and updating spike profiles and features. * Added pca and random projection Added the possibility to project data onto the plane of max variance as well as to a random projection plane. This is accessible from the HighDimsUi * Feature/new sd (#92) # New Spike Detection Added by: Federico (FB) ## Features ### Change in Structure `doSD` has been essentially re-written. * All the SD algorithms were moved to their own functions in the block private folder. * The same was done for feature extraction and artifact rejection. #### Motivation We want the pipeline to be more flexible, particularly in the spike detection step where many labs or users have their own preferred method. Reorganizing the structure of `doSD` in this way allows the specification of SD choices both in the `defaults.SD` file, as well as on an ad hoc basis through a user interface. #### To-Do * A short wiki entry will follow with the style specifications and naming convention of the files. ### Spike Detection Algorithm Updates * Revised and fixed `AdaptThresh`, `HardThresh`, `SNEO`, and `TIFCO`. * Fixed associated `SD` parameters. ### `configSD` Added a GUI to configure the SD parameters. #### Strategy It randomly extracts 1 minute of data from the recording, providing a graphical interface that displays the data and any spikes detected as the end-user tries different parameter configurations until they are satisfied that the current parameter settings can detect spikes well in the data. and you can try different parameters on the extracted data. * Artifact rejection parameters can be selected in the same way. * The "current channel" can be changed by the end-user. * noncritical bugfixes Fixed a potential bug in rereferencing Changed swtteo to compute wavelet coefficient at runtime instead of during the parameters loading. This makes it possible to try different wavelets in the configSD intereface * bugfix Fixed a couple of glitches in the configSD interface, Also, SD defaults are modified to some more sensible values * Fixed plotwaves * Explore data class added Added the explore data class to view snippets of data interactively. added datascrollerAxis class: mimics timescrolleraxis but with data instead of events and videos.It's used by exploredata class. added the possibility of plotting plotWaves into an existing axis * plotWaves reducees plots Now plotwaves reduces plots as well. Zooming is not suggested because of this. Inside the ExplroeData utility zooming is handled from the dataScrollerAxis dataScrollerAxis now uses the full lineplotreducer utility and zooming is therefore enabled Also fixed a small issue with nigelColors * bigfix Fixed small bugs and glitches for plotWaves and DataScrollerAxis * Spike Overlay Added the spike Overlay for plotWaved. It works by purgin the y-reduced signal by all non spike parts and plotting it on top of either CAR or Filt already plotted. It is selected by passing Spikes as first input * exploreData bugfix fixed small glitch * Updates to get RMS values of filtered data * not sure if change to analyzeRMS.m is technically correct, but it ran *updated block methods list to reflect changes in plotWaves * Visualization bugfix Added the empty case to ToString Added a fialsafe in case time is not set properly in block. fixed the usage of reduce_to_width Co-authored-by: Federico Co-authored-by: Page Hayley --- +nigeLab/+defaults/+SD/ART_HardThresh.m | 9 + +nigeLab/+defaults/+SD/FEAT_ica.m | 5 + +nigeLab/+defaults/+SD/FEAT_pca.m | 10 + +nigeLab/+defaults/+SD/FEAT_wavelet.m | 10 + +nigeLab/+defaults/+SD/SD_AdaptThresh.m | 10 + +nigeLab/+defaults/+SD/SD_HardThresh.m | 9 + +nigeLab/+defaults/+SD/SD_PTSD.m | 9 + +nigeLab/+defaults/+SD/SD_SNEO.m | 10 + +nigeLab/+defaults/+SD/SD_SWTTEO.m | 16 + +nigeLab/+defaults/+SD/SD_TIFCO.m | 15 + +nigeLab/+defaults/+SD/WIP_SD_MTEO.m | 8 + +nigeLab/+defaults/+SD/WIP_SD_WTEO.m | 8 + +nigeLab/+defaults/AutoClustering.m | 11 +- +nigeLab/+defaults/SD.m | 200 ++++--- +nigeLab/+defaults/Sort.m | 9 +- +nigeLab/+defaults/nigelColors.m | 10 +- +nigeLab/+evt/dataScrolled.m | 11 + +nigeLab/+evt/saveData.m | 32 ++ +nigeLab/+libs/@ChannelUI/ChannelUI.m | 3 +- .../@DataScrollerAxis/DataScrollerAxis.m | 197 +++++++ +nigeLab/+libs/@ExploreData/ExploreData.m | 154 ++++++ +nigeLab/+libs/@FeaturesUI/FeaturesUI.m | 87 ++- +nigeLab/+libs/@FeaturesUI/PlotFeatures.m | 28 +- +nigeLab/+libs/@FeaturesUI/PlotQuality.m | 4 +- +nigeLab/+libs/@FeaturesUI/ResetFeatureAxes.m | 3 - +nigeLab/+libs/@FeaturesUI/UpdateSil.m | 1 + +nigeLab/+libs/@HighDimsUI/HighDimsUI.m | 249 ++++++++- +nigeLab/+libs/@SortUI/SortUI.m | 2 +- +nigeLab/+libs/@SpikeImage/SpikeImage.m | 241 +++++---- +nigeLab/+libs/@configSD/configSD.m | 504 ++++++++++++++++++ +nigeLab/+libs/@nigelButton/nigelButton.m | 17 +- +nigeLab/+libs/@nigelPanel/nigelPanel.m | 1 - .../+utils/@LinePlotReducer/reduce_to_width.m | 113 ++++ +nigeLab/+utils/ToString.m | 28 + +nigeLab/+utils/buildWorkerConfigScript.m | 4 +- .../SNEOThreshold.m => +utils/fastsmooth.m} | 116 ---- +nigeLab/+utils/ginput.m | 2 +- +nigeLab/+utils/multiCallbackWrap.m | 27 +- +nigeLab/@Block/Block.m | 7 +- +nigeLab/@Block/analyzeRMS.m | 8 +- +nigeLab/@Block/doReReference.m | 9 +- +nigeLab/@Block/doSD.m | 337 ++++++------ +nigeLab/@Block/getClus.m | 131 +++-- +nigeLab/@Block/getSort.m | 84 +-- +nigeLab/@Block/getSpikeTimes.m | 2 +- +nigeLab/@Block/plotSpikes.m | 4 +- +nigeLab/@Block/plotWaves.m | 170 ++++-- ...rdArtifactRejection.m => ART_HardThresh.m} | 17 +- +nigeLab/@Block/private/FEAT_ica.m | 58 ++ +nigeLab/@Block/private/FEAT_pca.m | 48 ++ +nigeLab/@Block/private/FEAT_wavelet.m | 147 +++++ .../{AdaptiveThreshold.m => SD_AdaptThresh.m} | 45 +- .../{ThresholdDetection.m => SD_HardThresh.m} | 48 +- +nigeLab/@Block/private/SD_MTEO.m | 64 +++ +nigeLab/@Block/private/SD_PTSD.m | 39 ++ +nigeLab/@Block/private/SD_SNEO.m | 109 ++++ +nigeLab/@Block/private/SD_SWT2012.m | 133 +++++ +nigeLab/@Block/private/SD_SWTTEO.m | 155 ++++++ +nigeLab/@Block/private/SD_TIFCO.m | 304 +++++++++++ +nigeLab/@Block/private/SD_WTEO.m | 93 ++++ +nigeLab/@Sort/Sort.m | 4 +- +nigeLab/@Sort/saveData.m | 18 +- +nigeLab/@nigelObj/nigelObj.m | 51 +- 63 files changed, 3513 insertions(+), 745 deletions(-) create mode 100644 +nigeLab/+defaults/+SD/ART_HardThresh.m create mode 100644 +nigeLab/+defaults/+SD/FEAT_ica.m create mode 100644 +nigeLab/+defaults/+SD/FEAT_pca.m create mode 100644 +nigeLab/+defaults/+SD/FEAT_wavelet.m create mode 100644 +nigeLab/+defaults/+SD/SD_AdaptThresh.m create mode 100644 +nigeLab/+defaults/+SD/SD_HardThresh.m create mode 100644 +nigeLab/+defaults/+SD/SD_PTSD.m create mode 100644 +nigeLab/+defaults/+SD/SD_SNEO.m create mode 100644 +nigeLab/+defaults/+SD/SD_SWTTEO.m create mode 100644 +nigeLab/+defaults/+SD/SD_TIFCO.m create mode 100644 +nigeLab/+defaults/+SD/WIP_SD_MTEO.m create mode 100644 +nigeLab/+defaults/+SD/WIP_SD_WTEO.m create mode 100644 +nigeLab/+evt/dataScrolled.m create mode 100644 +nigeLab/+evt/saveData.m create mode 100644 +nigeLab/+libs/@DataScrollerAxis/DataScrollerAxis.m create mode 100644 +nigeLab/+libs/@ExploreData/ExploreData.m create mode 100644 +nigeLab/+libs/@configSD/configSD.m create mode 100644 +nigeLab/+utils/@LinePlotReducer/reduce_to_width.m create mode 100644 +nigeLab/+utils/ToString.m rename +nigeLab/{@Block/private/SNEOThreshold.m => +utils/fastsmooth.m} (70%) rename +nigeLab/@Block/private/{HardArtifactRejection.m => ART_HardThresh.m} (78%) create mode 100644 +nigeLab/@Block/private/FEAT_ica.m create mode 100644 +nigeLab/@Block/private/FEAT_pca.m create mode 100644 +nigeLab/@Block/private/FEAT_wavelet.m rename +nigeLab/@Block/private/{AdaptiveThreshold.m => SD_AdaptThresh.m} (62%) rename +nigeLab/@Block/private/{ThresholdDetection.m => SD_HardThresh.m} (66%) create mode 100644 +nigeLab/@Block/private/SD_MTEO.m create mode 100644 +nigeLab/@Block/private/SD_PTSD.m create mode 100644 +nigeLab/@Block/private/SD_SNEO.m create mode 100644 +nigeLab/@Block/private/SD_SWT2012.m create mode 100644 +nigeLab/@Block/private/SD_SWTTEO.m create mode 100644 +nigeLab/@Block/private/SD_TIFCO.m create mode 100644 +nigeLab/@Block/private/SD_WTEO.m 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