Skip to content

Commit

Permalink
DOC: Add raw fft comparison to Welch tutorial
Browse files Browse the repository at this point in the history
Changes to be committed:
	modified:   doc/source/_static/welch_basic_power.png
	new file:   doc/source/_static/welch_raw_fft_power.png
	modified:   doc/source/tutorials/welch.rst
  • Loading branch information
tensionhead committed Dec 23, 2022
1 parent 3f24615 commit d5209dc
Show file tree
Hide file tree
Showing 3 changed files with 31 additions and 7 deletions.
Binary file modified doc/source/_static/welch_basic_power.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/source/_static/welch_raw_fft_power.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
38 changes: 31 additions & 7 deletions doc/source/tutorials/welch.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@ is described in the following publication (`DOI link <https://doi.org/10.1109/TA
a method based on time averaging over short, modified periodograms.
IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73.`

In short, it splits the original time-domain signal into several, potentially overlapping, segments. A taper (window function) is then applied,
and the data is transferred into the frequency domain by computing the FFT. Computing the squared magnitude results in one periodogram per segment.
In short, it splits the original time-domain signal into several, potentially overlapping, segments. A taper (window function) is then applied on each segment individually,
and each segment is transferred into the frequency domain by computing the FFT. Computing the squared magnitude results in one periodogram per segment.
Welch refers to these as *modified periodograms* in the publication, because of the taper that has been applied. These
powers are averaged over the windows to obtain the final estimate of the power spectrum.

Expand Down Expand Up @@ -68,11 +68,10 @@ Let's inspect the resulting `SpectralData` instance by looking at its dimensions
The shape is as expected:

* The `time` axis contains two entries, one per trial, because by default there is no trial averaging (`keeptrials` is `True`). With trial averaging, there would only be a single entry here.
* The `taper` axis will always have size 1 for Welch, even for multi-tapering, as taper averaging must be active for Welch (`keeptapers` must be `False`), as explained in the function documentation.
* The size of the frequency axis (`freq`, 512 here), i.e., the frequency resolution, depends on the signal length of the input windows and is thus a function of the input signal, `t_ftimwin`, `toi`, and potentially other settings (like a `foilim`, i.e. a frequency selection).
* The `taper` axis will always have size 1 for Welch, even for multi-tapering, as tapers must be averaged for Welch's method (`keeptapers` must be `False`), as explained in the function documentation.
* The size of the frequency axis (`freq`, 251 here), i.e., the frequency resolution, depends on the signal length of the input windows and is thus a function of the input signal, `t_ftimwin` and `toi`, and potentially other settings (like a `foilim`, i.e. limiting the frequencies of interest).
* The channels are unchanged, as we receive one result per channel.


We can also visualize the power spectrum. Here we select the first of the two trials:

.. code-block:: python
Expand All @@ -96,8 +95,33 @@ Many settings affect the outcome of a Welch run, including:
* `t_ftimwin` : window length (a.k.a. segment length) in seconds.
* `toi` : overlap between windows, 0.5 = 50 percent overlap.
* `foilim` : frequencies of interest, a specific frequency range, e.g. set to ``[5, 100]`` to get results between 5 to 100 Hz.
* `taper` and `tapsmofrq` : for taper selection and multi-tapering. Note that in case of multi-tapering, the data in the windows will be averaged across the tapers first, then the Welch procedure will run.
* `keeptrials` : whether trials should be left as-is, or you want a trial-average. If `false`, and thus trial-averaging is requested, it will happen on the raw data in the time domain, before Welch is run.
* `taper` and `tapsmofrq` : for taper selection and multi-tapering. Note that in case of multi-tapering, the data in the windows will be averaged across the tapers first
* `keeptrials` : whether trials should be left as-is, or you want a trial-average. If ``False``, and thus trial-averaging is requested, it will happen on the raw data in the time domain, before Welch is run.

Comparison with raw FFT
------------------------

Let's compare Welch's result with the raw FFT estimate::

fft_psd = spy.freqanalysis(wn)
fft_psd.singlepanelplot(trials=0, channel=0, logscale=False)

.. image:: ../_static/welch_raw_fft_power.png

The power spectral esimtate is much more noisy, meaning the variance per frequency bin is considerably larger compared to Welch's estimate.

.. note::
We don't need any parameters for ``freqanalysis`` here, as ``method='mtmfft'`` and ``tapsmofrq=None`` are the defaults.

Note that the absolute power values are lower, as we have a lot more frequency bins when calculating the raw FFT at once for the unsegmented signal::

fft_psd.freq.shape # (10001,)
welch_psd.freq.shape # (251,)

Syncopy normalizes spectral power per 1Hz bin, meaning that the noise power gets diluted over many more frequency bins when using the raw FFT. We can check that by summing over the frequency axis for both estimates::

np.sum(fft_psd.show(trials=0, channel=0)) # will be ~1
np.sum(welch_psd.show(trials=0, channel=0)) # will also be ~1


Investigating the Effects of the Overlap Parameter as a Function of Signal Length
Expand Down

0 comments on commit d5209dc

Please sign in to comment.