Filter and Threshold Raw Data Into Snippets
Import raw streaming data into Matlab using TDTbin2mat
Digitally filter the single unit data using TDTdigitalfilter
Threshold and extract snippets using TDTthresh
Contents
Housekeeping
Clear workspace and close existing figures. Add SDK directories to Matlab path.
close all; clear all; clc; [MAINEXAMPLEPATH,name,ext] = fileparts(cd); % \TDTMatlabSDK\Examples DATAPATH = fullfile(MAINEXAMPLEPATH, 'ExampleData'); % \TDTMatlabSDK\Examples\ExampleData [SDKPATH,name,ext] = fileparts(MAINEXAMPLEPATH); % \TDTMatlabSDK addpath(genpath(SDKPATH));
Importing the Data
This example assumes you downloaded our example data sets (link) and extracted it into the \TDTMatlabSDK\Examples\ directory. To import your own data, replace BLOCKPATH with the path to your own data block.
In Synapse, you can find the block path in the database. Go to Menu > History. Find your block, then Right-Click > Copy path to clipboard.
BLOCKPATH = fullfile(DATAPATH,'Algernon-180308-130351');
Set up the variables for the data you want to extract. We will extract channel 1 from the Wav1 data store, created by the Stream Data Storage gizmo, and apply our threshold at -125uV to create 30 sample snippets.
STORE = 'Wav1'; CHANNEL = 1; THRESH = -125e-6; % define snippet when it crosses this threshold NPTS = 30; % size of the snippet OVERLAP = 0; % set to 1 to allow double crossings within same window to count as a different snippet
Now read the specified data from our block into a Matlab structure.
data = TDTbin2mat(BLOCKPATH, 'STORE', STORE, 'CHANNEL', CHANNEL);
read from t=0.00s to t=61.23s
Use TDTdigitalfilter to filter the streaming waveforms
We are interested in single unit activity in the 300Hz-5000Hz band.
su = TDTdigitalfilter(data, STORE, [300 5000]);
Use TDTthresh to extract snippet events
Mode is set to ‘manual’ in this example for extracting snippets around a hard threshold.
su = TDTthresh(su, STORE, 'MODE', 'manual', 'THRESH', THRESH, 'NPTS', NPTS, 'OVERLAP', OVERLAP); % You could also use 'auto' detection, which sets a moving threshold base % on a multiple of the RMS of a previous time period. For example, % set the threshold at 6*RMS of the previous 5 seconds of data %su = TDTthresh(su, STORE, 'MODE', 'auto', 'POLARITY', -1, 'STD', 6, 'TAU', 5);
window size is 30, without overlap. using absolute threshold method, set to -125.00uV.
And that’s it! Your data is now in Matlab and thresholded. TDTthresh adds a new store into the data structure called data.snips.Snip that we will now use for plotting. First, we’ll get some properties of the snippets that we’ll use later for scatter plots
maxvals = max(su.snips.Snip.data, [], 2)*1e6; minvals = min(su.snips.Snip.data, [], 2)*1e6;
Begin plotting
Plot the single unit data with snippet overlay.
subplot(2,3,[1 2 3]) ts = (1:numel(su.streams.(STORE).data(1,:)))/su.streams.(STORE).fs; plot(ts, su.streams.(STORE).data(1,:)*1e6, 'color', [.7 .7 .7]); axis([ts(1) ts(end) min(minvals) max(maxvals)]); ax1 = axis(); hold on; % Position the figure set(gcf, 'Position', [100 100 1000 700])
Plot the threshold line.
if numel(su.snips.Snip.thresh) == 1 thresh = su.snips.Snip.thresh*ones(size(ts)); else thresh = su.snips.Snip.thresh; end plot(ts, thresh*1e6, 'r--', 'LineWidth', 2) ylabel('uV'); xlabel('time (s)');
Define a color map based on the minimum values. The dots are not colored by any complex sorting algorithm, just by magnitude. We’ll use this coloring to track the snippets through the rest of the plots. This uses the vals2colormap function, which can be found at: https://github.com/vistalab/vistateach/blob/master/cogneuro/tutorial1_timeseries/vals2colormap.m.
colors = vals2colormap(minvals, 'spring', [min(minvals)*.9 max(minvals)]);
Overlay the min values as dots on the streaming plot.
scatter(su.snips.Snip.ts, minvals, 10, colors, 'filled'); legend({'stream', 'threshold', 'spike'}); % add a fancy latex title title(['\fontsize{14} Filtering and Re-thresholding Raw Data' char(10) ... '\fontsize{12} Filtered Single Unit Data'],'interpreter','tex'); % Alternatively, overlay entire snippets on time window. %npts = size(su.snips.Snip.data,2); %for ii = 1:size(su.snips.Snip.data, 1) % plot(su.snips.Snip.ts(ii)+((1:npts)-floor(npts/4))/su.snips.Snip.fs, su.snips.Snip.data(ii,:)*1e6, 'color', colors(ii,:)); %end
Plot the extracted snippets.
subplot(2,3,4) hold on; for ii = 1:size(su.snips.Snip.data, 1) plot(su.snips.Snip.data(ii,:)*1e6, 'color', colors(ii,:)); end axis tight; xlabel('samples') ylabel('uV') title(sprintf('Extracted snippets, N=%d', ii), 'FontSize', 12)
Do Min/Max k-Means and PCA analysis
Plot snippets in a simple feature space.
subplot(2,3,5) scatter(minvals, maxvals, 8, colors, 'filled') xlabel('Min Value (uV)') ylabel('Max Value (uV)') title('Min/Max k-means', 'FontSize', 12) hold on;
Do k-means for 3 clusters.
rand('seed', 1); % for reproducibility [idx, C] = kmeans([minvals, maxvals], 3);
Plot centroids and cluster circles
for ii = 1:max(idx) % plot centroid scatter(C(ii,1), C(ii,2), 30, '*k'); % find spike with max distance from center to use as radius for circle x_center = C(ii, 1); y_center = C(ii, 2); curr_mins = minvals(idx==ii); curr_maxs = maxvals(idx==ii); dist = sqrt((curr_mins - x_center).^2 + (curr_maxs - y_center).^2); r = max(dist); % plot circle th = 0:pi/50:2*pi; xunit = r * cos(th) + x_center; yunit = r * sin(th) + y_center; plot(xunit, yunit, '--'); end legend({'spikes','centroids'})
Plot snippets in PCA space.
subplot(2,3,6) hold on; try [coeff,score,latent] = pca(su.snips.Snip.data); % newer Matlab uses this catch [coeff,score,latent] = princomp(su.snips.Snip.data); end [idx3, C3] = kmeans(score(:,1:3), 3); for ii = 1:max(idx3) % get principle components and color for this cluster pca1 = score(idx3==ii,1); pca2 = score(idx3==ii,2); pca3 = score(idx3==ii,3); color = colors(idx3==ii,:); % plot dots scatter3(pca1, pca2, pca3, 15, color, 'filled') % plot centroid %scatter3(C3(ii,1), C3(ii,2), C3(ii,3), 75, 'k', 'filled'); end axis tight; view([-15 25]); % set initial view on 3-D plot title('PCA feature space', 'FontSize', 12)