Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Aperiodic Knee Fit #235

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Conversation

ryanhammonds
Copy link
Contributor

This updates the knee parameter to knee frequency, so that the parameter is more interpretable and can be more easily bounded to a frequency range of interest (#224). I've also included an new aperiodic knee mode which adds a constant to the current aperiodic fit (#234).

@TomDonoghue
Copy link
Member

Awesome! I still need to come back and review this properly, but the idea is great! This kind of tapering off at high frequencies is something I see fairly often, and capturing that like this seems like it could be really useful.

Practically, I wasn't planning a new release before the name shift (#205), so we should coordinate to merge those two together into the new major release, and then this can be the start of an increased set of fit_functions in the new release!

@ryanhammonds
Copy link
Contributor Author

Sounds good! Once #205 is ready, I'll resolve conflicts and rebase.

@fooof-tools fooof-tools deleted a comment from codecov-commenter May 28, 2022
@ryanhammonds
Copy link
Contributor Author

I have been using this model and wanted to make sure the knee / timescale I was getting wasn't incorrect since the knee location can visually appear to drift far from the knee frequency. I'm leaving this code here for myself to go back to add accuracy tests (and for when this PR is reviewed if it helps).

import numpy as np
import matplotlib.pyplot as plt
from fooof.core.funcs import expo_const_function

def lin_interp(x, y, i, half):
    """Interpolate at zerocrossing."""
    return x[i] + (x[i+1] - x[i]) * ((half - y[i]) / (y[i+1] - y[i]))

def compute_knee(x, y):
    """Compute knee as fwhm in linear space."""
    y = y - y.min()
    half = max(y)/2.0
    signs = np.sign(np.add(y, -half))
    zero_crossing = np.where((signs[0:-2] != signs[1:-1]))[0][0]
    knee_freq = lin_interp(x, y, zero_crossing, half)
    
    return knee_freq

# Settings
xs = np.arange(1, 200)
offset = 0
knee_freq = 20
exp = 2

# Check knee and plot
fig, axes = plt.subplots(ncols=2, figsize=(12, 5), sharex=True)

for ind, const in enumerate(np.linspace(1e-10, 1e-2, 40)):
    
    # Get ys
    ys = 10**expo_const_function(xs, offset, knee_freq, exp, const)
    
    # Plot
    if ind < 20:
        axes[0].loglog(xs, ys)
    else:
        axes[1].loglog(xs, ys)
     
    # Accuracy check
    knee = compute_knee(xs, ys)
    assert knee.round() == knee_freq
    
 
axes[0].axvline(knee_freq, color='k', ls='--')
axes[1].axvline(knee_freq, color='k', ls='--')

test

@TomDonoghue TomDonoghue added the 2.0 Targetted for the specparam 2.0 release. label Jun 28, 2023
@@ -172,7 +172,7 @@ def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_h
# Guess parameters for aperiodic fitting, [offset, knee, exponent]
# If offset guess is None, the first value of the power spectrum is used as offset guess
# If exponent guess is None, the abs(log-log slope) of first & last points is used
self._ap_guess = (None, 0, None)
self._ap_guess = (None, 1, None)
# Bounds for aperiodic fitting, as: ((offset_low_bound, knee_low_bound, exp_low_bound),
# (offset_high_bound, knee_high_bound, exp_high_bound))
# By default, aperiodic fitting is unbound, but can be restricted here, if desired
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if

self._ap_bounds = ((-np.inf, -np.inf, -np.inf), (np.inf, np.inf, np.inf))

should also be updated, as the knee parameters shouldn't be below 0 in any case?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
2.0 Targetted for the specparam 2.0 release.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants