✎ analysis:nsb2016:week3long []

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
===== Module 3: Visualizing neural data in MATLAB =====

Goals:

  * Master the basics of plotting data in MATLAB using the versatile ''plot()'' function
  * Understand how to use //handles// to fine-tune formatting
  * Meet some existing visualization tools: ''MultiRaster()'', ''PlotTSDfromIV()'' and ''ft_databrowser()''
  * Use interactive figure properties with callback functions for dynamic data viewing
  * Learn to use functions with variable numbers of input arguments (using ''varargin'' and ''cfg'' inputs) effectively
  * Understand which file format to export figures to
  * Advanced bonus section: GUIs, movies, and sonification of spiking data

Resources:

  * (reference) MATLAB documentation for ''plot'', [[ http://www.mathworks.com/help/matlab/ref/figure_props.html | Figure properties]], [[ http://www.mathworks.com/help/matlab/ref/axes_props.html | Axis properties]] and [[ http://www.mathworks.com/help/matlab/ref/line_props.html | Line Properties]]
  * (reference) [[http://www.mathworks.com/help/matlab/creating_plots/figures-plots-and-graphs.html | MATLAB plotting gallery ]]
  * (reference) [[http://www.vandermeerlab.org/BIOL377_assignment3_2012.pdf | Place cell assignment]] from %%MvdM%%'s BIOL377 course at Waterloo from a few years ago (if you want a gentle intro to what is going on in the data you are looking at; can be done in 30 min)
  * (optional) Tufte, [[ http://www.edwardtufte.com/tufte/books_vdqi | The Visual Display of Quantitative Information ]]

===== Introductory remarks =====

Visualization is an essential component of any data analysis workflow. You would never perform brain surgery or another intricate experimental procedure without a good view of what you are doing -- if conditions are met for a next step, whether a manipulation has the desired outcome -- and likewise, you should not do data analysis blind. This means that after every step you want to be able to see what changed. 

Also, for many analyses, a visualization is part of the final outcome. Most neuroscience papers have complex multipanel figures, such as this one from [[http://www.nature.com/nature/journal/v469/n7330/fig_tab/nature09633_F4.html | Dragoi and Tonegawa]] (Nature, 2011):

{{:analysis:course:nature09633-f4.2.jpg?nolink&600|}}

Notice that panel Aa, a highly selective and organized view of raw data, already contains LFP data (top) and spike data from 18 cells!

This sort of figure is constructed by first producing the individual panels in MATLAB, and then perhaps combining them in a graphic design package such as Adobe Illustrator/InDesign. We will do the MATLAB part here and aim to produce something similar to panel Aa, while introducing the fundamentals of plotting in MATLAB. Other types of plots, such as bar plots, errorbars, color bitmaps, filled shapes, and shading will be introduced as they are needed in later modules.

An important principle in graphical design for scientific communication is to only include those elements which contribute to the messages(s) you want to communicate. Remove and simplify as much as possible: this will make your figures more readable. The Tufte book linked to above is a beautiful statement of this (and other) guiding principles of graphic design.

Even if you do not go as far as reading it, do at least make the labels in your figures big enough to be readable! :-)

==== Loading the data ====

If you haven't done so already, this would be a good time for a ''git pull'' to make sure you have the latest version of the course codebase!

Now, make sure you have the folder ''R042-2013-08-18'' from the data share, and that this is placed in a sensible location (NOT in a %%GitHub%% or project folder!). 

This data set was collected from the dorsal CA1 region of the hippocampus of a rat performing a T-maze task, using a drive with 16 independently movable tetrodes (by [[https://github.com/aacarey|Alyssa Carey]] for her Master's thesis work). Spike and LFP data was recorded from each tetrode; possible spike events were detected online and stored for offline spike-sorting, and LFPs were sampled at 2kHz and bandpass filtered between 1-475Hz. (A quirk of this particular data set is that certain time intervals are cut out of the spike data, but not the LFP. So you may notice some odd looking gaps in the rasterplot later.)

As in previous modules, use the shortcut button created in the previous module to set your path. In your project folder, create a folder with today's date and ''cd'' to it. This will be your working directory with stuff that you don't necessarily intend to push to %%GitHub%%.

Create a ''sandbox.m'' file in your daily folder. In this file, use cell mode to load some spike trains, a LFP, and position data as follows (recall you can use Ctrl+Enter to execute the code in a cell):

<code matlab>
%% cd to data folder -- replace this with yours
fd = 'D:\data\DataAnalysisTutorial\R042-2013-08-18';
cd(fd);

%% load the data (note, may need to unzip position data first)
cfg = [];
S = LoadSpikes(cfg)
</code>

This should give a data structure containing spike train data:

<code matlab>
LoadSpikes: Loading 67 files...

S = 

     type: 'ts'
        t: {1x67 cell}
    label: {1x67 cell}
      cfg: [1x1 struct]
      usr: [1x1 struct]
</code>

The details of this ''ts'' structure are discussed in [[analysis:course-w16:week2|Module 2]]. In brief, each cell of ''S.t'' contains the spike times from a putative neuron. The qualifier "putative" is used because this is extracellular data and [[analysis:course-w16:week4|spike-sorting]] is not perfect, so it's likely there will be some spikes missing and some spikes included that are not from this neuron. Always remember this even if I will omit the "putative" from now on for short. You can see we have loaded 67 neurons for now.

The ''LoadSpikes()'' function has some features not discussed here, such as the contents of the ''usr'' field; you can type ''help LoadSpikes()'' to learn more.

Let's also load a LFP:

<code matlab>
cfg = [];
cfg.fc = {'R042-2013-08-18-CSC03a.ncs'};
csc = LoadCSC(cfg)
</code>

This gives:

<code matlab>
LoadCSC: Loading 1 file(s)...
LoadCSC: R042-2013-08-18-CSC03a.ncs 0/17193 bad blocks found (0.00%).

csc = 

     type: 'tsd'
     tvec: [8802816x1 double]
     data: [1x8802816 double]
    label: {'R042-2013-08-18-CSC03a.ncs'}
      cfg: [1x1 struct]
</code>

''csc'' is the common Neuralynx designation for "continuously sampled channel" and typically is an EEG or LFP type signal sampled and filtered so that high-frequency components such as spikes are not accessible. It is possible to have wide-band, 32kHz CSCs suitable for spike extraction, but these are not included in the current dataset. As discussed more extensively in [[analysis:course-w16:week2|Module 2]] a LFP is defined by matching sets of sample timestamps (''csc.tvec'') and sampled data (''csc.data'').

Before we proceed, let's restrict the data we have loaded (about 90 minutes' worth) to a more manageable size.

A useful function that works on ''ts'', ''tsd'' and ''iv'' objects is ''restrict()''.

☛ Use MATLAB's ''help'' function to look up the usage information for ''restrict()''. In a new editor cell, use this function to create ''S_r'' and ''csc_r'' variables, containing the spike and LFP data respectively, restricted between 5900 and 6000 seconds.

Hint: make a habit of using variables instead of hard-coding specific values. Applying this [[analysis:nsb2015:week0|good programming practice]] to the above means that you first define the interval of interest, and then use the same resulting variable to two calls to ''restrict()''. That way, if you ever want to change this interval, you can do so in a single, clearly visible location at the beginning of your code. This is much more robust than having to scan through the whole script to find and modify each instance of some hard-coded numbers!

☛ Verify that the restriction worked by inspecting the ''csc_r'' variable -- you should notice that the length of the ''tvec'' and ''data'' fields has been suitably reduced.

===== Do-it-yourself plotting =====

If you are already familiar with MATLAB's basic plotting functions and how to use object handles to set properties of drawing objects, you can skip this section.

==== Basic plot commands ====

Let's look at the spike times of the 17th cell in our data set (note the use of the curly brackets ''{}'', if this is confusing, review the MATLAB help on cell arrays):

<code matlab>
%%
iC= 17;
this_spk = S_r.t{iC}

this_spk =

   1.0e+03 *

    5.9046
    5.9046
    5.9046
    5.9046
</code>

That's not a very helpful format, it looks like all four spikes have the same spike time, but note the ''1.0e+03'' (meaning 10^3, or 1000)! 

☛ Type ''format bank'' and try again.

You should now see somewhat more clearly the spike times of the five spikes this neuron emitted (in our restriction interval).

An alternative to ''format'' is to use the ubiquitous ''fprintf()'' command:

<code matlab>
>> fprintf('%4.3f\n',this_spk)
5904.630
5904.635
5904.638
5904.643
</code>

This way, you have precise control over the formatting of command line output -- the cryptic ''%4.3f'' tag specifies that the contents of ''this_spk'' should be formatted as a floating-point number with 4 digits before, and 3 digits following, the decimal point; ''\n'' specifies a newline. (For reference: ''fprintf()'' has [[http://www.mathworks.com/help/matlab/ref/fprintf.html#input_argument_formatspec | many other formatting options]]).

Now that we know what we are dealing with, let's meet the plot() command:

<code matlab>
plot(this_spk)
</code>

You should get:
{{ :analysis:course-w16:plot_tutorial1.png?nolink&600 |}}

This may not be what you were expecting! If you give ''plot()'' only one input argument (''this_spk'' in this case) it will, by default, plot the **index** of each element against its **value**.

Let's unpack that statement so it is really clear: we have an array with five values here: the spike times of four spikes of neuron 17. ''plot()'' plots the value of the first element (''this_spk(1)'') at x-coordinate ''1'' (its index, i.e. position in the array), and so forth for the whole length of the array. It also connects the data points with a blue line.

Let's plot the same data differently:

<code matlab>
plot(this_spk,0,'.k');
</code>

This gives:

{{ :analysis:course-w16:plot_tutorial2.png?nolink&600 |}}

A pretty different result! Because we have given ''plot()'' multiple input arguments, it interprets the first (''this_spk'') as the x-coordinate, and the second (''0'') as the y-coordinate. Indeed you can see in the plot that the spike times now show up on the x-axis, in contrast to the first plot we made, and the y-coordinate is 0 for all spikes.

The final argument ('''.k''') specifies we want the data plotted as dots (''.'') and in black ('''k'''). 
 
Plotting spikes as dots is okay, but tickmarks (vertical lines, as in the example figure at the top) are better. Here is how to do it:

<code matlab>
h = plot([this_spk this_spk],[0.5 1.5],'k');
axis([xlim 0 5])
</code>

To understand what happened here, recall that ''plot()'' uses the first argument as x-coordinates, and the second as y-coordinates. So the above plot command used our spike times (''this_spk'') as the x-coordinates, drawing a black line at each spike time between the specified y-coordinates ''[0.5 1.5]''.

The ''axis'' command sets the axis limit according to the 4-dimensional array ''[xmin xmax ymin ymax]''; ''xlim'' returns the current ''xmin'' and ''xmax''. What the ''h'' output argument is for will be explained in the next section.

☛ (Test your MATLAB skills) Extend the example code above to make a rasterplot containing the spikes for all cells in ''S_r''. Each cell should have its spikes appear on the corresponding row, as in the example plot at the top. To facilitate re-use later, wrap your code in a function, %%PlotSpikeRaster%%, that plots the raster for any ''ts'' input. Thus, it should also work for task events, for instance, by virtue of following the common ''ts'' datatype specification. //Hint//: how will you handle the case where some cell doesn't have any spikes?

==== Using handles to fine-tune plotting ====

At the end of the previous section, I specified an //output argument// ''h'' with the ''plot()'' command. This output is a **handle**, a MATLAB data type that contains various properties associated with an object. In this case, ''h'' is a handle for the spikes we just plotted, and we can use it to change their appearance.

☛ Below are a few examples of how to use a plot handle. Try them -- one at a time -- and see what happens!

<code matlab>
set(h,'Color',[0 1 0]); % "set Color property of handle h to value [0 1 0] i.e. green in [red green blue] format
set(h,'LineWidth',2);
set(h,'Marker','.','MarkerSize',20); % note you can specify multiple properties in one set() function call
</code>

A list of all the properties of a handle can be obtained by typing ''get(h)''.

MATLAB has a number of built-in handles, such as ''gca'' ("get current axes"):

<code matlab>
set(gca,'LineWidth',2);
set(gca,'XTick',5904.62:0.01:5904.68,'YTickLabel',{});
set(gca,'FontSize',24);
</code>

☛ Look up what properties the ''gca'' handle has to find out how to change the color of the axes and labels, and change it to red.

Note that if you are making figures with multiple axes (perhaps by using ''subplot()'' or ''plotyy()'') each axes gets its own handle. One axes is active at any given time, which it can be by clicking on it, or by explicitly doing something like ''axis(h)''.

There is also a handle for the whole figure, ''gcf'':

<code matlab>
set(gcf,'Color',[0 0 0]);
</code>

Finally, MATLAB has a special handle, '0', for some defaults. For instance, a very useful default to change is the font size, ''set(0,'DefaultAxesFontSize',18)'' because the default is usually too small to be readable once exported (see the next section on this). A good place for such a default change is in a shortcut that also contains your path, or perhaps even in MATLAB's ''startup.m''.

==== Exporting figures ====

Now that our figure looks nice, let's save it as an image file:

<code matlab>
print(gcf,'-dpng','-r300','R042-2013-08-18-spikes.png');
</code>

Notice that the first argument of ''print'' is the current figure's handle. The other arguments specify that we want a ''PNG'' file, with 300dpi resolution, and the filename to write to. PNG format is a good choice for saving figures because it uses [[http://www.howtogeek.com/142174/what-lossless-file-formats-are-why-you-shouldnt-convert-lossy-to-lossless/ | lossless compression]], in contrast to JPEG images which use lossy compression and can have ugly artifacts as a result. Other useful save formats include ''-dill'' which saves ''.ai'' Illustrator files, and ''-deps'' which saves encapsulated PostScript; both of these are vector graphics formats.

☛ Look at the image file. Assuming you changed some of the figure and axis colors, you should find the colors in your image don't match those in MATLAB. Solve this problem by turning off the ''InvertHardCopy'' property of the figure, and save again. (MATLAB does this by default to facilitate printing images on white paper.)

You should get something like:

{{ :analysis:course-w16:plot_tutorial3.png?nolink&600 |}}

Obviously, I do not recommend formatting your rasterplots with this particular color scheme for publication!

==== Putting it all together ====

Let's add some fancy plot features together:

<code matlab>
%% restrict the data to interval of interest
this_iv = iv([5900 6000]);
S_r = restrict(S,this_iv);
csc_r = restrict(csc,this_iv);

%% plot spikes
SET_spkY = 0.4; % parameter to set spike height

figure; hold on;
for iC = 1:length(S_r.t)
    if ~isempty(S_r.t{iC})
        plot([S_r.t{iC} S_r.t{iC}],[iC-SET_spkY iC+SET_spkY],'k');
    end
end % of cells
ylabel('neuron #');

%% add a LFP
SET_cscY = [-5 0];
plot(csc_r.tvec,rescale(csc_r.data,SET_cscY(1),SET_cscY(2)),'r');
set(gca,'YTickLabel',{});

%% add multi-unit activity in separate axes
ax1 = gca; ax1_pos = get(ax1,'Position');
ax2 = axes('Position',ax1_pos); % new axes with same position as first
 
cfg = []; cfg.tvec = csc.tvec; cfg.sigma = 0.1;
mua = getMUA(cfg,S); % obtain multi-unit activity
 
xr = get(ax1,'XLim');
mua_r = restrict(mua,xr(1),xr(2)); % only keep the data we need
 
axes(ax2); % set current axes to the second set, so what follows happens there
mua_hdl = plot(mua_r.tvec,mua_r.data,'Color',[0.7 0.7 0.7]);
 
set(gca,'YAxisLocation','right','Box','off','XTick',[],'Color','none','YColor',[0.5 0.5 0.5])
ylabel('multi-unit activity (spk/s)');
linkaxes([ax1 ax2],'x'); 
</code>

A few things to note here:
  * Using the ''Position'' property of the current axes (ax1, containing the rasterplot) we created a new set of axes (ax2).
  * In these new axes, we plotted the multi-unit activity, obtaining an output argument. Like ''gca'' and ''gcf'' this is a **handle** to graphics object; in this case, a handle to the multi-unit activity signal we just drew.
  * We set the properties of the new multi-unit axes to have the y-axis on the right (plus a few other properties).
  * ''linkaxes()'' was used to link the x-axis of both axes so that both update when zooming in.

☛ Use the 'XLim' axes property to zoom in to 5965 to 5969 s.

You should now be looking at the synchronous activation of a substantial number of neurons, reflected in the MUA peak, and associated with a high-frequency oscillation (~150-250Hz) in the LFP. These are the neurophysiological signatures of a "sharp wave-ripple complex" (SWR for short), events which are thought to contribute to the consolidation and retrieval of episodic memories.

{{ :analysis:course-w16:plot_tutorial5.png?nolink&600 |}}

It would be nice to be able to update what time window of the data we are looking at, without having to type these ''XLim'' commands. To do this we need to use a special Figure property, introduced in the next section. 


===== Advanced topics =====

The topics that follow in this section are optional, in the sense that later modules do not assume you know how to do these things. Feel free to use this section as you see fit, but do make sure you go through the next section ("Using existing visualization tools")!

==== Interactive figures and callback functions ====

Figure windows in MATLAB have many properties (here is the complete [[http://www.mathworks.com/help/matlab/ref/figure-properties.html | list]]). A particularly useful one is the ''KeyPressFcn'' property, which specifies a function to be called when a key is pressed while the figure is active.

For example, we can write a simple function that enables us to scroll left and right using the arrow keys, as follows:

<code matlab>
function figscroll(src,event)

ax = get(src, 'CurrentAxes');
x_orig = get(ax, 'XLim');
x_step = (x_orig(1)-x_orig(2))/2;

switch event.Key 
    case 'leftarrow'
        x_new = x_orig + x_step;
    case 'rightarrow'
        x_new = x_orig - x_step;
end

set(ax,'XLim',x_new);
</code> 

You should be able to interpret what this function is doing, with the only potentially mysterious part the identity of the input arguments ''src'' and ''event''. These input arguments are common to all of MATLAB's [[http://www.mathworks.com/help/matlab/creating_plots/callback-definition.html | callback functions]], i.e. functions that are called in response to some user action like a key press:

  * The ''src'' ("source") argument contains details about where the callback came from; in this case this will be a figure, and as a result we can use this argument to find the handle to the figure's current axes, and make changes.
  * The ''evt'' ("event") argument tells us about the event that triggered the callback. As you can see, we use it here to find out which key was pressed, and react accordingly.

This ''figscroll'' function as shown above isn't going to do much by itself. We need to tie it to a figure.

☛ Create a rasterplot such as the one in the previous section. Then, set the ''KeyPressFcn'' property of the figure to this new function, like this: ''set(gcf,'KeyPressFcn',@figscroll)'' (the ''@'' indicates what follows is a [[http://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html | function handle]]).

Now you should be able to use the left and right arrow keys to scroll through the rasterplot. If nothing happens, make sure you saved the ''figscroll'' function in a place where MATLAB can find it (i.e. in the path, or MATLAB's current working directory).

If you haven't encountered function handles before, here is another example of a function handle (to an [[http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html | anonymous function]], i.e. a function that doesn't have a ''.m'' file):

<code matlab>
sqr_fn = @(x) x.^2;
sqr_fn(2)
</code>

The above example basically says "sqr_fn is a function of x, and should return x squared." (if you don't understand the ''.'' notation, it is important to look this up!)
==== Exporting to movies ====

As an alternative to scrolling manually through the data, we can use a loop:

<code matlab>
%% movie version
t = [5900 6000]; % start and end times (experiment time)
FPS = 30; % frame rate (per s)
twin = [-1 1]; % width of time window (in s)

tvec = t(1):1/FPS:t(2);
for iT = 1:length(tvec)
    
   set(gca,'XLim',twin+tvec(iT));
   drawnow; pause(1/FPS);
    
end
</code>

You may notice some annoying auto-scaling behavior for the MUA y-axis, you can fix this using the axis '''YLim''' property.

By making this MATLAB animation into a movie file, it is often easier to explore the results. To do this, we can run the animation code above, with a few small modifications. First, before entering the main plotting loop, set the figure to be used to a specific size:

<code matlab>
h = figure; set(h,'Position',[100 100 640 480]);
</code>

This is important first, to keep the size of the resulting movie file manageable (the above sets a 640x480 pixel figure size), and second, because many movie encoders (such as the excellent and free [[http://www.xvid.org/ | XVid]]) will only work with certain sizes.

Next, we need to store each frame into a variable that we can later write to file. Modify the last two lines inside the loop to:

<code matlab>
f(iT) = getframe(gcf); % store current frame
drawnow;
</code>

If you now run the code again, each frame gets stored in the ''f'' variable as the loop runs. Then, we can write the result to a file:

<code matlab>
fname = 'test.avi';
movie2avi(f,fname,'COMPRESSION','XVid','FPS',FPS);
</code>

The above will only work if you have the XVid codec installed: I highly recommend this because it creates movie files that are an order of magnitude smaller than uncompressed files. If you have trouble with XVid, you can of course still save an uncompressed file for now. For longer movies, it is often required to save a file, say, every 500 frames, to prevent the ''f'' variable getting too large. These segments can then be merged with a video editing program such as [[http://www.virtualdub.org/ | VirtualDub]] (Windows only AFAIK; please suggest OSX/Linux alternatives if you know any that work well!).

All the figure and axis properties we explored earlier can be used here, too -- so if you don't like the gray figure background, it's easily changed to e.g. black.
==== GUIs ====

An alternative to interactive keyboard input using Figure properties is to create a graphical user interface (GUI). MATLAB has a nice tool, ''guide'', which allows you to place various UI controls (such as buttons, text boxes, drop down menus, etc.) and display items (axes, text labels, etc.) using a graphical interface. The display items all have //handles// (whose names you can set in the ''guide'' tool), and the UI controls all have //callback functions// so that you can have UI control actions change what is displayed.

The GUIDE documentation has a nice [[http://www.mathworks.com/help/matlab/gui-building-basics.html | tutorial]], with this [[http://www.mathworks.com/help/matlab/creating_guis/about-the-simple-guide-gui-example.html | example]] as a good place to start.
==== Sonification ====

I'm a fan of silent movies such as [[http://www.imdb.com/title/tt0017136/ | Metropolis]], but if you plan to use movies of your data in any talks, you might want to add some sound. One way I've done this is to use Ken Schutte's [[https://github.com/kts/matlab-midi | matlab2midi]] toolbox to convert ensemble spiking data into a MIDI file, and then to synthesize the MIDI file into sound using something like Ableton Live.

To see how this works, first clone the matlab2midi repository and add the files to your MATLAB path.

A simple example is given by Ken:

<code matlab>
% initialize matrix:
N = 13;  % num notes
M = zeros(N,6);

M(:,1) = 1;         % all in track 1
M(:,2) = 1;         % all in channel 1
M(:,3) = (60:72)';      % note numbers: one octave starting at middle C (60)
M(:,4) = round(linspace(80,120,N))';  % lets have volume ramp up 80->120
M(:,5) = (.5:.5:6.5)';  % note on:  notes start every .5 seconds
M(:,6) = M(:,5) + .5;   % note off: each note has duration .5 seconds

midi_new = matrix2midi(M);
writemidi(midi_new, 'testout.mid');
</code>

From reading the comments, you can get an idea of how the MIDI file specification works -- it describes a piece of music using a number of different variables, here corresponding to the different rows of the ''M'' matrix. For our purposes, the most important are the "note on" and "note off" rows (5 and 6) and the pitch row (3). In converting spikes to sound, the spike times will determine what we put in rows 5 and 6, and we use the pitch row to assign a unique sound to each different cell.

So, we first need to create an array containing all spike times, with a corresponding array for each spike's cell number:

<code matlab>
t = [5900 6000]; % time window

all_spikes = []; all_ids = [];
for iC = 1:length(S_r.t)
   all_spikes = cat(1,all_spikes,S_r.t{iC});
   all_ids = cat(1,all_ids,repmat(iC,length(S_r.t{iC}),1)); % cell ID
end
[all_spikes,srt] = sort(all_spikes,'ascend');
all_ids = all_ids(srt); clear srt;

all_spikes = [all_spikes; t(end)]; % add spike at end time
all_ids = [all_ids; 0];

all_spikes = all_spikes - t(1); % start at time 0
</code>

Then, we can write the result:

<code matlab>
note_time = 0.05;
base = 20; % base pitch

nSpikes = length(all_spikes);
M = zeros(nSpikes,6); % output matrix

M(:,1) = 1; M(:,2) = 1; M(:,4) = 100;
M(:,3) = base+1.5*all_ids;
M(:,5) = all_spikes;
M(:,6) = all_spikes + note_time;

midi_new = matrix2midi(M);
writemidi(midi_new, 'r042.mid');
</code>

If you can't get this to work, the resulting file can be made to sound like this {{:analysis:course-w16:r042.mp3|here}} (with Ableton Live); Windows Media Player synthesizes it with a plain piano instrument, which sounds a bit different but is also effective at communicating the sudden, violently synchronous activity during SWRs!
===== Using existing visualization tools =====

==== MultiRaster ====

[[https://github.com/youkitan | Youki]] and [[https://github.com/aacarey | Alyssa]] in my lab wrote a fully-featured plotting function called ''MultiRaster()''. As an example of what it can do, try this:

<code matlab>
%% load data
S = LoadSpikes([]);
please = []; please.fc = {'R042-2013-08-18-CSC03a.ncs'};
csc = LoadCSC(please);

LoadMetadata; % load some experiment metadata such as trial start and end times

%
cfg_plot = [];
cfg_plot.lfp = csc;
cfg_plot.spkColor = 'jet';
cfg_plot.evt = metadata.taskvars.trial_iv_L; % "left" trials on the T-maze

h = MultiRaster(cfg_plot,S);
</code>

You should get:

{{ :analysis:course-w16:plot_tutorial6.png?nolink&600 |}}

Now, press the 'h' key to bring up a window explaining the various keyboard shortcuts available. These work through the ''KeyPressFcn'' method introduced above; you can examine MultiRaster's ''KeyPressFcn'' by typing ''edit navigate'' (it's called ''navigate.m'').

MultiRaster also uses a trick to speed up the plotting of spikes; it's many times faster than the method we used in this module so far, but it is more difficult to read the code -- take a look if you are interested in finding out how it works!

Notice the output argument ''h'', this is actually a struct with handles to the various plot objects. MultiRaster's help explains what they are.

Plotting time series data, such as the LFP used here, results in large numbers of data points which slow down MATLAB. A good way of further speeding up plotting is to decimate (i.e. downsample) time series data. So here, instead of plotting the full 2kHz signal, plotting it at 500Hz would be just as good and result in significant speedup, especially when scrolling/zooming a lot.

=== Diversion: handling variable numbers of function inputs ===

A standard function definition specifies the exact number of input and output arguments the function expects: ''function y = sqrt(x)'' has one input and one output argument. However, it is often useful to have a variable number of input arguments, so that you can override defaults or specify additional options as needed without complicating simple function calls. This section shows you two ways of doing that.

== Method 1: cfg input ==

As you can see in the above example, and in MultiRaster's help, using a ''cfg'' struct as an input is one way to handle variable numbers of inputs to a function. Inside the function, you can do things like, "if cfg has a field called lfp, then plot it". It is also a good way to deal with defaults which may be optionally overwritten. For instance, MultiRaster by default plots all spikes in black, but this can be overridden, as we did in the above example. Here is how that works (taken from inside MultiRaster, do not run this):

<code matlab>
cfg_def.spkColor = 'k';

cfg = ProcessConfig(cfg_def,cfg_in);
</code>

What happens is that inside MultiRaster, a default config struct (''cfg_def'') is defined. The function ''ProcessConfig()'' compares the default to the incoming config (''cfg_in''); if anything in ''cfg_in'' is different, then the default is overridden. That way, the function can be called simply without specifying each option, because it will just use the defaults. But, if you want, you can specify something else in the input config.

== Method 2: varargin ==

An alternative approach is to use the special input argument ''varargin''. Here is the idea:

<code matlab>
function test_fun(a,varargin)

b = 2; % set defaults for b and c
c = 3;

extract_varargin; % override b and c if specified

fprintf('a is %d, b is %d, c is %d\n',a,b,c);
</code>

Test it:

<code matlab>
>> test_fun(1)
a is 1, b is 2, c is 3
>> test_fun(1,'b',1)
a is 1, b is 1, c is 3
>> 
</code>

Note how specifying a "key-value" pair ('''b',1'') overwrote the default for ''b''. ''extract_varargin'' basically goes through the varargins and assigns each value to the corresponding key, similar to ''ProcessConfig()''.

☛ Verify your understanding by calling ''test_fun'' such that ''b'' and ''c'' are both set to ''0''. Also note what happens when you do ''test_fun(1,'B',1)''!

==== PlotTSDFromIV ====

''PlotTSDFromIV()'' is a simple, lightweight function designed to show time series data (TSD) for specific intervals (IV, hence the name). It has two different display modes, the first one ('''iv''') looks like this:

<code matlab>
cfg = [];
cfg.display = 'iv';
cfg.width = 0.1;
PlotTSDfromIV(cfg,metadata.SWRtimes,csc); % need to LoadMetadata to make this work!
</code>

{{ :analysis:course-w16:plot_tutorial7.png?nolink&600 |}}

(The intervals highlighted in red are manually selected, putative SWR events.)

☛ Try the other display mode, '''tsd'''.
==== FieldTrip's databrowser ====

As implied by the name, ''MultiRaster()'' is primarily intended for plotting spike data, although it can accept multiple tsd's for plotting as well. 

The %%FieldTrip%% toolbox, which will be covered in more detail in later modules, has a number of data visualization tools that are better suited for plotting time series data such as LFPs, EEG/MEG and single-trial fMRI data. 

To get our current data into the %%FieldTrip%% format, we first have to do a few steps:

<code matlab>
% convert to ft format
cfg = []; cfg.mode = 'resample';
csc_ft = TSDtoFT(cfg,csc);
t0 = csc_ft.time{1}(1); csc_ft.time{1} = csc_ft.time{1}-t0;

% create ft trials from iv
trl_cfg = [];
trl_cfg.t = IVcenters(metadata.SWRtimes)-t0;
trl_cfg.mode = 'neuralynx';
trl_cfg.hdr = csc_ft.hdr;
trl_cfg.twin = [-1 1];

trl = ft_maketrl(trl_cfg);

% use trials to create trialified data structure
temp_cfg = []; temp_cfg.trl = trl;

ft_in = ft_redefinetrial(temp_cfg,csc_ft);
</code> 

Now we can run the data browser by typing ''ft_databrowser([],ft_in);'', to get

{{ :analysis:course-w16:plot_tutorial8.png?nolink&600 |}}

Notice the various control buttons for moving between trials, channels (here we have only one loaded), and other functions, documented [[http://www.fieldtriptoolbox.org/faq/how_can_i_use_the_databrowser | here]]. %%FieldTrip%% is able to load many types of data files, see [[http://www.fieldtriptoolbox.org/reading_data | here]] for the list and how to use its data loaders.
User Tools

Site Tools

Page Tools
