# Event Classification Processor, Version 3.7, 02-MAY-2014
#
# To find spikes, spike bursts, and head shakes in mouse EEG with
# screw electrodes. Baseline power is dominated by white noise from
# the electronics, not the pink EEG noise. Thus we have event metric
# threshold at 3 x baseline power rather than our usual 5 x. This
# processor allows us to turn off metric calculation if we want,
# and baseline calibration as well. The metrics are event power,
# high-frequency power, asymmetry, spikiness, and intermittency. Our
# transient metric represents glitch power above the anti-aliasing
# filter cut-off, but does not represent low-frequency steps.
#
# Set the following variable to 1 to calculate event classification metrics.
set metric_calculation 1
# Set the following vareiable to 1 to calculate baseline power. If 0, we use
# the Neuroarchiver's existing baseline power value, which we assume has been
# read from the archive's metadata, or otherwise determined before-hand.
set baseline_calibration 1
# We set up the Neuroclassifier.
set config(classifier_types) "Spike red \
Spike_burst lightblue \
Glitch orange"
set config(classifier_metrics) "event_pwr \
hf_pwr \
spikiness \
asymmetry \
intermittency \
transient_pwr"
# We begin this channel's section of the characteristics line with the
# channel number.
append result "$info(channel_num) "
# Define the event band in Hertz.
set event_lo 2.0
set event_hi 160.0
# Define the limits for what we consider to be high-frequency power in Hertz.
set hf_lo 60.0
set hf_hi 160.0
# Define the limits for the intermittency envelope power in Hertz.
set int_lo 4.0
set int_hi 20.0
# Define the frequency band we use to detect the power of glitches that make
# it through our glitch filter.
set glitch_lo 200.0
set glitch_hi 250.0
# If we have too much signal loss, we ignore this interval.
set max_loss 20.0
# Calculate event power, which we always need.
set event_pwr [expr 0.001*[Neuroarchiver_band_power $event_lo $event_hi 0 0]]
# Set up the baseline_pwr variable to point to the Neuroarchiver's
# baseline power for the current channel.
upvar #0 Neuroarchiver_info(bp_$info(channel_num)) baseline_pwr
# Calculate baseline power if directed to do so. We set it to the
# minimum event power. When it appears that we have a new minimum
# power, we re-calculate the event power without the glitch threshold,
# in order to make sure that we are not looking at an interval that has
# been corrupted by an artifact of glitch filtering.
if $baseline_calibration {
if {($event_pwr < $baseline_pwr) && ($info(loss) < $max_loss) && ($event_pwr > 0)} {
set info(values) [Neuroarchiver_values $info(signal)]
set info(spectrum) [lwdaq_fft $info(values) -glitch 0 \
-window [expr round([llength $info(values)] * $config(window_fraction))]]
set pwr [expr 0.001*[Neuroarchiver_band_power $event_lo $event_hi 0 0]]
if {$pwr < $baseline_pwr} {set baseline_pwr [format %.3f $event_pwr]}
}
}
# Set the default values of the metrics.
foreach a "event_m hf_m spikiness_m asymmetry_m intermittency_m transient_m" {
set [set a] 0.0
}
# Calculate metrics if directed to do so.
if $metric_calculation {
# We choose our event band to include all activity that might contribute
# to the events we are looking for, but exclude both transient noise and
# very high frequency quantization noise. In order to calculate the spikiness
# of a signal, which we do below, we must look at the sample values of the
# signal itself rather than its spectrum. If we look at the unfiltered signal,
# we will find our calculation of spikiness is disturbed by transients. Here
# we instruct the band power routine to overwrite the info(values) string with
# the inverse transform of the event band spectrum.
set event_pwr [expr 0.001*[Neuroarchiver_band_power $event_lo $event_hi 0 1]]
# The spikiness of a signal is the ratio of its range to its standard deviation
# after filtering to the event frequency band.
scan [lwdaq ave_stdev $info(values)] %f%f%f%f ave stdev max min
if {$stdev > 0} {
set spikiness [expr ($max-$min)/($stdev)]
} {
set spikiness 1.0
}
# The asymmetry of a signal is a measure of how much it protrudes downwards
# or upwards from its average value. We construct a measure that is 0.5 whenprotrudes downwards, and greater
# than 0.5 for upwards.
set count 0
set upcount 0
set downcount 0
foreach v $info(values) {
if {($v - $ave) >= 2.0*$stdev} {incr upcount}
if {($ave - $v) >= 2.0*$stdev} {incr downcount}
if {abs($v - $ave) >= $stdev} {incr count}
}
if {$upcount > $downcount} {
set asymmetry [expr 0.5+1.0*abs($upcount-$downcount)/$count]
} {
set asymmetry [expr 0.5-1.0*abs($upcount-$downcount)/$count]
}
# Obtain the high-frequency power. We choose this band to include the
# power we see in high-frequency bursts, but to exclued the harmonics of
# eating and other rhythms. In order to calculate the intermittency
# of a signal, we overwrite the info(values) string with the inverse
# transform of the high-frequency band spectrum.
set hf_pwr [expr 0.001*[Neuroarchiver_band_power $hf_lo $hf_hi 0 1]]
# The intermittency of a signal is the low-frequency power of the rectified
# high-frequency signal. Thus we take the high-frequency signal, rectify it,
# take the transform, and calculate the power in the low-frequency band.
set rectified_values ""
foreach v $info(values) {
append rectified_values "[expr abs($v)] "
}
set spectrum [lwdaq_fft $rectified_values]
set f 0
set intermittency_pwr 0.0
foreach {a p} $spectrum {
if {$f >= $int_lo} {
set intermittency_pwr [expr $intermittency_pwr + 0.001*$a*$a]
}
set f [expr $f + $info(f_step)]
if {$f > $int_hi} {break}
}
# The transient power represents disturbances of the recording that
# would otherwise confuse our classification. Electrodes that are
# inadequately attached to the skull will give low-frequency transients.
# Poor solder joints will do the same. A bad message in the recording
# will usually result in a large glitch that we eliminate with our
# glitch filter, but sometimes the bad message value is close to our
# current signal value and produces a small, sharp spike. If we are
# having problems with such spikes, we include glitch power in the
# transient power measurement. Unlike artifact in the analog signal,
# glitches have power above the cut-off our our anti-aliasing filter.
set transient_pwr [expr 0.001*[Neuroarchiver_band_power $glitch_lo $glitch_hi 0 0]]
# We calculate the characteristics we use for determining similarity
# of events.
if {($info(loss) < $max_loss) && \
($event_pwr > 0)} {
set event_m [expr 1.0 / (1.0 + pow(3.0/($event_pwr/$baseline_pwr),1.0))]
set hf_m [expr 1.0 / (1.0 + pow(0.1/($hf_pwr/$event_pwr),1.0))]
set spikiness_m [expr 1.0 / (1.0 + pow(8.0/$spikiness,4.0))]
set asymmetry_m [expr 1.0 / (1.0 + pow(0.5/$asymmetry,2.0))]
set intermittency_m [expr 1.0 / (1.0 + pow(0.1/($intermittency_pwr/$hf_pwr),2.0))]
set transient_m [expr 1.0 / (1.0 + pow(0.5/($transient_pwr/$baseline_pwr),1.0))]
}
}
# We add a single-letter code that indicates the type of the interval, so
# far as this processor can deduce it. We have L for loss, G for glitch,
# U for unclassified and N for no event.
if {($info(loss) >= $max_loss) || ($event_pwr <= 0)} {
append result "L "
} elseif {$event_m>=0.5} {
append result "U "
} else {
append result "N "
}
# We add the baseline power to the characteristics. If we
# the signal is symmetric, less than 0.5 when it calculate
# the baseline power later, when looking through an archive, we could miss
# the defining low-power interval that sets the baseline.
append result "[format %.2f $baseline_pwr] "
# We record in the result line all the characteristics of the playback
# interval. The first one we record is the one we will use as a threshold
# for the occurance of an event, this threshold being 0.5 by default.
if $metric_calculation {
foreach a "event_m hf_m spikiness_m asymmetry_m intermittency_m transient_m" {
append result "[format %.3f [set $a]] "
}
}