Coverage for src / sdynpy / signal_processing / sdynpy_generator.py: 6%
155 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Functions for generating various excitation signals used in structural dynamics.
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import numpy as np
27def pseudorandom(fft_lines, f_nyq, signal_rms=1, min_freq=None, max_freq=None, filter_oversample=2,
28 integration_oversample=1, averages=1, shape_function=lambda freq: 1):
29 '''Creates a pseudorandom excitation given the specified signal parameters
31 This function creates a pseudorandom excitation from a sum of sine waves
32 with randomized phases. The function can return shaped input if the
33 optional shape_function argument is passed.
35 Parameters
36 ----------
37 fft_lines : int
38 The number of fft lines between 0 Hz and f_nyq, inclusive.
39 f_nyq : float
40 The maximum frequency. The final sampling frequency will be f_nyq x
41 filter_oversample x integration_oversample.
42 signal_rms : float
43 The RMS amplitude of the signal. Default is 1.
44 min_freq : float
45 The minimum frequency of the signal. If not specified, the minimum
46 frequency will be the first frequency line (no 0 Hz component)
47 max_freq : float
48 The maximum frequency of the signal. If not specified, the maximum
49 frequency will be the nyquist frequency f_nyq.
50 filter_oversample : float
51 The oversample factor that is used by some data acquisition systems to
52 provide bandwidth for antialiasing filters. For example, in the I-DEAS
53 data acquisition software, the true sampling frequency is 2.56x the
54 nyquist frequency, rather than the 2x required by the sampling theorem.
55 If not specified, the default factor is 2x.
56 integration_oversample : int
57 An oversampling factor that may be desired for integration of a signal.
58 Default is 1
59 averages : int
60 The number of averages, or frames of data to create. For a
61 pseudorandom input, each frame will be a replica of the first.
62 shape_function : function
63 A shaping that can be applied to the signal per frequency line. If
64 specified, it should be a function that takes in one argument (the
65 frequency in Hz) and returns a single argument (the amplitude of the
66 sine wave at that frequency).
68 Returns
69 -------
70 times : np.ndarray
71 A 1D array containing the time of each sample
72 signal : np.ndarray
73 A 1D array containing the specified signal samples.
74 '''
75 sampling_frequency = f_nyq * filter_oversample
76 oversample_frequency = sampling_frequency * integration_oversample
77 df = f_nyq / fft_lines
78 nbins = sampling_frequency / (2 * df)
79 nsamples = 2 * nbins
80 total_samples = nsamples * integration_oversample
81 oversample_dt = 1 / oversample_frequency
82 times = np.arange(total_samples) * oversample_dt
83 if min_freq is None:
84 min_freq = df
85 if max_freq is None:
86 max_freq = f_nyq
87 frequencies = np.arange(fft_lines) * df
88 max_freq_index = np.argmin(np.abs(frequencies - max_freq))
89 min_freq_index = np.argmin(np.abs(frequencies - min_freq))
90 frequencies = frequencies[min_freq_index:max_freq_index + 1]
91 phases = 2 * np.pi * np.random.rand(*frequencies.shape)
92 amplitudes = np.array([shape_function(frequency) for frequency in frequencies])
93 # Create the sine signal
94 signal = np.sum(amplitudes * np.sin(2 * np.pi * frequencies *
95 times[:, np.newaxis] + phases), -1)
96 # Create the number of averages
97 times = np.arange(total_samples * averages) * oversample_dt
98 signal = np.tile(signal, averages)
99 signal *= signal_rms / np.sqrt(np.mean(signal**2))
100 return times, signal
103def random(shape, n_samples, rms=1.0, dt=1.0,
104 low_frequency_cutoff=None, high_frequency_cutoff=None):
105 signal = np.random.randn(*shape, n_samples)
106 if (low_frequency_cutoff is not None) or (high_frequency_cutoff is not None):
107 freqs = np.fft.rfftfreq(n_samples, dt)
108 signal_fft = np.fft.rfft(signal, axis=-1)
109 if low_frequency_cutoff is not None:
110 zero_indices = freqs < low_frequency_cutoff
111 signal_fft[..., zero_indices] = 0
112 if high_frequency_cutoff is not None:
113 zero_indices = freqs > high_frequency_cutoff
114 signal_fft[..., zero_indices] = 0
115 signal = np.fft.irfft(signal_fft, axis=-1)
116 # Now set RMSs
117 current_rms = np.sqrt(np.mean(signal**2, axis=-1))
118 scale = rms / current_rms
119 signal *= scale[..., np.newaxis]
120 return signal
123def sine(frequencies, dt, num_samples, amplitudes=1, phases=0):
124 times = np.arange(num_samples) * dt
125 output = np.array(amplitudes)[..., np.newaxis] * np.sin(2 * np.pi *
126 np.array(frequencies)[..., np.newaxis] * times + np.array(phases)[..., np.newaxis])
127 return output
130def ramp_envelope(num_samples, ramp_samples, end_ramp_samples=None):
131 output = np.ones(num_samples)
132 output[:ramp_samples] = np.linspace(0, 1, ramp_samples)
133 end_ramp_samples = ramp_samples if end_ramp_samples is None else end_ramp_samples
134 output[-end_ramp_samples:] = np.linspace(1, 0, end_ramp_samples)
135 return output
138def burst_random(shape, n_samples, on_fraction, delay_fraction,
139 rms=1.0, dt=1.0, low_frequency_cutoff=None, high_frequency_cutoff=None,
140 ramp_fraction=0.05):
141 random_signal = random(shape, n_samples, rms, dt, low_frequency_cutoff, high_frequency_cutoff)
142 burst_window = np.zeros(random_signal.shape)
143 delay_samples = int(delay_fraction * n_samples)
144 ramp_samples = int(ramp_fraction * n_samples)
145 on_samples = int(on_fraction * n_samples)
146 burst_window[..., delay_samples:delay_samples +
147 on_samples] = ramp_envelope(on_samples, ramp_samples)
148 return random_signal * burst_window
151def chirp(frequency_min, frequency_max, signal_length, dt, force_integer_cycles=True):
152 n_times = int(signal_length / dt)
153 times = np.arange(n_times) * dt
154 if force_integer_cycles:
155 n_cycles = np.ceil(frequency_max * signal_length)
156 frequency_max = n_cycles / signal_length
157 # print(frequency_max)
158 frequency_slope = (frequency_max - frequency_min) / signal_length
159 argument = frequency_slope / 2 * times**2 + frequency_min * times
160 signal = np.sin(2 * np.pi * argument)
161 return signal
164def pulse(signal_length, pulse_time, pulse_width, pulse_peak=1, dt=1, sine_exponent=1):
165 signal = np.zeros(signal_length)
166 abscissa = np.arange(signal_length) * dt
167 pulse_time, pulse_width, pulse_peak = np.broadcast_arrays(pulse_time, pulse_width, pulse_peak)
168 for time, width, peak in zip(pulse_time.flatten(), pulse_width.flatten(), pulse_peak.flatten()):
169 period = width * 2
170 argument = 2 * np.pi / period * (abscissa - time)
171 signal += peak * np.cos(argument)**sine_exponent * (np.abs(argument) < np.pi / 2)
172 return signal
174def sine_sweep(dt, frequencies, sweep_rates, sweep_types, amplitudes = 1, phases = 0,
175 return_argument = False, return_frequency = False,
176 return_amplitude = False, return_phase = False):
177 """
178 Generates a sweeping sine wave with linear or logarithmic sweep rate
180 Parameters
181 ----------
182 dt : float
183 The time step of the output signal
184 frequencies : iterable
185 A list of frequency breakpoints for the sweep. Can be ascending or
186 decending or both. Frequencies are specified in Hz, not rad/s.
187 sweep_rates : iterable
188 A list of sweep rates between the breakpoints. This array should have
189 one fewer element than the `frequencies` array. The ith element of this
190 array specifies the sweep rate between `frequencies[i]` and
191 `frequencies[i+1]`. For a linear sweep,
192 the rate is in Hz/s. For a logarithmic sweep, the rate is in octave/s.
193 sweep_types : iterable or str
194 The type of sweep to perform between each frequency breakpoint. Can be
195 'lin' or 'log'. If a string is specified, it will be used for all
196 breakpoints. Otherwise it should be an array containing strings with
197 one fewer element than that of the `frequencies` array.
198 amplitudes : iterable or float, optional
199 Amplitude of the cosine wave at each of the frequency breakpoints. Can
200 be specified as a single floating point value, or as an array with a
201 value specified for each breakpoint. The default is 1.
202 phases : iterable or float, optional
203 Phases of the cosine wave at each of the frequency breakpoints. Can
204 be specified as a single floating point value, or as an array with a
205 value specified for each breakpoint. Be aware that modifying the phase
206 between breakpoints will effectively change the frequency of the signal,
207 because the phase will change over time. The default is 0.
208 return_argument : bool
209 If True, return cosine argument over time
210 return_frequency : bool
211 If True, return the instantaneous frequency over time
212 return_amplitude : bool
213 If True, return the instantaneous amplitude over time
214 return_phase : bool
215 If True, return the instantaneous phase over time
217 Raises
218 ------
219 ValueError
220 If the sweep rate and start and end frequency would result in a negative
221 sweep time, for example if the start frequency is above the end frequency
222 and a positive sweep rate is specified.
224 Returns
225 -------
226 ordinate : np.ndarray
227 A numpy array consisting of the generated sine sweep signal. The length
228 of the signal will be determined by the frequency breakpoints and sweep
229 rates.
230 arg_over_time : np.ndarray
231 A numpy array consisting of the argument to the cosine wave over time.
232 freq_over_time : np.ndarray
233 A numpy array consisting of the frequency of the cosine wave over time.
234 amp_over_time : np.ndarray
235 A numpy array consistsing of the amplitude of the cosine wave over time.
236 phs_over_time : np.ndarray
237 A numpy array consisting of the added phase of the cosine wave over time.
239 """
240 last_phase = 0
241 abscissa = []
242 ordinate = []
243 arg_over_time = []
244 freq_over_time = []
245 amp_over_time = []
246 phs_over_time = []
248 # Go through each section
249 for i in range(len(frequencies)-1):
250 # Extract the terms
251 start_frequency = frequencies[i]
252 end_frequency = frequencies[i+1]
253 omega_start = start_frequency*2*np.pi
254 try:
255 sweep_rate = sweep_rates[i]
256 except TypeError:
257 sweep_rate = sweep_rates
258 if isinstance(sweep_types,str):
259 sweep_type = sweep_types
260 else:
261 sweep_type = sweep_types[i]
262 try:
263 start_amplitude = amplitudes[i]
264 end_amplitude = amplitudes[i+1]
265 except TypeError:
266 start_amplitude = amplitudes
267 end_amplitude = amplitudes
268 try:
269 start_phase = phases[i]*np.pi/180
270 end_phase = phases[i+1]*np.pi/180
271 except TypeError:
272 start_phase = phases*np.pi/180
273 end_phase = phases*np.pi/180
274 # Compute the length of this portion of the signal
275 if sweep_type.lower() in ['lin','linear']:
276 sweep_time =+ (end_frequency-start_frequency)/sweep_rate
277 elif sweep_type.lower() in ['log','logarithmic']:
278 sweep_time = np.log(end_frequency/start_frequency)/(sweep_rate*np.log(2))
279 if sweep_time < 0:
280 raise ValueError('Sweep time for segment index {:} is negative. Check sweep rate.'.format(i))
281 sweep_samples = int(np.floor(sweep_time/dt))
282 # Construct the abscissa
283 this_abscissa = np.arange(sweep_samples+1)*dt
284 # Compute the phase over time
285 if sweep_type.lower() in ['lin','linear']:
286 this_argument = (1/2)*(sweep_rate*2*np.pi)*this_abscissa**2 + omega_start*this_abscissa
287 this_frequency = (sweep_rate)*this_abscissa + omega_start/(2*np.pi)
288 elif sweep_type.lower() in ['log','logarithmic']:
289 this_argument = 2**(sweep_rate*this_abscissa)*omega_start/(sweep_rate*np.log(2)) - omega_start/(sweep_rate*np.log(2))
290 this_frequency = 2**(sweep_rate*this_abscissa)*omega_start/(2*np.pi)
291 # Compute the phase at each time step
292 if end_frequency > start_frequency:
293 freq_interp = [start_frequency,end_frequency]
294 phase_interp = [start_phase,end_phase]
295 amp_interp = [start_amplitude,end_amplitude]
296 else:
297 freq_interp = [end_frequency, start_frequency]
298 phase_interp = [end_phase,start_phase]
299 amp_interp = [end_amplitude,start_amplitude]
300 this_phases = np.interp(this_frequency,freq_interp,phase_interp)
301 # if sweep_type in ['lin','linear']:
302 # this_phases = np.interp(this_frequency,freq_interp,phase_interp)
303 # elif sweep_type in ['log','logarithmic']:
304 # this_phase = np.interp(np.log(this_frequency),np.log(freq_interp),np.log(phase_interp))
305 # Compute the amplitude at each time step
306 this_amplitudes = np.interp(this_frequency,freq_interp,amp_interp)
307 this_ordinate = this_amplitudes*np.cos(this_argument+this_phases+last_phase)
308 arg_over_time.append(this_argument[:-1] + last_phase)
309 last_phase += this_argument[-1]
310 abscissa.append(this_abscissa[:-1])
311 ordinate.append(this_ordinate[:-1])
312 freq_over_time.append(this_frequency[:-1])
313 amp_over_time.append(this_amplitudes[:-1])
314 phs_over_time.append(this_phases[:-1])
315 ordinate = np.concatenate(ordinate)
316 return_vals = [ordinate]
317 if return_argument:
318 return_vals.append(np.concatenate(arg_over_time))
319 if return_frequency:
320 return_vals.append(np.concatenate(freq_over_time))
321 if return_amplitude:
322 return_vals.append(np.concatenate(amp_over_time))
323 if return_phase:
324 return_vals.append(np.concatenate(phs_over_time))
325 if len(return_vals) == 1:
326 return_vals = return_vals[0]
327 else:
328 return_vals = tuple(return_vals)
329 return return_vals