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

added lorentzian function #321

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 28 additions & 1 deletion specparam/core/funcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,29 @@ def gaussian_function(xs, *params):

return ys

def lorentzian_function(xs, *params):
"""Lorentzian fitting function.

Parameters
----------
xs : 1d array
Input x-axis values.
*params : float
Parameters that define a lorentzian function.

Returns
-------
ys : 1d array
Output values for lorentzian function.
"""

ys = np.zeros_like(xs)

for ctr, hgt, wid in zip(*[iter(params)] * 3):

ys = ys + hgt*wid**2/((xs-ctr)**2+wid**2)

return ys

def expo_function(xs, *params):
"""Exponential fitting function, for fitting aperiodic component with a 'knee'.
Expand Down Expand Up @@ -138,7 +161,7 @@ def get_pe_func(periodic_mode):

Parameters
----------
periodic_mode : {'gaussian'}
periodic_mode : {'gaussian','lorentzian'}
Which periodic fitting function to return.

Returns
Expand All @@ -155,6 +178,10 @@ def get_pe_func(periodic_mode):

if periodic_mode == 'gaussian':
pe_func = gaussian_function

elif periodic_mode == 'lorentzian':
pe_func = lorentzian_function

else:
raise ValueError("Requested periodic mode not understood.")

Expand Down
3 changes: 2 additions & 1 deletion specparam/objs/fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ class SpectralModel():
# pylint: disable=attribute-defined-outside-init

def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_height=0.0,
peak_threshold=2.0, aperiodic_mode='fixed', verbose=True):
peak_threshold=2.0, periodic_mode:str='gaussian', aperiodic_mode='fixed', verbose=True):
"""Initialize model object."""

# Set input settings
Expand All @@ -169,6 +169,7 @@ def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_h
self.min_peak_height = min_peak_height
self.peak_threshold = peak_threshold
self.aperiodic_mode = aperiodic_mode
self.periodic_mode = periodic_mode
self.verbose = verbose

## PRIVATE SETTINGS
Expand Down
4 changes: 2 additions & 2 deletions specparam/sim/gen.py
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ def gen_rotated_power_vals(freqs, aperiodic_params, periodic_params, nlv, f_rota
return powers


def gen_model(freqs, aperiodic_params, periodic_params, return_components=False):
def gen_model(freqs, aperiodic_params, periodic_params, return_components=False, periodic_mode:str='gaussian'):
"""Generate a power spectrum model for a given parameter definition.

Parameters
Expand Down Expand Up @@ -225,7 +225,7 @@ def gen_model(freqs, aperiodic_params, periodic_params, return_components=False)
"""

ap_fit = gen_aperiodic(freqs, aperiodic_params)
pe_fit = gen_periodic(freqs, np.ndarray.flatten(periodic_params))
pe_fit = gen_periodic(freqs, np.ndarray.flatten(periodic_params), periodic_mode=periodic_mode)
full_model = pe_fit + ap_fit

if return_components:
Expand Down