diff --git a/README.md b/README.md index bef19fb..87b2dda 100644 --- a/README.md +++ b/README.md @@ -108,7 +108,7 @@ Contributions to PyRASA are welcome! Whether it's raising issues, improving docu If you are using IRASA please cite the smart people who came up with the algorithm: -Wen, H., & Liu, Z. (2016). Separating fractal and oscillatory components in the power spectrum of neurophysiological signal. Brain topography, 29, 13-26. +Wen, H., & Liu, Z. (2016). Separating fractal and oscillatory components in the power spectrum of neurophysiological signal. Brain topography, 29, 13-26. https://doi.org/10.1007/s10548-015-0448-0 If you are using PyRASA it would be nice, if you could additionally cite us (whenever the paper is finally ready): diff --git a/examples/basic_functionality.ipynb b/examples/basic_functionality.ipynb index 8956d5d..9ee1028 100644 --- a/examples/basic_functionality.ipynb +++ b/examples/basic_functionality.ipynb @@ -630,7 +630,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 27, "metadata": {}, "outputs": [ { @@ -645,9 +645,6 @@ } ], "source": [ - "import seaborn as sns\n", - "sns.set_style('ticks')\n", - "sns.set_context('talk')\n", "f, axes = plt.subplots(ncols=2, figsize=(8, 4))\n", "axes[0].plot(irasa_out.freqs, irasa_out.periodic[0,:])\n", "axes[0].set_ylabel('Power (a.u.)')\n", @@ -661,51 +658,133 @@ "axes[1].set_xlabel('Frequency (Hz)')\n", "axes[1].set_title('Original + \\n Aperiodic Spectrum')\n", "\n", - "f.tight_layout()\n", - "f.savefig('../simulations/example_knee.png')\n" + "f.tight_layout()" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 26, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "| ch_name | cf | bw | pw |\n", - "|----------:|-----:|--------:|-------:|\n", - "| 0 | 9.5 | 1.44337 | 0.4146 |\n" - ] + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
ch_namecfbwpw
009.51.4433730.4146
\n", + "
" + ], + "text/plain": [ + " ch_name cf bw pw\n", + "0 0 9.5 1.443373 0.4146" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ "# %% get periodic stuff\n", - "peaks = irasa_out.get_peaks()\n", - "md = peaks.to_markdown(index=False)\n", - "print(md)" + "irasa_out.get_peaks()" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 25, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "| mse | r_squared | BIC | AIC | fit_type | ch_name |\n", - "|------------:|------------:|---------:|---------:|:-----------|----------:|\n", - "| 3.02402e-05 | 0.999894 | -2049.69 | -2062.86 | knee | 0 |\n" - ] + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
OffsetKneeExponent_1Exponent_2fit_typeKnee Frequency (Hz)ch_name
01.380983532.9097840.5119991.894478knee8.595540
\n", + "
" + ], + "text/plain": [ + " Offset Knee Exponent_1 Exponent_2 fit_type Knee Frequency (Hz) \\\n", + "0 1.380983 532.909784 0.511999 1.894478 knee 8.59554 \n", + "\n", + " ch_name \n", + "0 0 " + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "ap = irasa_out.fit_aperiodic_model(fit_func='knee').gof\n", - "md = ap.to_markdown(index=False)\n", - "print(md)" + "irasa_out.fit_aperiodic_model(fit_func='knee').aperiodic_params" ] }, { diff --git a/pyrasa/irasa_mne/mne_objs.py b/pyrasa/irasa_mne/mne_objs.py index c815349..b3e1c59 100644 --- a/pyrasa/irasa_mne/mne_objs.py +++ b/pyrasa/irasa_mne/mne_objs.py @@ -11,8 +11,6 @@ from pyrasa.utils.peak_utils import get_peak_params from pyrasa.utils.types import AperiodicFit -# FutureWarning: - class PeriodicSpectrumArray(SpectrumArray): """Subclass of SpectrumArray""" @@ -100,6 +98,7 @@ def plot_topo( def get_peaks( self: SpectrumArray, + smooth: bool = True, smoothing_window: float | int = 1, cut_spectrum: tuple[float, float] = (1, 40), peak_threshold: float = 2.5, @@ -108,26 +107,56 @@ def get_peaks( peak_width_limits: tuple[float, float] = (0.5, 6), ) -> pd.DataFrame: """ - This method can be used to extract peak parameters from the periodic spectrum extracted from IRASA. - The algorithm works by smoothing the spectrum, zeroing out negative values and - extracting peaks based on user specified parameters. - - Parameters: smoothing window : int, optional, default: 2 - Smoothing window in Hz handed over to the savitzky-golay filter. - cut_spectrum : tuple of (float, float), optional, default (1, 40) - Cut the periodic spectrum to limit peak finding to a sensible range - peak_threshold : float, optional, default: 1 - Relative threshold for detecting peaks. This threshold is defined in - relative units of the periodic spectrum - min_peak_height : float, optional, default: 0.01 - Absolute threshold for identifying peaks. The threhsold is defined in relative - units of the power spectrum. Setting this is somewhat necessary when a - "knee" is present in the data as it will carry over to the periodic spctrum in irasa. - peak_width_limits : tuple of (float, float), optional, default (.5, 12) - Limits on possible peak width, in Hz, as (lower_bound, upper_bound) - - Returns: df_peaks: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel + Extracts peak parameters from the periodic spectrum obtained via IRASA. + + This method identifies and extracts peak parameters such as center frequency (cf), bandwidth (bw), + and peak height (pw) from a periodic spectrum using scipy's find_peaks function. + The spectrum can be optionally smoothed prior to the peak detection. + + Parameters + ---------- + smooth : bool, optional + Whether to smooth the spectrum before peak extraction. Smoothing can help in reducing noise and + better identifying peaks. Default is True. + smoothing_window : int or float, optional + The size of the smoothing window in Hz, passed to the Savitzky-Golay filter. Default is 1 Hz. + polyorder : int, optional + The polynomial order for the Savitzky-Golay filter used in smoothing. The polynomial order must be + less than the window length. Default is 1. + cut_spectrum : tuple of (float, float) or None, optional + Tuple specifying the frequency range (lower_bound, upper_bound) to which the spectrum should be cut + before peak extraction. If None, peaks are detected across the full frequency range. Default is None. + peak_threshold : float, optional + Relative threshold for detecting peaks, defined as a multiple of the standard deviation of the + filtered spectrum. Default is 1.0. + min_peak_height : float, optional + The minimum peak height (in absolute units of the power spectrum) required for a peak to be recognized. + This can be useful for filtering out noise or insignificant peaks, especially when a "knee" is present + in the original data, which may persist in the periodic spectrum. Default is 0.01. + peak_width_limits : tuple of (float, float), optional + The lower and upper bounds for peak widths, in Hz. This helps in constraining the peak detection to + meaningful features. Default is (0.5, 12.0). + + Returns + ------- + pd.DataFrame + A DataFrame containing the detected peak parameters for each channel. The DataFrame includes the + following columns: + - 'ch_name': Channel name + - 'cf': Center frequency of the peak + - 'bw': Bandwidth of the peak + - 'pw': Peak height (power) + + + Notes + ----- + The function works by first optionally smoothing the periodic spectrum using a Savitzky-Golay filter. + Then, it performs peak detection using the `scipy.signal.find_peaks` function, taking into account the + specified peak thresholds and width limits. Peaks that do not meet the minimum height requirement are + filtered out. + + The `cut_spectrum` parameter can be used to focus peak detection on a specific frequency range, which is + particularly useful when the region of interest is known in advance. """ @@ -135,6 +164,7 @@ def get_peaks( self.get_data(), self.freqs, self.ch_names, + smooth=smooth, smoothing_window=smoothing_window, cut_spectrum=cut_spectrum, peak_threshold=peak_threshold, @@ -179,23 +209,53 @@ def fit_aperiodic_model( fit_bounds: tuple[float, float] | None = None, ) -> AperiodicFit: """ - This method can be used to extract aperiodic parameters from the aperiodic spectrum extracted from IRASA. - The algorithm works by applying one of two different curve fit functions and returns the associated parameters, - as well as the respective goodness of fit. - - Parameters: - fit_func : string - Can be either "fixed" or "knee". - fit_bounds : None, tuple - Lower and upper bound for the fit function, - should be None if the whole frequency range is desired. - Otherwise a tuple of (lower, upper) - - Returns: df_aps: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel - df_gof: DataFrame - DataFrame containing the goodness of fit of the specific fit function for each channel. - + Computes aperiodic parameters from the aperiodic spectrum using scipy's curve fitting function. + + This method can be used to model the aperiodic (1/f-like) component of the power spectrum. Per default, + users can choose between a fixed or knee model fit or specify their own fit method see examples + custom_fit_functions.ipynb for an example. + The method returns the fitted parameters for each channel along with some goodness of fit metrics. + + Parameters + ---------- + fit_func : str or type[AbstractFitFun], optional + The fitting function to use. Can be "fixed" for a linear fit or "knee" for a fit that includes a + knee (bend) in the spectrum or a class that is inherited from AbstractFitFun. The default is 'fixed'. + ch_names : Iterable or None, optional + Channel names corresponding to the aperiodic spectrum. If None, channels will be named numerically + in ascending order. Default is None. + scale : bool, optional + Whether to scale the data to improve fitting accuracy. This is useful in cases where + power values are very small (e.g., 1e-28), which may lead to numerical precision issues during fitting. + After fitting, the parameters are rescaled to match the original data scale. Default is False. + fit_bounds : tuple[float, float] or None, optional + Tuple specifying the lower and upper frequency bounds for the fit function. If None, the entire frequency + range is used. Otherwise, the spectrum is cropped to the specified bounds. Default is None. + + Returns + ------- + AperiodicFit + An object containing two pandas DataFrames: + - aperiodic_params : pd.DataFrame + A DataFrame containing the fitted aperiodic parameters for each channel. + - gof : pd.DataFrame + A DataFrame containing the goodness of fit metrics for each channel. + + Notes + ----- + This function fits the aperiodic component of the power spectrum using scipy's curve fitting function. + The fitting can be performed using either a simple linear model ('fixed') or a more complex model + that includes a "knee" point, where the spectrum bends. The resulting parameters can help in + understanding the underlying characteristics of the aperiodic component in the data. + + If the `fit_bounds` parameter is used, it ensures that only the specified frequency range is considered + for fitting, which can be important to avoid fitting artifacts outside the region of interest. + + The `scale` parameter can be crucial when dealing with data that have extremely small values, + as it helps to mitigate issues related to machine precision during the fitting process. + + The function asserts that the input data are of the correct type and shape, and raises warnings + if the first frequency value is zero, as this can cause issues during model fitting. """ return compute_aperiodic_model( @@ -306,6 +366,7 @@ def plot_topo( def get_peaks( self: EpochsSpectrumArray, + smooth: bool = True, smoothing_window: float | int = 1, cut_spectrum: tuple[float, float] = (1.0, 40.0), peak_threshold: float = 2.5, @@ -314,26 +375,56 @@ def get_peaks( peak_width_limits: tuple[float, float] = (0.5, 6.0), ) -> pd.DataFrame: """ - This method can be used to extract peak parameters from the periodic spectrum extracted from IRASA. - The algorithm works by smoothing the spectrum, zeroing out negative values and - extracting peaks based on user specified parameters. - - Parameters: smoothing window : int, optional, default: 2 - Smoothing window in Hz handed over to the savitzky-golay filter. - cut_spectrum : tuple of (float, float), optional, default (1, 40) - Cut the periodic spectrum to limit peak finding to a sensible range - peak_threshold : float, optional, default: 1 - Relative threshold for detecting peaks. This threshold is defined in - relative units of the periodic spectrum - min_peak_height : float, optional, default: 0.01 - Absolute threshold for identifying peaks. The threhsold is defined in relative - units of the power spectrum. Setting this is somewhat necessary when a - "knee" is present in the data as it will carry over to the periodic spctrum in irasa. - peak_width_limits : tuple of (float, float), optional, default (.5, 12) - Limits on possible peak width, in Hz, as (lower_bound, upper_bound) - - Returns: df_peaks: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel + Extracts peak parameters from the periodic spectrum obtained via IRASA. + + This method identifies and extracts peak parameters such as center frequency (cf), bandwidth (bw), + and peak height (pw) from a periodic spectrum using scipy's find_peaks function. + The spectrum can be optionally smoothed prior to the peak detection. + + Parameters + ---------- + smooth : bool, optional + Whether to smooth the spectrum before peak extraction. Smoothing can help in reducing noise and + better identifying peaks. Default is True. + smoothing_window : int or float, optional + The size of the smoothing window in Hz, passed to the Savitzky-Golay filter. Default is 1 Hz. + polyorder : int, optional + The polynomial order for the Savitzky-Golay filter used in smoothing. The polynomial order must be + less than the window length. Default is 1. + cut_spectrum : tuple of (float, float) or None, optional + Tuple specifying the frequency range (lower_bound, upper_bound) to which the spectrum should be cut + before peak extraction. If None, peaks are detected across the full frequency range. Default is None. + peak_threshold : float, optional + Relative threshold for detecting peaks, defined as a multiple of the standard deviation of the + filtered spectrum. Default is 1.0. + min_peak_height : float, optional + The minimum peak height (in absolute units of the power spectrum) required for a peak to be recognized. + This can be useful for filtering out noise or insignificant peaks, especially when a "knee" is present + in the original data, which may persist in the periodic spectrum. Default is 0.01. + peak_width_limits : tuple of (float, float), optional + The lower and upper bounds for peak widths, in Hz. This helps in constraining the peak detection to + meaningful features. Default is (0.5, 12.0). + + Returns + ------- + pd.DataFrame + A DataFrame containing the detected peak parameters for each channel. The DataFrame includes the + following columns: + - 'ch_name': Channel name + - 'cf': Center frequency of the peak + - 'bw': Bandwidth of the peak + - 'pw': Peak height (power) + + + Notes + ----- + The function works by first optionally smoothing the periodic spectrum using a Savitzky-Golay filter. + Then, it performs peak detection using the `scipy.signal.find_peaks` function, taking into account the + specified peak thresholds and width limits. Peaks that do not meet the minimum height requirement are + filtered out. + + The `cut_spectrum` parameter can be used to focus peak detection on a specific frequency range, which is + particularly useful when the region of interest is known in advance. """ @@ -346,6 +437,7 @@ def get_peaks( cur_epoch, self.freqs, self.ch_names, + smooth=smooth, smoothing_window=smoothing_window, cut_spectrum=cut_spectrum, peak_threshold=peak_threshold, @@ -405,23 +497,53 @@ def fit_aperiodic_model( fit_bounds: tuple[float, float] | None = None, ) -> AperiodicFit: """ - This method can be used to extract aperiodic parameters from the aperiodic spectrum extracted from IRASA. - The algorithm works by applying one of two different curve fit functions and returns the associated parameters, - as well as the respective goodness of fit. - - Parameters: - fit_func : string - Can be either "fixed" or "knee". - fit_bounds : None, tuple - Lower and upper bound for the fit function, - should be None if the whole frequency range is desired. - Otherwise a tuple of (lower, upper) - - Returns: df_aps: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel - df_gof: DataFrame - DataFrame containing the goodness of fit of the specific fit function for each channel. - + Computes aperiodic parameters from the aperiodic spectrum using scipy's curve fitting function. + + This method can be used to model the aperiodic (1/f-like) component of the power spectrum. Per default, + users can choose between a fixed or knee model fit or specify their own fit method see examples + custom_fit_functions.ipynb for an example. + The method returns the fitted parameters for each channel along with some goodness of fit metrics. + + Parameters + ---------- + fit_func : str or type[AbstractFitFun], optional + The fitting function to use. Can be "fixed" for a linear fit or "knee" for a fit that includes a + knee (bend) in the spectrum or a class that is inherited from AbstractFitFun. The default is 'fixed'. + ch_names : Iterable or None, optional + Channel names corresponding to the aperiodic spectrum. If None, channels will be named numerically + in ascending order. Default is None. + scale : bool, optional + Whether to scale the data to improve fitting accuracy. This is useful in cases where + power values are very small (e.g., 1e-28), which may lead to numerical precision issues during fitting. + After fitting, the parameters are rescaled to match the original data scale. Default is False. + fit_bounds : tuple[float, float] or None, optional + Tuple specifying the lower and upper frequency bounds for the fit function. If None, the entire frequency + range is used. Otherwise, the spectrum is cropped to the specified bounds. Default is None. + + Returns + ------- + AperiodicFit + An object containing two pandas DataFrames: + - aperiodic_params : pd.DataFrame + A DataFrame containing the fitted aperiodic parameters for each channel. + - gof : pd.DataFrame + A DataFrame containing the goodness of fit metrics for each channel. + + Notes + ----- + This function fits the aperiodic component of the power spectrum using scipy's curve fitting function. + The fitting can be performed using either a simple linear model ('fixed') or a more complex model + that includes a "knee" point, where the spectrum bends. The resulting parameters can help in + understanding the underlying characteristics of the aperiodic component in the data. + + If the `fit_bounds` parameter is used, it ensures that only the specified frequency range is considered + for fitting, which can be important to avoid fitting artifacts outside the region of interest. + + The `scale` parameter can be crucial when dealing with data that have extremely small values, + as it helps to mitigate issues related to machine precision during the fitting process. + + The function asserts that the input data are of the correct type and shape, and raises warnings + if the first frequency value is zero, as this can cause issues during model fitting. """ event_dict = {val: key for key, val in self.event_id.items()} diff --git a/pyrasa/utils/irasa_spectrum.py b/pyrasa/utils/irasa_spectrum.py index a7d79ae..c4abd0e 100644 --- a/pyrasa/utils/irasa_spectrum.py +++ b/pyrasa/utils/irasa_spectrum.py @@ -3,6 +3,7 @@ from attrs import define from pyrasa.utils.aperiodic_utils import compute_aperiodic_model +from pyrasa.utils.fit_funcs import AbstractFitFun from pyrasa.utils.peak_utils import get_peak_params from pyrasa.utils.types import AperiodicFit @@ -16,27 +17,59 @@ class IrasaSpectrum: ch_names: np.ndarray | None def fit_aperiodic_model( - self, fit_func: str = 'fixed', scale: bool = False, fit_bounds: tuple[float, float] | None = None + self, + fit_func: str | type[AbstractFitFun] = 'fixed', + scale: bool = False, + fit_bounds: tuple[float, float] | None = None, ) -> AperiodicFit: """ - This method can be used to extract aperiodic parameters from the aperiodic spectrum extracted from IRASA. - The algorithm works by applying one of two different curve fit functions and returns the associated parameters, - as well as the respective goodness of fit. - - Parameters: - fit_func : string - Can be either "fixed" or "knee". - fit_bounds : None, tuple - Lower and upper bound for the fit function, - should be None if the whole frequency range is desired. - Otherwise a tuple of (lower, upper) - - Returns: AperiodicFit - df_aps: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel - df_gof: DataFrame - DataFrame containing the goodness of fit of the specific fit function for each channel. + Computes aperiodic parameters from the aperiodic spectrum using scipy's curve fitting function. + + This method can be used to model the aperiodic (1/f-like) component of the power spectrum. Per default, + users can choose between a fixed or knee model fit or specify their own fit method see examples + custom_fit_functions.ipynb for an example. + The method returns the fitted parameters for each channel along with some goodness of fit metrics. + + Parameters + ---------- + fit_func : str or type[AbstractFitFun], optional + The fitting function to use. Can be "fixed" for a linear fit or "knee" for a fit that includes a + knee (bend) in the spectrum or a class that is inherited from AbstractFitFun. The default is 'fixed'. + ch_names : Iterable or None, optional + Channel names corresponding to the aperiodic spectrum. If None, channels will be named numerically + in ascending order. Default is None. + scale : bool, optional + Whether to scale the data to improve fitting accuracy. This is useful in cases where + power values are very small (e.g., 1e-28), which may lead to numerical precision issues during fitting. + After fitting, the parameters are rescaled to match the original data scale. Default is False. + fit_bounds : tuple[float, float] or None, optional + Tuple specifying the lower and upper frequency bounds for the fit function. If None, the entire frequency + range is used. Otherwise, the spectrum is cropped to the specified bounds. Default is None. + + Returns + ------- + AperiodicFit + An object containing two pandas DataFrames: + - aperiodic_params : pd.DataFrame + A DataFrame containing the fitted aperiodic parameters for each channel. + - gof : pd.DataFrame + A DataFrame containing the goodness of fit metrics for each channel. + + Notes + ----- + This function fits the aperiodic component of the power spectrum using scipy's curve fitting function. + The fitting can be performed using either a simple linear model ('fixed') or a more complex model + that includes a "knee" point, where the spectrum bends. The resulting parameters can help in + understanding the underlying characteristics of the aperiodic component in the data. + + If the `fit_bounds` parameter is used, it ensures that only the specified frequency range is considered + for fitting, which can be important to avoid fitting artifacts outside the region of interest. + + The `scale` parameter can be crucial when dealing with data that have extremely small values, + as it helps to mitigate issues related to machine precision during the fitting process. + The function asserts that the input data are of the correct type and shape, and raises warnings + if the first frequency value is zero, as this can cause issues during model fitting. """ return compute_aperiodic_model( aperiodic_spectrum=self.aperiodic, @@ -57,26 +90,56 @@ def get_peaks( peak_width_limits: tuple[float, float] = (0.5, 12), ) -> pd.DataFrame: """ - This method can be used to extract peak parameters from the periodic spectrum extracted from IRASA. - The algorithm works by smoothing the spectrum, zeroing out negative values and - extracting peaks based on user specified parameters. - - Parameters: smoothing window : int, optional, default: 2 - Smoothing window in Hz handed over to the savitzky-golay filter. - cut_spectrum : tuple of (float, float), optional, default (1, 40) - Cut the periodic spectrum to limit peak finding to a sensible range - peak_threshold : float, optional, default: 1 - Relative threshold for detecting peaks. This threshold is defined in - relative units of the periodic spectrum - min_peak_height : float, optional, default: 0.01 - Absolute threshold for identifying peaks. The threhsold is defined in relative - units of the power spectrum. Setting this is somewhat necessary when a - "knee" is present in the data as it will carry over to the periodic spctrum in irasa. - peak_width_limits : tuple of (float, float), optional, default (.5, 12) - Limits on possible peak width, in Hz, as (lower_bound, upper_bound) - - Returns: df_peaks: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel + Extracts peak parameters from the periodic spectrum obtained via IRASA. + + This method identifies and extracts peak parameters such as center frequency (cf), bandwidth (bw), + and peak height (pw) from a periodic spectrum using scipy's find_peaks function. + The spectrum can be optionally smoothed prior to the peak detection. + + Parameters + ---------- + smooth : bool, optional + Whether to smooth the spectrum before peak extraction. Smoothing can help in reducing noise and + better identifying peaks. Default is True. + smoothing_window : int or float, optional + The size of the smoothing window in Hz, passed to the Savitzky-Golay filter. Default is 1 Hz. + polyorder : int, optional + The polynomial order for the Savitzky-Golay filter used in smoothing. The polynomial order must be + less than the window length. Default is 1. + cut_spectrum : tuple of (float, float) or None, optional + Tuple specifying the frequency range (lower_bound, upper_bound) to which the spectrum should be cut + before peak extraction. If None, peaks are detected across the full frequency range. Default is None. + peak_threshold : float, optional + Relative threshold for detecting peaks, defined as a multiple of the standard deviation of the + filtered spectrum. Default is 1.0. + min_peak_height : float, optional + The minimum peak height (in absolute units of the power spectrum) required for a peak to be recognized. + This can be useful for filtering out noise or insignificant peaks, especially when a "knee" is present + in the original data, which may persist in the periodic spectrum. Default is 0.01. + peak_width_limits : tuple of (float, float), optional + The lower and upper bounds for peak widths, in Hz. This helps in constraining the peak detection to + meaningful features. Default is (0.5, 12.0). + + Returns + ------- + pd.DataFrame + A DataFrame containing the detected peak parameters for each channel. The DataFrame includes the + following columns: + - 'ch_name': Channel name + - 'cf': Center frequency of the peak + - 'bw': Bandwidth of the peak + - 'pw': Peak height (power) + + + Notes + ----- + The function works by first optionally smoothing the periodic spectrum using a Savitzky-Golay filter. + Then, it performs peak detection using the `scipy.signal.find_peaks` function, taking into account the + specified peak thresholds and width limits. Peaks that do not meet the minimum height requirement are + filtered out. + + The `cut_spectrum` parameter can be used to focus peak detection on a specific frequency range, which is + particularly useful when the region of interest is known in advance. """ diff --git a/pyrasa/utils/irasa_tf_spectrum.py b/pyrasa/utils/irasa_tf_spectrum.py index 7727da5..26de7b5 100644 --- a/pyrasa/utils/irasa_tf_spectrum.py +++ b/pyrasa/utils/irasa_tf_spectrum.py @@ -3,6 +3,7 @@ from attrs import define from pyrasa.utils.aperiodic_utils import compute_aperiodic_model_sprint +from pyrasa.utils.fit_funcs import AbstractFitFun from pyrasa.utils.peak_utils import get_peak_params_sprint from pyrasa.utils.types import AperiodicFit @@ -19,28 +20,56 @@ class IrasaTfSpectrum: ch_names: np.ndarray | None def fit_aperiodic_model( - self, fit_func: str = 'fixed', scale: bool = False, fit_bounds: tuple[float, float] | None = None + self, + fit_func: str | type[AbstractFitFun] = 'fixed', + scale: bool = False, + fit_bounds: tuple[float, float] | None = None, ) -> AperiodicFit: """ - This method can be used to extract aperiodic parameters from the aperiodic spectrum extracted from IRASA. - The algorithm works by applying one of two different curve fit functions and returns the associated parameters, - as well as the respective goodness of fit. - - Parameters: - fit_func : string - Can be either "fixed" or "knee". - fit_bounds : None, tuple - Lower and upper bound for the fit function, - should be None if the whole frequency range is desired. - Otherwise a tuple of (lower, upper) - - Returns: AperiodicFit - df_aps: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel - df_gof: DataFrame - DataFrame containing the goodness of fit of the specific fit function for each channel. + Extracts aperiodic parameters from the aperiodic spectrogram using scipy's curve fitting + function. + + This function computes aperiodic parameters for each time point in the spectrogram by applying either one of + two different curve fitting functions (`fixed` or `knee`) or a custom function specified by user to the data. + See examples custom_fit_functions.ipynb. The parameters, along with the goodness of + fit for each time point, are returned in a concatenated format. + + Parameters + ---------- + fit_func : str or type[AbstractFitFun], optional + The fitting function to use. Can be "fixed" for a linear fit or "knee" for a fit that includes a + knee (bend) in the spectrum or a class that is inherited from AbstractFitFun. The default is 'fixed'.. + scale : bool, optional + Whether to scale the data to improve fitting accuracy. This is useful when fitting a knee in cases where + power values are very small, leading to numerical precision issues. Default is False. + fit_bounds : tuple[float, float] or None, optional + Tuple specifying the lower and upper frequency bounds for the fit function. If None, the entire frequency + range is used. Otherwise, the spectrum is cropped to the specified bounds before fitting. Default is None. + + Returns + ------- + AperiodicFit + An object containing two pandas DataFrames: + - aperiodic_params : pd.DataFrame + A DataFrame containing the aperiodic parameters (e.g., center frequency, bandwidth, peak height) + for each channel and each time point. + - gof : pd.DataFrame + A DataFrame containing the goodness of fit metrics for each channel and each time point. + + Notes + ----- + This function iterates over each time point in the provided spectrogram to extract aperiodic parameters + using the specified fit function. It leverages the `compute_aperiodic_model` function for individual fits + at each time point, then combines the results across all time points into comprehensive DataFrames. + + The `fit_bounds` parameter allows for frequency range restrictions during fitting, which can help in focusing + the analysis on a particular frequency band of interest. + + Scaling the data using the `scale` parameter can be particularly important when dealing with very small power + values that might lead to poor fitting due to numerical precision limitations. """ + return compute_aperiodic_model_sprint( aperiodic_spectrum=self.aperiodic[np.newaxis, :, :] if self.aperiodic.ndim == min_ndim else self.aperiodic, freqs=self.freqs, @@ -62,26 +91,56 @@ def get_peaks( peak_width_limits: tuple[float, float] = (0.5, 12), ) -> pd.DataFrame: """ - This method can be used to extract peak parameters from the periodic spectrum extracted from IRASA. - The algorithm works by smoothing the spectrum, zeroing out negative values and - extracting peaks based on user specified parameters. - - Parameters: smoothing window : int, optional, default: 2 - Smoothing window in Hz handed over to the savitzky-golay filter. - cut_spectrum : tuple of (float, float), optional, default (1, 40) - Cut the periodic spectrum to limit peak finding to a sensible range - peak_threshold : float, optional, default: 1 - Relative threshold for detecting peaks. This threshold is defined in - relative units of the periodic spectrum - min_peak_height : float, optional, default: 0.01 - Absolute threshold for identifying peaks. The threhsold is defined in relative - units of the power spectrum. Setting this is somewhat necessary when a - "knee" is present in the data as it will carry over to the periodic spctrum in irasa. - peak_width_limits : tuple of (float, float), optional, default (.5, 12) - Limits on possible peak width, in Hz, as (lower_bound, upper_bound) - - Returns: df_peaks: DataFrame - DataFrame containing the center frequency, bandwidth and peak height for each channel + Extracts peak parameters from a periodic spectrogram obtained via IRASA. + + This method processes a time-resolved periodic spectrum to identify and extract peak parameters such as + center frequency (cf), bandwidth (bw), and peak height (pw) for each time point. It applies smoothing, + peak detection, and thresholding according to user-defined parameters + (see get_peak_params for additional Information). + + Parameters + ---------- + smooth : bool, optional + Whether to smooth the spectrum before peak extraction. Smoothing can help in reducing noise and better + identifying peaks. Default is True. + smoothing_window : int or float, optional + The size of the smoothing window in Hz, passed to the Savitzky-Golay filter. Default is 2 Hz. + polyorder : int, optional + The polynomial order for the Savitzky-Golay filter used in smoothing. The polynomial order must be less + than the window length. Default is 1. + cut_spectrum : tuple of (float, float) or None, optional + Tuple specifying the frequency range (lower_bound, upper_bound) to which the spectrum should be cut before + peak extraction. If None, the full frequency range is used. Default is (1, 40). + peak_threshold : int or float, optional + Relative threshold for detecting peaks, defined as a multiple of the standard deviation of the filtered + spectrum. Default is 1. + min_peak_height : float, optional + The minimum peak height (in absolute units of the power spectrum) required for a peak to be recognized. + This can be useful for filtering out noise or insignificant peaks, especially when a "knee" is present + in the data. Default is 0.01. + peak_width_limits : tuple of (float, float), optional + The lower and upper bounds for peak widths, in Hz. This helps in constraining the peak detection to + meaningful features. Default is (0.5, 12). + + Returns + ------- + pd.DataFrame + A DataFrame containing the detected peak parameters for each channel and time point. The DataFrame + includes the following columns: + - 'ch_name': Channel name + - 'cf': Center frequency of the peak + - 'bw': Bandwidth of the peak + - 'pw': Peak height (power) + - 'time': Corresponding time point for the peak + + Notes + ----- + This function iteratively processes each time point in the spectrogram, applying the `get_peak_params` + function to extract peak parameters at each time point. The resulting peak parameters are combined into + a single DataFrame. + + The function is particularly useful for analyzing time-varying spectral features, such as in dynamic or + non-stationary M/EEG data, where peaks may shift in frequency, bandwidth, or amplitude over time. """