**6. Cycle C1 CO2 frequency analysis**

As an illustration of frequency analysis we examine Cycle C1 CO2 data in details. Similar approach can be used for other data, such as temperature, methane, oxygen and dust. The first concern is to make the C1 CO2 data more homogeneous hence the sampling time can be more precise and useful.

#### **6.1 Data filling**

The original Cycle 1 Vostok data contains a total of 80 CO2 data points are very unevenly distributed in time. This poses issues with FFT or DFT harmonic analysis. Hence we corrected the situation somewhat by inserting certain number of additional "data points" using linear approximation between the original data points with the biggest time difference between the neighboring points, hence generating a total of 128 original and inserted data points. In particular, the very first point we added to the original Vostok data is the current value of CO2 which we stated at 350 [25] for the year 1990, as an example. Current levels are estimated at over 400 ppm. Our CPE model can be used for a variety of sensitivity analysis by playing "what if" and processing the data through CPE.

Next data is the first original Vostok data, 2362 years in the past. This will correspond to 64 FFT harmonics described next. Note that the time period of C1 CO2 data is 128,399 years (**Table 2**), time difference between the oldest C1 data and the current time, **Figure 9** above.

#### **6.2 Sampling time**

For sampling time, we proceeded as follows. As stated in Section 2 the cycle duration is measured between the maximum (CO2 and temperature) points in each cycle. Hence, we have an average of a bit over 1000 years time sampling for Cycle 1, i.e. *<sup>T</sup>* <sup>¼</sup> 128, 399 <sup>128</sup>�<sup>1</sup> = 1011 years, which is used as the CO2 data sampling time. The FFT analysis on 128 total points produces 64 harmonics with the highest frequency of *<sup>f</sup> max* <sup>¼</sup> *<sup>f</sup>* <sup>2</sup> <sup>¼</sup> <sup>1</sup> <sup>2</sup>*<sup>T</sup>* <sup>¼</sup> <sup>1</sup> <sup>2022</sup> <sup>10</sup>�<sup>3</sup> <sup>¼</sup> <sup>0</sup>*:*<sup>000495</sup> *<sup>x</sup>* <sup>10</sup>�<sup>3</sup> Hz. For the analysis we use *<sup>T</sup>* <sup>¼</sup> <sup>1</sup>*:*<sup>011</sup> with the understanding it is in 1000s of years. **Table 7** summarizes Cycle 1 harmonic content ordered by the amplitude, from the largest to the smallest to identify each harmonic energy contribution. First harmonic H1 is boldfaced. As shown in Section 7.1, for the discretization purposes, it is important to satisfy an additional sampling time requirement given in inequality (7) bellow which produces equivalent resonant frequency of continuous and discrete models, and in turn makes for more precise discrete KFHO and KFHB models as well. This increases required sampling frequency and reduces time sampling period bellow 1011 years, and we deal with it using energy analysis as described in **Table 8**. For other approaches to non-uniform data analysis see [26, 27] as well as general literature on non uniform and optimal Fourier analysis.

#### **6.3 Amplitude and energy analysis**

**Table 7** shows that the average signal value (DC or H0) is 231 and the rest of the energy content (reflected in a corresponding harmonic amplitude) of the harmonics follows the highest-to-lowest amplitude pattern H1, H2, H4, H3, H7, H10, H5, and so on.

The last (boldfaced) harmonic in this list is H27. For a bit smaller cumulative amplitude, say 88%+, we only need 19 harmonics, with H37 as the last one. **Table 8** further summarizes energy analysis. For example, a chosen sampling time of 1011 satisfies the requirement in inequality (7) bellow all the way to cumulative energy

**Cumulative amplitude %**

 0.293645 86.369 15 1.083996 1.011 0.293645 88.1983 19 1.083996 1.011 0.46365 90.2476 24 0.686531 **0.650** 0.486832 96.6477 44 0.653839 **0.650**

**Total H's**

**Sampling Ts <**

**Chosen Ts**

*Boldfaced entries correspond to the first harmonic H1.*

*Summary of cycle 1 CO2 harmonic analysis.*

*Cycle 1 CO2 amplitude and sampling time.*

**Maximum frequency**

*Boldfaced entries point to the lowest limit on Sampling Time Ts.*

**Table 7.**

**Table 8.**

**49**

**Highest frequency**

**H No Frequency f Period 1/f Phase Amplitude Cumulative %** 0 N/A 0 231.0068 58.83917 **0.007727498 129.408 0.754987 36.27814 68.07948** 0.015454995 64.704 0.555383 15.87177 72.12214 0.03090999 32.352 1.11877 9.409128 74.51872 0.023182493 43.136 0.478204 6.815518 76.25468 0.054092483 18.48686 1.435168 5.576594 77.67508 0.077274975 12.9408 0.316939 5.098453 78.97369 0.038637488 25.8816 1.14182 4.847574 80.20841 0.09272997 10.784 0.20146 4.780349 81.426 0.046364985 21.568 2.02131 3.49376 82.31588 0.131367458 7.612235 0.471006 3.462987 83.19793 0.108184965 9.243429 0.58073 3.138693 83.99738 0.115912463 8.6272 0.489932 2.959198 84.75111 0.200914936 4.977231 0.08324 2.30261 85.3376 0.193187438 5.17632 0.527865 2.031371 85.85501 0.293644906 3.405474 0.48454 2.018039 86.36902 0.247279921 4.044 0.171 1.841643 86.8381 0.085002473 11.76436 0.041037 1.827369 87.30354 0.262734916 3.806118 0.371439 1.794662 87.76066 **0.285917409 3.497514 0.402143 1.71814 88.19828** 0.463649852 2.1568 0.050809 1.694368 88.62985 0.162277448 6.162286 1.80194 1.617602 89.04186 0.177732443 5.626435 0.231643 1.60276 89.4501 0.216369931 4.621714 0.236801 1.570401 89.85009 **0.208642433 4.792889 0.21626 1.560695 90.24761**

*Kalman Filter Harmonic Bank for Vostok Ice Core Data Analysis and Climate Predictions*

*DOI: http://dx.doi.org/10.5772/intechopen.94263*

The cumulative column in **Table 7** indicates amplitude sum percentage of subsequent harmonics compared to the sum of all harmonic amplitudes. For example to reach 90% + of the total amplitude we need a total of 25 harmonics, including H0.


*Kalman Filter Harmonic Bank for Vostok Ice Core Data Analysis and Climate Predictions DOI: http://dx.doi.org/10.5772/intechopen.94263*

**Table 7.**

**6. Cycle C1 CO2 frequency analysis**

*Glaciers and the Polar Environment*

and processing the data through CPE.

the current time, **Figure 9** above.

**6.2 Sampling time**

i.e. *<sup>T</sup>* <sup>¼</sup> 128, 399

<sup>2</sup> <sup>¼</sup> <sup>1</sup>

<sup>2</sup>*<sup>T</sup>* <sup>¼</sup> <sup>1</sup>

and optimal Fourier analysis.

**6.3 Amplitude and energy analysis**

*<sup>f</sup> max* <sup>¼</sup> *<sup>f</sup>*

**48**

**6.1 Data filling**

the sampling time can be more precise and useful.

As an illustration of frequency analysis we examine Cycle C1 CO2 data in details. Similar approach can be used for other data, such as temperature, methane, oxygen and dust. The first concern is to make the C1 CO2 data more homogeneous hence

The original Cycle 1 Vostok data contains a total of 80 CO2 data points are very unevenly distributed in time. This poses issues with FFT or DFT harmonic analysis. Hence we corrected the situation somewhat by inserting certain number of additional "data points" using linear approximation between the original data points with the biggest time difference between the neighboring points, hence generating a total of 128 original and inserted data points. In particular, the very first point we added to the original Vostok data is the current value of CO2 which we stated at 350 [25] for the year 1990, as an example. Current levels are estimated at over 400 ppm. Our CPE model can be used for a variety of sensitivity analysis by playing "what if"

Next data is the first original Vostok data, 2362 years in the past. This will correspond to 64 FFT harmonics described next. Note that the time period of C1 CO2 data is 128,399 years (**Table 2**), time difference between the oldest C1 data and

For sampling time, we proceeded as follows. As stated in Section 2 the cycle duration is measured between the maximum (CO2 and temperature) points in each cycle. Hence, we have an average of a bit over 1000 years time sampling for Cycle 1,

analysis on 128 total points produces 64 harmonics with the highest frequency of

with the understanding it is in 1000s of years. **Table 7** summarizes Cycle 1 harmonic content ordered by the amplitude, from the largest to the smallest to identify each harmonic energy contribution. First harmonic H1 is boldfaced. As shown in Section 7.1, for the discretization purposes, it is important to satisfy an additional sampling time requirement given in inequality (7) bellow which produces equivalent resonant frequency of continuous and discrete models, and in turn makes for more precise discrete KFHO and KFHB models as well. This increases required sampling frequency and reduces time sampling period bellow 1011 years, and we deal with it using energy analysis as described in **Table 8**. For other approaches to non-uniform data analysis see [26, 27] as well as general literature on non uniform

**Table 7** shows that the average signal value (DC or H0) is 231 and the rest of the energy content (reflected in a corresponding harmonic amplitude) of the harmonics follows the highest-to-lowest amplitude pattern H1, H2, H4, H3, H7, H10, H5, and so on. The cumulative column in **Table 7** indicates amplitude sum percentage of subsequent harmonics compared to the sum of all harmonic amplitudes. For example to reach 90% + of the total amplitude we need a total of 25 harmonics, including H0.

<sup>128</sup>�<sup>1</sup> = 1011 years, which is used as the CO2 data sampling time. The FFT

<sup>2022</sup> <sup>10</sup>�<sup>3</sup> <sup>¼</sup> <sup>0</sup>*:*<sup>000495</sup> *<sup>x</sup>* <sup>10</sup>�<sup>3</sup> Hz. For the analysis we use *<sup>T</sup>* <sup>¼</sup> <sup>1</sup>*:*<sup>011</sup>

*Summary of cycle 1 CO2 harmonic analysis.*


## **Table 8.**

*Cycle 1 CO2 amplitude and sampling time.*

The last (boldfaced) harmonic in this list is H27. For a bit smaller cumulative amplitude, say 88%+, we only need 19 harmonics, with H37 as the last one. **Table 8** further summarizes energy analysis. For example, a chosen sampling time of 1011 satisfies the requirement in inequality (7) bellow all the way to cumulative energy


level of continuous-to-discrete approximation is based on one point derivative

*Kalman Filter Harmonic Bank for Vostok Ice Core Data Analysis and Climate Predictions*

<sup>¼</sup> <sup>1</sup> *<sup>T</sup>* �*ω*<sup>2</sup> <sup>0</sup>*T* 1 *x*1ð Þ*t*

with "*t* þ 1" standing for ð Þ *t* þ 1 *T*, and similarly ð Þ *t* þ 2 *T* for "ð Þ *t* þ 2 ", where *T* is dropped for simplicity. To the above discretized time model we can also add model uncertainty via additional stochastic zero mean Gaussian white inputs *r*<sup>1</sup> and *r*2, with certain variance values, V1 and V2 which can be fine-tuned. Hence

> <sup>¼</sup> <sup>1</sup> *<sup>T</sup>* �*ω*<sup>2</sup> <sup>0</sup>*T* 1 *x*1ð Þ*t*

The initial conditions are given as a transposed vector ð Þ *x*1ð Þj 0 *x*2ð Þ 0 with:

where *A*<sup>0</sup> and *θ*<sup>0</sup> are the amplitude and the phase of the harmonic *ω*. Better discrete approximation can be obtained by 2 point derivative approximation

where parameter *a*<sup>0</sup> is calculated to match discrete and continuous resonant

It is important to note that in this case one needs to choose sampling time

which is higher than the standard Nyquist frequency, *f* > 2 *f* <sup>0</sup>. We show in the next Section an example of this. Eq. (5) produces cosine function with the proper initial conditions *x*1ð Þ 0 in (4). The corresponding equation for *x*2ð Þ *t* þ 1 is equivalent to (5) with the initial condition *x*2ð Þ 0 in (4) to produce sin function. Instead of proceeding with two state model (3) we can produce another two state model

1 0 *<sup>x</sup>*1ð Þ*<sup>t</sup>*

where *w*1ð Þ*t* and *w*2ð Þ*t* are zero mean Gaussian white inputs with joint covariance

measurement error *v*0ð Þ*t* is simply defined in (9) and it corresponds to the harmonic *ω*0. Measurement error *v*0ð Þ*t* is zero mean stochastic process with variance *R*<sup>0</sup> assumed constant across all time samples, which is a reasonable assumption, unless

*x*1ð Þ *t* � 1 

, [2]. The measurement *<sup>y</sup>*0ð Þ*<sup>t</sup>* of *<sup>x</sup>*1ð Þ*<sup>t</sup>* with the

<sup>¼</sup> �*a*<sup>0</sup> �<sup>1</sup>

*Q*<sup>12</sup> *Q*<sup>22</sup>

*x*2ð Þ*t*

þ

*r*1ð Þ*t r*2ð Þ*t*

(3)

*x*2ð Þ*t* 

*x*1ð Þ¼ 0 *A*<sup>0</sup> *cos*ð Þ *θ*<sup>0</sup> , *x*2ð Þ¼� 0 *A*<sup>0</sup> *sin* ð Þ *θ*<sup>0</sup> (4)

*x*1ð Þ¼� *t* þ 1 *a*0*x*<sup>1</sup> ð Þ� t *x*1ð Þ *t* � 1 (5)

*a*<sup>0</sup> ¼ �2 cosð Þ *ω*0*T* (6)

*f* >*π f* <sup>0</sup>,*ω*<sup>0</sup> ¼ 2*π f* <sup>0</sup> (7)

þ

*y*0ðÞ¼ *t x*1ðÞþ*t v*0ð Þ*t* (9)

*w*1ð Þ*t w*2ð Þ*t*

(8)

(2)

approximation and sampling time T, whereas we obtain:

*DOI: http://dx.doi.org/10.5772/intechopen.94263*

*x*1ð Þ *t* þ 1 *x*2ð Þ *t* þ 2 

*x*1ð Þ *t* þ 1 *x*2ð Þ *t* þ 1 

we have:

whereas we obtain:

frequencies:

using (5):

**51**

*x*1ð Þ *t* þ 1 *x*1ð Þ*t* 

symmetric matrix Q = *<sup>Q</sup>*<sup>11</sup> *<sup>Q</sup>*<sup>12</sup>

*T* ¼ 1*=f* to satisfy:

**Table 9.**

*Kalman filter harmonic bank estimation errors.*

of 88.1983%. To reach the energy level of 90.2476%, much lower sampling time would be required, i.e. 650 years instead of 1011. In order to satisfy this requirement we would need to almost double the total number of data points and use more harmonics. Fortunately. not all the amplitudes are at their maximums all the time, hence the above percentages are only very conservative lower bounds. In reality we can calculate the errors by MAPE (Maximum Absolute Percentage Error) and these will be much more realistic, actually of order of 3–5% off of 100%, per **Table 9**. Also, energy of these first seven harmonics is of order of 99% + due to its calculation based on amplitude squares. This analysis gives us an idea of various issues at hand related to non-uniformity of Vostok data as well as number of harmonics to be used based on energy analysis. Our aim is to show that properly designed KFHB based CPE can deal with these issues in a very effective way. In Section 7 we illustrate various points raised here using first 7 harmonics, H1, H2, H4, H3, H7, H10, H5.

## **7. Kalman filter harmonic bank**

In an illustrative example in this paper we focus on C1 CO2 harmonic analysis and the corresponding KFHO for the strongest, amplitude wise, first harmonic H1 in **Table 7**. The same analysis can be repeated for any harmonic as we elaborate in Sections bellow. We proceed with the development of KFHB which models a series of harmonic components obtained by the energy analysis, which approximates to a reasonable degree (say 90% + or higher cumulative amplitude percentage) the original Vostok data. First step is to define a general Markov model [28] for a generic harmonic oscillator in discrete time described in the following Section. Other authors also used the Kalman filter approach to analyze ice core data, but with a different emphasis, [21, 22] as mentioned in Introduction.

#### **7.1 Discrete time harmonic oscillator**

To start we introduce a continuous model of a Harmonic Oscillator described in [29, 30]:

$$
\begin{pmatrix}
\dot{\boldsymbol{x}}\_1 \\
\dot{\boldsymbol{x}}\_2
\end{pmatrix} = \begin{pmatrix}
\mathbf{0} & \mathbf{1} \\
\end{pmatrix} \begin{pmatrix}
\boldsymbol{x}\_1 \\
\boldsymbol{x}\_2
\end{pmatrix} \tag{1}
$$

where *ω*<sup>0</sup> is radial frequency, *ω*<sup>0</sup> ¼ 2*π f* <sup>0</sup>, with *f* <sup>0</sup> frequency in Hz. The solution to the above state equation is sine or cosine function depending on the initial conditions. The state variables *x*<sup>1</sup> and *x*<sup>2</sup> are "position" and "velocity" of whatever variable we are dealing with, such as CO2 content and its rate of change. The first

*Kalman Filter Harmonic Bank for Vostok Ice Core Data Analysis and Climate Predictions DOI: http://dx.doi.org/10.5772/intechopen.94263*

level of continuous-to-discrete approximation is based on one point derivative approximation and sampling time T, whereas we obtain:

$$
\begin{pmatrix}
\varkappa\_1(t+1) \\
\varkappa\_2(t+2)
\end{pmatrix} = \begin{pmatrix}
1 & T \\
\end{pmatrix} \begin{pmatrix}
\varkappa\_1(t) \\
\varkappa\_2(t)
\end{pmatrix} \tag{2}
$$

with "*t* þ 1" standing for ð Þ *t* þ 1 *T*, and similarly ð Þ *t* þ 2 *T* for "ð Þ *t* þ 2 ", where *T* is dropped for simplicity. To the above discretized time model we can also add model uncertainty via additional stochastic zero mean Gaussian white inputs *r*<sup>1</sup> and *r*2, with certain variance values, V1 and V2 which can be fine-tuned. Hence we have:

$$
\begin{pmatrix}
\varkappa\_1(t+1) \\
\varkappa\_2(t+1)
\end{pmatrix} = \begin{pmatrix}
1 & T \\
\end{pmatrix} \begin{pmatrix}
\varkappa\_1(t) \\
\varkappa\_2(t)
\end{pmatrix} + \begin{pmatrix}
r\_1(t) \\
r\_2(t)
\end{pmatrix} \tag{3}
$$

The initial conditions are given as a transposed vector ð Þ *x*1ð Þj 0 *x*2ð Þ 0 with:

$$\mathbf{x}\_1(\mathbf{0}) = A\_0 \cos \left(\theta\_0\right), \mathbf{x}\_2(\mathbf{0}) = -A\_0 \sin \left(\theta\_0\right) \tag{4}$$

where *A*<sup>0</sup> and *θ*<sup>0</sup> are the amplitude and the phase of the harmonic *ω*. Better discrete approximation can be obtained by 2 point derivative approximation whereas we obtain:

$$\boldsymbol{\infty}\_{1}(t+1) = -\boldsymbol{a}\_{0}\boldsymbol{\infty}\_{1}\left(\mathbf{t}\right) - \boldsymbol{\infty}\_{1}(t-1) \tag{5}$$

where parameter *a*<sup>0</sup> is calculated to match discrete and continuous resonant frequencies:

$$a\_0 = -2\cos\left(a\_0 T\right) \tag{6}$$

It is important to note that in this case one needs to choose sampling time *T* ¼ 1*=f* to satisfy:

$$f > \pi f\_{\; 0}, \; \alpha\_0 = 2\pi f\_{\; 0} \tag{7}$$

which is higher than the standard Nyquist frequency, *f* > 2 *f* <sup>0</sup>. We show in the next Section an example of this. Eq. (5) produces cosine function with the proper initial conditions *x*1ð Þ 0 in (4). The corresponding equation for *x*2ð Þ *t* þ 1 is equivalent to (5) with the initial condition *x*2ð Þ 0 in (4) to produce sin function. Instead of proceeding with two state model (3) we can produce another two state model using (5):

$$
\begin{pmatrix} \varkappa\_1(t+1) \\ \varkappa\_1(t) \end{pmatrix} = \begin{pmatrix} -a\_0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \varkappa\_1(t) \\ \varkappa\_1(t-1) \end{pmatrix} + \begin{pmatrix} w\_1(t) \\ w\_2(t) \end{pmatrix} \tag{8}
$$

$$\mathcal{Y}\_0(t) = \varkappa\_1(t) + \nu\_0(t) \tag{9}$$

where *w*1ð Þ*t* and *w*2ð Þ*t* are zero mean Gaussian white inputs with joint covariance symmetric matrix Q = *<sup>Q</sup>*<sup>11</sup> *<sup>Q</sup>*<sup>12</sup> *Q*<sup>12</sup> *Q*<sup>22</sup> , [2]. The measurement *<sup>y</sup>*0ð Þ*<sup>t</sup>* of *<sup>x</sup>*1ð Þ*<sup>t</sup>* with the measurement error *v*0ð Þ*t* is simply defined in (9) and it corresponds to the harmonic *ω*0. Measurement error *v*0ð Þ*t* is zero mean stochastic process with variance *R*<sup>0</sup> assumed constant across all time samples, which is a reasonable assumption, unless

of 88.1983%. To reach the energy level of 90.2476%, much lower sampling time would be required, i.e. 650 years instead of 1011. In order to satisfy this requirement we would need to almost double the total number of data points and use more harmonics. Fortunately. not all the amplitudes are at their maximums all the time, hence the above percentages are only very conservative lower bounds. In reality we can calculate the errors by MAPE (Maximum Absolute Percentage Error) and these will be much more realistic, actually of order of 3–5% off of 100%, per **Table 9**. Also, energy of these first seven harmonics is of order of 99% + due to its calculation based on amplitude squares. This analysis gives us an idea of various issues at hand related to non-uniformity of Vostok data as well as number of harmonics to be used based on energy analysis. Our aim is to show that properly designed KFHB based CPE can deal with these issues in a very effective way. In Section 7 we illustrate various points raised here using first 7 harmonics, H1, H2, H4, H3, H7, H10, H5.

y 0.058 0.043 0.038 0.036 0.034 0.033 0.031

0.058 0.043 0.037 0.036 0.034 0.033 0.030

0.056 0.044 0.039 0.037 0.035 0.034 0.033

**1 1 + 2 1 + 2 + 4 1 + 2 + 4 + 3 1 + 2 + 4 + 3 + 7 1 + 2 + 4 + 3 + 7 + 10 1 + 2 + 4 + 3 + 7 + 10 + 5**

In an illustrative example in this paper we focus on C1 CO2 harmonic analysis and the corresponding KFHO for the strongest, amplitude wise, first harmonic H1 in **Table 7**. The same analysis can be repeated for any harmonic as we elaborate in Sections bellow. We proceed with the development of KFHB which models a series of harmonic components obtained by the energy analysis, which approximates to a reasonable degree (say 90% + or higher cumulative amplitude percentage) the original Vostok data. First step is to define a general Markov model [28] for a generic harmonic oscillator in discrete time described in the following Section. Other authors also used the Kalman filter approach to analyze ice core data, but

To start we introduce a continuous model of a Harmonic Oscillator described in

where *ω*<sup>0</sup> is radial frequency, *ω*<sup>0</sup> ¼ 2*π f* <sup>0</sup>, with *f* <sup>0</sup> frequency in Hz. The solution

*x*2 

(1)

<sup>¼</sup> 0 1 �*ω*<sup>2</sup> <sup>0</sup> 0 *x*<sup>1</sup>

to the above state equation is sine or cosine function depending on the initial conditions. The state variables *x*<sup>1</sup> and *x*<sup>2</sup> are "position" and "velocity" of whatever variable we are dealing with, such as CO2 content and its rate of change. The first

with a different emphasis, [21, 22] as mentioned in Introduction.

*x*\_ 1 *x*\_ 2 

**7. Kalman filter harmonic bank**

**MAPE Harmonics**

*Glaciers and the Polar Environment*

*Kalman filter harmonic bank estimation errors.*

x Corrected

x Predicted

**Table 9.**

**7.1 Discrete time harmonic oscillator**

[29, 30]:

**50**

there is a compelling reason to make it time varying. The model (8) and (9) above remains the same. Obviously equivalent model holds for *x*2ð Þ *t* þ 1 with the proper initial conditions. The model as given by (8) and (9) is our starting point for KFHO described next. The consideration holds for any harmonic *ω*.
