In [ ]:

    

import numpy as np
TSM = np.array([ 7.6,  7.4,  8.2,  9.2, 10.2, 11.5, 12.4, 13.4, 13.7, 11.8, 10.1, 9.0,
                 8.9,  9.5, 10.6, 11.4, 12.9, 12.7, 13.9, 14.2, 13.5, 11.4, 10.9, 8.1])

months = np.arange(0, 24) + 1


    

In [ ]:

    

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(months, TSM)
_ = ax.set_xlabel(r'Temperatura [$^\circ$C]')


    

In [ ]:

    

N = len(TSM)
A0, B0 = (2. / N) * sum(TSM), 0


    

In [ ]:

    

twopi = 2 * np.pi
y = TSM - TSM.mean()

A = np.zeros(12)
B = np.zeros(12)

for p in np.arange(1, N/2 + 1):
    for n in range(0, N):
        A[p-1] += 2. / N * y[n] * np.cos(twopi * p * (n+1) / N)

for p in np.arange(1, N/2 + 1):
    for n in range(0, N):
        B[p-1] += 2. / N * y[n] * np.sin(twopi * p * (n+1) / N)


    

In [ ]:

    

pmax = 1
tsm1 = np.zeros_like(TSM)
for n in range(0, N):
    tsm1[n] = A0 / 2
    for p in range(1, pmax+1):
        tsm1[n] += (A[p-1] * np.cos(twopi * p * n / N) +
                    B[p-1] * np.sin(twopi * p * n / N))

pmax = 2
tsm2 = np.zeros_like(TSM)
for n in range(0, N):
    tsm2[n] = A0 / 2
    for p in range(1, pmax+1):
        tsm2[n] += (A[p-1] * np.cos(twopi * p * n / N) +
                    B[p-1] * np.sin(twopi * p * n / N))
        
pmax = 3
tsm3 = np.zeros_like(TSM)
for n in range(0, N):
    tsm3[n] = A0 / 2
    for p in range(1, pmax+1):
        tsm3[n] += (A[p-1] * np.cos(twopi * p * n / N) +
                    B[p-1] * np.sin(twopi * p * n / N))


    

In [ ]:

    

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(months, tsm1, label='First harmonic')
ax.plot(months, tsm2, label='Second harmonic')
ax.plot(months, tsm3, label='Third harmonic')
ax.set_xlabel(r'Temperatura [$^\circ$C]')
_ = ax.legend(numpoints=1, loc='lower right')


    

In [ ]:

    

ynew = np.zeros_like(TSM)
pmax = 12
for n in range(0, N):
    ynew[n] = A0 / 2
    for p in range(1, pmax+1):
        ynew[n] += (A[p-1] * np.cos(twopi * p * n / N) +
                    B[p-1] * np.sin(twopi * p * n / N))


    

In [ ]:

    

fig, ax = plt.subplots()
ax.plot(months, TSM, label='original')
ax.plot(months, ynew, label='reconstructed')
ax.set_xlabel(r'Temperatura [$^\circ$C]')
_ = ax.legend(numpoints=1)


    

In [ ]:

    

p = np.arange(1, N / 2 + 1)
Cp = np.sqrt(A**2 + B**2)
fp = np.ones_like(Cp) / p[::-1]

fig, ax = plt.subplots()
#ax.bar(fp, Cp)
ax.plot(fp, Cp, '--k')
ax.axis([0.05, 1, 0, 3])
ax.set_title('Spectra')
ax.set_xlabel('Frequencies [cpm]')
_ = ax.set_ylabel('Amplitude')


    

In [ ]:

    

fig, ax = plt.subplots()
_ = ax.fill_between(p, Cp, color='0.5')