Consistency with results from literature

[joeybarreto]

I've been trying to understand the output of harminv on a specific time series, and can't seem to understand why the output is worse than what I can verify from the literature. I've installed it and can verify that the sines example works. Following the example at the end of this paper, there is specific test signal to process, with the following frequencies and decay rates

ws = np.array([32.85, 36.51, 38.29, 48.82, 105.63, 146.96, 183.24, 191.10, 217.34, 241.39, 253.59, 283.60, 298.20, 299.85, 304.95])
gs = np.array([4.61, 9.44, 12.95, 0.09, 9.24, 14.88, 10.19, 7.46, 11.21, 2.18, 11.03, 9.64, 4.42, 12.25, 9.20])

which are both given in rad/s. The amplitudes are all 1, and the phases are all 0. I can reproduce the exact time series provided in Fig. 1 (by inspection), and also its spectral form in Fig. 2. With

ts = np.linspace(0,1,128)
dt = ts[1]-ts[0]
fmin = np.min(ws)/(2*np.pi)
fmax = np.max(ws)/(2*np.pi)

I run harminv, but the output shows no converged modes, contrary to the results of that paper which achieve numerically exact results with those parameter settings (they don't provide the code that they used). The Nyquist frequency is 63.5 Hz and the max frequency is 48.53, so the criterion holds, and the "information uncertainty principle" $N > 4 \pi \rho(\omega) / dt \equiv N_c$ also holds, where $N = 128$ and $N_c \approx 87.97$ based on the data and values above. I computed $\rho(\omega) = \text{num of modes} / \omega_{max} - \omega_{min} \approx 0.055$ modes per rad/s. If increase $T$ and $N$ by some large multiplicative factor , then $dt$ remains unchanged but the uncertainty principle holds by a much larger margin and I start seeing some converged modes, but many are still missed. Can anyone verify this behavior on their end for me? I'm trying to understand if there is some bug on my end or if harminv is behaving differently than the Fortran implementation alluded to in that paper, which seems to work very well.