Knee Parameter

[ryanhammonds]

@TomDonoghue do you think it would be useful to add an option to fit a knee frequency parameter? I have a fork where I've converted y = 10^offset * (1/(knee + x^exp)) to y = 10^offset * (1/(knee_freq^exp + x^exp)). This is based on the knee to knee_freq conversion: knee_freq = knee**(1./exponent). I think this makes the knee parameter a more interpretable feature of spectra, and it allows passing bounds like (1, 40 Hz) to curve_fit if one wants to limit the range of knee freq locations.

[TomDonoghue]

Ooops - my apologies, I apparently missed this issue / question!

I think updating the knee from the 'parameter' to the 'frequency' generally makes sense. This is a significant change to what we report, but could be bundled and described into the 2.0 major release, which is a breaking change anyways.

I wonder if there is any particular difference / benefit to fitting the knee frequency directly, versus continuing to fit the parameter and then converting post hoc to frequency to update the value we return? We could do either - I'm just not sure if there is any particular reason to prefer one or the other.

Also - tagging in knee expert @rdgao for any comments!

[rdgao]

chiming in here:

the knee parameter is technically a "well-defined" quantity that just falls out of the "generalized Lorentzian". Like, if you plug that number back into the equation, you get what you want out.

knee frequency, on the other hand, is sort of an "interpretation", if you will. It's "roughly where the spectrum falls off from the plateau", but it's a rather vague definition, as opposed to something like "knee frequency is where the power is half max in the Lorentzian (or whatever)". In other words, if you fit knee_freq^exp instead, the relative power at knee_freq will probably vary.

In the case exp=2, there is a correspondence, but for any other exponents, it's more like a guideline. I just kinda hacked it and extrapolated from exp=2 to in general take the inverse power.

All that being said, I think if you just defined knee_freq = knee**exp formally, then there's no difference whether you take the inverse power before or after, theoretically. I guess it remains to be tested whether doing it that way will then impact your exponent fit in some weird way, or perhaps make the knee fitting more stable bc now it's a pre-exponentiated number so perhaps more normally distributed.

(clearing my inbox and just found this, sorry!)

[rdgao]

actually to be more precise: the correspondence I meant for exp=2 refers to the fact that the PSD of an exponential decay with a time constant of T is:

PSD( exp(-(t/T)) ) ~ 1/((1/T)**2 + f**2). So technically speaking, this allows you to interpret the knee as the timescale exactly, by inserting this dummy "knee frequency" variable because both (1/T) and f are squared, so they have the same units (and 1/T is inverse of time which is frequency, so all conveniently works out even though there's no formal definition of "knee frequency").

when the fitted exponent is not 2, you can't go the inverse direction and find T in an exact way (or at least I don't know the analytical correspondence). But this doesn't stop you from taking that first thing as f_k**exp = knee still and take the inverse power of that as the "knee frequency" since it will then still have the same unit as f. So afaik, there's nothing theoretically wrong with this, and I'd probably make the decision of which to fit on the numerical stability of fitting knee vs. f_k**exp

[ryanhammonds]

Thanks for the follow up @Rdago. All of this makes sense. I ran tests and don't think one model is consistently more numerically stable than the other. There's tradeoff in stability depends on combinations of parameters (knee, exp, f_range, oscillations, curve_fit bounds, etc) and the difference in knee and exponent error is on the order of 1e-16.

[chapochn]

Hello,
I would also have the need to limit the position of the knee, and that would work well with the proposed definition of the knee parameter above:

y = 10^offset * (1/(knee_freq^exp + x^exp))

Would anyone have an implementation/fork of that already?

[chapochn]

Seems I can just include the changes from here in my code?
48b7cb7

[ryanhammonds]

@chapochn A fork for this is here. It's old and needs updating, but should work.

[chapochn]

Thanks! I got it to work with the code. Can't wait for this to be integrated in a new version.