wiki

~~DISCUSSION~~ ===== Time-frequency analysis: spectrograms ===== Goals: * Understand the idea of the spectrogram and the inputs and outputs of typical implementations * Apply your knowledge of spectral leakage to the selection of spectrogram parameters * Cultivate an awareness of common pitfalls in spectrogram-based analysis * Get started with using the fieldtrip toolbox Deliverables: * Code to construct a event-triggered LFP plot Resources: * (for reference) [[http://www.chronux.org/tutorials | Chronux toolbox slides]] * (for reference) [[http://fieldtrip.fcdonders.nl/ | Fieldtrip documentation]] and [[http://fieldtrip.fcdonders.nl/tutorial/timefrequencyanalysis | spectrogram tutorial]] ==== Introduction ==== The spectral estimation methods introduced in the previous module average spectral content over a number of different windows. This means that we get an improved estimate of the //overall// spectral content of the signal, but such an average cannot establish if the spectral content is changing over time. Neural signals rarely have stationary spectral properties, as illustrated by this snippet of a LFP from the rat ventral striatum (source: [[http://www.frontiersin.org/integrative_neuroscience/10.3389/neuro.07.009.2009/abstract | van der Meer and Redish, 2009]]) {{ http://c431376.r76.cf2.rackcdn.com/544/fnint-03-009/image_m/fnint-03-009-g004.jpg?600 |}} The top trace illustrates some striking properties of the ventral striatal LFP: under certain conditions, gamma oscillations appear in "events" of a relatively stereotyped length, with frequencies alternating between ~80 Hz ("gamma-80" or "high gamma") and ~50 Hz ("gamma-50" or "low gamma"). How can we characterize the spectral properties of such a rapidly changing signal? One approach is illustrated in the colorful panel directly below the raw LFP trace -- this is a **spectrogram**. A spectrogram estimates the spectral content of a signal over time. It has time and frequency axes, with power shown in the color dimension (note the colorbar: in this plot, red is high power, blue is low power). The spectrogram is constructed by sliding a window along the signal and computing the power spectrum at each position of the window. As with the PSD estimators in the previous module such as Welch's method, the window is characterized by its size and shape (rectangular, Hamming, etc.) which control the tradeoffs resulting from spectral leakage. Unlike these PSD estimators however, to obtain a spectrogram we do not average across windows. Instead, we keep the spectrum for each position of the window, filling in each "time" column of the spectrogram as it moves along. (As an aside, spectrograms are a widely used tool for time-frequency analysis; for instance, they are often the first step in voice recognition software. Plotting a spectrogram of your own voice can be done with just a few lines of MATLAB code.) ==== Constructing a basic spectrogram ==== === Understanding the spectrogram() function === MATLAB's signal processing toolbox includes basic functionality for constructing spectrograms. ☛ Look at the documentation for the ''spectrogram()'' function. Notice that the input arguments it takes are familiar to you from the last module: window, overlap, nfft, Fs. The usage we are looking for is the last one: <code matlab> [S,F,T,P] = spectrogram(...) P is a matrix representing the Power Spectral Density (PSD) of each segment. </code> Let's give it a spin: <code matlab> %% load the data % remember to cd to the correct folder here, may need to get this file from the lab database fname = 'R016-2012-10-03-CSC04a.Ncs'; csc = myLoadCSC(fname); hdr = getHeader(csc); Fs = hdr.SamplingFrequency; %% restrict data cscR = Restrict(csc,3282,3286); % if you don't have this, implement it (it's one line of code!) plot(cscR); % note there are various gamma oscillations present, as well as a large negative-going transient %% construct and plot the spectrogram [S,F,T,P] = spectrogram(Data(cscR),hanning(512),256,1:200,Fs); imagesc(T,F,10*log10(P)); % converting to dB as usual set(gca,'FontSize',20); axis xy; xlabel('time (s)'); ylabel('Frequency (Hz)'); </code> You should see: {{ :analysis:course:week6_fig1.png?600 |}} As in the example, color indicates signal power at a particular point (bin), corresponding to a specific frequency and a specific time. Some power is apparent at low frequencies (<10Hz), as well as in the gamma range (~50-100Hz). Let's deconstruct this plot a little further. Since we specified 512-sample windows with a 256-sample overlap, we expect roughly 8 estimates per second. This is because our Fs is 2000, and therefore 8 "steps" of our moving window, moving at 256 samples per step (the window size minus the overlap equals the step size) fit into one second. Indeed we can discern approximately 8 "columns" (time bins) per second in the plot. For the frequency axis, we specified ''1:200'' as the bins of interest, so our frequency axis has a resolution of 1 Hz per bin. The centers of the time and frequency bins are returned as the ''T'' and ''F'' outputs respectively. An important point about the time axis is this: <code matlab> >> T(1) ans = 0.1280 </code> That is, the time corresponding to the first column (time bin) of the spectrogram is //not// 0. It is also not the time of the first sample in our data, as is readily verified: <code matlab> >> tvec = Range(csc_postRec); >> tvec(1) ans = 3.2820e+03 </code> In fact, the time of the first bin returned by ''spectrogram()'' is the //centre of the first window for which there is enough data//. The data passed to ''spectrogram()'' is assumed to start at time 0. This makes sense when we remember that the input arguments do not contain timestamps (only voltage samples) so this function has no way of knowing the "experiment time" (''Range(csc)''). Thus, the time of the first spectrogram bin (''0.128'') is because this is the centre of the first 512-sample window for which there is enough data (''256/2000'' is ''0.128''), and the first sample of the data is assigned time 0. === Aligning the spectrogram with a LFP === Before we start investigating the effects of changing the spectrogram parameters, it is helpful to be able to see the raw LFP overlaid onto the spectrogram: <code matlab> hold on; tvec = Range(cscR); data = Data(cscR); lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting tvec0 = tvec - tvec(1); % align LFP with spectrogram data = rescale(data,-lfp_minmax,lfp_minmax); data = data+lfp_cent; lfp_h = plot(tvec0,data,'k'); </code> Note that to make the LFP and the spectrogram line up, the time of the first LFP sample needs to be zero. The rescaling of the LFP is just to make it easier to see when plotted overlaid onto the spectrogram; of course, if you prefer, you can plot it below or above without overlap by changing ''lfp_cent''. Now that we can see the LFP and the spectrogram together in one plot, you can observe the correspondence between them. The sharp transient early in the data shows up as a dark red in the low frequency range. This is perhaps not obvious, because such a transient can hardly be described as an oscillation -- it does not even complete a full period (cycle). So beware: if you see strong power in a spectrogram (or a PSD for that matter) this does not imply that there is a true oscillation in the signal. In contrast to this low-frequency transient, you can see that there are various gamma-range oscillations which extend over multiple cycles, and can therefore be regarded as true oscillations rather than transients. For instance, at the start of this data segment there is some low gamma around 50-60Hz, and later some higher frequencies appear. Another thing you will have noticed is that the spectrogram looks pretty coarse. Let's see if we can improve this by changing the input parameters. ==== Spectrogram parameters ==== === Frequency and time resolution === In the example above, we used a 512-sample Hanning window with an overlap of 256; that is, the window moves in steps of 256 samples, and as a result we get a new power estimate every 0.128 seconds. In addition, we requested power estimates at integer frequencies from 1 to 200. ☛ Change the code above so that ''spectrogram()'' uses an overlap of 384 samples (instead of 256) and returns power with 0.25 Hz resolution (instead of 1). The spectrogram should now look somewhat smoother. Because of the larger overlap, the spectrogram now contains an estimate every 512-384 = 128 samples, rather than every 256. However, it is important to appreciate that increasing the frequency and time resolution in this way does not change the underlying estimate. This can be seen, for instance, by setting the frequency resolution to something really small (e.g. 0.01). This will take a while to compute, but if you compare the spectrograms side-by-side, it should be clear that the frequency estimates do not become more precise. Similarly, you can decrease the overlap to ''window_length - 1'', resulting in a spectrogram that is ultra-smooth in appearance, but this cannot escape the underlying limitation that spectral content is estimated over a window that is 512 samples wide. An important parameter that we know can affect our estimate is the shape of the window used. ☛ Change the window used to a rectangular one (''rectwin''). What happens to the spectrogram? Tracking changes over time using a spectrogram generally means that we don't want to use very large time windows. A large window would obscure rapid changes in spectral content, such as the tendency for a low-gamma to follow high-gamma events illustrated in the first figure, above. However, as you know from the previous module, the effects of spectral leakage are more prominent in smaller windows: there is fewer data to make up for the edge effects. This means that spectrograms tend to be particularly sensitive to the choice of window. A Hanning window is often a good default choice. As we will see below, however, there are methods we can use to improve further on this. === Frequency and time tradeoffs === The most important choice in constructing spectrograms is the size of the time window. This controls the fundamental tradeoff between our ability to resolve //which frequencies// are present in the signal, and at //what time// those frequencies occur. This tradeoff arises because of spectral leakage. With a large window, the edge effects are relatively small, and frequency content can be estimated relatively accurately. However, we cannot know at what time in this large window a given oscillatory event occurs. In contrast, using a small window, we may be able to localize more precisely in time when an event occurs. The price we pay is that we have increased spectral leakage, and we are now limited in our ability to determine the frequency of the event. This idea is illustrated in panel (b) of the figure below: {{ http://fieldtrip.fcdonders.nl/_media/tutorial/timefrequencyanalysis/tfrtiles.png?600 |}} Rectangular "time-frequency tiles", with different widths and heights, correspond to the area in the spectrogram that would be highlighted by given oscillatory event of a fixed frequency for different choices of time window. For a small window, we can know with more precision when the event occurred, but we are uncertain about the frequency. In contrast, a large window constrains the frequency estimate, but now we are uncertain about time. ☛ For the piece of data used here, compare side by side (using ''subplot()'') the spectrogram obtained with 512- and 1024-sample windows. For a fair comparison, set the overlap such that you end up with a comparable number of time bins. You should get something like this: {{ :analysis:course:week6_fig2.png?600 |}} As expected, the larger window has tighter frequency estimates, but the spectrogram is now more smeared out in time. These observations raise the question of if there are ways of picking the window size that are somehow optimal. In general, there is no substitute for looking at the raw data and tuning the parameters so that features of interest, such as the gamma events here, are well resolved in the spectrogram. However, there are some sophisticated methods available for which we have to move to the FieldTrip toolbox, which we will get to later in this module. ==== Pitfalls ==== === Mind the gaps === As noted above, the input to ''spectrogram()'' is not a tsd containing samples with accompanying timestamps; rather, it is simply a list of samples assumed to be acquired at a constant sampling frequency ''Fs''. If the input in fact does not have a constant ''Fs'', as can occur for instance if there are gaps in the data, there will be trouble: <code matlab> cscR = Restrict(csc,3300,3340); % section of data with a gap [S,F,T,P] = spectrogram(Data(cscR),rectwin(256),128,1:200,Fs); imagesc(T,F,10*log10(P)); set(gca,'FontSize',20); axis xy; xlabel('time (s)'); ylabel('Frequency (Hz)'); hold on; tvec = Range(cscR); data = Data(cscR); lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting tvec0 = tvec - tvec(1); % align LFP with spectrogram data = rescale(data,-lfp_minmax,lfp_minmax); data = data+lfp_cent; lfp_h = plot(tvec0,data,'k'); xlim([tvec0(1) tvec0(end)]); </code> {{ :analysis:course:week6_fig3.png?600 |}} As you can see, ''spectrogram()'' ignores the gap present in the LFP, causing the spectrogram and the LFP to become misaligned. In general, always verify that your data does not have gaps when computing spectrograms! === Sampling frequency mismatches === Another way for the spectrogram and LFP to become misaligned is when an incorrect sampling frequency is specified. ☛ Return to the original data segment (''cscR = Restrict(csc,3282,3286);'') and compute the spectrogram using ''Fs-100'' instead of ''Fs''. Plot the spectrogram and overlay the LFP. What do you notice? ==== Loading task events ==== Computing a spectrogram enables a number of further analyses, such as characterizing the relationship between power at different frequency bands. An example of this is shown by the correlation matrix in panel C of the first figure on this page, which reveals that ventral striatal gamma-50 and gamma-80 tend to be slightly anticorrelated in time. Spectrograms can also be used to obtain a time series of power in a certain band over time, similar to the filtering-based approach in the previous module. Such a time series can then for instance be used to determine if power in a certain frequency band is related to a behavioral or task variable such as running speed. Another common spectrogram-based analysis involves aligning the spectrogram to task events of interest, such as a button press, stimulus onset, or receipt of a reinforcer. To do this, we first need to extract the times of these events from the data files. === Processing Neuralynx events files (.nev) === Apart from neural data stored in ''.ncs'' and ''.ntt'' files, the Neuralynx system also stores //events// in a file called ''Events.Nev''. Let's load it to see what is inside: <code matlab> fn = FindFile('*Events.nev'); [EVTimeStamps, EventIDs, TTLs, EVExtras, EventStrings, EVHeader] = Nlx2MatEV(fn,[1 1 1 1 1],1,1,[]); </code> If you look at the ''EventStrings'' variable, you will see something like: <code matlab> 'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).' '1 or 5 pellet cue' 'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0006).' 'Feeder 1 nosepoke' 'TTL Output on AcqSystem1_0 board 0 port 0 value (0x0002).' '1 pellet dispensed' 'TTL Output on AcqSystem1_0 board 0 port 0 value (0x0000).' 'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).' 'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0006).' 'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).' </code> Each //cell// of ''EventStrings'' contains a string that corresponds to an event. Events are stored based on a number of different triggers, including: * Any change that is detected on the input/output (I/O) ports of the system; these ports interface with components of the experimental setup such as reward pellet dispensers, light cues, levers, photosensors, and others * System events, such as "Start Recording", "Data lost", et cetera * Manual keyboard input into the Cheetah software * Events sent from MATLAB In the above example, the "TTL Input..." and "TTL Output..." entries are all automatically generated because an input was detected (e.g. a photobeam break or lever press) or an output was sent (e.g. dispense a reward pellet). Messages such as ''Feeder 1 nosepoke'' were sent from the MATLAB script that controls the experiment. If you look through the ''EventStrings'' you will also find some system messages. The ''EventStrings'' variable has the same size as ''EVTimeStamps'' which contains the corresponding times for each event. This way, we can line up the events with the neural data. === An example getEvents.m file === Typically we are interested in a specific set of events, such as all times of reward receipt, and want to extract these from the Events file. For the current task, the ''getEvents_value.m'' function extracts a number of events of interest: <code matlab> >> evt = getEvents_value evt = c1: [1x35 double] c3: [1x19 double] c5: [1x34 double] d1: [1x40 double] d3: [1x25 double] d5: [1x42 double] n0: [1x58 double] n1: [1x57 double] </code> As you can see, this function returns a struct with a number of fields. ''c_x'' are the times of three different audio cues, indicating the availability of 1, 3 and 5 food pellets; ''d_x'' are the times at which the corresponding number of pellets were dispensed, and ''n_x'' are the times at which the rat nosepoked in receptacle 0 (on the left end of a linear track) or 1 (on the right). The events returned by ''getEvents_value()'' are already present in the Events file; this function merely presents them in a more usable format. However, for many tasks, some more processing is needed to obtain event times of interest, for instance "only those nosepokes that follow presentation of a cue predicting 5 pellets, and for which the rat does not un-poke for at least one second". Creating appropriate "trial functions" that return usable event times is unglamorous, but a key analysis step for any data set. (Remember: "Garbage in, garbage out"!) In any case, armed with some event times, we are now in a position to ask questions of the sort, "does reward delivery elicit changes in the ventral striatal LFP?" and "is the ventral striatal LFP sensitive to reward-predictive cues?" To do so, we will use the FieldTrip toolbox, which massively expands MATLAB's suite of built-in functions for spectral analysis. ==== Event-triggered spectrograms using FieldTrip ==== [[http://fieldtrip.fcdonders.nl/ | FieldTrip]] is a toolbox originally developed for the analysis of human EEG and MEG data. Because these settings have historically faced complex signal processing challenges, such as eyeblink artifact removal, source reconstruction, and mapping across many electrodes, the methods developed to address these challenges are mature and have much to offer to rodent electrophysiologists. The FieldTrip toolbox continues to be actively developed and maintained, primarily at the [[http://www.ru.nl/donders/ | Donders Institute for Brain, Cognition and Behaviour]] in Nijmegen. It is supported by a popular [[http://fieldtrip.fcdonders.nl/discussion_list | mailing list]], extensive documentation with [[http://fieldtrip.fcdonders.nl/tutorial | tutorials]], and various workshops organized across the globe (such as [[http://fieldtrip.fcdonders.nl/workshop/toronto | this one]] in Toronto). === FieldTrip setup === ☛ Download FieldTrip from the lab database using the BIOL680 account. For this course, do not use the official version from the FieldTrip website: we currently run a custom version modified to allow the processing of Neuralynx data. Unzip the archive in a suitable code subfolder -- do not mix code with your data! ☛ Next, create a new MATLAB shortcut that is a copy of your existing course shortcut, but with the addition of the fieldtrip folder (and all its subfolders) to the path in addition to the GitHub codebase. I have named this shortcut ''mvdmlab-ft''. Click this shortcut to add FieldTrip to your path. If you can type ''ft_defaults'' and not get errors, the installation is successful. === Loading Neuralynx data into FieldTrip === Let's load our familiar ventral striatal LFP into FieldTrip: <code matlab> %% remember to cd to your data folder fc = {'R016-2012-10-03-CSC04a.ncs'}; data = ft_read_neuralynx_interp(fc); </code> You will get a number of warnings about discontinuous recording and expected sample mismatches. This is because FieldTrip expects an exactly constant ''Fs'' in the data, which we know Neuralynx typically does not provide. Next, let's look at the contents of the ''data'' structure: <code matlab> >> data data = hdr: [1x1 struct] label: {'R016-2012-10-03-CSC04a'} time: {[1x5256982 double]} trial: {[1x5256982 double]} fsample: 2000 sampleinfo: [1 5256982] cfg: [1x1 struct] </code> A quirk of FieldTrip is that it treats all data as organized into //trials//, a legacy of its event-related potential roots in human work. However, a quick ''plot(data.time{1},data.trial{1})'' should assure you that all the data is there. At this point, FieldTrip thinks of this full recording session as one single "trial", the first cell of ''data.label'', ''data.time'', and ''data.trial''. === Creating trial data in FieldTrip === To convert this single "trial" containing the entire recording into proper trials surrounding events of interest, FieldTrip uses a "trialfun": a function that returns the times of the events of interest. These times, along with some parameters like how much time around the events to include, are then passed to the ''ft_redefinetrial()'' function, which cuts up the data into trials: <code matlab> cfg = []; cfg.trialfun = 'ft_trialfun_lineartracktone2'; cfg.trialdef.hdr = data.hdr; cfg.trialdef.pre = 2.5; cfg.trialdef.post = 5; cfg.trialdef.eventtype = 'nosepoke'; % could be 'nosepoke', 'reward', 'cue' cfg.trialdef.location = 'both'; % could be 'left', 'right', 'both' cfg.trialdef.block = 'both'; % could be 'value', 'risk' cfg.trialdef.cue = {'c1','c3','c5'}; % cell array with choice of elements {'c1','c3','c5','lo','hi'} [trl, event] = ft_trialfun_lineartracktone2(cfg); cfg.trl = trl; data_trl = ft_redefinetrial(cfg,data); </code> The above reorganizes the data into 7-second windows (from -2 to 5 seconds) around nosepoke times following the 1- 3- and 5-pellet cues. The function ''ft_trialfun_lineartracktone2()'' is a more powerful version of ''getEvents_value()'' we encountered above: it can accept a variety of requests for particular event times. For instance, if we only wanted nosepokes on the left side of the maze that followed the 3-pellet cue, we can specify this by changing 'both' to 'left' and only including 'c3'. ☛ Inspect the ''data_trl'' structure. Notice that there are now 88 trials instead of 1, and that each trial has its own time and data vectors. Plot the 47th trial as an example to verify this. === Constructing the event-triggered spectrogram === Now we are ready to construct the event-triggered spectrogram: <code matlab> cfg = []; % start with empty cfg cfg.output = 'pow'; cfg.channel = 'R016-2012-10-03-CSC04a'; cfg.method = 'mtmconvol'; cfg.taper = 'hanning'; cfg.foi = 10:2:100; % frequencies of interest cfg.t_ftimwin = ones(size(cfg.foi)).*0.5; % window size: fixed at 0.5s cfg.toi = -1:0.05:4; % times of interest TFR = ft_freqanalysis(cfg, data_trl); figure cfg = []; cfg.channel = 'R016-2012-10-03-CSC04a'; ft_singleplotTFR(cfg, TFR); </code> That is quite a mouthful. Almost any FieldTrip function requires a ''cfg'' struct as input that specifies the parameters of the task to be completed. To avoid unexpected results, it is good practice to start with an empty cfg (''cfg = [];'') when moving on to a new type of operation, such as from trial segmentation to spectrogram construction in the above example. The function ''ft_freqanalysis()'' computes the spectrogram based on the data and the parameters included in the ''cfg'' file. ''ft_singleplotTFR()'' plots the result, again based on the information in the ''cfg'' file and the output of ''ft_freqanalysis()''. You should get: {{ :analysis:course:week6_fig4.png?600 |}} This plot is an //average// spectrogram over 88 trials, with time 0 corresponding to the time of nosepoking into the reward receptacle. As with the examples using ''spectrogram()'', we can control the smoothness of the spectrogram by specifying the time and frequency steps. For the frequency steps this is accomplished the same way in FieldTrip, but for the time axis we can now simply specify a vector with the window centers instead of a window size and overlap. Looking at the results, it appears there are some beta (~20Hz) and gamma events visible in the average spectrogram. ☛ Change the ''cfg.toi'' and ''cfg.foi'' parameters to verify you can change the resolution of the time and frequency axes. What happens if you try to extend the time axis from -2 to 5? === Options and parameters for FieldTrip spectrograms === FieldTrip provides many powerful tools for spectral analysis. This section highlights a few useful ones, but the [[http://fieldtrip.fcdonders.nl/tutorial/timefrequencyanalysis | tutorials]] and documentation of the functions involved (e.g. [[http://fieldtrip.fcdonders.nl/reference/ft_freqanalysis | ft_freqanalysis]]) give a more complete overview. == Baseline correction == Often we are not interested in absolute power levels, as provided per default by a spectrogram. Rather, we wish to compare spectral content following some event of interest to a baseline. This is easily accomplished: <code matlab> figure cfg = []; cfg.baseline = [-2 0]; cfg.baselinetype = 'relative'; cfg.channel = 'R016-2012-10-03-CSC04a'; ft_singleplotTFR(cfg, TFR); </code> Note that the beta and gamma events now stand out clearly. == Frequency-dependent windowing == Up until now, we have used a //fixed// window size (number of samples, or equivalently, time) in computing spectrograms. The above example specified a 0.5s window. ☛ Recompute and plot the event-triggered spectrogram using a //frequency-dependent// window of 20 cycles by setting ''cfg.t_ftimwin = 20./cfg.foi;'' and a frequency range of 1:100. You should now have a beautiful image that balances time and frequency resolution: {{ :analysis:course:week6_fig5.png?600 |}} ☛ What causes the whitespace at the bottom of the image? This frequency-dependent windowing can be thought of as computing many spectrograms with different windows in parallel, and then keeping different pieces of each so that low frequencies are estimated with good precision in the frequency domain (but poor temporal precision) and high frequencies are resolved from good temporal precision (but relatively poor frequency precision). This tradeoff is illustrated by the time-frequency tiles in the "B" panel of the figure earlier in this module. (As a technical note, this frequency-dependent windowing is similar to what can be obtained with wavelet transforms.) == Statistical tests == To set up a statistical comparison of the pre-nosepoke baseline and the oscillation patterns following the nosepoke, we perform two different frequency analyses, with different times of interest: <code matlab> cfg = []; cfg.output = 'pow'; cfg.channel = 'R016-2012-10-03-CSC04a'; cfg.method = 'mtmconvol'; cfg.taper = 'hanning'; cfg.foi = 1:1:100; cfg.keeptrials = 'yes'; % need this for stats later cfg.t_ftimwin = 20./cfg.foi; % 20 cycles per time window cfg.toi = -2:0.05:0; % pre-nosepoke baseline TFR_pre = ft_freqanalysis(cfg, data_trl); cfg.toi = 0:0.05:2; % post-nosepoke TFR_post = ft_freqanalysis(cfg, data_trl); TFR_pre.time = TFR_post.time; % time axes should be identical for comparison </code> Now we run ''ft_freqstatistics()'' to do the comparison between the two frequency analyses: <code matlab> %% t-test cfg = []; cfg.channel = 'R016-2012-10-03-CSC04a'; cfg.latency = 'all'; cfg.trials = 'all'; cfg.frequency = 'all'; cfg.avgoverchan = 'no'; cfg.avgovertime = 'no'; cfg.avgoverfreq = 'no'; cfg.parameter = 'powspctrm'; cfg.method = 'stats'; cfg.statistic = 'ttest2'; cfg.alpha = 0.05; nTrials1 = size(TFR_pre.powspctrm,1); nTrials2 = size(TFR_post.powspctrm,1); cfg.design = cat(2,ones(1,nTrials1),2*ones(1,nTrials2)); % two conditions cfg.ivar = 1; % dimension of design var which contains the independent variable (group) stat = ft_freqstatistics(cfg,TFR_post,TFR_pre); cfg.parameter = 'stat'; ft_singleplotTFR(cfg,stat); % plot the t-statistic </code> The result shows the t-statistic of the bin-by-bin comparison between pre- and post-nosepoke power. Depending on what correction for multiple comparisons is to be applied, some time-frequency bins may approach significance. This is however only a single recording session. The Statistics section of [[http://fieldtrip.fcdonders.nl/walkthrough | this walkthrough]] has some background on the different steps performed here; the [[http://fieldtrip.fcdonders.nl/tutorial/eventrelatedstatistics || event-related statistics]] tutorial has further examples of different statistics. == Multichannel data == FieldTrip excels at processing many data channels in parallel. Processing all 16 channels from this session is as simple for the user as processing 1 channel. (Be aware however, that this will be mild torture for older computers; the following took about 5 minutes to run on my 3-year old laptop.) <code matlab> fc = FindFiles('*.ncs'); % get filenames of all LFPs recorded data = ft_read_neuralynx_interp(fc); % load them all -- this will take a while % define layout cfg = []; cfg.layout = 'ordered'; cfg.channel = data.label; layout = ft_prepare_layout(cfg, data); </code> Defining a channel layout is the only addition when dealing with multichannel data. Here, we are using one of the default layouts ('ordered'). However, layouts specific to particular EEG/MEG and intracranial probes are also available, so that results can be plotted topographically (i.e. matching the physical sensor layout). Remarkably, the trial definition and spectrogram computations are identical to those for a single channel: <code matlab> cfg = []; cfg.trialfun = 'ft_trialfun_lineartracktone2'; cfg.trialdef.hdr = data.hdr; cfg.trialdef.pre = 2.5; cfg.trialdef.post = 5; cfg.trialdef.eventtype = 'nosepoke'; % could be 'nosepoke', 'reward', 'cue' cfg.trialdef.location = 'both'; % could be 'left', 'right', 'both' cfg.trialdef.block = 'both'; % could be 'value', 'risk' cfg.trialdef.cue = {'c1','c3','c5'}; % cell array with choice of elements {'c1','c3','c5','lo','hi'} [trl, event] = ft_trialfun_lineartracktone2(cfg); cfg.trl = trl; data_trl = ft_redefinetrial(cfg,data); %% cfg = []; cfg.output = 'pow'; %cfg.channel = 'R016-2012-10-03-CSC04a'; cfg.method = 'mtmconvol'; cfg.taper = 'hanning'; cfg.foi = 1:1:100; cfg.keeptrials = 'yes'; % should be needed for subsequent statistics... cfg.t_ftimwin = 20./cfg.foi; % 20 cycles per time window cfg.toi = -1:0.05:4; TFR = ft_freqanalysis(cfg, data_trl); </code> Only the final plot command is slightly different: <code matlab> figure cfg = []; cfg.baseline = [-2 0]; cfg.baselinetype = 'relative'; cfg.layout = layout; ft_multiplotTFR(cfg, TFR); % note this is now multiplot rather than singleplot </code> You should get: {{ :analysis:course:week6_fig6.png?600 |}} This shows the baseline-corrected, event-aligned spectrograms for all 16 channels in this session. There is also an average shown (the subplot on the lower right). Note that a subset of the channels, presumably located in the ventral striatum, show a very similar time-frequency pattern, whereas another subset do not show this at all. The channels that do not show this same pattern are likely located in the hippocampus. ==== Assignment ==== Because the event-triggered spectrograms above show an average over trials, it is possible to lose touch with the properties of the raw data that generated it. For instance, the average may arise from a distribution of very similar looking individual trials, or it may be dominated by one atypical but extreme trial. For this reason it is important to compare the spectrogram to the raw data. An example of a study where this was done well is [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/ | Leventhal et al]]. Neuron 2012. [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/figure/F3/ | Figure 3]] in this paper shows a number of event-triggered LFP traces, one for each trial. ☛ Implement a function eventLFPplot(csc,event_times,varargin) to make such a plot, as follows: <code matlab> % function eventLFPplot(csc,event_times,varargin) % % INPUTS % % csc: [1 x 1] mytsd, LFP signal to be plotted % event_times: [nEvents x 1] double with event times to align LFP on % % varargins (with defaults): % % t_window: [2 x 1] double indicating time window to use, e.g. [-1 3] for 1 second before to 3 seconds after event times </code> You can use a skeleton like the following: <code matlab> % for each event: % *extract the corresponding piece of LFP % *replace the original time axis with a new one based on the time window asked for % *optionally, rescale the LFP % *add a y-offset to the LFP to plot one above the other % *plot the current piece of LFP % add a vertical line to the plot indicating time zero % further niceties: add an option to decimate, to color, to filter before plotting </code> Test your code with the same events used to generate the event-triggered spectrogram above. Do you think the average is a reasonable charaterization of what is happening in the LFP for these events?