傅里叶分析与最小二乘谱分析（LSSA、Lomb-Scargle）之间的关系

Question

我将最小二乘谱分析（LSSA、Lomb-Scargle）应用于信号 f(t)，该信号既不是以等距时间 t 采样的，也不是完全无噪声的。

t > 0 的理想信号是

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

使用拟合从 f(t) 中提取 beta 和 Omega 效果很好。

我还想从 f(t) 的 LSSA 变换中提取这些参数。

问题：使用Python和scipy.signal.lombscargle：上面定义的理想信号的LSSA变换F(omega)的exakt形式是什么？

我期待某种洛伦兹

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

Answer 1

您定义的理想信号的 LSSA 变换确实是洛伦兹变换，其形式如下：

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

其中

omega_0

是信号频率，

beta

是阻尼常数。

要使用 Python 和 scipy.signal.lombscargle 从 LSSA 转换中提取这些参数，您可以使用以下步骤：

1. 使用

   lombscargle()

   函数计算信号的 LSSA 周期图。
2. 使用

   np.argmax()

   函数查找周期图中的峰值频率。
3. 使用

   scipy.optimize.curve_fit()

   函数将洛伦兹曲线拟合到峰值频率周围的周期图。
4. 提取的参数

   omega_0

   和

   beta

   将是洛伦兹拟合的系数。

这里是如何使用 Python 和 scipy.signal.lombscargle 从 LSSA 变换中提取阻尼余弦信号参数的示例：

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)

输出：

Extracted omega_0: 6.283185307179586
Extracted beta: 0.5

如您所见，提取的参数非常接近用于生成信号的

omega_0

和

beta

的真实值。