Frequency spectrum

UML diagrams presenting the flow of the analysis for each module are available here: https://github.com/ANCPLabOldenburg/MEG-QC-code/tree/main/diagrams

meg_qc.calculation.metrics.PSD_meg_qc.PSD_meg_qc(psd_params: dict, channels: dict, chs_by_lobe: dict, raw_orig: Raw, m_or_g_chosen: list, verbose_plots: bool, helperplots: bool)[source]

Main psd function. Calculates:

  • PSD for each channel

  • amplitudes (area under the curve) of functionally distinct frequency bands, such as

    delta (0.5-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), beta (12-30 Hz), and gamma (30-100 Hz) for each channel and average amplitude of band over all channels

  • average psd over all channels

  • noise frequencies for average psd + creates a band around them

  • noise frequencies for each channel + creates a band around them

  • noise amplitudes (area under the curve) for each noisy frequency band for average psd

  • noise amplitudes (area under the curve) for each noisy frequency band for each channel.

Frequency spectrum peaks we can often see:

  • Hz 50, 100, 150 - powerline EU

  • Hz 60, 120, 180 - powerline US

  • Hz 6 - noise of shielding chambers

  • Hz 44 - MEG device noise

  • Hz 17 - train station

  • Hz 10 - specific for MEG device in Nessy

  • Hz 1 - highpass filter.

  • flat frequency spectrum is white noise process. Has same energy in every frequency (starts around 50Hz or even below)

Parameters:

  • psd_params (dict) – dictionary with psd parameters originating from config file

  • channels (dict) – dictionary with channel names for each channel type: ‘mag’ or/and ‘grad’

  • chs_by_lobe (dict) – dictionary with channel objects sorted by ch type and lobe

  • raw_orig (mne.io.Raw) – raw data

  • m_or_g_chosen (list) – list with chosen channel types: ‘mag’ or/and ‘grad’

  • verbose_plots (bool) – True for showing plot in notebook.

  • helperplots (bool) – if True, plots with noisy freq bands for average PSD + for 3 different channels will be created (but not added to report).

Returns:

  • derivs_psd (list) – list with the psd derivatives as QC_derivative objects (figures)

  • simple_metric (dict) – dictionary with the psd simple metrics

  • psd_str (str) – string with notes about PSD for report

  • noisy_freqs_global (dict) – dictionary with noisy frequencies for average psd - used in Muscle artifact detection

meg_qc.calculation.metrics.PSD_meg_qc.Plot_psd_old(m_or_g: str, freqs: ndarray, psds: ndarray, channels: list, chs_by_lobe: dict, method: str, verbose_plots: bool)[source]

Plotting Power Spectral Density for all channels. OLD version

Parameters:

  • m_or_g (str) – ‘mag’ or ‘grad’

  • freqs (np.ndarray) – frequencies

  • psds (np.ndarray) – power spectral density for each channel

  • channels (list) – list of channel names

  • chs_by_lobe (dict) – dictionary with channel objects sorted by lobe

  • method (str) – ‘welch’ or ‘multitaper’ or other method

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

QC_derivative object with plotly figure as content

Return type:

QC_derivative

meg_qc.calculation.metrics.PSD_meg_qc.add_log_buttons(fig: Figure)[source]

Add buttons to switch scale between log and linear. For some reason only swithcing the Y scale works so far.

Parameters:

fig (go.Figure) – The figure to be modified withot buttons

Returns:

fig – The modified figure with the buttons

Return type:

go.Figure

meg_qc.calculation.metrics.PSD_meg_qc.assign_psds_to_channels(chs_by_lobe, freqs, psds)[source]

TODO: fix docstrings TODO: move as a part of MEG_channels in plots

Assign std or ptp values of each epoch as list to each channel. This is done for easier plotting when need to plot epochs per channel and also color coded by lobes.

Parameters:

  • what_data (str) – ‘peaks’ for peak-to-peak amplitudes or ‘stds’

  • chs_by_lobe (dict) – dictionary with channel objects sorted by lobe.

  • df_std_ptp (pd.DataFrame) – Data Frame containing std or ptp value for each chnnel and each epoch

Returns:

chs_by_lobe – updated dictionary with channel objects sorted by lobe - with info about std or ptp of epochs.

Return type:

dict

meg_qc.calculation.metrics.PSD_meg_qc.cut_the_noise_from_psd(noisy_bands_indexes: List[dict], freqs: list, one_psd: list, helper_plots: bool, ch_name: str = '', noisy_freqs_indexes: list = [], unit: str = '', verbose_plots: bool = True)[source]

Cut the noise peaks out of PSD curve. By default, it is not used, but can be turned on. If turned on, in the next steps, the area under the curve will be calculated only for the cut out peaks.

By default, the area under the curve is calculated under the whole peak, uncluding the ‘main brain signal’ psd area + peak area. This is done, because in reality we can not define, which part of the ‘noisy’ frequency is signal and which is noise. In case later, during preprocessing this noise will be filtered out, it will be done completely: both the peak and the main psd area.

Process:

  1. Find the height of the noise peaks. For this take the average between the height of the start and end of this noise bend.

  2. Cut the noise peaks out of PSD curve at the found height.

  3. Baseline the peaks: all the peaks are brought to 0 level.

Function also can prodece helper plot to demonstrate the process.

Parameters:

  • noisy_bands_indexes (list[list]) – list of lists with indexes of noisy bands. Indexes! Not frequency bands themselves. Index is defined by fequency/freq_resolution.

  • freqs (list) – vector of frequencies

  • one_psd (list) – vector if psd values for 1 channel (or 1 average over all channels)

  • helper_plots (bool) – if True, helper plots will be produced

  • ch_name (str, optional) – channel name, by default ‘’, used for plot display

  • noisy_freqs_indexes (list, optional) – list of indexes of noisy frequencies. Indexes! not frequencies themselves. Index is defined by fequency/freq_resolution., by default [] because we might have no noisy frequencies at all. Used for plot display.

  • unit (str, optional) – unit of the psd values, by default ‘’, used for plot display

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

psd_only_peaks_baselined – vector of psd values for 1 channel (or 1 average over all channels) with the noise peaks cut out and baselined to 0 level. Later used to calculate area under the curve for the noise peaks only.

Return type:

list

meg_qc.calculation.metrics.PSD_meg_qc.find_noisy_freq_bands_complex(ch_name: str, freqs: list, one_psd: list, helper_plots: bool, m_or_g: str, prominence_lvl_pos: int, verbose_plots: bool)[source]

Detect the frequency band around the noise peaks. Complex approach: This function is trying to detect the actual start and end of peaks.

  1. Bands around the noise frequencies are created based on detected peak_width.

  2. If the found bands overlap, they are cut at the lowest point between 2 neighbouring noise peaks pn PSD curve.

This function is not used by default, becausesuch a complex approach, even though can accurately find start and end of the noise bands, is not very robust. It can sometimes take too much of the area arouund the noise peak, leading to a large part of the signel folsely counted as noise. By default, the more simple approach is used. See find_noisy_freq_bands_simple() function.

Parameters:

  • ch_name (str) – channel name like ‘MEG1234’ or just ‘Average’. For plotting purposes only.

  • freqs (list) – list of frequencies

  • one_psd (list) – list of psd values for one channels or for the average over multiple channels

  • helper_plots (bool) – if True, helper plots will be shown

  • m_or_g (str) – ‘mag’ or ‘grad’ - for plotting purposes only - to get the unit of the psd values

  • prominence_lvl_pos (int) – prominence level for peak detection. The higher the value, the more peaks will be detected.

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

  • noisy_freqs (list) – list of noisy frequencies

  • noisy_freqs_indexes (list) – list of indexes of noisy frequencies in the psd

  • noisy_bands_final (list[list]) – list of lists of noisy frequency bands. Each list contains 2 values: start and end of the band.

  • noisy_bands_final_indexes (list[list]) – list of lists of indexes of noisy frequency bands. Each list contains 2 values: start and end of the band.

  • split_indexes (list) – list of indexes of the split points in the psd

meg_qc.calculation.metrics.PSD_meg_qc.find_noisy_freq_bands_simple(ch_name: str, freqs: list, one_psd: list, helper_plots: bool, m_or_g: str, prominence_lvl_pos: int, band_half_length: float, verbose_plots: bool)[source]

Form a frequency band around the noise peaks. Simple approach: used by default.

  1. Create frequency band around central noise frequency just by adding -x…+x Hz around.

  2. If the found bands overlap, they are cut at the lowest point between 2 neighbouring noise peaks pn PSD curve.

Parameters:

  • ch_name (str) – channel name like ‘MEG1234’ or just ‘Average’. For plotting purposes only.

  • freqs (list) – list of frequencies

  • one_psd (list) – list of psd values for one channels or for the average over multiple channels

  • helper_plots (bool) – if True, helper plots will be shown

  • m_or_g (str) – ‘mag’ or ‘grad’ - for plotting purposes only - to get the unit of the psd values

  • prominence_lvl_pos (int) – prominence level for peak detection. The higher the value, the more peaks will be detected.

  • band_half_length (float) – length of the frequency band before and after the noise peak in Hz. The band will be created by adding -band_half_length…+band_half_length Hz around the noise peak.

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

  • noisy_freqs (list) – list of noisy frequencies

  • noisy_freqs_indexes (list) – list of indexes of noisy frequencies in the psd

  • noisy_bands_final (list[list]) – list of lists of noisy frequency bands. Each list contains 2 values: start and end of the band.

  • noisy_bands_final_indexes (list[list]) – list of lists of indexes of noisy frequency bands. Each list contains 2 values: start and end of the band.

  • split_indexes (list) – list of indexes of the split points in the psd

meg_qc.calculation.metrics.PSD_meg_qc.find_number_and_ampl_of_noise_freqs(ch_name: str, freqs: list, one_psd: list, pie_plotflag: bool, helper_plots: bool, m_or_g: str, cut_noise_from_psd: bool, prominence_lvl_pos: int, simple_or_complex: str = 'simple', verbose_plots: bool = True)[source]

The function finds the number and amplitude of noisy frequencies in PSD function in these steps:

  1. Calculate average psd curve over all channels

  2. Run peak detection on it -> get number of noise freqs. Create the bands around them. Split blended freqs.

  3. (Optional) Fit a curve to the general psd OR cut the noise peaks at the point they start and baseline them to 0. Optional. By default not used

  4. Calculate area under the curve for each noisy peak (amplitude of the noise)):

    • If 3 was done: area is limited to where noise band crosses the fitted curve. - count from there.

    • If not (default): area is limited to the whole area under the noise band, including the psd of the signal.

  5. Calculate what part of the whole psd is the noise (noise amplitude) and what part is the signal (signal amplitude) + plot as pie chart

Parameters:

  • ch_name (str) – name of the channel or ‘average’

  • freqs (list) – list of frequencies

  • one_psd (list) – list of psd values for one channel or average psd

  • pie_plotflag (bool) – if True, plot the pie chart

  • helper_plots (bool) – if True, plot the helper plots (will show the noise bands, how they are split and how the peaks are cut from the psd if this is activated).

  • m_or_g (str) – ‘mag’ or ‘grad’

  • cut_noise_from_psd (bool) – if True, cut the noise peaks at the point they start and baseline them to 0. Optional. By default not used

  • prominence_lvl_pos (int) – prominence level for peak detection (central frequencies of noise bands). The higher the value, the more peaks will be detected. prominence_lvl will be different for average psd and psd of 1 channel, because average has small peaks smoothed.

  • simple_or_complex (str) – ‘simple’ or ‘complex’ approach to create the bands around the noise peaks. Simple by default. See functions above for details.

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

  • noise_pie_derivative (list) – list with QC_derivative object containing the pie chart with the noise amplitude and signal amplitude

  • noise_ampl (list) – list of noise amplitudes for each noisy frequency band

  • noise_ampl_relative_to_signal (list) – list of noise amplitudes relative to the signal amplitude for each noisy frequency band

  • noisy_freqs (list) – list of noisy frequencies

meg_qc.calculation.metrics.PSD_meg_qc.get_ampl_of_brain_waves(channels: list, m_or_g: str, freqs: ndarray, psds: ndarray, avg_psd: ndarray, plotflag: bool, verbose_plots: bool)[source]

Amplitude of frequencies calculation for all channels. If desired: creating a pie chart of mean amplitude of every band over the entire data.

Parameters:

  • channels (list) – List of channel names

  • m_or_g (str) – ‘mag’ or ‘grad’ - to choose which channels to calculate amplitude for.

  • freqs (np.ndarray) – numpy array of frequencies for mag or grad

  • psds (np.ndarray) – numpy array of power spectrum dencities for mag or grad

  • avg_psd (np.ndarray) – numpy array of average power spectrum dencities for mag or grad

  • plotflag (bool) – need to plot pie chart or not

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

  • psd_pie_derivative (QC_derivative object or empty list) – If plotflag is True, returns one QC_derivative object, which is a plotly piechart figure. If plotflag is False, returns empty list.

  • dfs_with_name (list) – List of dataframes with amplitude of each frequency band in each channel (dfs: absolute amplitude, relative amplitude, amplitude divided by number of frequencies in the band)

  • mean_brain_waves_dict (dict) – Dictionary of mean amplitude of each frequency band (used for simple metric json)

meg_qc.calculation.metrics.PSD_meg_qc.get_ampl_of_noisy_freqs(channels, freqs, avg_psd, psds, m_or_g, pie_plotflag=True, helperplots=True, cut_noise_from_psd=False, prominence_lvl_pos_avg=50, prominence_lvl_pos_channels=15, simple_or_complex='simple', verbose_plots: bool = True)[source]

Find noisy frequencies, their absolute and relative amplitude for averages over all channel (mag or grad) PSD and for each separate channel.

Parameters:

  • channels (list) – list of channel names

  • freqs (list) – list of frequencies

  • avg_psd (list) – list of average PSD values over all channels

  • psds (list) – list of PSD values for each channel

  • m_or_g (str) – ‘mag’ or ‘grad’

  • pie_plotflag (bool) – if True, plot pie chart of SNR

  • helperplots (bool) – if True, plot helper plots

  • cut_noise_from_psd (bool) – if True, cut the noise peaks at the point they start and baseline them to 0.

  • prominence_lvl_pos_avg (int) – prominence level of peak detection for finding noisy frequencies in the average PSD

  • prominence_lvl_pos_channels (int) – prominence level of peak detection for finding noisy frequencies in the PSD of each channel

  • simple_or_complex (str) – ‘simple’ or ‘complex’ - method of finding noisy frequencies. see find_number_and_ampl_of_noise_freqs() for details

  • verbose_plots (bool) – True for showing plot in notebook.

Returns:

  • noise_pie_derivative (QC_derivative object or empty list if pie_plotflag is False) – QC_derivative containig a pie chart of SNR

  • noise_ampl_global (dict) – dictionary for simple metric with info about noisy frequencies in the average PSD, absolute values for bands

  • noise_ampl_relative_to_all_signal_global (dict) – dictionary for simple metric with info about noisy frequencies in the average PSD, relative values for bands

  • noisy_freqs_global (dict) – dictionary for simple metric with info about noisy frequencies in the average PSD, central frequencies

  • noise_ampl_local_all_ch (dict) – dictionary for simple metric with info about noisy frequencies in the PSD of each channel, absolute values for bands

  • noise_ampl_relative_to_all_signal_local_all_ch (dict) – dictionary for simple metric with info about noisy frequencies in the PSD of each channel, relative values for bands

  • noisy_freqs_local_all_ch (dict) – dictionary for simple metric with info about noisy frequencies in the PSD of each channel, central frequencies

meg_qc.calculation.metrics.PSD_meg_qc.get_mean_bands_amplitude(freq_bands: list, freqs: list, psds: ndarray, channels: list, bands_names: list | None = None)[source]

Calculate the area under the curve of frequency bands (e.g. alpha, beta, gamma, delta, …) for mag or grad. Adopted from: https://raphaelvallat.com/bandpower.html Take mean over all channels (mags/grads).

Parameters:

  • freq_bands (list) – List of lists of frequencies. Expects list of lists: [[f_low, f_high], [f_low, f_high], …]

  • freqs (list) – List of frequencies.

  • psds (np.ndarray or list) – numpy array (or list) of power spectrum dencities. Expects array of arrays: channels*psds. (or list of lists) Will not work properly if 1 dimentional array given. In this case do: np.array([your_1d_array])

  • channels (list) – List of channel names. Expects list of strings: [‘MEG 0111’, ‘MEG 0112’, …] If only one channel is given, it should be a list of one string: [‘Average]

  • bands_names (list, optional) – List of names of the bands. If None, the names will be generated automatically, by default None Important! taht these names are in the same order as in freq_bands.

Returns:

  • band_ampl_df (pd.DataFrame) – Dataframe of amplitudes of each frequency band like: [abs_power_of_delta, abs_power_of_gamma, etc…] - in absolute values

  • band_ampl_relative_to_signal_df (pd.DataFrame) – Dataframe of amplitudes of each frequency band divided by the total amplitudeof the signal for this channel. Shows how much amplitude this particular band takes in the entire signal.

  • ampl_by_Nfreq_per_ch_list_df (pd.DataFrame) – Dataframe of amplitudes of each frequency band divided by the number of frequencies in the band. (This is done to compare with RMSE later. But not used any more).

  • total_signal_amplitude (list) – List of total signal amplitude for each channel.

meg_qc.calculation.metrics.PSD_meg_qc.get_nfft_nperseg(raw: Raw, psd_step_size: float)[source]

Get nfft and nperseg parameters for Welch psd function. Allowes to always have the step size in psd which is chosen by the user. Recommended 0.5 Hz.

Parameters:

  • raw (mne.io.Raw) – raw data

  • psd_step_size (float) – step size for PSD chosen by user, recommended 0.5 Hz

Returns:

  • nfft (int) – Number of points for fft. Used in welch psd function from mne. The length of FFT used, must be >= n_per_seg (default: 256). The segments will be zero-padded if n_fft > n_per_seg. If n_per_seg is None, n_fft must be <= number of time points in the data.

  • nperseg (int) – Number of points for each segment. Used in welch psd function from mne. Length of each Welch segment (windowed with a Hamming window). Defaults to None, which sets n_per_seg equal to n_fft.

meg_qc.calculation.metrics.PSD_meg_qc.make_dict_global_psd(mean_brain_waves_dict: dict, noisy_freqs_global: list, noise_ampl_global: list, noise_ampl_relative_to_all_signal_global: list)[source]

Create a dictionary for the global part of psd simple metrics. Global: overall part of noise in the signal (all channels averaged).

Parameters:

  • mean_brain_waves_dict (dict) – dictionary with the mean brain waves (alpha, beta, etc) metrics in the form: {wave band name: {mean realtive value, mean absolute value}, …}

  • noisy_freqs_global (list) – list of noisy frequencies

  • noise_ampl_global (list) – list of noise amplitudes for each noisy frequency band

  • noise_ampl_relative_to_all_signal_global (list) – list of noise amplitudes relative to the total signal amplitude for each noisy frequency band

Returns:

dict_global – dictionary with the global part of psd simple metrics

Return type:

dict

meg_qc.calculation.metrics.PSD_meg_qc.make_dict_local_psd(noisy_freqs_local: dict, noise_ampl_local: dict, noise_ampl_relative_to_all_signal_local: dict, channels: list)[source]

Create a dictionary for the local part of psd simple metrics. Local: part of noise in the signal for each channel separately.

Parameters:

  • noisy_freqs_local (dict) – dictionary with noisy frequencies for each channel

  • noise_ampl_local (dict) – dictionary with noise amplitudes for each noisy frequency band for each channel

  • noise_ampl_relative_to_all_signal_local (dict) – dictionary with noise amplitudes relative to the total signal amplitude for each noisy frequency band for each channel

Returns:

dict_local – dictionary with the local part of psd simple metrics

Return type:

dict

meg_qc.calculation.metrics.PSD_meg_qc.make_simple_metric_psd(mean_brain_waves_dict: dict, noise_ampl_global: dict, noise_ampl_relative_to_all_signal_global: dict, noisy_freqs_global: dict, noise_ampl_local: dict, noise_ampl_relative_to_all_signal_local: dict, noisy_freqs_local: dict, m_or_g_chosen: list, freqs: dict, channels: dict)[source]

Create a dictionary for the psd simple metrics.

Parameters:

  • mean_brain_waves_dict (dict) – dictionary with mean brain waves (alpha, beta, etc) for each channel type. Inside each channel type: dictionary in the form: {wave band name: {mean realtive value, mean absolute value}, …}

  • noise_ampl_global (dict) – dictionary with noise amplitudes for each noisy frequency band

  • noise_ampl_relative_to_all_signal_global (dict) – dictionary with noise amplitudes relative to the total signal amplitude for each noisy frequency band

  • noisy_freqs_global (dict) – dictionary with noisy frequencies

  • noise_ampl_local (dict) – dictionary with noise amplitudes for each noisy frequency band for each channel

  • noise_ampl_relative_to_all_signal_local (dict) – dictionary with noise amplitudes relative to the total signal amplitude for each noisy frequency band for each channel

  • noisy_freqs_local (dict) – dictionary with noisy frequencies for each channel

  • m_or_g_chosen (list) – list with chosen channel types: ‘mag’ or/and ‘grad’

Returns:

simple_metric – dictionary with the psd simple metrics

Return type:

dict

meg_qc.calculation.metrics.PSD_meg_qc.plot_one_psd(ch_name: str, freqs: List, one_psd: List, peak_indexes: List, noisy_freq_bands_indexes: List[list], unit: str)[source]

Plot PSD for one channels or for the average over multiple channels with noise peaks and split points using plotly.

Parameters:

  • ch_name (str) – channel name like ‘MEG1234’ or just ‘Average’

  • freqs (list) – list of frequencies

  • one_psd (list) – list of psd values for one channels or for the average over multiple channels

  • peak_indexes (list) – list of indexes of the noise peaks in the psd

  • noisy_freq_bands_indexes (list[list]) – list of lists of indexes of the noisy frequency bands in the psd. Indexes! Not frequency bands themselves. Index is defined by fequency/freq_resolution.

  • unit (str) – unit of the psd values. For example ‘T/Hz’

Returns:

plotly figure of the psd with noise peaks and bands around them.

Return type:

fig

meg_qc.calculation.metrics.PSD_meg_qc.split_blended_freqs_at_the_lowest_point(noisy_bands_indexes: List[list], one_psd: List[dict], noisy_freqs_indexes: List[dict])[source]

If there are 2 bands that are blended together, split them at the lowest point between 2 central noise frequencies.

Parameters:

  • noisy_bands_indexes (list[list]) – list of lists with indexes of noisy bands. Indexes! not frequency bands themselves. Index is defined by fequency/freq_resolution.

  • one_psd (list) – vector if psd values for 1 channel (or 1 average over all channels)

  • noisy_freqs_indexes (list) – list of indexes of noisy frequencies. Indexes! not frequencies themselves. Index is defined by fequency/freq_resolution.

Returns:

  • noisy_bands_final_indexes (list[list]) – list of lists with indexes of noisy bands After the split. Indexes! not frequency bands themselves. Index is defined by fequency/freq_resolution.

  • split_indexes (list) – list of indexes at which the bands were split (used later for plotting only).