Coverage for src / thunderlab / fourier.py: 98%
52 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-12 18:00 +0000
1"""Fourier series
3Extract and normalize Fourier coefficients of a Fourier series of
4periodic functions, or synthesize periodic functions from Fourier
5coefficients.
7## Functions
9- `fourier_coeffs()`: extract Fourier coefficients from data.
10- `normalize_fourier_coeffs()`: set phase of first harmonics and offset to zero.
11- `fourier_synthesis()`: compute waveform from Fourier coefficients.
13"""
15import numpy as np
18def fourier_coeffs(data, ratetime, freq, max_harmonics):
19 """ Extract Fourier coefficients from data.
21 Decompose a periodic, real-valued signal \\(x(t)\\) with known fundamental frequency
22 \\(f_1\\) into a Fourier series:
23 \\[ x(t) \\approx \\Re \\sum_{k=0}^n c_k e^{i 2 \\pi k f_1 t} \\]
24 with the Fourier coefficients
25 \\[ c_0 = \\frac{1}{jT} \\int_{0}^{jT} x(t) \\, dt , \\quad k = 0 \\]
26 and
27 \\[ c_k = \\frac{2}{jT} \\int_{0}^{jT} x(t) e^{-i 2 \\pi k f_1 t} \\, dt , \\quad k > 0 \\]
28 integrated over integer multiples of the period \\(T=1/f_1\\).
30 Parameters
31 ----------
32 data: 1D array of float
33 Time series of data.
34 ratetime: float or 1-D array of float
35 Times corresponding to `data`.
36 If single float, then sampling rate of the data.
37 freq: float
38 Fundamental frequency of Fourier series.
39 max_harmonics: int
40 The highest harmonics for which to compute the Fourier coefficient.
41 The number of coefficients returned is one more than `max_harmonics`,
42 because the first coefficient is the zeroth harmonics.
43 For example, if `max_harmonics` is set to three,
44 then the zeroth harmonics (offset), the fundamental
45 frequency (first harmonics), and the second harmonics are returned.
47 Returns
48 -------
49 coeffs: 1D array of complex
50 For each harmonics the complex valued Fourier coefficient.
51 The first one is the offset and the second one is the coefficient
52 of the fundamental. If the number of data samples is less than
53 a single period, then a zero-sized array is returned.
54 """
55 if isinstance(ratetime, (list, tuple, np.ndarray)):
56 time = ratetime
57 else:
58 time = np.arange(len(data))/ratetime
59 # integrate over full periods:
60 n_periods = int(np.floor((time[-1] - time[0])*freq))
61 if n_periods < 1 or max_harmonics < 0:
62 return np.zeros(0, dtype=complex)
63 n_max = np.argmax(time >= time[0] + n_periods/freq)
64 data = data[:n_max]
65 time = time[:n_max]
66 # Fourier projections:
67 iomega = -2j*np.pi*freq*time
68 fac = 2/len(data) # = 2*deltat/T
69 coeffs = np.zeros(max_harmonics + 1, dtype=complex)
70 coeffs[0] = np.sum(data)*fac/2
71 for k in range(1, max_harmonics + 1):
72 coeffs[k] = np.sum(data*np.exp(iomega*k))*fac
73 return coeffs
76def normalize_fourier_coeffs(coeffs):
77 """Set phase of first harmonics and offset to zero.
79 Parameters
80 ----------
81 coeffs: 1D array of complex
82 For each harmonic the complex valued Fourier coefficient
83 as, for example, returned by `fourier_coeffs()`.
84 The first one is the offset and the second one is the coefficient
85 of the fundamental.
87 Returns
88 -------
89 coeffs: 1D array of complex
90 The normalized Fourier coefficients.
91 """
92 phi1 = np.angle(coeffs[1])
93 for k in range(1, len(coeffs)):
94 coeffs[k] *= np.exp(-1j*k*phi1)
95 coeffs[0] = 0 + 0j
96 return coeffs
99def fourier_synthesis(freq, coeffs, ratetime, n=None):
100 """ Compute periodic waveform from Fourier coefficients.
102 Given the Fourier coefficients
103 \\[ c_0 = \\frac{1}{jT} \\int_{0}^{jT} x(t) \\, dt , \\quad k = 0 \\]
104 and
105 \\[ c_k = \\frac{2}{jT} \\int_{0}^{jT} x(t) e^{-i 2 \\pi k f_1 t} \\, dt , \\quad k > 0 \\]
106 integrated over integer multiples of the period \\(T=1/f_1\\) of a
107 periodic, real-valued signal \\(x(t)\\) with fundamental frequency \\(f_1\\),
108 compute the Fourier series
109 \\[ x(t) \\approx \\Re \\sum_{k=0}^n c_k e^{i 2 \\pi k f_1 t} \\]
111 Parameters
112 ----------
113 freq: float
114 Fundamental frequency.
115 coeffs: 1D array of complex
116 For each harmonics the complex valued Fourier coefficient
117 as, for example, returned by `fourier_coeffs()`.
118 The first one is the offset.
119 ratetime: float or 1-D array of float
120 Time points for which the waveform is calculated.
121 If single float, then sampling rate of the computed waveform.
122 n: int
123 Number of samples if `ratetime` is float.
125 Returns
126 -------
127 wave: 1D array of float
128 Real-valued waveform computed from Fourier series with fundamental frequency
129 `freq` and Fourier coefficients `coeffs` for each harmonic.
130 The waveform is computed for a sampling rate `rate` and contains
131 `n` samples.
132 """
133 if isinstance(ratetime, (list, tuple, np.ndarray)):
134 time = ratetime
135 else:
136 time = np.arange(n)/ratetime
137 iomega = 2j*np.pi*freq*time
138 wave = np.zeros(len(time))
139 if len(coeffs) == 0:
140 return wave
141 for k, c in enumerate(coeffs):
142 wave += np.real(coeffs[k]*np.exp(iomega*k))
143 return wave
146def main():
147 """Demonstrate the fourier module.
148 """
149 import matplotlib.pyplot as plt
151 f = 1/1.5
152 t = np.linspace(0, 3.0, 200)
153 x = 1.7 + 2.3*np.cos(2*np.pi*f*t) + 0.8*np.cos(2*np.pi*2*f*t + 0.5*np.pi) + 0.4*np.cos(2*np.pi*3*f*t - 1.6*np.pi)
154 coeffs = fourier_coeffs(x, t, f, 6)
155 y = fourier_synthesis(f, coeffs, t)
156 fig, (axt, axc) = plt.subplots(1, 2, layout='constrained')
157 axt.plot(t, x, label='original', lw=6)
158 axt.plot(t, y, label='synthesized', lw=2)
159 axt.legend()
160 axc.plot(np.real(coeffs), '-o', label='real')
161 axc.plot(np.imag(coeffs), '-o', label='imag')
162 axc.legend()
163 plt.show()
166if __name__ == '__main__':
167 main()