Calculating time-frequency representations over power using multitapers

From chbhwiki
Jump to: navigation, search

Developed by Ole Jensen August 19, 2018


The following Matlab exercise covers how to calculate time-frequency representations of power using a simple time-window FFT approach in combination with multitapers. The consequences of using multitapers will be discussed in terms of controlling the spectral smoothing. Please go through the tutorial step-by-step and answer the question. The exercises will help you to work through and adapt the scripts. If this is an assignment please return the answers to the exercise with relevant figure

The test signal

The sampling frequency is set to 1024 Hz:

Fs = 1024;

Define the test signal to be analyzed. Following a flat signal of 0.5 s a broadband ‘gamma band’ signal at 70 – 90 Hz turning on. This signal somehow mimics human visual gamma band activity in response to a grating. The gamma signal is created by filtering white noise in the 70 – 90 Hz band using a simple Butterworth filter:

t = (1:1024)/Fs;
s = [zeros(1,512) gammas(1,256:767)];

Add a 50 Hz line noise signal and some Gaussian noise:

s = s + 0.5*sin(2*pi*50*t) + 0.2*randn(1,1024) ;

Then plot the signal:


The sliding time-window

The first step it to drag a time window over the data. This is done by selecting 512 points epocs of data around (0.5 s) each sample points. For instance, for sample point 300 we will select the epoch from 300-256 = 44 to 300+256-1 = 555. The boundaries create problems. In this case we will set the existing value to zero (zero-padding). E.g. for sample 1, the 255 values prior to the sample will be considered 0. The following script with cut out the epochs; one epoch for each sample point (note that the script uses several ‘for loops’. In general this inefficient Matlab code; however, it is done for illustrative purposes).

Each epoch (time window) will be 512 points long, i.e. 0.5 s long:

Nwindow = 512;

Make a matrix to save computation time:

swindow = zeros(Nwindow,length(s));

Now loop over the sample points and cut out the epochs of interest.

for k=1:length(s) % loop over all the sample points
    nstart = k - Nwindow/2; % Define the start of time window k
    nstop = k + Nwindow/2 - 1; % Define the end of time window k
    for m=nstart:nstop
        n = n + 1;
        if m < 1 || m > length(s)
            swindow(n,k) = 0; % Zero-pad at the boundary
            swindow(n,k) = s(m); % Copy the data to the time window

Time-frequency representations of power using conventional Hanning tapers

First perform a ‘conventional’ time-frequency presentation analysis of power using Hanning tapers. This is done by multiplying each epoch to a Hanning taper of the same length prior to performing the Fourier transform.

Prior to using Matlabs definition of the Hanning taper it needs to be normalized such that ∑hannk2 = 1 :

hann = hanning(Nfft); 
hann = hann/ sqrt(sum(hann.^2));

Nfft = Nwindow; 
for k=1:length(s)
    Sfft = fft(hann.*swindow(:,k),Nfft)/(Nfft/2);
    Sfft = Sfft(1:Nfft/2+1);
    Phann(:,k) = abs(Sfft).^2;
f = (Fs/Nfft)*(0:Nfft/2);
axis xy; % flip the axis for convenience
ylim([0 100]) % show from 0 to 100 Hz

This procedure generated the time-frequency representation of power for the trial.

Exercise 1:

  • Is the 60 – 90 Hz ‘gamma band’ signal well represented?
  • How can the gamma power representation be improved?
  • Is the 50 Hz ‘line noise’ well represented?

Applying multitapers when calculating time-frequency representations of power

The short-coming of the Hanning tapers can be alleviated using multitapers. Instead of using one taper (e.g. one Hanning taper), multiple orthogonal tapers are applied. We will make use of discrete prolate spheroidal sequences (DPSS) for creating the tapers. The Slepian tapers will serve to smooth the power estimates in the frequency domain.

To generate 6 Slepian tapers write

Nfft = 512; 
[dps,lambda] = dpss(Nfft,3); 

The 6 tapers will increase the spectal smoothing according to K = 2*ΔT*ΔFs-1

Since ΔT = 0.5 and K = 6 the smoothing is about ±7 Hz.

Exercise 2:

  • Plot the tapers one by one
  • Check that the tapers are orthogonal

The tapers are applied by multiplying them one by one to each segment. After the multiplication, the Fourier transform is calculated. The power is then derived.

for k=1:length(s)
    for l = 1:size(dps,2)  % Loop over the 6 tapers
        Pt(:,l) = (abs(fft(dps(:,l).*swindow(:,k),Nfft)/(Nfft/2))).^2;
    Pmult(:,k) = mean(Pt(1:Nfft/2+1,:),2);
f = (Fs/Nfft)*(0:Nfft/2);
axis xy; % flip the axis for conveniency
ylim([0 100]) % only show 0 to 100 Hz

Exercise 3:

  • How did the multi-tapers change the estimate of the 60 – 90 Hz gamma?
  • What happened to the 50 Hz line noise?
  • Interpret that maximum power estimate and compared to the Hanning taper estimate. Why is it reduced?

Exercise 4:

  • Increase the number of tapers to 8.
  • Estimate the spectral smoothing (ΔF)
  • How does is change the spectral estimate?

Visualizing the spectral smoothing from tapering

To study the effect of the tapers we will consider their effects in the frequency domain. This can be done qualitatively by multiplying the function cos(2πf t) + i sin(2πf t) to the taper(s) and then study the power spectra.

We will investigate of the spectral influence of the Hanning taper at 80 Hz (note that 80 Hz is arbitrary and a similar calculation could be done for any frequency producing similar results in terms spectral smoothing):

Nfft = 512; 
tw = (1:Nwindow)/Fs;
sw = cos(2*pi*80*tw) + i*sin(2*pi*80*tw);
hann = hanning(Nfft); 
hann = hann/ sqrt(sum(hann.^2));

The power spectrum of the above function will provide insight into spectral influence of the Hanning taper. We will calculate the spectrum with a relatively high frequency resolution in order to get a better estimate. To do so zero-pad the signal to 10240 points. This result in 0.1 Hz frequency steps:

Swfft = fft(hann'.*sw,10240)/(10240/2);
Swfft = Swfft(1:10240/2+1);
Pw = abs(Swfft).^2;
fw = (Fs/10240)*(0:10240/2);
plot(fw,Pw);xlim([40 100]);
plot(fw,20*log10(Pw));xlim([40 100]);ylim([-350 -80])
From the power spectrum we can estimate the spectral smoothing with resulting from the Hanning taper. The log- transform allows us to assess the filter properties in decibel in terms of attenuation.  

Exercise 5:

  • What is the spectral smoothing resulting from the 0.5 s Hanning taper?

The next step is to calculate the spectral smoothing resulting from the application of the 6 multitapers. This is done as above by calculating power spectra for the individual tapers multiplied to a sinusoid. This spectra are then averaged.

Nfft = 512; 
[dps,lambda] = dpss(Nfft,3); 
clear Pmw
for k = 1 : size(dps,2)
    Swfft = fft(dps(:,k)'.*sw,10240)/(10240/2);
    Swfft = Swfft(1:10240/2+1);
    Pmw(:,k) = abs(Swfft).^2;
fw = (Fs/10240)*(0:10240/2);
plot(fw,mean(Pmw,2));xlim([40 100]);
plot(fw,20*log10(mean(Pmw,2)));xlim([40 100]);ylim([-200 -80])

Exercise 6:

  • Compare the spectral smoothing between the Hanning taper and the 6 Slepians tapers.
  • What are the advantages and disadvantages of the increased smoothing?

Exercise 7:

  • Shorten the time-window to 256 points and compared the Hanning and multi-taper (6 tapers) approaches (6 tapers).
  • What are the consequences for the temporal resolution?
  • What are the consequences the frequency smearing? – consider both the 60 – 90 Hz and 50 Hz components

Exercise 8:

  • Increase the time-window to 256 points and compare the Hanning and multi-taper (6 tapers).
  • What are the consequences for the temporal resolution?
  • What are the consequences the frequency smearing? – consider both the 60 – 90 Hz and 50 Hz components