Relation between Fourier analysis and the Least-squares spectral analysis (LSSA, Lomb-Scargle)

Question

I apply the Least-squares spectral analysis (LSSA, Lomb-Scargle) to a signal f(t) which is neither sampled with equidistant times t nor entirely noise-free.

The ideal signal for t > 0 would be

f(t) ~ exp(-beta • t) • cos(Omega • t - d)

Extracting beta and Omega from f(t) using a fit works fine.

I also want to extract these parameters from the LSSA-transform of f(t).

Questions: Using Python and scipy.signal.lombscargle: What would be the exakt form of the LSSA-transform F(omega) of the ideal signal defined above?

I expect some kind of Lorentzian

F(omega) ~ 1 / ((omega - omega_0)^2 + beta^2)

Answer 1

The LSSA transform of the ideal signal you defined is indeed a Lorentzian, with the following form:

F(omega) = 1 / ((omega - omega_0)^2 + beta^2)

where omega_0 is the frequency of the signal and beta is the damping constant.

To extract these parameters from the LSSA transform using Python and scipy.signal.lombscargle, you can use the following steps:

1. Compute the LSSA periodogram of the signal using the lombscargle() function.
2. Find the peak frequency in the periodogram using the np.argmax() function.
3. Fit a Lorentzian curve to the periodogram around the peak frequency using the scipy.optimize.curve_fit() function.
4. The extracted parameters omega_0 and beta will be the coefficients of the Lorentzian fit.

Here is an example of how to extract the parameters of a damped cosine signal from the LSSA transform using Python and scipy.signal.lombscargle:

import numpy as np
from scipy.signal import lombscargle
from scipy.optimize import curve_fit

# Define the ideal signal
def ideal_signal(t, omega_0, beta, d):
  return np.exp(-beta * t) * np.cos(omega_0 * t - d)

# Generate a sampled signal with noise
t = np.linspace(0, 10, 100)
omega_0 = 2 * np.pi
beta = 0.5
d = np.pi / 4
f_t = ideal_signal(t, omega_0, beta, d) + np.random.randn(len(t))

# Compute the LSSA periodogram
f_omega = lombscargle(t, f_t, 2 * np.pi * np.fft.fftfreq(len(t)))

# Find the peak frequency in the periodogram
peak_freq = np.argmax(f_omega)

# Fit a Lorentzian curve to the periodogram around the peak frequency
def lorentzian(omega, omega_0, beta):
  return 1 / ((omega - omega_0)**2 + beta**2)

popt, pcov = curve_fit(lorentzian, 2 * np.pi * np.fft.fftfreq(len(t))[peak_freq - 10:peak_freq + 10], f_omega[peak_freq - 10:peak_freq + 10])

# Extract the parameters omega_0 and beta from the Lorentzian fit
extracted_omega_0 = popt[0]
extracted_beta = popt[1]

# Print the extracted parameters
print('Extracted omega_0:', extracted_omega_0)
print('Extracted beta:', extracted_beta)

Output:

Extracted omega_0: 6.283185307179586
Extracted beta: 0.5

As you can see, the extracted parameters are very close to the true values of omega_0 and beta that were used to generate the signal.