# 傅里叶模型
# 关于傅里叶级数模型
傅里叶级数将周期函数描述为正弦函数和余弦函数之和。您可以使用傅里叶级数将任意周期函数分解为简单的分量。这些分量易于积分、微分和分析。因此,傅里叶级数通常用于近似周期信号。
傅里叶级数有多种表示形式。曲线拟合工具箱使用三角傅里叶级数形式
其中
有关傅里叶级数的更多信息,请参见傅里叶分析和滤波。
# 在曲线拟合器中以交互方式拟合傅里叶模型
此示例展示如何使用曲线拟合器将傅里叶模型拟合到数据。
加载声音信号样本数据。
using TyCurveFitting
using TyPlot
using TyBase
filename = pkgdir(TyCurveFitting) * "/examples/docs/gong.mat"
load(filename)
变量 y 和 Fs 分别包含锣铃的声音信号和频率数据。通过将 y 的前 1000 个元素存储在名为 gongClip 的向量中来创建声音片段。
gongClip = y[1:1000];
要计算 gongClip 中每个元素对应的时间,将元素的索引除以 Fs。
t = (1:1000)./Fs;
打开曲线拟合器。
using TyCurveFitter
curveFitter()
在曲线拟合器中,选择拟合的数据变量。在曲线拟合器选项卡的数据部分,单击选择数据。在选择拟合数据对话框中,选择 t 作为 X 数据值,选择 gongClip 作为 Y 数据值。
应用程序会在您选择变量时绘制数据点。默认情况下,应用程序会将多项式拟合到数据中。要拟合傅里叶模型,请在曲线拟合器选项卡的拟合类型部分中单击"傅里叶"。
曲线拟合器拟合单项傅里叶模型。
拟合的单项傅里叶模型是具有简单振荡行为的周期函数。结果面板显示模型的一般方程、具有 95% 置信区间的拟合系数估计值、基频和拟合优度统计量。
拟合的单项傅里叶模型的均方根误差 (RMSE) 为 0.1996。要将单项傅里叶模型与具有四个项的傅里叶模型进行比较,请在拟合选项面板中将项数选择为 4。该应用程序将具有四个项的傅里叶模型拟合到数据中。
单击导出部分中的导出,然后选择导出到工作区,将拟合的四项傅里叶模型导出到工作区。
您可以使用函数 sound 来收听 gongClip 中的声音数据。
using TyDSPSystem
using TyBase
sound(gongClip,Fs)
pause(2) # 允许 gongClip 在执行下一行之前播放
要获取 gongClip 的傅里叶模型近似的声音数据,在 t 时刻评估 gongFourierModel。播放近似的声音数据。
gongClipApprox = gongFourierModel(t);
sound(gongClipApprox,Fs)
两个剪辑的近似平均音调相同。但是,近似声音数据的音调波动没有 gongClip 中的声音数据那么多。
# 在命令行中拟合傅里叶模型
此示例显示如何使用拟合函数将傅里叶模型拟合到数据。
拟合二项傅里叶模型
加载 El Niño-Southern Oscillation (ENSO) 数据。
using TyCurveFitting
using TyPlot
include(pkgdir(TyCurveFitting) * "/examples/docs/enso.jl")
变量 pressure 包含智利复活节岛和澳大利亚达尔文之间的平均气压差数据。变量 month 包含发生每个气压差的月份的数据。
绘制气压与月份的关系图。
plot(month,pressure)
压力数据在 0 到 18 之间波动,这表明可以用傅里叶级数来描述。
使用傅里叶库模型拟合两项傅里叶模型。指定模型类型为傅里叶,后跟项数。保存拟合优度统计数据以供以后比较。
f2 = fit("fourier2",month,pressure)
gof2 = f2.s_data
f2 =
常规模型 fourier2:
y(a0,a1,b1,a2,b2,w,x) = a0 + a1*cos(x*w) + b1*sin(x*w) + a2*cos(2*x*w) + b2*sin(2*x*w)
系数(置信边界为 95%):
a0 = 10.629238239654885 (10.227603403547397, 11.030873075762374)
a1 = 2.9231726647295493 (2.2704754102861098, 3.575869919172989)
b1 = 1.0586921519420507 (0.015925548515366872, 2.1014587553687347)
a2 = -0.5052309430300345 (-1.0857797958737512, 0.07531790981368225)
b2 = 0.21874179615470057 (-0.42014929006021184, 0.857632882369613)
w = 0.5258041509261316 (0.5222289410227758, 0.5293793608294874)
gof2 =
sse: 1122.998540270975
rsquare: 0.4279254551746283
dfe: 162.0
adjrsquare: 0.4102688334207588
rmse: 2.632886202394457
f2 是一个结构体,包含通用公式、具有 95% 置信区间的系数估计值以及拟合 w 的基频。a2 和 b2 的置信区间越过零,因此没有足够的证据来得出它们不同于零或拟合模型不同于单项傅里叶模型的结论。2.6329 的均方根误差 (RMSE) 可用于将 f2 的精度与其他拟合的精度进行比较。
要根据基频计算周期,请使用公式 T = 2*pi/w。
w = f2.params[end]
w = 0.5258041509261316
T = 2*pi/w
T = 11.94966851462207
拟合的二项傅里叶模型的周期约为 12 个月,即一年。
用数据散点图绘制 f2。
plotfit(f2,month,pressure)
f2 的形状与单项傅里叶模型的形状相似,振荡峰值大约每 12 个月出现一次。
拟合七项傅里叶模型
将七项傅里叶模型拟合到数据。保存拟合优度统计数据。
f7 = fit("fourier7",month,pressure)
gof7 = f7.s_data
f7 =
常规模型 fourier7:
y(a0,a1,b1,a2,b2,a3,b3,a4,b4,a5,b5,a6,b6,a7,b7,w,x) = a0 + a1*cos(x*w) + b1*sin(x*w) + a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w) + a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w) + a6*cos(6*x*w) + b6*sin(6*x*w) + a7*cos(7*x*w) + b7*sin(7*x*w)
系数(置信边界为 95%):
a0 = 10.626175900322497 (10.28312712810244, 10.969224672542556)
a1 = 0.566907981109863 (0.08284736189608055, 1.0509686003236454)
b1 = 0.1968777205762401 (-0.2900397779692975, 0.6837952191217777)
a2 = -1.203140600336192 (-1.6873783363047883, -0.7189028643675957)
b2 = -0.8084946163884884 (-1.3069496668200422, -0.31003956595693466)
a3 = 0.9322913113507187 (0.4324807907673119, 1.4321018319341257)
b3 = 0.7599040785837621 (0.26220357610276895, 1.2576045810647551)
a4 = -0.6653270807580326 (-1.1489756367147226, -0.18167852480134256)
b4 = -0.20375093826463636 (-0.6995043393129357, 0.2920024627836629)
a5 = -0.02913333093452589 (-0.5129370182753944, 0.45467035640634257)
b5 = -0.370111329342404 (-0.8566049383243401, 0.11638227963953213)
a6 = -0.04841078362788808 (-0.5437201844873867, 0.44689861723161056)
b6 = -0.13670328350175978 (-0.6286010201939334, 0.3551944531904139)
a7 = 2.811957474941592 (2.1904659741906034, 3.4334489756925803)
b7 = 1.3329899378027341 (0.40166526272427827, 2.26431461288119)
w = 0.07527069918816204 (0.07478387088796279, 0.07575752748836129)
gof7 =
sse: 768.3656006879204
rsquare: 0.608581502548569
dfe: 152.0
adjrsquare: 0.5699546771421777
rmse: 2.2483409097890505
f7 包含几个置信区间超过零的系数,因此没有足够的证据来得出结论,即它们的对应项会提高拟合傅里叶模型的准确性。2.2483 的 RMSE 小于 f2 的 RMSE 误差,这证实了七项傅里叶模型比二项傅里叶模型更准确地预测压力。
要根据基频计算周期,请使用公式 T = 2*pi/w 来计算周期。
w = f7.params[end]
w = 0.07527069918816204
T = (2*pi)/w
T = 83.4745176402952
拟合的七项傅里叶模型的周期约为 83 个月,或大约七年。拟合系数的幅度决定了哪些项对压力差的预测值贡献最大。
形式为 sin(Ax) 或 cos(Ax) 的正弦曲线的周期由公式 T = 2*pi/|A| 给出。a7 和 b7 是最大系数。
T = 2*pi/(w*7)
T = 11.924931091470743
a7 和 b7 对应的项的周期约为 12 个月,表明年周期最强。
使用相同的公式计算以下项的周期:
项 a1 和 b1 的周期均为 7 年;
项 a2 和 b2 的周期均为 3.5(7/2)年。a2 和 b2 系数的幅度大于 a1 和 b1,因此 3.5 年周期对压力差预测值的贡献大于 7 年周期;
项 a3 和 b3 较强,表明 2.3 年(7/3)周期。
较小的项(例如 a6、b6、a5 和 b5)对于拟合不太重要。
使用数据的散点图绘制 f7。
plotfit(f7,month,pressure)
与单项傅里叶模型相比,七项傅里叶模型以更复杂的模式振荡,并能捕捉到更广泛的压力差值。该周期大约每 84 个月或 7 年重复一次。通常,厄尔尼诺变暖以两到七年的不规则间隔发生,持续九个月到两年。平均周期长度为五年。模型结果反映了其中一些时期。
设置起点
拟合函数使用数据输入参数来计算系数和基频计算的优化起点。傅里叶级数模型对起点特别敏感,优化值可能仅对相关方程中的几个项准确。您可以通过指定 startpt 名称值参数来覆盖优化的起点。
数据散点图中的极值表明可能存在四年周期。要确认此建议,请将基频的起点设置为对应于八年或 96 个月周期的值。拟合的傅里叶模型的八年周期将项 a2 和 b2 的周期从 3.5 增加到 4。
w_8 = (2*pi)/96
w_8 = 0.06544984694978735
使用 coeffnames 函数在 f7 系数名称的单元格向量中查找基频的索引。
coeffnames(f7)
ans = 16-element Vector{String}:
"a0"
"a1"
"b1"
"a2"
"b2"
"a3"
"b3"
"a4"
"b4"
"a5"
"b5"
"a6"
"b6"
"a7"
"b7"
"w"
基频位于系数名称向量的最后一个条目中。根据 f7 的系数创建一个系数值向量,并将基频的值替换为对应于八年周期的值。
coeffs = coeffvalues(f7);
coeffs[end] = w_8;
coeffs
coeffs = 16-element Vector{Float64}:
10.626175900322497
0.566907981109863
0.1968777205762401
-1.203140600336192
-0.8084946163884884
0.9322913113507187
0.7599040785837621
-0.6653270807580326
-0.20375093826463636
-0.02913333093452589
-0.370111329342404
-0.04841078362788808
-0.13670328350175978
2.811957474941592
1.3329899378027341
0.06544984694978735
使用具有新基频值的系数作为起点,对压力差数据进行七项傅里叶模型拟合。保存拟合优度统计数据。
f7_8 = fit("fourier7",month,pressure,startpt=coeffs)
gof7_8 = f7_8.s_data
f7_8 =
常规模型 fourier7:
y(a0,a1,b1,a2,b2,a3,b3,a4,b4,a5,b5,a6,b6,a7,b7,w,x) = a0 + a1*cos(x*w) + b1*sin(x*w) + a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w) + a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w) + a6*cos(6*x*w) + b6*sin(6*x*w) + a7*cos(7*x*w) + b7*sin(7*x*w)
系数(置信边界为 95%):
a0 = 10.577077952243792 (10.052149016817237, 11.102006887670347)
a1 = 0.3285994158695221 (-0.43390523008621773, 1.0911040618252619)
b1 = -0.05917359822972016 (-0.7884162629091295, 0.6700690664496892)
a2 = -0.8667416510075757 (-1.7377440219709746, 0.004260719955823111)
b2 = 1.094012542633506 (0.28187448928510306, 1.9061505959819087)
a3 = -0.45238031477194096 (-1.2319777812970956, 0.32721715175321353)
b3 = -0.3117432158256887 (-1.0987962709441939, 0.4753098392928165)
a4 = 0.1810015125997952 (-0.7948764778577488, 1.156879503057339)
b4 = 0.5805870676286697 (-0.17957007445943196, 1.3407442097167714)
a5 = 0.032626029186113945 (-0.7174414205083444, 0.7826934788805722)
b5 = -0.22990658390721652 (-0.9766931057715873, 0.5168799379571543)
a6 = 0.3725835225469138 (-0.3899987065763351, 1.1351657516701628)
b6 = -0.274492997106749 (-1.1650960945353424, 0.6161101003218442)
a7 = 0.43093851207758027 (-0.49106022421793105, 1.3529372483730917)
b7 = -0.35466610679191746 (-1.315497042687062, 0.6061648291032271)
w = 0.06794680668718348 (0.06519252635634366, 0.07070108701802329)
gof7_8 =
sse: 1685.1430522148953
rsquare: 0.1415594856170891
dfe: 152.0
adjrsquare: 0.05684496117140703
rmse: 3.3296347320362645
f7_8 的系数与 f7 的系数略有不同。f7_8 的 RMSE 较高,表明 f7 更适合数据。绘制两个拟合图以直观地比较模型。
plotfit(f7_8,month,pressure)
hold("on")
plotfit(f7, "b")
hold("off")
legend(["Data","f7_8","f7"])
该图显示,f7 比 f7_8 更准确地捕捉了压力差数据的变化。
为了进一步研究傅里叶模型拟合,您可以尝试指定 NonlinearLeastSquares 拟合算法可用的不同选项。有关更多信息,请参见 fitoptions。
# 另请参阅
# APP
# 函数
fit | fittype | fitoptions