Frequency Domain

Frequency

Let’s create a signal with amplitude of 1 V and a frequency of 10 Hz:

    # generate data
    data = create_data(n_chan=2, signal='sine', amplitude=1)

    traces = [
        go.Scatter(
            x=data.time[0],
            y=data(trial=0, chan='chan00'))
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Time (s)'),
        yaxis=dict(
            title='Amplitude (V)'),

We can transform it to the frequency-domain, by computing the power spectral density (default option):

    save_plotly_fig(fig, 'freq_01_data')

    # default options
    freq = frequency(data, detrend=None)

    traces = [
        go.Scatter(
            x=freq.freq[0],
            y=freq(trial=0, chan='chan00'))
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Frequency (Hz)',
            range=(0, 20)),
        yaxis=dict(

Note that the Parseval’s theorem holds:

    p_time = math(data, operator_name=('square', 'sum'), axis='time')
    p_freq = math(freq, operator_name='sum', axis='freq')
    assert_array_almost_equal(p_time(trial=0), p_freq(trial=0) * data.s_freq)

Welch’s Method

If you want to apply the Welch’s method (compute PSD on 1s long, 50% overlapping windows and Hann window) use:

    # generate very long data
    data = create_data(n_chan=1, signal='sine', time=(0, 100))
    freq = frequency(data, taper='hann', duration=1, overlap=0.5)

    traces = [
        go.Scatter(
            x=freq.freq[0],
            y=freq(trial=0, chan='chan00')),
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Frequency (Hz)',
            range=(0, 20)),
        yaxis=dict(
            title='Amplitude (V<sup>2</sup>/Hz)'),
        )

Multitapers

A common approach is to use multitaper to suppress the activity outside the frequency band of interest and to smooth the frequency band of interest:

    # dpss
    data = create_data(n_chan=1, signal='sine')
    freq = frequency(data, taper='dpss', halfbandwidth=5)

    traces = [
        go.Scatter(
            x=freq.freq[0],
            y=freq(trial=0, chan='chan00')),
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Frequency (Hz)',
            range=(0, 20)),
        yaxis=dict(
            title='Amplitude (V<sup>2</sup>/Hz)'),
        )

Note

You can either specify the half-bandwidth smoothing in the frequency domain (halfbandwidth) or the normalized halfbandwidth (NW), where:

NW = halfbandwidth * duration

Energy Spectral Density

All the above examples were transformed to the Power Spectral Density (PSD). PSD is appropriate for signals that are (roughly) periodic. However, if your signal is limited in the time domain, it makes sense to compute the Energy Spectral Density (ESD, see wikipedia for discussion).

So, for a signal that has a clear start and end:

    DURATION = 2
    data = create_data(n_chan=1, signal='sine', time=(0, DURATION))
    data.data[0][0, :] *= hann(data.data[0].shape[1])

    traces = [
        go.Scatter(
            x=data.time[0],
            y=data(trial=0, chan='chan00'))
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Time (s)'),
        yaxis=dict(
            title='Amplitude (V)'),
        )
    fig = go.Figure(data=traces, layout=layout)

You can compute the ESD in this way:

    freq = frequency(data, detrend=None, scaling='energy')

    traces = [
        go.Scatter(
            x=freq.freq[0],
            y=freq(trial=0, chan='chan00'))
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Frequency (Hz)',
            range=(0, 20)),
        yaxis=dict(
            title='Amplitude (V<sup>2</sup>)'),
        )
    fig = go.Figure(data=traces, layout=layout)

Note

The units for the PSD are V²/ Hz while those for the ESD are V².

The Parseval’s theorem holds in this case as well, but we need to make sure to include the duration as well:

    p_time = math(data, operator_name=('square', 'sum'), axis='time')
    p_freq = math(freq, operator_name='sum', axis='freq')
    assert_array_almost_equal(p_time(trial=0), p_freq(trial=0) * data.s_freq * DURATION)

Complex Output

To get the complex output (so, not the PSD or ESD), you should pass the argument output='complex'.

    data = create_data(n_chan=1, signal='sine')
    freq = frequency(data, output='complex', sides='two', scaling='energy')

    traces = [
        go.Scatter(
            x=freq.freq[0],
            y=abs(freq(trial=0, chan='chan00', taper=0)))
        ]
    layout = go.Layout(
        xaxis=dict(
            title='Frequency (Hz)'
            ),
        yaxis=dict(
            title='Amplitude (V)'),
        )
    fig = go.Figure(data=traces, layout=layout)

Note

The complex frequency domain returns positive and negative frequencies.

Note

The complex output has an additional dimension, called taper. This is necessary in the case of taper='dpss' because it does not make sense to average in the complex domain over multiple tapers. For consistency, there is also a taper dimension also when using single tapers (such as hann or boxcar).

Comparison to Scipy

The code reflects quite some ideas from scipy, but the nomenclature in scipy might give rise to some confusion, especially if compared to other sources (well summarized by wikipedia).

Name	Units	scipy parameters	wonambi parameters
Power Spectral Density	V²/ Hz	`output='psd'` `scaling='density'`	`output='spectraldensity'` `scaling='power'`
Energy Spectral Density	V²	`output='psd'` `scaling='spectrum'`	`output='spectraldensity'` `scaling='energy'`
Complex Fourier	V	`output='complex'` `scaling='spectrum'`	`output='complex'` `scaling='energy'`

Time-Frequency

There are two main approaches to the time-frequency analysis:

Spectrogram / short-time Fourier transform

Morlet wavelets

Spectrogram / STFT

The first approach is identical to computing frequency() on small epochs. The duration of the epochs is defined by duration, and you can specify either the overlap (between 0, no overlap, and 1, complete overlap) or the step (distance between epochs, in seconds). The output of timefrequency() has a different name than the output of frequency() for consistency with the literature, so:

timefrequency()	frequency()
`output='spectrogram'`	`output='spectraldensity'`
`output='stft'`	`output='complex'`

The other arguments of timefrequency() are the same as frequency().

Wavelet

Morlet Wavelets can be used for time-frequency analysis and have an intuitive tradeoff between time and frequency resolution. The catch is there is no straightforward way to normalize the output of the wavelets. The Parseval’s theorem does not hold because wavelets are not orthogonal.