Consistency in power spectra computed in STUDY by std_spec?

[amisepa]

Description

Unclear what scales/units are returned by some STUDY power spectrum modes. Is it either absolute power or dB? Units in ylabel don't seem to correspond (see images below). It always indicates dB even though values and settings are different.
And outputs should potentially be checked to make sure the correct values are returned. There are a lot of 'if' conditions, conversions, and methods employed, so it's hard to keep track (e.g., psd computed in dB by spectcomp in spectopo with either pwelch or spec method depending on presence/absence of signal processing toolbox, then converted to dB in spectopo, and converted again to abs power in std_spec).
Also, calculating psd is extremly long compared to when I compute it with pwelch directly in Matlab. Not sure why.
Also, for FFT, there doesn't seem to be a difference if you chose the option 'output' 'power' or not.

• FFT with 'logtrials', 'off' (default):
  

• FFT with 'logtrials', 'on':
  

• FFT with 'logtrials', 'off', 'output', 'power':
  

• FFT with 'logtrials', 'on', 'output', 'power':
  

• PSD with 'logtrials, 'off':
  

• PSD with 'logtrials, 'on':
  

---

Using the 'test' feature in the precompute window, different plots are generated when using the log 'on' option (compared to above after precomputation).

• FFT with 'logtrials', 'on':
  

• PSD with 'logtrials, 'on':
  

#### Steps to Reproduce

1. Load a study
2. precompute channel power spectra:

 [STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','fft','logtrials','off'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);


[STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','fft','logtrials','on'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);

[STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','fft','logtrials','off', 'output', 'power'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);

[STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','fft','logtrials','on', 'output', 'power'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);

[STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','psd','logtrials','off'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);

[STUDY, ALLEEG] = std_precomp(STUDY, ALLEEG, {},'savetrials','on','recompute','on','spec','on','specparams',...
    {'specmode','psd','logtrials','on'});
STUDY = pop_specparams(STUDY, 'averagechan','on', 'freqrange',[1 50]);
STUDY = std_specplot(STUDY,ALLEEG,'channels',{'Fp1','AF7','AF3','F1','F3','F5','F7','FT7','FC5','FC3','FC1','C1','C3','C5','T7','TP7','CP5','CP3','CP1','P1','P3','P5','P7','P9','PO7','PO3','O1','Iz','Oz','POz','Pz','CPz','Fpz','Fp2','AF8','AF4','AFz','Fz','F2','F4','F6','F8','FT8','FC6','FC4','FC2','FCz','Cz','C2','C4','C6','T8','TP8','CP6','CP4','CP2','P2','P4','P6','P8','P10','PO8','PO4','O2'}, 'design', 1);

3. see code in std_spec and spectopo

[arnodelorme]

@cll008 what would you recommend

[arnodelorme]

@amisepa what change would your recommend?

[arnodelorme]

@MakotoMiyakoshi what would you recommend

[MakotoMiyakoshi]

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

no-option results must be >= logtrial results because of triangle inequality. If it does not change the result, it must be a bug.

[amisepa]

Based on this present issue, and past ones (#175 and #172), in a world where someone had time to do it, it would be great to review all the functions that underly power spectra computation and plotting, correct them, and make a simpler version. Currently, if I understand correctly, PSD is first computed in dB by spectcomp in spectopo, using either pwelch or spec method depending on presence or absence of signal processing toolbox. Then it is converted to dB in spectopo, and converted again to abs power in std_spec.

Looks like a quick fix might start with the discrepancies due to eeg_checkset which scales power to average RMS microVolt (see issue 172).

[MakotoMiyakoshi]

Is it true that first 'computed in dB' 'Then it is converted to dB' i.e. dB-conversion of dB is used? Or am I reading it wrong?

[amisepa]

Sorry, that is what I said but it's a mistake.

In std_spec:
if 'logtrials' is 'on' (line 303:)
line 311: compute PSD in dB by calling spectopo, which call spectcomp, which uses either pwelch or spec functions (here not sure what the units are. Watts per Hz?)
line 321: if 'logtrials' is 'off', converts to abs, else keeps in dB
else
line 328 and 329: compute PSD in dB and conversion to abs
end

So it looks correct. No dB-conversion of dB. But could use some simplification, and corrections in the ylabels of the plots (see above).

I think it's the "test" feature in STUDY that has some weird things (see above).

Note there is also some normalization line 982 of spectopo that might explain some differences in scales as well (not sure why this is done)?

[MakotoMiyakoshi]

Dear Cedric,

line 311: compute PSD in dB by calling *spectopo*, which call *spectcomp*,

which uses either *pwelch* or *spec* functions (here not sure what the units are. Watts per Hz?) Just reading the code, but spectopo() line 882 shows 'eegspecdB' as the first output from spectcomp() function, but the variable eegspec is a sum of pwelch() results. So eegspec seems to have a unit of microV^2/Hz, not dB-converted yet. Later in line 372 for example, I found eegspecdB = 10*log10(eegspecdB); This is where the output is dB-convereted. It's just a confusing naming of the variables and I understand that it is misleading indeed. Makoto

…

On Thu, Sep 2, 2021 at 4:24 PM Cedric Cannard ***@***.***> wrote: Sorry, that is what I said but it's a mistake. In std_spec: if 'logtrials' is 'on' (line 303:) line 311: compute PSD in dB by calling *spectopo*, which call *spectcomp*, which uses either *pwelch* or *spec* functions (here not sure what the units are. Watts per Hz?) line 321: if 'logtrials' is 'off', converts to abs, else keeps in dB else line 328 and 329: compute PSD in dB and conversion to abs end So it looks correct. No dB-conversion of dB. But could use some simplification, and corrections in the ylabels of the plots (see above). I think it's the "test" feature in STUDY that has some weird things (see above). Note there is also some normalization line 982 of spectopo that might explain some differences in scales as well (not sure why this is done)? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_sccn_eeglab_issues_364-23issuecomment-2D912127128&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=XpgQqocIrCV6TVqnYh-LdkOfTZUwl7XPmRMdz3vpxyQ&s=KRieI4_cZVtl93J9lHL7X8_YgJAR5vt-xNk8nAeCYQE&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AK7GVNMVWJBIVG2JF2FGMFDUAABRDANCNFSM5BXOEB2A&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=XpgQqocIrCV6TVqnYh-LdkOfTZUwl7XPmRMdz3vpxyQ&s=pMZiFslUDpvhOGskVMzkCzdW0WmINabLnP42Bqjtc4Y&e=> . Triage notifications on the go with GitHub Mobile for iOS <https://urldefense.proofpoint.com/v2/url?u=https-3A__apps.apple.com_app_apple-2Dstore_id1477376905-3Fct-3Dnotification-2Demail-26mt-3D8-26pt-3D524675&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=XpgQqocIrCV6TVqnYh-LdkOfTZUwl7XPmRMdz3vpxyQ&s=Agf1LE_XyjGld-ShEGqSCQXV5jRUq41WD-OFn0sIC5E&e=> or Android <https://urldefense.proofpoint.com/v2/url?u=https-3A__play.google.com_store_apps_details-3Fid-3Dcom.github.android-26referrer-3Dutm-5Fcampaign-253Dnotification-2Demail-2526utm-5Fmedium-253Demail-2526utm-5Fsource-253Dgithub&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=XpgQqocIrCV6TVqnYh-LdkOfTZUwl7XPmRMdz3vpxyQ&s=Q-ChMNsiITOHtNOhk6087451D7vRPDOgXvUPGDRbxhs&e=>.

[amisepa]

@MakotoMiyakoshi

I see, indeed, thanks!

[arnodelorme]

Can we close this report?

[MakotoMiyakoshi]

Sure.

…

On Mon, Mar 28, 2022 at 11:24 AM Arnaud Delorme ***@***.***> wrote: Can we close this report? — Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_sccn_eeglab_issues_364-23issuecomment-2D1080992876&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=AH-TraQo_7aPV5X8Y6GWiIouJRn_3xa0dQxjpI0NwnIbqUFniANWmhuBBA0UOvc_&s=8azlEqaSAJRli-ZEbmxepkXml52eDHDg_ufjiyHNFUw&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AK7GVNOZXLYH657VQ46UCF3VCH2LJANCNFSM5BXOEB2A&d=DwMCaQ&c=-35OiAkTchMrZOngvJPOeA&r=pyiMpJA6aQ3IKcfd-jIW1kWlr8b1b2ssGmoavJHHJ7Q&m=AH-TraQo_7aPV5X8Y6GWiIouJRn_3xa0dQxjpI0NwnIbqUFniANWmhuBBA0UOvc_&s=P2dvUIc2Rz5i7CkakI47ZM_f0dTnn3Ph7W8_-mmRyRY&e=> . You are receiving this because you were mentioned.Message ID: ***@***.***>