这一切都在Python完成。 我有一个带有二氧化碳数据的数据框架，相关列是“年”（CE年）和“ CO2”（在PPM中）。数据范围从1959年到2024年，这是我观察到的数据。我

Question

# Build model based on equation
def co2_model(x, a, b, A, phi, c):
  # Model containing linear, quadratic and sinusoidal component 
    return (a*x) + (b*(x**2)) + (A*np.sin((2*np.pi*x)+ phi)) + c 

x = (co2['time']).values - co2['time'].values[0]
y = co2['co2'].values

# Split data into training and validation periods 
co2_train = co2[co2['year'] <= 2004]
co2_val = co2[co2['year'] > 2004]

# Rescaling time periods and specifying dependent variable 
x_train = co2_train['year'].values - co2_train['year'].values[0]
y_train = co2_train['co2'].values

x_val = co2_val['year'].values - co2_val['year'].values[0]
y_val = co2_val['co2'].values

# Fit the co2 model to the training data
params_train, covar_train = curve_fit(co2_model, x_train, y_train)
a, b, A, phi, c = params_train  # Extract fitted parameters

#Calculate predicted values for the training and validation data
y_train_predict = co2_model(x_train, a, b, A, phi, c)
y_val_predict = co2_model(x_val, a, b, A, phi, c)

# Plot the original data with the predictions 
plt.figure(figsize=(15, 6))
plt.scatter(x, y, color='grey', label='Observed temperature', alpha=0.5)

# Plot the training period predictions
plt.plot(x_train, y_train_predict, label='Training Period', color='blue', linestyle='--')

# Plot the validation period predictions
plt.plot(x_val, y_val_predict, label='Validation Period', color='red', linestyle='--')

plt.xlabel('Year (CE)')
plt.ylabel('CO$_2$ (ppm)') 

# Show the plot
plt.show()

这就是我得到的

我得到的图

，但是我期望的是类似的东西，但一半的红色和另一个蓝色不同的曲线

请让我知道我是否只是犯了一个愚蠢的错误

ok，我相信您可能正在绘制每月的均值数据来自

Https：//gml.noaa.gov/ccgg/trends/data.html 我自由下载链接“ mauna loa二氧化碳月度均值数据（文本）”，然后将其删除到桌子上（加上两个标头线）。我认为您的数据位于“小数日期”和“每月平均”列中。 （这里有一个很小的缺陷，用于随后的傅立叶分析 - 几个月的长度完全相同。您可以纠正它 - 请参阅下面的代码 - 但它可以忽略不计。）

）

然后，我进行了FFT，并查看了频谱（下图）。您会收到年度信号（每年F = 1次）和适度的季节性谐波（每年F = 2次）。

Answer 1

I使用FFT从数据中剥离了年度和季节性的谐波正弦，然后对剩余数据进行了标准二次拟合。

curve_fit

输出中给出了完全拟合。 （如果您只想使用一年的正弦体，则不必添加谐波。）

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import rfft
from scipy.optimize import curve_fit

PI = math.pi

#---------------------------------------------------------------------

def amplitudePhaseAngle( p, q ):
    '''
    Convert   p.cos(h)+q.sin(h)   to    A.sin(h+phi)
    '''
    A = math.sqrt( p ** 2 + q ** 2 )
    phi = math.atan2( p, q )
    return A, phi

#---------------------------------------------------------------------

def plotSpectrum( freq, absft ):
    '''
    Plot the absolute values of the Fourier coefficients
    '''
    plt.plot( freq[:len(absft)], absft, '-o' )
    plt.grid()
    plt.show()

#---------------------------------------------------------------------

def fitSeries( N, nfreq, ft, w, t ):
   '''
   Fit a (truncated) Fourier series, using only the frequency indices in nfreq[]
   '''
   fit = np.zeros_like( t )
   for n in nfreq:
       if n == 0 or ( N % 2 == 0 and n == N // 2 ):
           A, phi = amplitudePhaseAngle( (1/N) * ft[n].real, 0.0 )
       else:
           A, phi = amplitudePhaseAngle( (2/N) * ft[n].real, -(2/N) * ft[n].imag )
       print( "Adding amplitude, circular frequency, phase =", A, w[n], phi )
       fit += A * np.sin( w[n] * t + phi )
   return fit

#---------------------------------------------------------------------

def quadfit( t, a, b, c ):
    return a + b * t + c * t ** 2

#---------------------------------------------------------------------


############
# Get data #
############
t, CO2 = np.loadtxt( 'co2_mm_mlo.txt', skiprows=2, usecols=[2,3], unpack=True )
t -= t[0]
N = len( t )
#t = np.linspace( 0, t[-1], N )         # use this if you want to correct for the fact that months have different length!
dt = ( t[-1] - t[0] ) / ( N - 1 )


########################################
# Remove the annual periodic component #
########################################
ft = rfft( CO2 )
ft[0] = 0
absft = np.abs( ft )

w = np.linspace( 0, 2 * PI / dt, N, endpoint=False )
dw = w[1] - w[0]
annual   = int( 2 * PI     / dw + 0.5 )          # index of nearest frequency to 1 time per year
harmonic = int( 2 * PI * 2 / dw + 0.5 )          # index of nearest frequency to 2 times per year
Fourier = fitSeries( N, [annual, harmonic], ft, w, t )


##########################################
# Now do a quadratic fit to what is left #
##########################################
CO2rev = CO2 - Fourier
p0 = ( CO2rev[0], 0.0, 0.0 )
popt, pcov = curve_fit( quadfit, t, CO2rev, p0=p0 )
print( 'Quadratic fit a + b.t + c.t^2;  a, b, c = ', *popt )


##################################
# Compare fit with original data #
##################################
plt.plot( t, CO2, 'o', color='gray', markersize=3, label='Observed' )
plt.plot( t, Fourier + quadfit( t, *popt ), 'b-', label='Fitted' )
plt.legend()
plt.show()


##################
# Check spectrum #
##################
plotSpectrum( w / ( 2 * PI ), absft )                      # plot against annual frequency
plotSpectrum( range( len(absft) ), absft )                 # plot against index of frequency