The current knee model is slightly incorrect translation of an exponentially decaying ACF, the difference is visually subtle but lies in the the power at the nyquist, e.g. the ACF(1) coeff. Below is a picture of a spectrum that matches what is simulated in OU processes. The power at the nyquist slightly tapers off but the current knee model assumes it should be linearly decreasing.

Below is the model we could be add. It's described in this paper

Where $w$ is the AR(1) or ACF(1) coefficient, which directly translates to the knee frequency, $f_k$. This math is also nice because it lets us compute the knee directly from the AR(1), which could be used as the guess knee param when oscillations aren't present:

# Get knee freq directly from signal
s1 = np.linalg.norm(sig) # singular value
e1 = s1**2 # eigenvalue
w = (1/e1) * sig[:-i] @ sig[i:] # AR(1) weight
fk = -(np.log(w) * fs) / (2 * np.pi) # knee freq

This is the full simulation/model/implementation:

import matplotlib.pyplot as plt
import numpy as np

fs = 1000

fk = 20 # Knee frequency
w = np.exp(-(fk * 2 * np.pi) / fs) # AR(1) coefficient
freqs = np.linspace(.1, 500, 1000)

# Update knee model
num = 1-w**2
den = 1 + w**2 - 2*w * np.cos(2*np.pi*freqs/fs)
powers = num / den

plt.loglog(freqs, powers);

In practice, this won't change much unless fitting up to the nyquist. It could be useful to help validate sims since sim_knee in neurodsp doesn't produce this updated PSD, and the mapping between knee freq and ACF(1) doesn't hold like it does in other timescale sims.