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
« prev ^ index » next coverage.py v7.6.1, created at 2025-05-06 04:57 +0000
1import numpy as np
3_TOL = 14
6def get_sampling_frequency(time_array):
7 """
8 Calculate sampling frequency from a time series
10 Attributes
11 ==========
12 time_array: array_like
13 Time array to get sampling_frequency from
15 Returns
16 =======
17 Sampling frequency of the time series: float
19 Raises
20 ======
21 ValueError: If the time series is not evenly sampled.
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)
31def get_sampling_frequency_and_duration_from_time_array(time_array):
32 """
33 Calculate sampling frequency and duration from a time array
35 Attributes
36 ==========
37 time_array: array_like
38 Time array to get sampling_frequency/duration from: array_like
40 Returns
41 =======
42 sampling_frequency, duration: float, float
44 Raises
45 ======
46 ValueError: If the time_array is not evenly sampled.
48 """
50 sampling_frequency = get_sampling_frequency(time_array)
51 duration = len(time_array) / sampling_frequency
52 return sampling_frequency, duration
55def get_sampling_frequency_and_duration_from_frequency_array(frequency_array):
56 """
57 Calculate sampling frequency and duration from a frequency array
59 Attributes
60 ==========
61 frequency_array: array_like
62 Frequency array to get sampling_frequency/duration from: array_like
64 Returns
65 =======
66 sampling_frequency, duration: float, float
68 Raises
69 ======
70 ValueError: If the frequency_array is not evenly sampled.
72 """
74 tol = 1e-10
75 if np.ptp(np.diff(frequency_array)) > tol:
76 raise ValueError("Your frequency series was not evenly sampled")
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)
82 sampling_frequency = np.round(2 * (number_of_frequencies - 1) / duration, decimals=14)
83 return sampling_frequency, duration
86def create_time_series(sampling_frequency, duration, starting_time=0.):
87 """
89 Parameters
90 ==========
91 sampling_frequency: float
92 duration: float
93 starting_time: float, optional
95 Returns
96 =======
97 float: An equidistant time series given the parameters
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)
107def create_frequency_series(sampling_frequency, duration):
108 """ Create a frequency series with the correct length and spacing.
110 Parameters
111 ==========
112 sampling_frequency: float
113 duration: float
115 Returns
116 =======
117 array_like: frequency series
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)
124 return np.linspace(start=0,
125 stop=sampling_frequency / 2,
126 num=number_of_frequencies)
129def _check_legal_sampling_frequency_and_duration(sampling_frequency, duration):
130 """ By convention, sampling_frequency and duration have to multiply to an integer
132 This will check if the product of both parameters multiplies reasonably close
133 to an integer.
135 Parameters
136 ==========
137 sampling_frequency: float
138 duration: float
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 )
153def create_white_noise(sampling_frequency, duration):
154 """ Create white_noise which is then coloured by a given PSD
156 Parameters
157 ==========
158 sampling_frequency: float
159 duration: float
160 duration of the data
162 Returns
163 =======
164 array_like: white noise
165 array_like: frequency array
166 """
167 from .random import rng
169 number_of_samples = duration * sampling_frequency
170 number_of_samples = int(np.round(number_of_samples))
172 frequencies = create_frequency_series(sampling_frequency, duration)
174 norm1 = 0.5 * duration**0.5
175 re1, im1 = rng.normal(0, norm1, (2, len(frequencies)))
176 white_noise = re1 + 1j * im1
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
184 # python: transpose for use with infft
185 white_noise = np.transpose(white_noise)
186 frequencies = np.transpose(frequencies)
188 return white_noise, frequencies
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).
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.
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.
208 """
209 frequency_domain_strain = np.fft.rfft(time_domain_strain)
210 frequency_domain_strain /= sampling_frequency
212 frequency_array = np.linspace(
213 0, sampling_frequency / 2, len(frequency_domain_strain))
215 return frequency_domain_strain, frequency_array
218def infft(frequency_domain_strain, sampling_frequency):
219 """ Inverse FFT for use in conjunction with nfft.
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.
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
239class IllegalDurationAndSamplingFrequencyException(Exception):
240 pass