Skip to content

Fourier

coef_conversion(amp_coef=None, phase_coef=None, cos_coef=None, sin_coef=None)

The term for the \(n^{th}\) harmonic of a Fourier expansion can be written in two ways:

  1. \(F_n\cos(2n\pi ft - \Phi_n)\)
  2. \(F_{n, cos}\cos(2n\pi ft) + F_{n, sin}\sin(2n\pi ft)\)

Given the coefficients of one form, this returns the coefficients in the other form. Note \(\sin(x) = \cos(x - \pi /2)\).

Parameters:

Name Type Description Default
amp_coef Optional[Union[float, ndarray]]

float [n_coefs]
\(F_n\) coefficients. If provided, will return \(F_{n, cos}\) and \(F_{n, sin}\).

 None
phase_coef Optional[Union[float, ndarray]]

float [n_coefs]
\(\Phi_n\) coefficients. If provided, will return \(F_{n, cos}\) and \(F_{n, sin}\).

 None
cos_coef Optional[Union[float, ndarray]]

float [n_coefs]
\(F_{n, cos}\) coefficients. If provided, will return \(F_{n}\) and \(\Phi_n\).

 None
sin_coef Optional[Union[float, ndarray]]

float [n_coefs]
\(F_{n, sin}\) coefficients. If provided, will return \(F_{n}\) and \(\Phi_n\).

 None

Returns:

Name Type Description
coef1 Union[float, ndarray]

float [n_coefs]
Either \(F_n\) or \(F_{n, cos}\) depending on input.
coef2 Union[float, ndarray]

float [n_coefs]
Either \(\Phi_n\) or \(F_{n, sin}\) depending on input.
Source code in isca_tools/utils/fourier.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def coef_conversion(amp_coef: Optional[Union[float, np.ndarray]] = None,
                    phase_coef: Optional[Union[float, np.ndarray]] = None,
                    cos_coef: Optional[Union[float, np.ndarray]] = None,
                    sin_coef: Optional[Union[float, np.ndarray]] = None) -> Tuple[Union[float, np.ndarray],
Union[float, np.ndarray]]:
    """
    The term for the $n^{th}$ harmonic of a Fourier expansion can be written in two ways:

    1. $F_n\\cos(2n\\pi ft - \\Phi_n)$
    2. $F_{n, cos}\\cos(2n\\pi ft) + F_{n, sin}\\sin(2n\\pi ft)$

    Given the coefficients of one form, this returns the coefficients in the other form.
    Note $\\sin(x) = \\cos(x - \\pi /2)$.

    Args:
        amp_coef: `float [n_coefs]`</br>
            $F_n$ coefficients. If provided, will return $F_{n, cos}$ and $F_{n, sin}$.
        phase_coef: `float [n_coefs]`</br>
            $\Phi_n$ coefficients. If provided, will return $F_{n, cos}$ and $F_{n, sin}$.
        cos_coef: `float [n_coefs]`</br>
            $F_{n, cos}$ coefficients. If provided, will return $F_{n}$ and $\Phi_n$.
        sin_coef: `float [n_coefs]`</br>
            $F_{n, sin}$ coefficients. If provided, will return $F_{n}$ and $\Phi_n$.

    Returns:
        coef1: `float [n_coefs]`</br>
            Either $F_n$ or $F_{n, cos}$ depending on input.
        coef2: `float [n_coefs]`</br>
            Either $\Phi_n$ or $F_{n, sin}$ depending on input.
    """
    if amp_coef is not None:
        cos_coef = amp_coef * np.cos(phase_coef)
        sin_coef = amp_coef * np.sin(phase_coef)
        return  cos_coef, sin_coef
    else:
        if cos_coef == 0:
            phase_coef = np.pi/2
            amp_coef = sin_coef
        else:
            phase_coef = np.arctan(sin_coef/cos_coef)
            amp_coef = cos_coef / np.cos(phase_coef)
        return amp_coef, phase_coef

fourier_series(time, coefs_amp, coefs_phase)

For \(N\) harmonics, the fourier series with frequency \(f\) is:

\[F(t) \approx \frac{F_0}{2} + \sum_{n=1}^{N} F_n\cos(2n\pi ft - \Phi_n)\]

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)) covering entire period such that period is time[-1] - time[0] + 1.

 required
coefs_amp Union[List[float], ndarray]

float [N+1]
The amplitude coefficients, \(F_n\)

 required
coefs_phase Union[List[float], ndarray]

float [N]
The phase coefficients in radians, \(\Phi_n\)

 required

Returns:

Type Description
ndarray

float [n_time]
Value of Fourier series solution at each time
Source code in isca_tools/utils/fourier.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def fourier_series(time: np.ndarray, coefs_amp: Union[List[float], np.ndarray],
                   coefs_phase: Union[List[float], np.ndarray]) -> np.ndarray:
    """
    For $N$ harmonics, the fourier series with frequency $f$ is:

    $$F(t) \\approx \\frac{F_0}{2} + \\sum_{n=1}^{N} F_n\\cos(2n\\pi ft - \\Phi_n)$$

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`) covering entire period
            such that period is `time[-1] - time[0] + 1`.
        coefs_amp: `float [N+1]`</br>
            The amplitude coefficients, $F_n$
        coefs_phase: `float [N]`</br>
            The phase coefficients in radians, $\Phi_n$

    Returns:
        `float [n_time]`</br>
            Value of Fourier series solution at each time
    """
    period = float(time[-1] - time[0] + 1)
    n_harmonics = len(coefs_amp)
    ans = 0.5 * coefs_amp[0]
    for n in range(1, n_harmonics):
        ans += coefs_amp[n] * np.cos(2*n*np.pi*time/period - coefs_phase[n-1])
    return ans

fourier_series_deriv(time, coefs_amp, coefs_phase, day_seconds=86400)

For \(N\) harmonics, the derivative of a fourier series with frequency \(f\) is:

\[\frac{dF}{dt} \approx -\sum_{n=1}^{N} 2n\pi fF_n\sin(2n\pi ft - \Phi_n)\]

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)) covering entire period such that period is time[-1] - time[0] + 1.

 required
coefs_amp Union[ndarray, List[float]]

float [N+1]
The amplitude coefficients, \(F_n\). \(F_0\) needs to be provided even though it is not used.

 required
coefs_phase Union[ndarray, List[float]]

float [N]
The phase coefficients in radians, \(\Phi_n\)

 required
day_seconds float

Duration of a day in seconds.

 86400

Returns:

Type Description
ndarray

float [n_time]
Value of derivative to Fourier series solution at each time. Units is units of \(F\) divided by seconds.
Source code in isca_tools/utils/fourier.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def fourier_series_deriv(time: np.ndarray, coefs_amp: Union[np.ndarray, List[float]],
                         coefs_phase: Union[np.ndarray, List[float]], day_seconds: float = 86400) -> np.ndarray:
    """
    For $N$ harmonics, the derivative of a fourier series with frequency $f$ is:

    $$\\frac{dF}{dt} \\approx -\\sum_{n=1}^{N} 2n\\pi fF_n\\sin(2n\\pi ft - \\Phi_n)$$

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`) covering entire period
            such that period is `time[-1] - time[0] + 1`.
        coefs_amp: `float [N+1]`</br>
            The amplitude coefficients, $F_n$. $F_0$ needs to be provided even though it is not used.
        coefs_phase: `float [N]`</br>
            The phase coefficients in radians, $\Phi_n$
        day_seconds: Duration of a day in seconds.

    Returns:
        `float [n_time]`</br>
            Value of derivative to Fourier series solution at each time. Units is units of $F$ divided by seconds.
    """
    period = float(time[-1] - time[0] + 1)
    n_harmonics = len(coefs_amp)
    ans = np.zeros_like(time, dtype=float)
    for n in range(1, n_harmonics):
        ans -= coefs_amp[n] * np.sin(2*n*np.pi*time/period - coefs_phase[n-1]) * (2*n*np.pi/period)
    return ans / day_seconds     # convert units from per day to per second

get_fourier_coef(time, var, n, integ_method='spline')

This calculates the analytic solution for the amplitude and phase coefficients for the nth harmonic 🔗

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)) covering entire period such that period is time[-1] - time[0] + 1.

 required
var ndarray

float [n_time]
Variable to fit fourier series to.

 required
n int

Harmonic to find coefficients for, if 0, will just return amplitude coefficient. Otherwise, will return an amplitude and phase coefficient.

 required
integ_method str

How to perform the integration.
If spline, will fit a spline and then integrate the spline, otherwise will use scipy.integrate.simpson.

 'spline'

Returns: amp_coef: The amplitude fourier coefficient \(F_n\). phase_coef: The phase fourier coefficient \(\Phi_n\). Will not return if \(n=0\).

Source code in isca_tools/utils/fourier.py 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def get_fourier_coef(time: np.ndarray, var: np.ndarray, n: int,
                     integ_method: str = 'spline') -> [Union[float, Tuple[float, float]]]:
    """
    This calculates the analytic solution for the amplitude and phase coefficients for the `n`th harmonic
    [🔗](https://www.bragitoff.com/2021/05/fourier-series-coefficients-and-visualization-python-program/)

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`) covering entire period
            such that period is `time[-1] - time[0] + 1`.
        var: `float [n_time]`</br>
            Variable to fit fourier series to.
        n: Harmonic to find coefficients for, if 0, will just return amplitude coefficient.
            Otherwise, will return an amplitude and phase coefficient.
        integ_method: How to perform the integration.</br>
            If `spline`, will fit a spline and then integrate the spline, otherwise will use `scipy.integrate.simpson`.
    Returns:
        `amp_coef`: The amplitude fourier coefficient $F_n$.
        `phase_coef`: The phase fourier coefficient $\\Phi_n$. Will not return if $n=0$.
    """
    # Computes the analytical fourier coefficients for the n harmonic of a given function
    # With integrate method = spline works very well i.e. fit spline then use spline.integrate functionality
    # Otherwise, there are problems with the integration especially at the limits e.g. t=0 and t=T.
    period = float(time[-1] - time[0] + 1)
    if integ_method == 'spline':
        var = np.append(var, var[0])
        time = np.append(time, time[-1]+1)
    if n == 0:
        if integ_method == 'spline':
            spline = CubicSpline(time, var, bc_type='periodic')
            return 2/period * spline.integrate(0, period)
        else:
            return 2/period * scipy.integrate.simpson(var, time)
    else:
        # constants for acos(t) + bsin(t) form
        if integ_method == 'spline':
            spline = CubicSpline(time,var * np.cos(2*n*np.pi*time/period), bc_type='periodic')
            cos_coef = 2/period * spline.integrate(0, period)
            sin_curve = var * np.sin(2*n*np.pi*time/period)
            sin_curve[-1] = 0       # Need first and last value to be the same to be periodic spline
                                    # usually have last value equal 1e-10 so not equal
            spline = CubicSpline(time,sin_curve, bc_type='periodic')
            sin_coef = 2/period * spline.integrate(0, period)
        else:
            cos_coef = 2/period * scipy.integrate.simpson(var * np.cos(2*n*np.pi*time/period), time)
            sin_coef = 2/period * scipy.integrate.simpson(var * np.sin(2*n*np.pi*time/period), time)
        # constants for Acos(t-phi) form
        phase_coef = np.arctan(sin_coef/cos_coef)
        amp_coef = cos_coef / np.cos(phase_coef)
        return amp_coef, phase_coef

get_fourier_fit(time, var, n_harmonics, integ_method='spline')

Obtains the Fourier series solution for \(F=\)var, using \(N=\)n_harmonics:

\[F(t) \approx \frac{F_0}{2} + \sum_{n=1}^{N} F_n\cos(2n\pi ft - \Phi_n)\]

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)) covering entire period such that period (\(1/f\)) is time[-1] - time[0] + 1.

 required
var ndarray

float [n_time]
Variable to fit fourier series to.

 required
n_harmonics int

Number of harmonics to use to fit fourier series, \(N\).

 required
integ_method str

How to perform the integration when obtaining Fourier coefficients.
If spline, will fit a spline and then integrate the spline, otherwise will use scipy.integrate.simpson.

 'spline'

Returns:

Type Description
ndarray

fourier_solution: float [n_time]
The Fourier series solution that was fit to var.
ndarray

amp_coef: float [n_harmonics+1]
The amplitude Fourier coefficients \(F_n\).
ndarray

phase_coef: float [n_harmonics]
The phase Fourier coefficients \(\Phi_n\).
Source code in isca_tools/utils/fourier.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def get_fourier_fit(time: np.ndarray, var: np.ndarray, n_harmonics: int,
                    integ_method: str = 'spline') -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Obtains the Fourier series solution for $F=$`var`, using $N=$`n_harmonics`:

    $$F(t) \\approx \\frac{F_0}{2} + \\sum_{n=1}^{N} F_n\\cos(2n\\pi ft - \\Phi_n)$$

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`) covering entire period
            such that period ($1/f$) is `time[-1] - time[0] + 1`.
        var: `float [n_time]`</br>
            Variable to fit fourier series to.
        n_harmonics: Number of harmonics to use to fit fourier series, $N$.
        integ_method: How to perform the integration when obtaining Fourier coefficients.</br>
            If `spline`, will fit a spline and then integrate the spline, otherwise will use `scipy.integrate.simpson`.

    Returns:
        `fourier_solution`: `float [n_time]`</br>
            The Fourier series solution that was fit to `var`.
        `amp_coef`: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients $F_n$.
        `phase_coef`: `float [n_harmonics]`</br>
            The phase Fourier coefficients $\\Phi_n$.
    """
    # Returns the fourier fit of a function using a given number of harmonics
    amp_coefs = np.zeros(n_harmonics+1)
    phase_coefs = np.zeros(n_harmonics)
    amp_coefs[0] = get_fourier_coef(time, var, 0, integ_method)
    for i in range(1, n_harmonics+1):
        amp_coefs[i], phase_coefs[i-1] = get_fourier_coef(time, var, i, integ_method)
    return fourier_series(time, amp_coefs, phase_coefs), amp_coefs, phase_coefs