Python改进傅里叶级数参数求解：周期变星数据相位提取问题

AI 小北

2026-10-5

问题描述

编辑说明：为明确问题，我需找到可对周期数据做傅里叶分解、精准计算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，需使用正弦形式的傅里叶级数。

本站有大量关于傅里叶级数拟合的问题，但没有符合我需求的答案。我见过使用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型模型。但查看LS模型定义可知，通过cos(a)=sin(a+pi/2)可以将模型改写为和Deb和 Singh论文中公式(2)近似的形式，只是偶数项的phi参数固定为0、奇数项的phi参数固定为pi/2，相当于对我需要求解的参数做了先验约束。

因此我有三个问题：

是否我的计算有误，LS模型的结果其实可以转换为我需要的参数形式？
如果不能，是否有办法修改LS模型中的级数，使其适配更通用的2N+1项形式？
如果也不行，有没有其他方法可以对这些数据拟合通用2N+1项傅里叶级数，从而提取phi1、phi3等参数？

问题解答

首先澄清你对Lomb-Scargle模型的误解：你提到的参数数量问题是对astropy实现的理解偏差，LombScargle类中nterms=N对应的模型形式为：

y(t) = theta0 + sum_{n=1}^N [ theta_{2n-1} * sin(2π n f t) + theta_{2n} * cos(2π n f t) ]

你输出的theta数组长度是2*N +1，正好对应你需要的2N+1个参数，不存在参数缺失的问题。

问题1解答：LS模型参数可以直接转换为你需要的正弦形式参数

转换逻辑非常简单，利用三角恒等式a*sinx + b*cosx = A*sin(x+phi)，其中：

振幅A_n = sqrt(theta_{2n-1}^2 + theta_{2n}^2)
相位phi_n = arctan2(theta_{2n}, theta_{2n-1})
你需要的phi_31 = phi_3 - 3*phi_1直接用转换后的相位计算即可，不需要额外修改模型。

补充一段转换代码，直接接在你现有代码后面运行即可：

def convert_ls_params(theta, nterms):
    A0 = theta[0]
    A_list = []
    phi_list = []
    for n in range(1, nterms+1):
        sin_coeff = theta[2*n -1]
        cos_coeff = theta[2*n]
        A_n = np.sqrt(sin_coeff**2 + cos_coeff**2)
        phi_n = np.arctan2(cos_coeff, sin_coeff)
        A_list.append(A_n)
        phi_list.append(phi_n)
    return A0, A_list, phi_list

A0, A_list, phi_list = convert_ls_params(theta, nterms=5)
phi1 = phi_list[0]
phi3 = phi_list[2]
phi_31 = phi3 - 3*phi1
print(f"phi_31 = {phi_31}")

问题2解答：不需要修改LS模型

你之前误以为LS模型固定了相位，是对模型形式的误解，实际上它的正弦和余弦项的系数都是自由拟合的，没有任何相位先验约束，完全适配你需要的通用级数形式。

问题3解答：如果要换其他方案，可以用scipy的curve_fit配合LS给出的初始值

如果担心LS的结果有偏差，可以用LS得到的参数作为初始值输入非线性拟合，避免局部最优：

from scipy.optimize import curve_fit

def fourier_model(t, A0, *args):
    nterms = len(args)//2
    res = A0
    for n in range(nterms):
        A = args[2*n]
        phi = args[2*n+1]
        res += A * np.sin(2*np.pi*(n+1)*best_frequency*t + phi)
    return res

# 用LS转换后的参数作为初始值
p0 = [A0] + [val for pair in zip(A_list, phi_list) for val in pair]
popt, pcov = curve_fit(fourier_model, time, signal, p0=p0, sigma=signalerror, absolute_sigma=True)
# popt里直接就是你要的A0、A1、phi1、A2、phi2...参数

内容的提问来源于stack exchange，提问作者Daryl Janzen