Power spectrum density estimation using Welch

Bugs
1

Hi Francois,

If you drag any raw file to the process window and select process->frequency->power spectrum density (Welch), you compute the power spectrum using a sliding window of a given length. Say:
Time window: 54 - 397sec
Window length: 1000ms

If I keep the window length constant, but run the same process with modified time window (say I use half the raw data), I get almost the same spectrum, but with very different scaling. (some times orders of magnitude).

Can you please check the scaling when averaging spectra? Perhaps you are not dividing by the number of window segments?

Also, a minor thing: if you first divide the raw file into small segments (say 10sec each), Brainstorm does not let you compute the power spectrum for each of them simultaneously. Even if raw segments are placed the process window, only the first segment will get a power spectrum.

Thank you,
Dimitrios

2

Hi Dimitrios,

If I keep the window length constant, but run the same process with modified time window (say I use half the raw data), I get almost the same spectrum, but with very different scaling.

You can check the code of bst_psd, it's a short and simple function.

Also, a minor thing: if you first divide the raw file into small segments (say 10sec each), Brainstorm does not let you compute the power spectrum for each of them simultaneously.

You need to set the same time vector for all the blocks first: process "Standardize > Uniform epoch time"
I know this is not that intuitive...

Francois

3

Hi Francois,

I believe I have spotted the mistake. On line 79 of bst_psd.m, you divide with the time length of the entire raw file (nTime, ~300,000msec in my case), but you should divide with the time length of the window (length(iTimes), 1000ms).

TFwin = 2 * Ffft(:,1:NFFT/2+1) ./ nTime;
should be
TFwin = 2 * Ffft(:,1:NFFT/2+1) ./ length(iTimes);

Best,
Dimitrios

4

Yes you’re right, I fixed the code.
The PSD window length is constant, hence the number of frequency bins is constant. The total length of recordings should not impact the range of values of the PSD results, just its accuracy.

Thank you for spotting this one!
Francois

5

Francois, I just updated my Brainstorm version and I am having issues with the PSD of my estimated sources. Every time I run the Welch algorithm, I get these warnings:

BST> Warning: Frequency band “delta” is empty.
BST> Warning: Frequency band “theta” is empty.
BST> Warning: Frequency band “alpha” is empty.

I tried both the linear and TF bands on sources and raw recordings. For the linear, the lowest frequency is 15.8Hz. It’s odd because I have previously computed a PSD on these files with Brainstorm.

Thanks,
YP

6

Your frequency resolution is extremely low, most likely because the length of the moving window used in the PSD is not defined properly.
The PSD is the average of the power of the FFT of small time window taken form the original file.
Increasing the length of this time window in the PSD options will lead the FFT to use more frequency bins, giving a better frequency resolution.

Does that solve your problem?
Francois

7

Yes, that worked. My mistake. I accidentally made the time window 10 ms instead of 1000ms. Fixed it now!

-YP