function Data = spikePlot(data, Fs)
%% spikePlot plots the given data and provides an interactive window discrimination function
% two horizontal cursors can be dragged with the mouse pointer to discriminate
% the spikes, which will be displayed in the upper part of the figure.
% The signal can be flipped on the y-axis around 0 [V, mV, ...]. Two vertical
% cursors can be added to select a range of spikes.
%
% Function parameters: The waveform 'data' sampled with the frequency
% 'Fs'
%
% Output value: Data, a struct containing all the information of the
% waveform and the calculated values, like thresholds, times of
% spike occurrence, interspike intervals and the reciprocal instantanious rate.
% Calculations are executed in the nested functions
% Data = calcSpikeTimes(Data) and Data = calcFreq(Data) which are executed after
% pressing the 'Update' button.
% Here is a good place to add basic custom calculations.
% In the callback 'buttonYour_callback' of the yourButton 'func' you
% can call your own functions.
%
% 'Data' can be assigned to the base workspace, so it is available
% for further calculations
%
% Example: generate noisy signal (from Matlab Help: FFT)
% Fs = 1000; % Sampling frequency
% T = 1/Fs; % Sample time
% L = 1000; % Length of signal
% t = (0:L-1)*T; % Time vector
% x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
% y = x + 2*randn(size(t)); % Sinusoids plus noise
%
% Data = spikePlot(y(1:50), Fs)
%
% Kai Voges
% kai.voges@gmx.net
% 09.10.2009
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add your own function calls here
function buttonYour_callback(src, evnt) %#ok<INUSD,INUSD>
if ~isfield(Data.user, 'freq')
buttonUpdate_callback
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here you can call your own functions for fast evaluations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Calculate values
%if ~isfield(Data.user, 'times')
Data.user.times.sampleFreq = Fs; % Sampling frequency [Hz]
Data.user.times.timeVector = 0:1/Data.user.times.sampleFreq:length(data)/Fs; % Time points for recording [s]
Data.user.times.timeVector = Data.user.times.timeVector(1:end-1);
%end
% Look for stored data
%if ~isfield(Data.user, 'threshold')
% if there is no data stored use these default values
% window discrimination (high and low thresholds)
Data.user.threshold.high = 0.02;
Data.user.threshold.low = 0.01;
Data.user.threshold.reversed = 'normal';
Data.user.threshold.cut = 'off';
%end
Data.chan1 = data;
hThreshold = Data.user.threshold.high;
lThreshold = Data.user.threshold.low;
reversed = Data.user.threshold.reversed;
cut = Data.user.threshold.cut;
% Plot commands
h_figure = figure;
scnSize = get(0,'ScreenSize'); % screen size [pixel]
figPos = [0 0 scnSize(3) scnSize(4) * 0.45];
set(h_figure, 'Position', figPos, 'Toolbar', 'figure', 'WindowButtonMotionFcn',@mouse_motion_figure);
% Subplot for raster diagram
h_subplot(1) = subplot(2,1,1);
plotRaster('initial') % call the plot function for the raster the first time
% Subplot for original recording
h_axes = axes('Drawmode','fast');
h_subplot(2) = subplot(2,1,2, h_axes);
newplot(h_subplot(2))
plot(Data.user.times.timeVector, Data.chan1) % plot original data
xlim([0 Data.user.times.timeVector(end)])
hold on
% Plot thresholds @subbplot 2
h_hLine = line('XData', [0 Data.user.times.timeVector(end)], 'YData', [hThreshold hThreshold],...
'Color', 'r','LineStyle','--','ButtonDownFcn',@button_down_high); % plot high threshold
h_lLine = line('XData', [0 Data.user.times.timeVector(end)], 'YData', [lThreshold lThreshold],...
'Color', 'g','LineStyle','--','ButtonDownFcn',@button_down_low); % plot low threshold
set(h_subplot(2), 'Position',[0.12 0.1 0.85 0.75])
xlabel('Time [sec]')
ylabel('Amplitude [V]')
hold off
linkaxes(h_subplot, 'x')
%% Ui controlls
flipToggle = uicontrol('Style','Togglebutton',...
'Position',[10 140 65 20],...
'Callback', @toggleFlip_callback,...
'String','Flip');
cutToggle = uicontrol('Style','Togglebutton',...
'Position',[10 110 65 20],...
'Callback', @toggleCut_callback,...
'String','Cut');
updateButton = uicontrol('Style','Pushbutton',...
'Position',[10 80 65 20],...
'Callback', @buttonUpdate_callback,...
'String','Update');
yourButton = uicontrol('Style','Pushbutton',...
'Position',[10 180 65 20],...
'Callback', @buttonYour_callback,...
'String','func'); %#ok<NASGU>
assignBox = uicontrol('Style','Pushbutton',...
'Position',[10 40 65 20],...
'Callback', @boxAssign_callback,...
'String','Assign in'); %#ok<NASGU>
% set value for y axis orientation
if strcmp(reversed, 'normal')
set(flipToggle, 'Value', 0)
set(gca,'YDir','normal')
else
set(flipToggle, 'Value', 1)
set(gca,'YDir','reverse')
end
% set value for cut or not
if strcmp(cut, 'off')
set(cutToggle, 'Value', 0)
else
set(cutToggle, 'Value', 1)
toggleCut_callback
end
colorButton = get(cutToggle,'BackgroundColor');
set(updateButton, 'BackgroundColor',colorButton)
%% Callback functions
% graphics controll for toggle buttons 'Flip' and 'Cut'
function toggleFlip_callback(varargin)
state = get(flipToggle, 'Value');
% if flipping: reverse y axis, updata 'Data', reset thresholds
if state
set(h_subplot(2),'YDir','reverse')
Data.user.threshold.reversed = 'reverse';
set(h_hLine,'YData',[-0.02 -0.02]);
set(h_lLine,'YData',[-0.01 -0.01]);
Data.user.threshold.high = -0.02;
Data.user.threshold.low = -0.01;
elseif ~state
set(h_subplot(2),'YDir','normal')
Data.user.threshold.reversed = 'normal';
set(h_hLine,'YData',[0.02 0.02]);
set(h_lLine,'YData',[0.01 0.01]);
Data.user.threshold.high = 0.02;
Data.user.threshold.low = 0.01;
end
set(updateButton, 'BackgroundColor','r')
end
function toggleCut_callback(varargin)
state = get(cutToggle, 'Value');
if strcmp(cut, 'off')
% if there are no stored values, use these for the initial plot
yLab = get(h_axes, 'YTick');
xLab = get(h_axes, 'XTick');
yLab = [min(yLab) max(yLab)];
cut = [xLab(2) xLab(3)];
else
yLab = get(h_axes, 'YTick');
yLab = [min(yLab) max(yLab)];
end
if state
axes(h_subplot(2))
handles.h_clLine(1) = line('XData', [cut(1) cut(1)], 'YData', [yLab(1) yLab(2)],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_left); % plot left cut
handles.h_crLine(1) = line('XData', [cut(2) cut(2)], 'YData', [yLab(1) yLab(2)],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_right); % plot right cut
% need to save handles, because plot command is inside of function
% and handles aren't visible outside the function
setappdata(h_figure, 'handles', handles)
% this part is for the cut lines in the raster plot; these lines
% have their own handle !!
if ~isappdata(h_figure, 'handles_raster')
axes(h_subplot(1))
handles_raster.h_clLine = line('XData', [cut(1) cut(1)], 'YData', [0 1],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_left); % plot left cut in raster
handles_raster.h_crLine = line('XData', [cut(2) cut(2)], 'YData', [0 1],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_left); % plot right cut in raster
% need to save handles, because plot command is inside of function
% and handles aren't visible outside the function
setappdata(h_figure, 'handles_raster', handles_raster)
end
elseif ~state
% delete handles and set 'cut' to 'off' if button is released
handles = getappdata(h_figure, 'handles');
handles_raster = getappdata(h_figure, 'handles_raster');
delete(handles.h_crLine)
delete(handles.h_clLine)
delete(handles_raster.h_crLine)
delete(handles_raster.h_clLine)
rmappdata(h_figure, 'handles_raster')
rmappdata(h_figure, 'handles')
cut = 'off';
end
Data.user.threshold.cut = cut;
set(updateButton, 'BackgroundColor','r')
end
% graphics controll for threshold lines
function button_down_high(src,evnt)%#ok<INUSD,INUSD>
% src - the object that is the source of the event
% evnt - empty for this property
set(h_figure,'WindowButtonMotionFcn',@mouse_motion)
set(h_figure,'WindowButtonUpFcn',@button_up)
beep on
% if mouse click on line, get motion of pointer and redraw line
function mouse_motion(src,evnt)%#ok<INUSD,INUSD>
cp = get(h_axes,'CurrentPoint');
% constraint, that threshold can't cross zero; for both cases:
% 'YDir', 'normal' or 'reversed'
if (cp(1,2) <= 0 && strcmp(Data.user.threshold.reversed, 'normal'))
set(h_hLine,'YData',[0 0]);
beep
beep off
return
elseif (cp(1,2) >= 0 && strcmp(Data.user.threshold.reversed, 'reverse'))
set(h_hLine,'YData',[0 0]);
beep
beep off
return
end
ydat = [cp(1,2) cp(1,2)];
set(h_hLine,'YData',ydat);
drawnow
end
function button_up(src,evnt)%#ok<INUSD,INUSD>
if strcmp(get(src,'SelectionType'),'normal')
% reset the 'WindowButtonMotionFcn' to figures mouse_motion
% function
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
set(src,'WindowButtonUpFcn','')
cp = get(h_axes,'CurrentPoint');
% if pointer crossed zero, set threshold to 0; if not, set
% actual value
if (cp(1,2) <= 0 && strcmp(Data.user.threshold.reversed, 'normal'))
Data.user.threshold.high = 0;
elseif (cp(1,2) >= 0 && strcmp(Data.user.threshold.reversed, 'reverse'))
Data.user.threshold.low = 0;
else
Data.user.threshold.high = cp(1,2);
end
set(updateButton, 'BackgroundColor','r')
else
return
end
end
end
function button_down_low(src,evnt)%#ok<INUSD,INUSD>
% src - the object that is the source of the event
% evnt - empty for this property
set(h_figure,'WindowButtonMotionFcn',@mouse_motion)
set(h_figure,'WindowButtonUpFcn',@button_up)
beep on
% if mouse click on line, get motion of pointer and redraw line
function mouse_motion(src,evnt)%#ok<INUSD,INUSD>
cp = get(h_axes,'CurrentPoint');
% constraint, that threshold can't cross zero; for both cases:
% 'YDir', 'normal' or 'reversed'
if (cp(1,2) <= 0 && strcmp(Data.user.threshold.reversed, 'normal'))
set(h_lLine,'YData',[0 0]);
beep
beep off
return
elseif (cp(1,2) >= 0 && strcmp(Data.user.threshold.reversed, 'reverse'))
set(h_lLine,'YData',[0 0]);
beep
beep off
return
end
ydat = [cp(1,2) cp(1,2)];
set(h_lLine,'YData',ydat);
drawnow
end
function button_up(src,evnt)%#ok<INUSD,INUSD>
if strcmp(get(src,'SelectionType'),'normal')
% reset the 'WindowButtonMotionFcn' to figures mouse_motion
% function
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
set(src,'WindowButtonUpFcn','')
cp = get(h_axes,'CurrentPoint');
% if pointer crossed zero, set threshold to 0; if not, set
% actual value
if (cp(1,2) <= 0 && strcmp(Data.user.threshold.reversed, 'normal'))
Data.user.threshold.low = 0;
elseif (cp(1,2) >= 0 && strcmp(Data.user.threshold.reversed, 'reverse'))
Data.user.threshold.low = 0;
else
Data.user.threshold.low = cp(1,2);
end
set(updateButton, 'BackgroundColor','r')
else
return
end
end
end
% graphics controll for cut lines
function button_down_left(src,evnt)%#ok<INUSD,INUSD>
% src - the object that is the source of the event
% evnt - empty for this property
set(h_figure,'WindowButtonMotionFcn',@mouse_motion)
set(h_figure,'WindowButtonUpFcn',@button_up)
handles = getappdata(h_figure, 'handles');
handles_raster = getappdata(h_figure, 'handles_raster');
% if mouse click on line, get motion of pointer and redraw line
function mouse_motion(src,evnt)%#ok<INUSD,INUSD>
cp = get(h_axes,'CurrentPoint');
xdat = [cp(1,1) cp(1,1)];
set(handles.h_clLine,'XData',xdat);
set(handles_raster.h_clLine,'XData',xdat);
drawnow
end
function button_up(src,evnt)%#ok<INUSD,INUSD>
if strcmp(get(src,'SelectionType'),'normal')
% reset the 'WindowButtonMotionFcn' to figures mouse_motion
% function
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
set(src,'WindowButtonUpFcn','')
cp = get(h_axes,'CurrentPoint');
cut(1) = cp(1,1);
Data.user.threshold.cut = cut;
set(updateButton, 'BackgroundColor','r')
else
return
end
end
end
function button_down_right(src,evnt)%#ok<INUSD,INUSD>
% src - the object that is the source of the event
% evnt - empty for this property
set(h_figure,'WindowButtonMotionFcn',@mouse_motion)
set(h_figure,'WindowButtonUpFcn',@button_up)
handles = getappdata(h_figure, 'handles');
handles_raster = getappdata(h_figure, 'handles_raster');
% if mouse click on line, get motion of pointer and redraw line
function mouse_motion(src,evnt)%#ok<INUSD,INUSD>
cp = get(h_axes,'CurrentPoint');
xdat = [cp(1,1) cp(1,1)];
set(handles.h_crLine,'XData',xdat);
set(handles_raster.h_crLine,'XData',xdat);
drawnow
end
function button_up(src,evnt)%#ok<INUSD,INUSD>
if strcmp(get(src,'SelectionType'),'normal')
% reset the 'WindowButtonMotionFcn' to figures mouse_motion
% function
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
set(src,'WindowButtonUpFcn','')
cp = get(h_axes,'CurrentPoint');
cut(2) = cp(1,1);
Data.user.threshold.cut = cut;
set(updateButton, 'BackgroundColor','r')
else
return
end
end
end
% graphics controll for mouse pointer
function mouse_motion_figure(src,evnt)%#ok<INUSD,INUSD>
ptrOff = 0.01; % offset for mouse pointer sensitivity
ptrOffc = 0.1; % offset for mouse pointer sensitivity
cp = get(h_axes,'CurrentPoint');
lThreshold = Data.user.threshold.low;
hThreshold = Data.user.threshold.high;
% if pointer is over threshold +/- offset and change pointer to crosshair
if ((lThreshold-ptrOff < cp(1,2) && cp(1,2) < lThreshold+ptrOff) ||...
hThreshold-ptrOff < cp(1,2) && cp(1,2) < hThreshold+ptrOff)
set(h_figure, 'Pointer', 'crosshair')
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
return
else
set(h_figure, 'Pointer', 'arrow')
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
end
% if there are 'cut' lines, check if pointer is over line +/-
% offset and change pointer to crosshair
if ~strcmp(cut, 'off')
if ((cut(1)-ptrOffc < cp(1,1) && cp(1,1) < cut(1)+ptrOffc) ||...
cut(2)-ptrOffc < cp(1,1) && cp(1,1) < cut(2)+ptrOffc)
set(h_figure, 'Pointer', 'crosshair')
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
return
else
set(h_figure, 'Pointer', 'arrow')
set(src,'WindowButtonMotionFcn',@mouse_motion_figure)
end
end
end
% graphics controll for push button 'Update'
function buttonUpdate_callback(src,evnt)%#ok<INUSD,INUSD>
Data = calcSpikeTimes(Data);
Data = calcFreq(Data);
plotRaster('update')
% addSaveData(Data.user.file, Data.user)
if ~strcmp(cut, 'off')
axes(h_subplot(1))
handles_raster.h_clLine = line('XData', [cut(1) cut(1)], 'YData', [0 1],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_left);
handles_raster.h_crLine = line('XData', [cut(2) cut(2)], 'YData', [0 1],...
'Color', 'k','LineStyle','--','ButtonDownFcn',@button_down_left);
% need to save handles, because plot command is inside of function
% and handles aren't visible outside the function
setappdata(h_figure, 'handles_raster', handles_raster)
end
%user = Data.user;
%fileLocation = [Data.user.file.pathName, Data.user.file.fileName];
%save(fileLocation,'-append', 'user')
colorButton = get(cutToggle,'BackgroundColor');
set(updateButton, 'BackgroundColor',colorButton)
end
% graphics controll for rasterplot
function plotRaster(str)
if strcmp(str, 'initial')
xlimits = [0 Data.user.times.timeVector(end)];
else
xlimits = get(h_subplot(2), 'XLim');
end
if isfield(Data.user, 'spikes') && isfield(Data.user.spikes, 'spikeTime') && ~isempty(Data.user.spikes.spikeTime)
if length(Data.user.spikes.spikeTime)<=2 % this is necessary to prevent connection of 2 datapoints
for i=1:length(Data.user.spikes.spikeTime)
plot(h_subplot(1), [Data.user.spikes.spikeTime(i) Data.user.spikes.spikeTime(i)], [0 1], 'k');
end
else % this should be the usual case
plot(h_subplot(1), [Data.user.spikes.spikeTime Data.user.spikes.spikeTime], [0 1], 'k');
end
elseif ~isfield(Data.user, 'spikes') || ~isfield(Data.user.spikes, 'spikeTime') || isempty(Data.user.spikes.spikeTime)
cla(h_subplot(1))
end
%title(h_subplot(1), ['File: ', Data.user.file.fileName],'Interpreter', 'None', 'FontWeight','Bold')
xlim(xlimits)
set(h_subplot(1), 'Position',[0.12 0.88 0.85 0.05], 'XTick', [], 'YTick', [], 'XTickLabel', '', 'YTickLabel','')
end
function boxAssign_callback(src, evnt)%#ok<INUSD,INUSD>
str = sprintf('Data_%i',Data.user.file.fileNum);
assignin('base', str, Data);
end
% calculating spiketimes and frequencies
function Data = calcSpikeTimes(Data)
pattern = [1 -1];
% check if 'YDir' is reversed, find spikes
if strcmp(Data.user.threshold.reversed, 'normal')
% check if low is 'lower' than 'high', if not exchange values
if (Data.user.threshold.high > Data.user.threshold.low)
hThresh = Data.user.threshold.high;
lThresh = Data.user.threshold.low;
else
lThresh = Data.user.threshold.high;
hThresh = Data.user.threshold.low;
end
lSpikes = Data.chan1 > lThresh; % all data points bigger than lower threshold; logical
changes = conv(pattern, lSpikes); % find begin and end of each spike; logical
% convert logical into original datapoints and get max
spikeLength = arrayfun(@(x,y) Data.chan1(x: y), find(changes==1), find(changes==-1),'uni',false);
spikeMax = cellfun(@max, spikeLength);
% find spikes smaller than border
spikesAll = [find(changes==1), find(changes==-1)];
spikes = spikesAll(find(spikeMax < hThresh),1); %#ok<FNDSB>
Data.user.spikes.spikeTime = spikes / Data.user.times.sampleFreq;
else
% check if low is 'lower' than 'high', if not exchange values
if (Data.user.threshold.high < Data.user.threshold.low)
hThresh = Data.user.threshold.high;
lThresh = Data.user.threshold.low;
else
lThresh = Data.user.threshold.high;
hThresh = Data.user.threshold.low;
end
lSpikes = Data.chan1 < lThresh; % all data points bigger than lower threshold; logical
changes = conv(pattern, lSpikes); % find begin and end of each spike; logical
% convert logical into original datapoints and get max
spikeLength = arrayfun(@(x,y) Data.chan1(x: y), find(changes==1), find(changes==-1),'uni',false);
spikeMin = cellfun(@min, spikeLength);
% find spikes smaller than border
spikesAll = [find(changes==1), find(changes==-1)];
spikes = spikesAll(find(spikeMin > hThresh),1); %#ok<FNDSB>
Data.user.spikes.spikeTime = spikes / Data.user.times.sampleFreq;
end
end
function Data = calcFreq(Data)
%% Instantanious frequency
cut = Data.user.threshold.cut;
spikes = Data.user.spikes.spikeTime;
if ~strcmp(cut, 'off')
cutMax = max(cut); % bigger cut value [sec]
cutMin = min(cut); % bigger cut value [sec]
spikesIndex = find(spikes > cutMin);
spikes = spikes(spikesIndex); %#ok<FNDSB>
spikesIndex = find(spikes < cutMax);
spikes = spikes(spikesIndex); %#ok<FNDSB> % time spikes occured after cutting [sec]
end
interSpike = zeros(length(spikes) - 1, 1);
for i = 2:length(spikes)
interSpike(i-1) = spikes(i) - spikes(i-1); % time between two succesiv spikes [sec]
end
instantFreq = 1./interSpike; % instantaneous frequency = reciprocal of interspike time
numFreq = length(instantFreq);
meanInstFreq = mean(instantFreq);
firingVariability = sqrt((sum((instantFreq - meanInstFreq).^2))/numFreq);
Data.user.freq.instantFreq = instantFreq;
Data.user.freq.meanInstFreq = meanInstFreq;
Data.user.freq.firingVariability = firingVariability;
end
end
