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