fooof icon indicating copy to clipboard operation
fooof copied to clipboard

[ENH] Aperiodic Knee Fit

Open ryanhammonds opened this issue 3 years ago • 3 comments

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).

ryanhammonds avatar Apr 18 '22 22:04 ryanhammonds

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!

TomDonoghue avatar Apr 19 '22 14:04 TomDonoghue

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

ryanhammonds avatar Apr 19 '22 18:04 ryanhammonds

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

ryanhammonds avatar Jun 16 '22 22:06 ryanhammonds