Coverage for src / thunderlab / fourier.py: 97%

1"""Fourier series 

2 

3Extract and normalize Fourier coefficients of a Fourier series of 

4periodic functions, or synthesize periodic functions from Fourier 

5coefficients. 

6 

7## Functions 

8 

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. 

12 

13""" 

14 

15import numpy as np 

16 

17 

18def fourier_coeffs(data, ratetime, freq, n_harmonics): 

19 """ Extract Fourier coefficients from data. 

20 

21 Decompose a periodic 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_k = \\frac{2}{jT} \\int_{0}^{jT} x(t) e^{-i 2 \\pi k f_1 t} \\, dt \\] 

26 integrated over integer multiples of the period \\(T=1/f_1\\). 

27 

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 n_harmonics: int 

38 Number of harmonics inclusively the zeroth one. 

39 That is, if `n_harmonics` is set to two, 

40 then the zeroth harmonics (offset) and the fundamental 

41 frequency (first harmonics) are returned. 

42 

43 Returns 

44 ------- 

45 coeffs: 1D array of complex 

46 For each harmonics the complex valued Fourier coefficient. 

47 The first one is the offset and the second one is the coefficient 

48 of the fundamental. If the number ofdata samples is less than 

49 a single period, then a zero-sized array is returned. 

50 """ 

51 if isinstance(ratetime, (list, tuple, np.ndarray)): 

52 time = ratetime 

53 else: 

54 time = np.arange(len(data))/ratetime 

55 # integrate over full periods: 

56 n_periods = int(np.floor((time[-1] - time[0])*freq)) 

57 if n_periods < 1: 

58 return np.zeros(0, dtype=complex) 

59 n_max = np.argmax(time > time[0] + n_periods/freq) 

60 data = data[:n_max] 

61 time = time[:n_max] 

62 # Fourier projections: 

63 iomega = -2j*np.pi*freq*time 

64 fac = 2/len(data) # = 2*deltat/T 

65 coeffs = np.zeros(n_harmonics, dtype=complex) 

66 for k in range(n_harmonics): 

67 coeffs[k] = np.sum(data*np.exp(iomega*k))*fac 

68 return coeffs 

69 

70 

71def normalize_fourier_coeffs(coeffs): 

72 """Set phase of first harmonics and offset to zero. 

73 

74 Parameters 

75 ---------- 

76 coeffs: 1D array of complex 

77 For each harmonic the complex valued Fourier coefficient 

78 as, for example, returned by `fourier_coeffs()`. 

79 The first one is the offset and the second one is the coefficient 

80 of the fundamental. 

81 

82 Returns 

83 ------- 

84 coeffs: 1D array of complex 

85 The normalized Fourier coefficients. 

86 """ 

87 phi1 = np.angle(coeffs[1]) 

88 for k in range(1, len(coeffs)): 

89 coeffs[k] *= np.exp(-1j*k*phi1) 

90 coeffs[0] = 0 + 0j 

91 return coeffs 

92 

93 

94def fourier_synthesis(freq, coeffs, ratetime, n=None): 

95 """ Compute periodic waveform from Fourier coefficients. 

96 

97 Given the Fourier coefficients 

98 \\[ c_k = \\frac{2}{jT} \\int_{0}^{jT} x(t) e^{-i 2 \\pi k f_1 t} \\, dt \\] 

99 integrated over integer multiples of the period \\(T=1/f_1\\) of a signal 

100 \\(x(t)\\) with fundamental frequency \\(f_1\\), compute 

101 the Fourier series 

102 \\[ x(t) \\approx \\Re \\sum_{k=0}^n c_k e^{i 2 \\pi k f_1 t} \\] 

103 

104 Parameters 

105 ---------- 

106 freq: float 

107 Fundamental frequency. 

108 coeffs: 1D array of complex 

109 For each harmonics the complex valued Fourier coefficient 

110 as, for example, returned by `fourier_coeffs()`. 

111 The first one is the offset. 

112 ratetime: float or 1-D array of float 

113 Time points for which the waveform is calculated. 

114 If single float, then sampling rate of the computed waveform. 

115 n: int 

116 Number of samples if `ratetime` is float. 

117 

118 Returns 

119 ------- 

120 wave: 1D array of float 

121 Waveform computed from Fourier series with fundamental frequency 

122 `freq` and Fourier coefficients `coeffs` for each harmonic. 

123 The waveform is computed for a sampling rate `rate` and contains 

124 `n` samples. 

125 """ 

126 if isinstance(ratetime, (list, tuple, np.ndarray)): 

127 time = ratetime 

128 else: 

129 time = np.arange(n)/ratetime 

130 iomega = 2j*np.pi*freq*time 

131 wave = np.zeros(len(time)) 

132 for k in range(len(coeffs)): 

133 wave += np.real(coeffs[k]*np.exp(iomega*k)) 

134 return wave 

135