Skip to content

Instantly share code, notes, and snippets.

@rcassani
Created March 27, 2025 14:26
Show Gist options
  • Select an option

  • Save rcassani/1262dde3bf677228a7ed8ecd3febadd0 to your computer and use it in GitHub Desktop.

Select an option

Save rcassani/1262dde3bf677228a7ed8ecd3febadd0 to your computer and use it in GitHub Desktop.
Understanding FOOOF parameters in Brainstorm
%% Understanding FOOOF parameters in Brainstorm, using one channel as example
% This is based on the files from the FOOOF tutorial
% https://neuroimage.usc.edu/brainstorm/Tutorials/Fooof
%
% Raymundo Cassani, 2025
%% FIND PSD AND FOOOF FILES
% Process: Select timefreq files in: Subject01/@rawS01_AEF_20131218_01_600Hz_notch
sFiles = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
'subjectname', 'Subject01', ...
'condition', '@rawS01_AEF_20131218_01_600Hz_notch', ...
'tag', '', ...
'includebad', 0, ...
'includeintra', 0, ...
'includecommon', 0, ...
'outprocesstab', 'no'); % No
iPsd = find(~cellfun(@isempty, strfind({sFiles.Comment}, 'PSD')));
iFooof = find(~cellfun(@isempty, strfind({sFiles.Comment}, 'specparam')));
% PSD and FOOOF files
psdFileName = sFiles(iPsd).FileName;
fooofFileName = sFiles(iFooof).FileName;
%% LOAD FILES
sPsd = in_bst_timefreq(psdFileName);
sFooof = in_bst_timefreq(fooofFileName);
% Select one channel
ix = 1 ; % MLC11
%% GET PSD AND FOOOF DATA
% PSD
freqs_psd = sPsd.Freqs;
y_psd_u2_hz = squeeze(sPsd.TF(ix,:,:));
y_psd_db = 10*log10(y_psd_u2_hz);
% FOOOF parameters
fooof_params = sFooof.Options.FOOOF;
freqs_foo = fooof_params.freqs;
[~, iFreqs_psd] = ismember(freqs_foo, freqs_psd);
fooof_data = fooof_params.data(ix);
% Aperiodic
y_ap = 10*fooof_data.aperiodic_params(1) + 10*log10(freqs_foo.^-fooof_data.aperiodic_params(2));
% Peaks
y_peak1 = 10*fooof_data.peak_params(1,2) .* exp((-((freqs_foo - fooof_data.peak_params(1,1)) ./ fooof_data.peak_params(1,3)) .^2)./2);
y_peak2 = 10*fooof_data.peak_params(2,2) .* exp((-((freqs_foo - fooof_data.peak_params(2,1)) ./ fooof_data.peak_params(2,3)) .^2)./2);
%% PLOTS
figure()
% PSD
plot(freqs_foo, y_psd_db(iFreqs_psd));
hold on
% Aperiodic
plot(freqs_foo, y_ap);
% Aperiodic+Peaks
plot(freqs_foo, y_ap + y_peak1 + y_peak2);
% Peaks
ylims = ylim();
plot(freqs_foo, y_peak1 + y_peak2 + ylims(1));
xlabel('frequency (Hz)')
ylabel('power (dB)')
legend({'Original PSD', 'Aperiodic', 'FOOOF model', 'FOOOF peaks'})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment