Skip to content

Fourier

coef_conversion(amp_coef=None, phase_coef=None, cos_coef=None, sin_coef=None, pos_amp=False, neg_amp=False, take_cos_sign=False)

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, DataArray]]

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

 None
phase_coef Optional[Union[float, ndarray, DataArray]]

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

 None
cos_coef Optional[Union[float, ndarray, DataArray]]

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

 None
sin_coef Optional[Union[float, ndarray, DataArray]]

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

 None
pos_amp bool

If pos_amp=True: amp_coef >= 0, phase in (\([-\pi, \pi]\) (via atan2). If pos_amp=False: allow signed amp_coef and shift phase by +/-pi so that phase is constrained to \([-\pi/2, \pi/2]\) (minimize |phase|).

 False
neg_amp bool

Force amplitude to be negative.

 False
take_cos_sign bool

Will make amp_coef take sign of provided cos_coef.

 False

Returns:

Name Type Description
coef1 Union[float, ndarray, DataArray]

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

float [n_coefs]
Either \(\Phi_n\) or \(F_{n, sin}\) depending on input.
Source code in isca_tools/utils/fourier.py 
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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def coef_conversion(amp_coef: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                    phase_coef: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                    cos_coef: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                    sin_coef: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                    pos_amp: bool = False,
                    neg_amp: bool = False,
                    take_cos_sign: bool = False
                    ) -> Tuple[Union[float, np.ndarray, xr.DataArray], Union[float, np.ndarray, xr.DataArray]]:
    """
    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$.
        pos_amp: If pos_amp=True: amp_coef >= 0, phase in ($[-\pi, \pi]$ (via atan2).
            If pos_amp=False: allow signed amp_coef and shift phase by +/-pi so that
            phase is constrained to $[-\pi/2, \pi/2]$ (minimize |phase|).
        neg_amp: Force amplitude to be negative.
        take_cos_sign:
            Will make amp_coef take sign of provided cos_coef.

    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:
        # amp/phase -> cos/sin (sign convention is already encoded in amp/phase)
        cos_coef = amp_coef * np.cos(phase_coef)
        sin_coef = amp_coef * np.sin(phase_coef)
        return  cos_coef, sin_coef
    else:
        phase_coef = np.arctan2(sin_coef, cos_coef)     # robust quadrant-aware phase
        amp_coef = np.sqrt(sin_coef**2 + cos_coef**2)   # nonnegative amplitude
        if take_cos_sign:
            # make amplitude signed to match cos_coef
            sgn = np.where(cos_coef < 0, -1.0, 1.0)
            amp_coef = amp_coef * sgn

            # compensate phase where we flipped the sign
            phase_coef = phase_coef + (sgn < 0) * np.pi

            # optional: wrap to [-pi, pi)
            phase_coef = (phase_coef + np.pi) % (2 * np.pi) - np.pi
        elif neg_amp:
            amp_coef = -amp_coef  # force negative amplitude

            # keep the same modeled signal by shifting phase by pi
            phase_coef = phase_coef + np.pi

            # optional: wrap back to [-pi, pi)
            phase_coef = (phase_coef + np.pi) % (2 * np.pi) - np.pi
        elif not pos_amp:
            # Flip sign to keep phase in [-pi/2, pi/2]
            # (works for scalars, numpy arrays, and xarray DataArray via broadcasting)
            over = phase_coef > (0.5 * np.pi)
            under = phase_coef < (-0.5 * np.pi)

            phase_coef = phase_coef - np.pi * over + np.pi * under
            amp_coef = amp_coef * (1.0 - 2.0 * (over | under))
        return amp_coef, phase_coef

fourier_series(time, coefs_amp, coefs_phase, pad_coefs_phase=False)

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
pad_coefs_phase bool

If True, expect coefs_fourier_phase to be of length n_harmonics+1 with first value equal to zero.

 False

Returns:

Type Description
ndarray

float [n_time]
Value of Fourier series solution at each time
Source code in isca_tools/utils/fourier.py 
 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
33
34
35
def fourier_series(time: np.ndarray, coefs_amp: Union[List[float], np.ndarray],
                   coefs_phase: Union[List[float], np.ndarray], pad_coefs_phase: bool = False) -> 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$
        pad_coefs_phase: If `True`, expect `coefs_fourier_phase` to be of length `n_harmonics+1` with first value
            equal to zero.

    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 if pad_coefs_phase else 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 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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', pos_amp=False)

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'
pos_amp bool

If True, amplitude coefficient will always be positive with phase_coef in \([-\pi, \pi]\), Otherwise will choose the sign of amp_coef to minimize the magnitude of phase_coef i.e., keep phase_coef in range between \([-\pi/2, \pi/2]\).

 False

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 
 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
114
115
116
117
118
def get_fourier_coef(time: np.ndarray, var: np.ndarray, n: int,
                     integ_method: str = 'spline', pos_amp: bool = False) -> [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`.
        pos_amp: If `True`, amplitude coefficient will always be positive with `phase_coef` in $[-\pi, \pi]$,
            Otherwise will choose the sign of `amp_coef` to minimize the magnitude of `phase_coef`
            i.e., keep `phase_coef` in range between $[-\pi/2, \pi/2]$.
    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
        # use arctan2 not arctan as more robust (quadrant aware), always keeps amp positive
        return coef_conversion(cos_coef=cos_coef, sin_coef=sin_coef, pos_amp=pos_amp)

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

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'
pad_coefs_phase bool

If True, will set coefs_phase to length n_harmonics+1, with the first value as zero. Otherwise will be size n_harmonics.

 False

Returns:

Type Description
ndarray

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

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

coefs_phase: float [n_harmonics]
The phase Fourier coefficients \(\Phi_n\). If pad_coefs_phase, will pad at start with a zero so of length n_harmonics+1.
Source code in isca_tools/utils/fourier.py 
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
148
149
150
151
152
153
154
155
156
157
158
def get_fourier_fit(time: np.ndarray, var: np.ndarray, n_harmonics: int,
                    integ_method: str = 'spline',
                    pad_coefs_phase: bool = False) -> 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`.
        pad_coefs_phase: If `True`, will set `coefs_phase` to length `n_harmonics+1`, with
            the first value as zero. Otherwise will be size `n_harmonics`.

    Returns:
        `fourier_solution`: `float [n_time]`</br>
            The Fourier series solution that was fit to `var`.
        `coefs_amp`: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients $F_n$.
        `coefs_phase`: `float [n_harmonics]`</br>
            The phase Fourier coefficients $\\Phi_n$.
            If `pad_coefs_phase`, will pad at start with a zero so of length `n_harmonics+1`.
    """
    # 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)
    if pad_coefs_phase:
        phase_coefs = np.hstack((np.zeros(1), phase_coefs))
    return fourier_series(time, amp_coefs, phase_coefs, pad_coefs_phase), amp_coefs, phase_coefs