Coverage for bilby/core/utils/series.py: 100%

58 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2025-05-06 04:57 +0000

1import numpy as np 

2 

3_TOL = 14 

4 

5 

6def get_sampling_frequency(time_array): 

7 """ 

8 Calculate sampling frequency from a time series 

9 

10 Attributes 

11 ========== 

12 time_array: array_like 

13 Time array to get sampling_frequency from 

14 

15 Returns 

16 ======= 

17 Sampling frequency of the time series: float 

18 

19 Raises 

20 ====== 

21 ValueError: If the time series is not evenly sampled. 

22 

23 """ 

24 tol = 1e-10 

25 if np.ptp(np.diff(time_array)) > tol: 

26 raise ValueError("Your time series was not evenly sampled") 

27 else: 

28 return np.round(1. / (time_array[1] - time_array[0]), decimals=_TOL) 

29 

30 

31def get_sampling_frequency_and_duration_from_time_array(time_array): 

32 """ 

33 Calculate sampling frequency and duration from a time array 

34 

35 Attributes 

36 ========== 

37 time_array: array_like 

38 Time array to get sampling_frequency/duration from: array_like 

39 

40 Returns 

41 ======= 

42 sampling_frequency, duration: float, float 

43 

44 Raises 

45 ====== 

46 ValueError: If the time_array is not evenly sampled. 

47 

48 """ 

49 

50 sampling_frequency = get_sampling_frequency(time_array) 

51 duration = len(time_array) / sampling_frequency 

52 return sampling_frequency, duration 

53 

54 

55def get_sampling_frequency_and_duration_from_frequency_array(frequency_array): 

56 """ 

57 Calculate sampling frequency and duration from a frequency array 

58 

59 Attributes 

60 ========== 

61 frequency_array: array_like 

62 Frequency array to get sampling_frequency/duration from: array_like 

63 

64 Returns 

65 ======= 

66 sampling_frequency, duration: float, float 

67 

68 Raises 

69 ====== 

70 ValueError: If the frequency_array is not evenly sampled. 

71 

72 """ 

73 

74 tol = 1e-10 

75 if np.ptp(np.diff(frequency_array)) > tol: 

76 raise ValueError("Your frequency series was not evenly sampled") 

77 

78 number_of_frequencies = len(frequency_array) 

79 delta_freq = frequency_array[1] - frequency_array[0] 

80 duration = np.round(1 / delta_freq, decimals=_TOL) 

81 

82 sampling_frequency = np.round(2 * (number_of_frequencies - 1) / duration, decimals=14) 

83 return sampling_frequency, duration 

84 

85 

86def create_time_series(sampling_frequency, duration, starting_time=0.): 

87 """ 

88 

89 Parameters 

90 ========== 

91 sampling_frequency: float 

92 duration: float 

93 starting_time: float, optional 

94 

95 Returns 

96 ======= 

97 float: An equidistant time series given the parameters 

98 

99 """ 

100 _check_legal_sampling_frequency_and_duration(sampling_frequency, duration) 

101 number_of_samples = int(duration * sampling_frequency) 

102 return np.linspace(start=starting_time, 

103 stop=duration + starting_time - 1 / sampling_frequency, 

104 num=number_of_samples) 

105 

106 

107def create_frequency_series(sampling_frequency, duration): 

108 """ Create a frequency series with the correct length and spacing. 

109 

110 Parameters 

111 ========== 

112 sampling_frequency: float 

113 duration: float 

114 

115 Returns 

116 ======= 

117 array_like: frequency series 

118 

119 """ 

120 _check_legal_sampling_frequency_and_duration(sampling_frequency, duration) 

121 number_of_samples = int(np.round(duration * sampling_frequency)) 

122 number_of_frequencies = int(np.round(number_of_samples / 2) + 1) 

123 

124 return np.linspace(start=0, 

125 stop=sampling_frequency / 2, 

126 num=number_of_frequencies) 

127 

128 

129def _check_legal_sampling_frequency_and_duration(sampling_frequency, duration): 

130 """ By convention, sampling_frequency and duration have to multiply to an integer 

131 

132 This will check if the product of both parameters multiplies reasonably close 

133 to an integer. 

134 

135 Parameters 

136 ========== 

137 sampling_frequency: float 

138 duration: float 

139 

140 """ 

141 num = sampling_frequency * duration 

142 if np.abs(num - np.round(num)) > 10**(-_TOL): 

143 raise IllegalDurationAndSamplingFrequencyException( 

144 '\nYour sampling frequency and duration must multiply to a number' 

145 'up to (tol = {}) decimals close to an integer number. ' 

146 '\nBut sampling_frequency={} and duration={} multiply to {}'.format( 

147 _TOL, sampling_frequency, duration, 

148 sampling_frequency * duration 

149 ) 

150 ) 

151 

152 

153def create_white_noise(sampling_frequency, duration): 

154 """ Create white_noise which is then coloured by a given PSD 

155 

156 Parameters 

157 ========== 

158 sampling_frequency: float 

159 duration: float 

160 duration of the data 

161 

162 Returns 

163 ======= 

164 array_like: white noise 

165 array_like: frequency array 

166 """ 

167 from .random import rng 

168 

169 number_of_samples = duration * sampling_frequency 

170 number_of_samples = int(np.round(number_of_samples)) 

171 

172 frequencies = create_frequency_series(sampling_frequency, duration) 

173 

174 norm1 = 0.5 * duration**0.5 

175 re1, im1 = rng.normal(0, norm1, (2, len(frequencies))) 

176 white_noise = re1 + 1j * im1 

177 

178 # set DC and Nyquist = 0 

179 white_noise[0] = 0 

180 # no Nyquist frequency when N=odd 

181 if np.mod(number_of_samples, 2) == 0: 

182 white_noise[-1] = 0 

183 

184 # python: transpose for use with infft 

185 white_noise = np.transpose(white_noise) 

186 frequencies = np.transpose(frequencies) 

187 

188 return white_noise, frequencies 

189 

190 

191def nfft(time_domain_strain, sampling_frequency): 

192 """ Perform an FFT while keeping track of the frequency bins. Assumes input 

193 time series is real (positive frequencies only). 

194 

195 Parameters 

196 ========== 

197 time_domain_strain: array_like 

198 Time series of strain data. 

199 sampling_frequency: float 

200 Sampling frequency of the data. 

201 

202 Returns 

203 ======= 

204 frequency_domain_strain, frequency_array: (array_like, array_like) 

205 Single-sided FFT of time domain strain normalised to units of 

206 strain / Hz, and the associated frequency_array. 

207 

208 """ 

209 frequency_domain_strain = np.fft.rfft(time_domain_strain) 

210 frequency_domain_strain /= sampling_frequency 

211 

212 frequency_array = np.linspace( 

213 0, sampling_frequency / 2, len(frequency_domain_strain)) 

214 

215 return frequency_domain_strain, frequency_array 

216 

217 

218def infft(frequency_domain_strain, sampling_frequency): 

219 """ Inverse FFT for use in conjunction with nfft. 

220 

221 Parameters 

222 ========== 

223 frequency_domain_strain: array_like 

224 Single-sided, normalised FFT of the time-domain strain data (in units 

225 of strain / Hz). 

226 sampling_frequency: int, float 

227 Sampling frequency of the data. 

228 

229 Returns 

230 ======= 

231 time_domain_strain: array_like 

232 An array of the time domain strain 

233 """ 

234 time_domain_strain_norm = np.fft.irfft(frequency_domain_strain) 

235 time_domain_strain = time_domain_strain_norm * sampling_frequency 

236 return time_domain_strain 

237 

238 

239class IllegalDurationAndSamplingFrequencyException(Exception): 

240 pass