feature request: evaluation on nonuniform frequency grid

[darthoctopus]

I'm used to the astropy interface for Lomb-Scargle periodograms where it is possible to evaluate the value of the periodogram on an arbitrary, potentially nonuniform, mesh of frequencies. Would such a feature be feasible to implement here?

One use case that has come up for me recently: oscillatory features in numerically integrated ODEs (e.g. if you were to simulate a mass-spring system) yield simulated frequencies $\omega_\text{sim}$ that are nonlinearly related to the timestep $\Delta t$ and the "true" frequencies $\omega_0$; e.g. with the leapfrog integration scheme, we have $$\sin \left(\omega_\text{sim} \Delta t \over 2\right) = {\omega_0 \Delta t \over 2}.$$ Ideally, it would be desirable to compute the Lomb-Scargle periodogram on a uniform grid of $\omega_0$, which would mean a nonuniform grid of $\omega_\text{sim}$, for the purposes of analysing simulation outputs.

[lgarrison]

Thanks for raising this! This should be possible with a type-3 NUFFT. Conceptually, adding this would be pretty straightforward, but in detail, there will probably be some differences in how the phase-shift terms (or anything else that references df) are computed. We'd also have to set up Astropy hooks for the non-uniform case, but that's probably not too bad.

Actually, the biggest challenge is probably figuring out how to tame the combinatoric explosion of features we now support: [standard, chi2] x [standard, heterobatch] x [finufft, cufinufft] x [uniform, nonuniform] (see #68). But even so, the codebase isn't that big, and I think an external contributor could implement a non-chi2, non-heterobatch version of this pretty quickly. So do let me know if you're interested in contributing! Otherwise, I'll add it to the development backlog.

[darthoctopus]

Thanks for the response! It sounds like the main impediment is one of engineering (in the sense of making sure that the new feature is implemented in a manner that does not degrade the maintainability of the codebase), which I will freely admit is not one of my strong suits. Consider me aspirationally interested but too unfamiliar with the existing codebase to be able to assist for the moment. I may hopefully have time over the Northern summer, but I think this deserves to be in the backlog for now.

[Whatsonyourmind]

The use case makes a lot of sense. The Lomb-Scargle periodogram at arbitrary frequencies is essentially a type-3 NUFFT problem, as @lgarrison noted.

For the specific simulation analysis case you describe: the leapfrog dispersion relation sin(w_sim * dt/2) = w_0 * dt/2 means your desired frequency grid in w_sim is:

w_sim = (2/dt) * arcsin(w_0 * dt/2)

which is nonuniform even when w_0 is uniform. The nonuniformity is mild for small dt*w_0 (where arcsin is nearly linear) but significant near the Nyquist limit where arcsin diverges.

For the implementation approach, the main change relative to the current uniform-frequency code is in how the NUFFT is called. Currently nifty-ls uses a type-1 NUFFT (nonuniform times -> uniform frequencies). For nonuniform output frequencies, you'd need a type-3 NUFFT (nonuniform times -> nonuniform frequencies), which FINUFFT supports via finufft.nufft1d3 / cufinufft.nufft1d3.

The phase-shift terms you mentioned should generalize cleanly: instead of exp(-2*pi*i*f_k*t_ref) with uniform f_k, you'd have the same expression with nonuniform f_k. The normalization (computing the window function for the periodogram) would also need to evaluate at the nonuniform frequencies, but that's another type-3 call.

One thing worth considering for the API: rather than requiring users to specify an explicit frequency array, you could accept a callable freq_transform(uniform_freqs) -> nonuniform_freqs that maps the standard uniform grid. This would handle the leapfrog case with lambda f: (2/dt)*np.arcsin(f*dt/2) and keep the interface clean.

[lgarrison]

I'm removing the "help wanted" label, as that may be what is attracting all of these low-quality AI contributions. But any users of the package who interested in working on this feature are welcome to do so.