0

编辑：为了澄清问题，我专门寻找一个例程，它可以在周期性数据的傅里叶分解中准确计算 phi1、phi2、phi3 等，给定一系列形式

m(t) = A0 + sum_{n=1}^{N} A_i * sin(2*pi*n*f + phi_n)

正如下面所讨论的，我已经探索了一些关于 Stack Overflow 的相关讨论，但似乎都没有产生准确的结果（除非我在我的代码或数学中犯了错误）。

我正在尝试将傅里叶级数拟合到一些周期性变星数据，我想要从中获得特定的模型参数。此处的方程式 (1) 或 (2)或任何等效项都是有用的（我为使用链接而道歉；我不允许发布图片，并且除了纯文本之外，看不到如何排版方程式）。例如，这些结果将用于计算 RRab Lyrae 变星的金属丰度，如本文所述。本质上，这需要准确确定 phi_31=phi_3-3*phi_1，在这种情况下使用傅里叶级数的正弦形式（等式 2）。

这个网站上有很多关于拟合傅立叶级数的问题，但我认为没有一个答案能完全满足我的需要。这里的答案实际上使用 symfit 对上面链接的 Deb 和 Singh 论文中的方程 (1) 进行建模。使用这种方法，我得到了一个非常适合 RRab Lyrae 数据的模型，但经验结果似乎与此相去甚远，表明该模型没有找到正确的参数。我还将该代码修改为正弦傅立叶级数，最佳拟合模型会产生不同的结果，所以我认为 symfit 做得不够好（我的代码对于这篇文章来说有点长，但我会如果有人愿意，很乐意分享）。

jakevdp在这里给出的答案表明 symfit 陷入了局部最小值，这对我来说是有道理的。他讨论的 Lomb-Scargle 方法做得很好，由于他的答案中包含的代码有一些被贬低的功能，我将对其进行更新。首先，我们需要一些数据来处理。这是来自 Catalina Surveys DR2 的 RRab Lyrae 变量 SV Hya 的一些光度数据。下面的代码大致按照 jakevdp 的代码拉入数据，围绕恒星的实际周期（0.47855 天）绘制 LS 周期图，并拉出最适合的 LS 模型参数：

import numpy as np
data=np.genfromtxt('result_web_filePAevp8.csv', delimiter=',')
time=data[1:281,5]
signal=data[1:281,1]
signalerror=data[1:281,2]
import matplotlib.pyplot as plt
plt.figure()
from astropy.timeseries import LombScargle
ls = LombScargle(time, signal, signalerror, nterms=5)
freq, power = ls.autopower(minimum_frequency=2.05,maximum_frequency=2.15,samples_per_peak=10)
plt.plot(freq, power);
plt.show()
best_frequency = freq[np.argmax(power)]
theta=ls.model_parameters(best_frequency)
print(theta)

这真的很酷，但我的问题是我不知道如何，正如 jakevdp 所说，

这些线性正弦/余弦幅度可以通过一点三角函数转换回非线性幅度和相位。

我看不出如何将此处描述的 Lomb-Scargle 模型重写为 Deb 和 Singh 的等式 (1) 或 (2)。LS 模型有 N+1 项，而不是 2N+1，所以我想需要 2*nterms LS 模型来产生 N 项 Deb 和 Singh 型模型。但是，如果您查看链接中的模型，很明显，设置 cos(a)=sin(a+pi/2) 使我们能够以类似于 Deb 和 Singh 中的等式 (2) 的形式重写模型，仅phi 参数固定为偶数项为 0，奇数项为 pi/2。因此，该模型似乎先验地约束了我想要确定的实际参数。

所以，我的问题是：（1）我的数学错了，我能把 LS 模型结果转换成能满足我需要的形式吗？(2) 如果没有，有没有办法将 LS 模型中的级数编辑为更一般的 2N+1 项形式？(3) 如果没有，我可以使用另一种方法将一般 2N+1 项傅里叶级数拟合到这些数据中，以便提取 phi1、phi3 等吗？

4

1 回答 1

0

标记为“编辑”的问题是光谱分析的一个众所周知的结果。这是一个简单的例程，将计算 phi1、phi2 等。（请注意，下面的函数计算“完整”幅度谱，因此一半功率在负频率，但相位与您的兴趣相关，并且不受影响。）

import numpy as np
import matplotlib.pyplot as plt

def fft(t, x):
    fftspec = np.fft.fft(x)
    freqs = np.fft.fftfreq(x.size, np.mean(np.ediff1d(t)))
    amp = np.abs(fftspec)/len(fftspec)
    phase = np.arctan2(np.imag(fftspec), np.real(fftspec))
    #unfuck the phases to only show up
    #when amplitude is above some threshold
    phase = phase*(amp>1e-12)
    return freqs, amp, phase

if __name__ == "__main__":
    fs = 200.0 #Hz
    ampl1, ampl2   = 5, 7 #Volts
    freq1, freq2   = 15, 22 #Hz
    phase1, phase2 = 1, 2 #radians
    t = np.arange(0,10,1/fs) #timesteps, 10 seconds
    #2-component signal
    x1 = ampl1*np.cos(2*np.pi*freq1*t+phase1) #first component
    x2 = ampl2*np.cos(2*np.pi*freq2*t+phase2) #second component
    x = x1 + x2
    freqs, amp, phase = fft(t, x)
    fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3,
                                        figsize = (15,6))
    ax0.plot(t,x)
    ax0.set_title('Input Data')
    ax0.set_ylabel('Volts')
    ax0.set_xlabel('Time[s]')

    ax1.plot(freqs, amp)
    ax1.set_title('Full Amplitude Spectrum')
    ax1.set_ylabel('Volts')
    ax1.set_xlabel('Freq[Hz]')

    ax2.plot(freqs, phase)
    ax2.set_title('Phase spectrum')
    ax2.set_ylabel('Phase')
    ax2.set_xlabel('Freq[Hz]')

    plt.show()
于 2021-09-06T00:01:01.240 回答