# Customizing common M/EEG plots part 2: the time-frequency representation (TFR)

This is the second installment (out of probably three) of blog posts based on the data visualization workshop I gave for the Language in Interaction consortium. (Part 1 was on plotting ERP/Fs, or 1D data over time with uncertainty.). This time I’ll go over some things you can do to improve the appearance of time-frequency representations (TFRs), often of oscillatory power. The customizations I’m describing here will be a bit more concerned with aesthetics and personal preference than with actually “showing the data” (like in part 1), but I still think/hope this will be useful to some.

For the same experiment as in part 1, we found a nice gamma band response after visual stimulation. Set up the path and load the results; available for download here (note that only some interesting channels are present in the result file, to speed up loading it in):

ft_defaults();

% Also note that trials are still present in the freq structure; these are
% needed later on.

Let’s start by plotting it in the FieldTrip standard way:

figure();
cfg = [];
ft_singleplotTFR(cfg, freq); Clearly we need to make some changes. While we can do most of these things by simply changing FieldTrip configuration options (and using FieldTrip functions is great for interactively exploring the data), for publication-quality figures it can be useful to use default Matlab plotting options. (Just as we did for the ERF example.) We will let FieldTrip do a baseline correction on the data, though:

% First compute the average over trials:
cfg = [];
freq_avg = ft_freqdescriptives(cfg, freq);

% And baseline-correct the average:
cfg = [];
cfg.baseline = [-0.8 0];
cfg.baselinetype = 'db'; % Use decibel contrast here
freq_avg_bsl = ft_freqbaseline(cfg, freq_avg);

Compute the mean over channels of the baseline-corrected powerspectrum (recall that we only have channels of interest here, so we don’t need to select those first):

meanpow = squeeze(mean(freq_avg_bsl.powspctrm, 1));

And plot it:

figure();
imagesc(freq.time, freq.freq, meanpow);
xlim([-0.7 0.8]); % to hide the NaN values at the edges of the TFR
axis xy; % to have the y-axis extend upwards

Label the principal axes:

xlabel('Time (s)');
ylabel('Frequency (Hz)'); Set colour limits symmetric around zero:

clim = max(abs(meanpow(:)));
caxis([-clim clim]); The default colourmap that was introduced in Matlab R2014b (parula) is infinitely better than the old default (jet), but we can do better than this by using a properly diverging colour map. Parula was created such that it would be reasonably OK for both diverging and linearly increasing data, but since we know our data is diverging (positive and negative deflections of comparable amplitude, at least in principle), we can choose a colour map accordingly. We will use one of the excellent ColorBrewer colour maps, available from GitHub or the Matlab FileExchange.

Use a fancy diverging colour map, add a colour bar, and label it:

addpath('thirdparty\brewermap'); % or wherever you put brewermap

colormap(brewermap(256, '*RdYlBu'));
h = colorbar();
ylabel(h, 'Power vs baseline (dB)'); Looks pretty decent, but still a bit bland, right? In particular, it might be nice to have it look less ‘blocky’. To achieve this, we could have used more finely spaced time and frequency axes in the estimation of our TFR. However, this increases computational time quite a bit. Moreover, the intrinsic time/frequency resolution is determined by the data, so oversampling at the `ft_freqanalysis` step does not buy us any extra information. For visualization purposes only, therefore, we can simply interpolate the TFR onto a finer grid.

% The finer time and frequency axes:
tim_interp = linspace(-0.7, 0.8, 512);
freq_interp = linspace(20, 100, 512);

% We need to make a full time/frequency grid of both the original and
% interpolated coordinates. Matlab's meshgrid() does this for us:
[tim_grid_orig, freq_grid_orig] = meshgrid(freq.time, freq.freq);
[tim_grid_interp, freq_grid_interp] = meshgrid(tim_interp, freq_interp);

% And interpolate:
pow_interp = interp2(tim_grid_orig, freq_grid_orig, meanpow,...
tim_grid_interp, freq_grid_interp, 'spline');

(You can ignore any warnings about NaNs being ignored here.)

Plot our interpolated TFR:

figure();
imagesc(tim_interp, freq_interp, pow_interp);
xlim([-0.7 0.8]);
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
clim = max(abs(meanpow(:)));
caxis([-clim clim]);
colormap(brewermap(256, '*RdYlBu'));
h = colorbar();
ylabel(h, 'Power vs baseline (dB)');

% Let's also add a line at t = 0s for clarity:
hold on;
plot(zeros(size(freq_interp)), freq_interp, 'k:'); Now things are starting to look nice. At least I quite like the smoother appearance of the interpolated TFR, better than the blocky original. I realize that this is definitely a matter of personal preference; not everyone will agree with me here. The interpolation emphasizes that power is evolving smoothly over time; the discreteness suggested by the original TFR is merely a consequence of how we chose to estimate power (i.e., with (partially overlapping) bins sliding over the data). But like I said, feel free to disagree here and continue to use the non-interpolated plots!

We’ll add an extra improvement: marginal plots for the TFR, so that power as a function of a time and frequency of interest is emphasized. This requires some pretty heavy customization.

First, get the data we actually want to draw in the marginal plots:

time_of_interest = 0.3;
freq_of_interest = 50;
timind = nearest(tim_interp, time_of_interest);
freqind = nearest(freq_interp, freq_of_interest);
pow_at_toi = pow_interp(:,timind);
pow_at_foi = pow_interp(freqind,:);

Initialize the axes for our plots before we start drawing:

figure();
ax_main = axes('Position', [0.1 0.2 0.55 0.55]);
ax_right = axes('Position', [0.7 0.2 0.1 0.55]);
ax_top = axes('Position', [0.1 0.8 0.55 0.1]); Draw the TFR in the main axes:

axes(ax_main);
im_main = imagesc(tim_interp, freq_interp, pow_interp);
% note we're storing a handle to the image im_main, needed later on
xlim([-0.7 0.8]);
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
clim = max(abs(meanpow(:)));
caxis([-clim clim]);
colormap(brewermap(256, '*RdYlBu'));
hold on;
plot(zeros(size(freq_interp)), freq_interp, 'k:'); Draw the power over time:

axes(ax_top);
area(tim_interp, pow_at_foi,...
'EdgeColor', 'none', 'FaceColor', [0.5 0.5 0.5]);
xlim([-0.7 0.8]);
ylim([0 clim]);
box off;
ax_top.XTickLabel = [];
ylabel('Power (dB)');
hold on;
plot([0 0], [0 clim], 'k:'); Draw the power over frequency:

axes(ax_right);
area(freq_interp, pow_at_toi,...
'EdgeColor', 'none', 'FaceColor', [0.5 0.5 0.5]);
view([270 90]); % this rotates the plot
ax_right.YDir = 'reverse';
ylim([0 clim]);
box off;
ax_right.XTickLabel = [];
ylabel('Power (dB)'); Add a colour bar at a pre-specified position:

h = colorbar(ax_main, 'manual', 'Position', [0.85 0.2 0.05 0.55]);
ylabel(h, 'Power vs baseline (dB)'); It might be nice to additionally indicate which times and frequencies of interest the marginal plots were created from. We can do this via lines in the main and marginal plots:

% Main plot:
axes(ax_main);
plot(ones(size(freq_interp))*time_of_interest, freq_interp,...
'Color', [0 0 0 0.1], 'LineWidth', 3);
plot(tim_interp, ones(size(tim_interp))*freq_of_interest,...
'Color', [0 0 0 0.1], 'LineWidth', 3);

% Marginals:
axes(ax_top);
plot([time_of_interest time_of_interest], [0 clim],...
'Color', [0 0 0 0.1], 'LineWidth', 3);
axes(ax_right);
hold on;
plot([freq_of_interest freq_of_interest], [0 clim],...
'Color', [0 0 0 0.1], 'LineWidth', 3); Now it’s time to add some uncertainty to the plot. You have by now probably gotten used to seeing error bars on line and bar graphs (1D data), and are surprised when they are missing (hopefully). But for 2D data, indications of uncertainty are only rarely presented. One option is to use colour saturation or transparency as an indication of spread. We’ll be using the latter here, encoding t-scores into the transparency channel. T-scores are a combination of effect strength (versus baseline, in this case) and data spread, which is what we want in this case.

(Question: why would using spread only, and no effect strength, be inappropriate for determining the transparency? This is in opposition to plotting error bars around curves.)

Compute single-trial power versus baseline first:

cfg = [];
cfg.baseline = [-0.8 0];
cfg.baselinetype = 'absolute';
freq_bsl = ft_freqbaseline(cfg, freq);

Then, compute mean and standard error over trials:

cfg = [];
cfg.variance = 'yes';
freq_sem = ft_freqdescriptives(cfg, freq_bsl);

And compute a t-score from that:

tscore = freq_sem.powspctrm ./ freq_sem.powspctrmsem;

% Average the t-score over our channels:
tscore = squeeze(mean(tscore, 1));

Plot the t-scores, to get an idea of what we’re dealing with:

figure();
imagesc(freq.time, freq.freq, tscore);
axis xy;
colorbar(); Since we want to keep using our nicely interpolated plot, we will have to interpolate the t-scores as well:

tscore_interp = interp2(tim_grid_orig, freq_grid_orig, tscore,...
tim_grid_interp, freq_grid_interp, 'spline');

Now convert our map of t-values into transparency scores. 0 will be fully transparent, 1 will be fully opaque. As a cutoff point for fully opaque, we will use a critical t-value based on a 99% confidence interval:

alpha = 0.01;
tcrit = tinv(1-alpha/2, size(freq.powspctrm, 1)-1);

And linearly ramp up opacity with absolute t-score:

opacity = abs(tscore_interp) / tcrit;
opacity(opacity > 1) = 1;

Set the opacity map like so: One final note on the marginal plots I chose. Although I like the aesthetic appearance of the area curves, they have some drawbacks. First of all, I would only use area curves for data that is strictly positive or negative. This is more or less the case in this example (positive), but may in general not be the case. Second, it’s difficult to add measures of spread or uncertainty. You can fix this by using `boundedline()` instead, as I did in the post about plotting 1D data over time with uncertainty.