using 'psd' or 'fft' in std_spec creates 60dB difference

[neuromechanist]

Hi, I use std_precomp to create the power spectra of the ICA components in a source cluster. I tried computing the spectra using psd or fft and results have almost 60 dB difference (I am aware that it is not exactly dB, since it does not have a reference). From experience and also the time-frequency data, I know that the results of the fft mode is correct. Also, the fft results appeared to be very noisy. Do you have any comments for adjusting parameters to get smoother results? A friend of mine at another university also observed this discrepancy using a totally different dataset but a similar (yet independent) pipeline.

the full script for the computing spectra is as follows:

% PSD
[f_STUDY, f_ALLEEG] = std_precomp(f_STUDY, f_ALLEEG, 'components', 'recompute' ,'on','spec','on',...
                'specparams', {'specmode', 'psd', 'logtrials', 'off', 'freqrange',[3 80], 'freqfac', 10}, ...
                'recompute', 'on', 'savetrials', 'on');

% FFT
[f_STUDY, f_ALLEEG] = std_precomp(f_STUDY, f_ALLEEG, 'components', 'recompute' ,'on','spec','on',...
                'specparams', {'specmode', 'fft', 'logtrials', 'off', 'freqrange',[3 80]}, ...
                'recompute', 'on', 'savetrials', 'on');

BTW, changing logtrials and/or freqfac values do not seem to make any change.

EEGLAB version: 2019 MATLAB version: 2018b

[arnodelorme]

One is power spectral density (psd) and one is power spectrum (fft). The FFT scale is wrong as it should not be per Hz since it is absolute power.

[neuromechanist]

Thank you for your reply. I believe in many contexts, power spectral density PSD, and power spectrum PS are used interchangeably, please see the Sciencedirect description for power spectrum:

Also, assuming that there is indeed a difference between PS and PSD, how we can draw a relation between those? A solution might be to multiply PSD by the Equivalent Noise Bandwidth (ENBW) of the Hamming window. ENBW is in the order of ~0.1-1 for epoch sizes of ~1 seconds at 512 Hz, which actually makes PS values smaller than PSD. Still, I don't think the 60 dB difference can be justified.

[arnodelorme]

In my understanding, power spectral density is expressed per Hz. Power spectrum is simply the square of the amplitude of the FFT - often log transform for scaling, so it is not the spectral density per Hz.

You are right that FFT should be normalized by window length. This is the formula we have in newtimef. We could try to use it in std_spec.m

alltfX = 2/0.375*alltfX/g.winsize; % TF and MC (12/11/2006): normalization, divide by g.winsize
P  = alltfX.*conj(alltfX); % power    
% TF and MC (12/14/2006): multiply by 2 account for negative frequencies,
% and ounteract the reduction by a factor 0.375 that occurs as a result of 
% cosine (Hann) tapering. Refer to Bug 446
% Modified again 04/29/2011 due to comment in bug 1032

See also

https://sccn.ucsd.edu/bugzilla/show_bug.cgi?id=1032 https://sccn.ucsd.edu/bugzilla/show_bug.cgi?id=446

Copying @widmann who is more expert than me on these issues.

[neuromechanist]

To provide an update for this issue:

1. As Arno mentioned, fft option for std_precomp which calls std_spec, only does the Fast Fourier Transform (FFT), Not computing power spectral density (PSD) using FFT (see Matlab example for PSD using FFT)
2. ICA weights, by default, are scaled to average RMS microvolt when eeg_checkset is called during loading data. This is particularly useful for a better time-frequency (TF) analysis. Since TF analysis is w.r.t. a baseline, scaling would not hurt the results. But computing PSD seems to be different. The difference between PSD using the original ICA weights and the PSD using RMS-microvolt-scaled ICA is very telling (note the y scale in the figure below). Again, RMS scaling happens only during loading the dataset, or whenever eeg_checkset() is called. So, if you save your datasets to process it later, your ERSPs should be fine, but your PSDs may not be Ok. (Look below for a possible workaround)

3. Original weights are conserved in EEG.etc. Substituting EEG.icaweights with EEG.etc.icaweights_beforerms and EEG.etc.icasphere matrix with EEG.etc.icasphere_beforerms, and of course, clearing EEG.icawinv to have EEGLAB create it again before running std_spec seem to solve this issue. Just note that you might still need to use the RMS-corrected ICA weights when running std_precomp for computing TF. So, I would use std_precomp once with RMS-corrected ICA for TF and another time with the original ICA weights for PSD.

[arnodelorme]

See also #175

[arnodelorme]

PP(:,j) = 2/0.375*abs(tmpX).^2; % power % TF and MC (12/14/2006): multiply by 2 account for negative frequencies, % Counteract the reduction by a factor 0.375 % that occurs as a result of cosine (Hann) tapering. Refer to Bug 446

@cll008 to add density correction line 390 of std_spec.m and also apply correction above, then compare with the other method.