Equalizing Sampling Rate for 1-D Correlations in MATLAB

We often have two signals with different sampling rates that we want to correlate. One example is when we want to correlate local field potential (LFP) activity of a neural signal to some average spike rate of a neuron. Our LFP signal is recorded at 24 kHz, but the spike density estimate (SDE) we create from individual spike timestamps is given an artificial sampling rate of 1 kHz. End-to-end, the LFP and SDE are potentially correlated, but to operate on those as variables in MATLAB they need to be the same length.

In this function, we use 1-D interpolation to force d1 data into the length of d2 and call that d1new. Ideally, your d2 data has the higher sampling rate so that d1 is being ‘upsampled’ instead of ‘downsampled’.

Consider two signals (d1 and d2) where f is either sin() or cos() and the time variable t is 0-100s, but each t has a different sampling rate.

t1 = linspace(0,100,100); % 0-100s, 100 points, 'low' sampling rate
t2 = linspace(0,100,400); % 0-100s, 400 points, 'high' sampling rate
d1 = sin(t1);
d2 = cos(t2);
d1new = equalVectors(d1,d2);

The top plot shows the signals overlaid in Time with each point from t marked in black. The middle plot highlights that when plotted by Samples, the two signals are unequal in size. The bottom plot shows how 1-D interpolation through equalVectors and the d1new variable equalize the vector lengths, thereby making a correlation possible.


Infinite Highlight Color Marks in Ruby (modulo based)


Consider the case where you have many sentences pieced together from different sources and want to highlight each one of those sentences based on the source number. In HTML, the standard for highlighting text utilizes the <mark> tag with the background-color attribute set to a light yellow color.

My background color is set using style="background-color: hsl(50,90%,90%);"

The first value in hsl() is a value in degrees for hue which ranges from 0 to 360, spanning from red to green to blue and back to red in circular form. The goal of the following function (i.e., line of code) is to return maximally different hues based on the total amount of colors you want to use.

If the first color is a standard yellow (used above) then startHue is 50. nColors is the total number of colors you want to use and ii is the current color you want to select. Consider where nColors = 6:

This is from the first source (where ii = 1).1 The second source has a different color.2 So does the third!3 All of the colors are meant to contrast.4 Unless nColors is very large.5 This is the where ii = 6; that's all of them.6

Pro Tip: If you want to selectively highlight text based on a word or regex pattern in Ruby on Rails there is a great Text Helper called highlight.

MATLAB Line Colors for The University of Michigan

One maize. One blue. One brand. This function makes the University of Michigan colors available in an array (just like how lines() works in MATLAB). To get all the primary colors:

colors = linesUM;

To limit the amount of colors returned (for instance, to use the colors as a colormap):

colors = linesUM(3);

Finally, to get all the secondary colors just pass in Inf for the amount of colors and set onlySecondary to true:

colors = linesUM(Inf,true);

A Better Way to Operate: The Surgical Placemat

In early 2018, I met an expert on checklists—yes, they exist. She works closely with NASA to develop checklists that keep astronauts safe and productive. I asked her, "what is one way to use a checklist that I wouldn't have thought of?" After all, I bought into the fact that someone could be good at making checklists, but I wanted to know what else I was missing. She asked me to give her a situation or problem that might lend itself to a checklist. I said, "surgery," followed with, "we have protocols, and logs, but we don't use a checklist and it burns us a lot." Not that our animals were dying, but I've seen the antiseptic skipped over, lidocaine forgotten, and I have showed up to empty oxygen bottles because the valve wasn't closed at the end of the day.

My new checklist expert friend went on to ask, "does each step have an item, or tool associated with it?". "Yes," I followed, "just about all of them." She then said, "well here's one thing to think about doing," as if this was one of thirty tricks up her sleeve, "you could make a large placemat with a bunch of pouches, and each pouch is preloaded with your tool and the order of pouches walks you through the surgery." How simple, I thought; and brilliant. For the rest of the afternoon all I could think about was sewing custom, ordered pouches into canvas for everything I do in a day.

Atul Gwande's, The Checklist Manifesto, is a very readable exposé on the utility of checklists not only in surgery, but almost all personal and professional endeavors. Checklists assist pilots in starting a plane as well as any and all emergency situations in the air. "FLY THE PLANE" is often number one on the emergency checklist, testament to how simple and direct a good checklist should be. A checklist can help an investor stay true to a method and keep them from letting their gut get the best of them (and their clients). However, what was missing from the book is what my checklist expert honed in on: tactics to increase compliance.

A checklist is only worth its weight if it's followed; step by step, line by line. That's where a checklist—a good checklist—beats a protocol. Protocols are developed as a "Do-Confirm" style list that is better at informing an auditor to the procedure than being something an operator can efficiently utilize. The "Read-Do" style checklist is more useful, more tactful, and built for the experienced operator. Still, the Read-Do method suffers from compliance, as it is easy for anyone with experience to think they know the next three steps, skip ahead, and ink the check-check-check after skipping ahead.

That's where the physical checklist comes into play. The next step in the list is not possible or revealed until a real-world action satisfying the previous step has taken place. For me, the surgical placemat was a perfect solution. It outlines all the items and tools we need in surgery and sits directly on the operating table. It can be preloaded before the operation and thrown out at the end. Although we perform many types of operations in our lab, the setup and teardown is mostly the same.

This is just a prototype—a start—as it is a constant challenge to balance too much with too little, and succinctly separate what needs to be checked-off from what simply needs to be performed. But it seems like the right step forward.


When Movement Breaks - Lecture & Commentary

When Movement Breaks - Lecture & Commentary

Central Nervous System (CNS) Aspects of Motor Control II

In this lecture I review how the basal ganglia, thalamus and cerebellum work with the motor cortex to influence, and produce movement. I use historic cases of motor deficits and disease, including Parkinson's disease, to reflect on how knowledge has progressed. Included are discussions on how the MPTP model of Parkinson's disease was founded, the standard "rate model" developed by Albin et al., and modern theories concerning firing rates and firing patterns than are employed by the thalamocortical system.

Simpler Subplots: Converting Grid Coordinates to Axes Position

When you're working with a subplot matrix where each row is some new analysis, using the standard "position" or p input to MATLAB subplot is not intuitive. Rather, you want to specify your subplot by rows and columns. This tiny function, prc, does that conversion.

The function requires the number of columns in your figure (cols) and the two-element position of your desired subplot (rc) given as an array in the format, [row, col].



The figure above was generated with the following code.

Making your own colormap (cmap) in MATLAB

This function will use the first row of an image and allow you to use that as a colormap, or as individual colors for your own application. It just extracts colors and puts them into an array. For instance, I designed this "traffic light" color scheme and want to apply it to a plot where each color represents some reaction time.


You can use any image, but here is how to design the gradient in Adobe Photoshop:

  1. Open a new document with a width greater or equal to the amount of elements you will be coloring in MATLAB (for instance, width: 1500 pixels, height: 200 pixels). The height is only so you can visualize the color scheme; the function will only require a single row.

  2. Use the Shape tool to draw a rectangle the size of the canvas.

  3. Double click the rectangle layer and apply a Gradient Overlay. Use the color tools to make your gradient.

  4. Save your document as a Photoshop file first if you want to archive the original, but then export it to a JPG of the original size either using Save As, or Save for Web.

Screen Shot 2018-03-09 at 9.09.24 AM.png

The mycmap function has one mandatory input: the image filename. There is one optional input: the number of elements to be returned. For instance, if you are using the colors to set cmap, you probably don't need to set the number of elements (it will just return a colors array the width of your image). On the other hand, if you are using the colors to set the color of say, 55 different lines on a plot, you will pass that number as the second input to the function, and the returned array colors will have 55 rows.


Below are phase plots from an electrophysiology experiment where each line represents a single trial. Each line has been colored based on reaction time (low-to-high) using the colors from the mycmap function.


Finding Consecutive Numbers that Exceed a Threshold in MATLAB


The goal is to determine if there exists a span of n p-values that exceed threshold t

pVals = [0.91 0.66 0.23 0.96 0.99 0.97 0.99 0.83 0.12 0.96 0.99 0.97 0.76];

n = 3;

t = 0.95;

The total amount of values exceeding the threshold can be found simply:

sum(pVals > t);

However, we are also requiring that those p-values be consecutive. We can using a moving sum centered on the current element plus n on the binary index of p-value threshold crossings to do this:

ntpIdx = movsum(pVals > t,[0 n-1]) == n;

ntpIdx =

  1×13 logical array

   0   0   0   1   1   0   0   0   0   1   0   0   0

If you only want the initial crossing indexes:

ntpIdx_init = logical(diff([0 ntpIdx]) == 1);

ntpIdx_init =

  1×13 logical array

   0   0   0   1   0   0   0   0   0   1   0   0   0

Characterizing Tremor from Video using Frequency Analysis in MATLAB

Previously I shared a function useful for Creating an Actogram from Video in MATLAB. However, in cases where movement is rhythmic, like in cases of Parkinson's disease and essential tremor, understanding movement in the frequency domain is helpful. Using similar video analysis principles (assessing the change in pixel values from frame-to-frame) extracting frequency information is rather straight forward.

[allFrames,pos] = videoFreqAnalysis(videoFile,resizeFactor,ROItimestamp,freqList);

This method allows you to analyze a video for the frequencies provided in freqList given that the maximum frequency is less than the video frame rate divided by two (Shannon's sampling theorem).

Find these functions on Github: MoveAlgorithms


Exporting your original video to a smaller size will speed processing. This sample video of my hand is 15 seconds, taken with an iPhone, then exported as 540p with no audio.

The video processing follows these simple steps:

  1. Prompt the user to select a region of interest (ROI). The resulting ROI is used for analysis; anything outside of the selected ROI is not used. This helps focus the analysis and also speeds the analysis.
  2. Process all video frames. Every frame is converted to black and white, slightly contrasted, and then put into the matrix allFrames.
  3. Threshold the data so only pixels that show large movement (in the 80th percentile) are analyzed.
  4. Convert time domain data in the frequency domain using a complex scalogram (tip: we also use this with low-pass filtered electrophysiology data).
  Top:  Based on the  fpass  and resulting  freqList , a scalogram heatmap is presented that shows the power of each frequency band over the duration of the video.   Bottom:  This plot removes all temporal information and represents the mean frequency transform of the scalogram, useful for finding peaks in the frequency content of the video.

Top: Based on the fpass and resulting freqList, a scalogram heatmap is presented that shows the power of each frequency band over the duration of the video.

Bottom: This plot removes all temporal information and represents the mean frequency transform of the scalogram, useful for finding peaks in the frequency content of the video.

Together, the tremor from the video is present from the 5-10 second mark and peaks in the ~3 Hz band. The one major pitfall of this method is the lack of a standard amplitude measure. Depending on the video resolution and lighting conditions the amplitude will change. Controlling for these variables would make for a more meaningful comparisons between videos or subjects.

Frame-by-frame Heatmap

One visualization of the frequency domain is to take the pixels from the ROI and analyze them over a ±1 second window, where the red colors represent greater amplitude of the fpass band. The ROI is very large in the video below, but shows where movement is identified using this method.

Creating an Actogram from Video in MATLAB

An actogram is a way of quantifying movement, useful in behavioral studies in order to understand how an animal moves over the duration of an experiment. An actogram is usually based on movement sensors (e.g., IR beams), but because video is often being recorded anyways, it's a natural method of extracting movement information. The method I present here produces actogram data based on the changes in pixel values from frame-to-frame:

frameData = videoActogram(videoFile,frameInterval,resizePx);

For large videos you may also use a frameInterval to skip frames between each data point. The result is a frameData array consisting of frame number, frame timestamp, and movement data. More details are included within the file comments.

Find these functions on Github: MoveAlgorithms


Exporting your original video to a smaller size will speed processing. This sample video of my hand is 15 seconds, taken with an iPhone, then exported as 540p with no audio. It takes about 10 seconds to get the actogram data using every frame.

Using the test_videoActogram.m function you can easily plot the actogram data:


For a full view, you can use the overlayActogram.m function to create a video with a running actogram below the original video. Note, if you used a frameInterval your resulting video will only contain those frames (i.e. your video will be abbreviated).

Using Compose in MATLAB for Pretty Tick Labels

Creating descriptive and well formatted text labels for x- or y-ticks in MATLAB is essential to the presentation of your data. I have always struggled to remember the best way to replace tick marks and default numerical labels with text labels. I can not claim that the compose function in MATLAB is hidden, but I don't see it used often and it performs precisely what is required in many of cases.

For example, you have sampled a reaction time (RT) distribution using 10 quantiles. In RT experiments the RT histogram will generally have some type of long-tailed distribution, therefore each quantile contains a different range of RTs.


Let's say you want to compare some neuronal firing rate associated with those quantiles of RTs.

% RTs = reaction times
% Zscores = Z score of neuronal firing associated with RT
RTs = [0.10758143999999 0.113643519999982 0.119705600000088 0.12435456000005 0.127467520000003 0.131112959999996 0.134963199999959 0.138751999999897 0.14278655999999 0.147066880000125 0.151367680000021 0.155709439999896 0.160829439999929 0.166154239999969 0.173260799999866 0.184340479999889 0.197160960000019 0.215386720000026 0.24209264000001 0.293253120000031];
Zscores = [1.91375996684532 2.11512451719412 2.31635992789693 1.97677278262753 1.83302457695781 1.76635184611673 2.58036012006912 1.40779739402557 2.29017317616969 1.44587687689824 1.76950209157802 1.60579505116933 1.69767216641066 1.87616271740106 1.75912513938938 1.41276581199477 1.23089211422587 1.03332795801901 0.863934150049046 0.683887429578394];

colors = cool(numel(Zscores));
markerSize = 50;
for iiRT = 1:numel(RTs)
    lns(iiRT) = plot(RTs(iiRT),Zscores(iiRT),'.','markerSize',markerSize,'color',colors(iiRT,:));
    hold on;

ylabel('Z score');
title('Z score vs. RT');

However, you may not want to all of those data points to cluster where the quantile spacing is tight. Therefore, you can plot the x-axis on a linear scale, but now your x-ticks do not represent RTs, rather, they are a linear scale.

colors = cool(numel(Zscores));
markerSize = 50;
for iiRT = 1:numel(RTs)
    lns(iiRT) = plot(iiRT,Zscores(iiRT),'.','markerSize',markerSize,'color',colors(iiRT,:));
    hold on;
ylabel('Z score');
title('Z score vs. RT');

To correct this situation, you would want to replace those x-ticks with RT values. The primary issue when it comes to data presentation is formatting those 'floats' into something manageable and readable. Although it appears that xtickformat would accomplish this, it doesn't appear to retain values set by xticklabels. In addition, using compose gives you access to a cell string that can be used in a legend.

RTs_labels = compose('%1.3f',RTs);

* As a final note, although scatter would work instead of a FOR loop, it does not return independent handles for each point, so a colorized legend is not as easy to achieve.

Rat Brain Atlas API for MATLAB

I've made a Rat Brain Atlas API available on MATLAB Central File Exchange (linked to the GitHub Repo) based on my Rat Brain Atlas Webtool. The API wraps a web request to a JSON-formatted API endpoint, accessible in raw form using the following URL structure:


The function ratBrainAtlas(ap,ml,dv) reads the web request, downloads the atlas images, and outputs a struct in S. If the Image Processing Toolbox is installed S will include marked images (red dot).

A Dark MATLAB Theme for the University of Michigan

This theme—Dark Mich—is based on the University of Michigan Style Guide and uses Scott Lowe's Schemer tool to load in MATLAB. The colors have been slightly modified to improve readability and contrast.


  1. Download Schemer and unzip it to your MATLAB directory.
  2. Download darkmich.prf to the Schemer schemes folder and then load it by running:
>> schemer_import

While you're at it, treat yourself to the free font Source Code Pro to beautify things a bit (change in Preferences > Fonts).

Pro Tips

  1. f.lux will automatically adjust your screen color at night and remove those pesky blues that have a devastating effect on your circadian rhythm.
  2. Buy a pair of J+S Vision Computer Glasses to filter out 90% of your ambient blue light. I was so impressed with their product (e.g. filter quality and price) that I reached out to J+S Vision and they now sponsor my neuroscience research.

Surgical Tool for Cranial Drill Alignment in Rodents

Patent Pending #62369132

This enables the surgeon to drill a hole that travels longitudinally along the lateral bone structures, thereby increasing surface area for a bone screw to purchase to the skull. Often, these bone screws are used along with some type of epoxy or cement to attach larger implantable devices. This method of bone screw placement improves on many other methods that only drill into the cranial plates which are very thin and preclude proper threading and purchase of a screw. Developing the surgical skills to properly traverse the cranial ridge may take a surgeon a very long time and failures in placement can result in brain damage or improper purchase of the bone screw/s, leading to an eventual failure of the implant, often resulting in the premature sacrificing of the animal. Therefore, proper bone screw placement is both a technical challenge, and furthermore, a concern for animal well-being.

Other methods for placing bone screws along the cranial ridge involve stereotactic tools that are cumbersome and time consuming to employ. Minimizing surgical complexity and the time an animal spends anesthetized are crucial factors that contribute to operational success and animal recovery. This surgical tool has been designed to quickly align to anatomical features of the rodent cranium and guides a drill and bit at the proper angle into the cranial ridge bone. The dimensions of the tool are such that a surgeon may place multiple screws on each side of the skull, following the cranial ridge (circled) from the anterior/frontal aspect (F), to the posterior/parietal aspect (P). The cranial ridge described is conserved in several rodent species and is thick enough, in the described areas, to drill holes that do not enter the brain cavity on the medial aspect, or protrude out of the skull on the lateral aspect. The “angle of attack” for the drill bit has been determined through consultation of anatomical drawings as well as trial in the operating room with anesthetized animals.

Finding Nearest Value in MATLAB Using min()

I'm often converting between samples or video frames and time. For example, I'll have some video of animal behavior and an electrophysiological recording that each run on their own clocks (e.g. 30 frames-per-second for video and 30,000 samples-per-second for ephys). Therefore, it makes a lot of sense to use some absolute unit (like seconds) in the software so frame or sampling rates can disappear. This does create some syncing issues, when say, you want to find a video frame that occurs at t = 2.3539 seconds based on ephys events and your frame rate only resolved frames at t = 2.3333 and t = 2.3666. Given an array of values and a target value, this function returns the nearest index within that array, and the value which is nearest the target value.

Let's use my example and create an array of thirty values between 2 and 3:

frameTimes = [2, 2.0345, 2.069, 2.1034, 2.1379, 2.1724, 2.2069, 2.2414, 2.2759, 2.3103, 2.3448, 2.3793, 2.4138, 2.4483, 2.4828, 2.5172, 2.5517, 2.5862, 2.6207, 2.6552, 2.6897, 2.7241, 2.7586, 2.7931, 2.8276, 2.8621, 2.8966, 2.931, 2.9655, 3];

And now we can use that stupid-simple function to find a nearest index and value:

[idx, val] = closest(frameTimes,2.3539);

Here's the output and a figure to show what's going on:

idx =
val =

Now, I would know that whatever important thing happened at 2.3539 seconds can probably be revealed by looking at frame 11, or 2.3448 seconds into my video.

The [Simpler] Tetrode Spinner

We wanted a simpler, off-the-shelf tetrode twister that could be bought, built, and assembled by anyone. We took some tips from the "big brother" Twister that was already on Open Ephys, and made a few feature cuts. The simple twister is based on Arduino microcontroller hooked to a stepper motor, and we hope that our parts list helps you build your own semi-automated twister for about $300.


  1. [$8.94] Hammond Manufacturing Clear Polycarbonate Enclosure
  2. [$4.95] Momentary Pushbutton Switch
  3. [$4.99] 10kOHm Resistors 1/4 Watt (100pcs)


  1. [1$85.00] Aluminum Breadboard 6"x12"x1/2" with 1/4"-20 Taps
  2. [$26.00] Sorbothane Feet 1/4"-20 Thread
  3. [$10.82] Ø1/2"x12" Optical Post 1/4"-20 Tap
  4. [$6.77]  Ø1/2"x6" Optical Post 1/4"-20 Tap
  5. [$12.64] Ø1/2"x6" Post Holder, Spring-Loaded Hex-Locking Thumbscrew
  6. [$9.76] Right-Angle Post Clamp Fixed 90° Adapter
  7. [$25.70] Slim Right-Angle Bracket with Counterbored & M6 Tapped Holes
  8. [$9.50] 1/4"-20 Stainless Steel Cap Screw, 1" Long, Pack of 25


  1. [$24.95] Arduino Uno R3
  2. [$19.95] Adafruit Motor/Stepper/Servo Shield for Arduino v2 Kit - v2.3
  3. [$14.00] Stepper Motor 200 steps/rev, 12V 350mA
  4. [8.95] Stepper Motor Mount with Hardware


You mat need a few other parts, like hookup wire, a small needle, and some standoffs to mount the Arduino, but you might already have those. Let's get started with the Arduino.

  1. Head over to Adafruit's wonderful documentation on how to install the motor shield.
  2. Install a 10kOhm resistor between IO pin 12 and ground (see photo below), and install one end of the switch directly to pin 12, and the other side to 3.3V+.
  3. Download the Arduino Software.
  4. Clone our Tetrode Arduino repository on GitHub (or download the ZIP) and compile it onto the Arduino.
  5. Test the stepper motor by pressing the switch.

If you get the Arduino to run the stepper motor, all you have left to do is mount the components and stick a bent (and blunted) needle onto the 8/32" screw of the horizontal Ø1/2" post. You'll notice that we are leaving a few things up to you, which may change based on your abilities:

  • We used the right-angle mount from Thor Labs to place the motor, but we had to drill two small holes in the steel motor mount for it to work

  • The momentary switch shown is different than the one recommended in the parts list, it's all we had lying around and can be found at All Electronics

  • You will need to decide how you wan the tetrode clip to interface with the motor, and you can reference the original Open Ephys Twister for 3D parts, or do what we did, and scrap one together


The momentary switch will spin the motor 80-turns in one direction, and 20-turns in the other, and then stop. These parameters can be adjusted in the INO file provided in Step 4 (above). We found that the stepper motor introduces a slight vibration into our plastic "goal posts" used to hold the tetrode clamp, although we were able to dampen this using some rubber heat-shrink tubing. It is recommended that you unplug the Arduino from the AC charger when the spinner is not in use.

We hope this serves as a starting point for making your own spinner, and please share any innovations you have made in your quest for the perfect tetrode.

Extracting Spikes from Neural Electrophysiology in MATLAB

 Neural spikes extracted using this method.

Neural spikes extracted using this method.

One hour of neural recordings amasses to nearly forty-three gigabytes of raw data for me. This is streamed through fiber optics onto our storage system and accounts for 128-channels, sampled at just over 24 kHz. If it was a good day, these files contain hundreds of thousands of spikes, so how to extract them?

Extracting spikes is just one step in our lab's multi-step protocol to analyze animal behavior. With the flux of everyone from undergraduates to post-docs working with the code base, extracting spikes can’t be something only one person understands. This imposes a few constraints on the algorithm we implement—we want something that is reliable, but more importantly, simple.

The first step is the same no matter what extraction you use: we want to exclude the low frequency content of the signal. This will “flatten” the signal, and hopefully begin to highlight the high frequency neural spikes. I use a butterworth bandpass filter, with cutoff frequencies (Fc) between 244 Hz and 6.104 kHz. Depending on your sampling frequency (Fs), you can easily calculate your own Wn values (Wn = Fc/(Fs/2)).

>> [b,a] = butter(4, [0.02 0.5]);
>> filteredData = filtfilt(b,a,double(data));

At this point, you might think about removing high amplitude artifacts. Movement-related potentials will often soar above your spike amplitude, so everything above (and below) a certain level can be removed. It makes most sense to me to just apply a zero-amplitude segment in place of the artifact. You can find my code for artifactThresh.m, and here’s what it’s doing:

  1. Identify peaks above given threshold.
  2. Move forward and backward in time and identify when the signal reaches reasonable amplitude (when it "resets").
  3. Replace the artifact spans with zeros.

Finally, it’s time to do the detection itself. What we really want to know is, at what times (or sample numbers) are there spikes? Bestel et al. reviewed some of the detection methods in use by others [1], they include:

If you have a great signal to noise ratio, just drawing a threshold is by far the simplest method. You can even get pretty darn close just using peakseek.m by Peter O’Connor. Although we’ve done our best to eliminate any signal not within the “spiking” frequency band, unless your recordings are perfect you will still find that threshold method will catch a lot of noise. Being constrained to simple ideas, the Nonlinear Energy Operator (NEO) is the next obvious candidate, because as Sudipta & Ray show [2], you are basically squaring your signal but also subtracting the amplitude of neighboring samples. This ensures that the waveform is, indeed, spikey.

 Where  ψ  is nonlinear energy,&nbsp; x &nbsp;is your data, and  n &nbsp;is your sample number.

Where ψ is nonlinear energy, x is your data, and n is your sample number.

This serves to first, make big signals bigger, and small signals smaller (squaring), but also acts as a more selective high pass filter, as signals that are broad will have large amplitudes for neighboring samples. Now, you perform a threshold on the nonlinear energy value. Choosing a threshold is still a subjective process, although I have had success using a multiple of the median of nonlinear energy values. You can find all these operations in my function, getSpikeLocations.m (with two dependencies: peekseak.m and snle.m). The next step, which I don't cover here, will be extracting the waveforms based on the timestamps and spike sorting them using either more MATLAB software, or a commercial product, like Plexon's Offline Sorter.

1.   Bestel, R., Daus, A. W. & Thielemann, C. A novel automated spike sorting algorithm with adaptable feature extraction. J. Neurosci. Methods 211, 168–78 (2012).

2.   Mukhopadhyay, S. & Ray, G. C. A new interpretation of nonlinear energy operator and its efficacy in spike detection. IEEE Trans. Biomed. Eng. 45, 180–187 (1998).


Mouse and Rat Brain Atlas: An Interactive Online Tool

I've created online tools for the Rat Brain Atlas and Mouse Brain Atlas based on the Paxinos et al. work, in stereotactic coordinates. All the sections have been extracted from the original PDF documents. My software accounts for the subtle differences in page formatting and finds the nearest section match to your desired coordinate. Once you plot a point, you can send or save the URL—the tool is purposefully very stateful. Finally, I have included print style sheets, so you can print out your sections, should you need them for surgery (this is why I built the tools).