Coverage for src / sdynpy / signal_processing / sdynpy_srs.py: 3%
610 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Created on Fri Nov 3 09:23:59 2023
5@author: dprohe
6"""
8import numpy as np
9from scipy.signal import lfilter
10import matplotlib.pyplot as plt
11from scipy.signal import oaconvolve
12from scipy.optimize import minimize, NonlinearConstraint, nnls
15def srs(signal, dt, frequencies=None, damping=0.05, spectrum_type=9,
16 b_filter_weights=None, a_filter_weights=None):
17 """
18 Computes shock response spectrum of a signal
20 Computes the shock response spectrum using a ramp invariant digital filter
21 simulation of the single degree-of-freedom system.
23 Parameters
24 ----------
25 signal : np.ndarray
26 A shape (..., num_samples) ndarray containing the sampled signals to
27 compute SRSs from.
28 dt : float
29 The time between samples (1/sample rate)
30 frequencies : np.ndarray
31 An iterable of frequency lines at which the shock response spectrum
32 will be calculated. If not specified, it will be computed from the
33 sample rate/10000 to the sample rate/4 over 50 lines.
34 damping : float
35 Fraction of critical damping to use in the SRS calculation (e.g. you
36 should specify 0.03 to represent 3%, not 3). The default is 0.03.
37 spectrum_type : int
38 The type of spectrum desired:
39 If `spectrum_type` > 0 (pos) then the SRS will be a base
40 acceleration-absolute acceleration model
41 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
42 displacement model (expressed in equivalent static acceleration units).
43 If abs(`spectrum_type`) is:
44 1--positive primary, 2--negative primary, 3--absolute maximum primary
45 4--positive residual, 5--negative residual, 6--absolute maximum residual
46 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
47 9 -- maximax, the largest absolute value of 1-8
48 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
49 b_filter_weights : np.ndarray, optional
50 Optional filter weights with shape (frequencies.size,3).
51 The default is to automatically compute filter weights using
52 `sdof_ramp_invariant_filter_weights`.
53 a_filter_weights : np.ndarray, optional
54 Optional filter weights with shape (frequencies.size,3).
55 The default is to automatically compute filter weights using
56 `sdof_ramp_invariant_filter_weights`.
58 Returns
59 -------
60 srs : np.ndarray
61 The shock response spectrum at the frequencies specified. If
62 `spectrum_type` == 10 or -10, then this will have shape
63 (... x 9 x len(frequencies)) where each row is a different type of SRS
64 frequencies : np.ndarray
65 The frequencies at which the SRSs were computed.
66 """
67 # Compute default parameters
68 sample_rate = 1/dt
69 if frequencies is None:
70 frequencies = np.logspace(np.log10(sample_rate/1e4),
71 np.log10(sample_rate/4),
72 50)
73 else:
74 frequencies = np.array(frequencies).flatten()
75 if frequencies.size == 0:
76 raise ValueError('Frequencies must have nonzero size')
77 signal = np.array(signal)
78 num_samples = signal.shape[-1]
80 # Compute bad parameters
81 abs_type = abs(spectrum_type)
82 if abs_type > 10:
83 raise ValueError("`spectrum_type` must be in the range [-10,10]")
84 if abs_type == 0:
85 raise ValueError('`spectrum_type` must not be 0')
86 if sample_rate <= 0:
87 raise ValueError('`dt` must be positive')
88 if frequencies.min() < 0:
89 raise ValueError('`frequencies` must all be positive')
90 if damping < 0:
91 raise ValueError('`damping` must be positive')
92 if damping == 0:
93 print('Warning: Damping is Zero in SRS calculation!')
95 # Get the filter coefficients if not specified
96 if b_filter_weights is None or a_filter_weights is None:
97 b_filter_weights, a_filter_weights = sdof_ramp_invariant_filter_weights(
98 frequencies, sample_rate, damping, spectrum_type)
99 else:
100 b_filter_weights = np.atleast_2d(np.array(b_filter_weights))
101 a_filter_weights = np.atleast_2d(np.array(a_filter_weights))
103 if b_filter_weights.shape != (frequencies.size, 3):
104 raise ValueError('`b_filter_weights` and `a_filter_weights` must have shape (frequencies.size,3)')
106 srs_output = np.zeros(signal.shape[:-1]
107 + ((9,) if np.abs(spectrum_type) == 10 else ())
108 + (frequencies.size,))
109 for i_freq, (frequency, b, a) in enumerate(zip(
110 frequencies, b_filter_weights, a_filter_weights)):
111 # Append enough zeros to get 1/100 of the residual response
112 zeros_to_append = int(np.max([2, np.ceil(0.01*sample_rate/frequency)]))
113 signals_extended = np.concatenate((
114 signal, np.zeros(signal.shape[:-1] + (zeros_to_append,))), axis=-1)
116 # Filter the signal
117 response = sdof_filter(b, a, signals_extended)
119 # Now compute the SRS parameters
120 # Primary response
121 primary_maximum = np.array(np.max(response[..., :num_samples], axis=-1))
122 primary_maximum[primary_maximum < 0] = 0
123 primary_minimum = np.array(np.abs(np.min(response[..., :num_samples], axis=-1)))
124 primary_minimum[primary_minimum < 0] = 0
125 primary_abs = np.array(np.max([primary_maximum, primary_minimum], axis=0))
127 # Now compute the residual response. we need time steps and
128 # responses at two points
129 residual_times = [0, (zeros_to_append-1)/sample_rate]
130 residual_responses = np.concatenate([response[..., num_samples, np.newaxis],
131 response[..., num_samples+zeros_to_append-1, np.newaxis]],
132 axis=-1)
134 # Compute the peak response after the structure has stopped responding
135 peak_residual_responses = sdof_free_decay_peak_response(
136 residual_responses, residual_times, frequency, damping)
138 # Pull out the residual peaks
139 residual_maximum = np.array(np.max(peak_residual_responses, axis=-1))
140 residual_maximum[residual_maximum < 0] = 0
141 residual_minimum = np.array(np.abs(np.min(peak_residual_responses, axis=-1)))
142 residual_minimum[residual_minimum < 0] = 0
143 residual_abs = np.array(np.max([residual_maximum, residual_minimum], axis=0))
145 all_srs_this_frequency_line = np.array([
146 primary_maximum, primary_minimum, primary_abs,
147 residual_maximum, residual_minimum, residual_abs,
148 np.max([primary_maximum, residual_maximum], axis=0),
149 np.max([primary_minimum, residual_minimum], axis=0),
150 np.max([primary_abs, residual_abs], axis=0)])
152 if abs_type < 10:
153 srs_output[..., i_freq] = all_srs_this_frequency_line[abs_type-1]
154 else:
155 srs_output[..., i_freq] = np.moveaxis(all_srs_this_frequency_line, 0, -1)
156 return srs_output, frequencies
159def sdof_ramp_invariant_filter_weights(frequencies, sample_rate, damping,
160 spectrum_type):
161 """
162 Computes filter weights for SDOF resonators using a ramp-invariant filter.
164 The weights are used in conjunction with the function
165 `sdof_filter` and `srs` to calculate the shock response
166 spectrum of an acceleration time history using a ramp
167 invarient simulation.
169 Parameters
170 ----------
171 frequencies : np.ndarray
172 The frequencies at which filters should be computed
173 sample_rate : float
174 The sample rate of the measurement
175 damping : float
176 The fraction of critical damping for the system (e.g. 0.03, not 3 for 3%)
177 spectrum_type : int
178 See `spectrum_type` argument for compute_srs
180 Returns
181 -------
182 b : np.ndarray
183 Filter coefficients with shape (frequencies.shape,3)
184 a : np.ndarray
185 Filter coefficients with shape (frequencies.shape,3)
186 """
187 # Make sure it's a numpy array
188 frequencies = np.array(frequencies)
189 # Preallocate arrays
190 a = np.ones(frequencies.shape+(3,))
191 b = np.ones(frequencies.shape+(3,))
193 normalized_frequencies = 2*np.pi*frequencies/sample_rate
195 small_freq_indices = normalized_frequencies < 0.001
197 # Small frequencies
198 small_freqs = normalized_frequencies[small_freq_indices]
200 # 2*z*w + w*w*(1-2*z*z) where z is damping and w is normalized frequencies
201 a[small_freq_indices, 1] = (
202 2*damping*small_freqs
203 + small_freqs**2*(1-2*damping**2)
204 )
205 # -2*z*w + 2*z*z*w*w where z is damping and w is normalized frequencies
206 a[small_freq_indices, 2] = (
207 -2*damping*small_freqs
208 + 2*damping**2*small_freqs**2
209 )
211 if spectrum_type > 0:
212 # Absolute Acceleration Model
213 # z*w + (w*w)*( (1/6) - 2*z*z/3 );
214 b[small_freq_indices, 0] = (
215 damping*small_freqs
216 + small_freqs**2
217 * (1/6 - 2*damping**2/3)
218 )
219 # 2*w*w*(1-z*z)/3
220 b[small_freq_indices, 1] = (
221 2*small_freqs**2*(1-damping**2)/3
222 )
223 # -z*w + w*w*( (1/6) - 4*z*z/3 );
224 b[small_freq_indices, 2] = (
225 - damping*small_freqs
226 + small_freqs**2
227 * (1/6 - 4*damping**2/3)
228 )
229 else:
230 # Relative Displacement Model
231 # -w*w/6;
232 b[small_freq_indices, 0] = -small_freqs**2/6
233 # -2*w*w/3;
234 b[small_freq_indices, 1] = -2*small_freqs**2/3
235 # -w*w/6;
236 b[small_freq_indices, 2] = -small_freqs**2/6
238 # Large frequencies
239 large_freqs = normalized_frequencies[~small_freq_indices]
241 # Define some helper variables
242 sq = np.sqrt(1-damping**2)
243 e = np.exp(-damping*large_freqs)
244 wd = large_freqs*sq
245 sp = e*np.sin(wd)
246 fact = (2*damping**2 - 1)*sp/sq
247 c = e*np.cos(wd)
249 a[~small_freq_indices, 1] = 2*(1-c)
250 a[~small_freq_indices, 2] = -1 + e**2
252 if spectrum_type > 0:
253 # Absolute Acceleration Model
254 spwd = sp/wd
255 b[~small_freq_indices, 0] = 1-spwd
256 b[~small_freq_indices, 1] = 2*(spwd - c)
257 b[~small_freq_indices, 2] = e**2 - spwd
258 else:
259 # Relative Displacement Model
260 b[~small_freq_indices, 0] = -(2*damping*(c-1) + fact + large_freqs)/large_freqs
261 b[~small_freq_indices, 1] = -(-2*c*large_freqs + 2*damping*(1-e**2) - 2*fact)/large_freqs
262 b[~small_freq_indices, 2] = -(e**2*(large_freqs+2*damping) - 2*damping*c + fact)/large_freqs
264 return b, a
267def sdof_filter(b, a, signal, zi=None):
268 """
269 Applies a filter to simulate a single degree of freedom system
271 Parameters
272 ----------
273 b : np.ndarray
274 Size 3 array representing filter coefficients used by scipy.signal.lfilter
275 a : np.array
276 Size 3 array representing filter coefficients used by scipy.signal.lfilter
277 signal : np.ndarray
278 The signals that are to be filtered
279 zi : np.ndarray, optional
280 Optional initial state for the filters, having length
281 `max(len(a), len(b)) - 1`. If not specified, zero initial conditions
282 are assumed.
284 Returns
285 -------
286 filtered_signal : np.ndarray
287 The filtered signal
289 """
290 aa = np.array([a[0], a[1]-2, a[2] + 1])
291 return lfilter(b, aa, signal, axis=-1, zi=zi)
294def sdof_free_decay_peak_response(responses, times_at_responses, frequency, damping):
295 """
296 Calculates peak response of a freely-decaying sdof system.
298 The residual response is the peak response of the sdof system
299 as it decays after the input has ended, i.e., the input is zero.
300 The first two peaks in the decaying response will be the largest.
301 One peak will be negative and one positive. We don't know before
302 the calculation if the first and largest peak in amplitude will be
303 positive or negative.
305 Parameters
306 ----------
307 responses : np.array
308 A (...,2) shape array giving the response at two time steps
309 times_at_responses : np.array
310 A length-2 array giving the time values that the responses occur at
311 frequency : float
312 The frequency at which the computation is performed
313 damping : float
314 The fraction of critical damping for the system (e.g. 0.03, not 3 for 3%)
316 Returns
317 -------
318 response_peaks : np.ndarray
319 A shape (...,2) array giving the amplitude at the two peaks after the
320 response has decayed
322 Notes
323 -----
324 The response has the form
325 a(t) = exp(-zeta wn t)[z[0]sin(wd t) + z[1]cos(wd t)].
326 If I know a(t) at two values of time, t, I can calculate the constants
327 z[0] and z[1]. The general form of the response can then be solved for
328 the maximum by finding the time of the first maximum by setting the
329 derivative to zero. Then substituting the time of the maximum response
330 back into the general equation to find the maximum response.
331 The second peak will occur half a cycle later.
333 """
334 responses = np.array(responses)
335 times_at_responses = np.array(times_at_responses)
336 # Set up some initial helper variables
337 fd = frequency*np.sqrt(1-damping*damping)
338 wd = 2*np.pi*fd
339 wn = 2*np.pi*frequency
340 wdt = wd*times_at_responses
341 e = np.exp(-wn*damping*times_at_responses)
342 sd = np.sin(wdt)
343 cd = np.cos(wdt)
344 c = e[:, np.newaxis]*np.array((sd, cd)).T
345 z = np.linalg.solve(c, responses[..., np.newaxis])[..., 0]
347 # the response has the form
348 # a(t) = exp(-zeta wn t)[z[0]sin(wd t) + z[1]cos(wd t)]
349 # the first maximum response will occur at the time tmax
350 tmax = np.array((1/wd) * np.arctan((wd*z[..., 0] - damping*wn*z[..., 1])/(wd*z[..., 1] - damping*wn*z[..., 0])))
351 tmax[tmax < 0] = tmax[tmax < 0] + np.pi/wd
352 tmax[tmax > np.pi/wd] = tmax[tmax > np.pi/wd] - np.pi/wd
354 # We can now plug it into the equation to find the maximum
355 response_peaks = []
356 response_peaks.append(np.exp(-damping*wn*tmax)*(z[..., 0]*np.sin(wd*tmax) + z[..., 1]*np.cos(wd*tmax)))
357 tmin = tmax + np.pi/wd
358 response_peaks.append(np.exp(-damping*wn*tmin)*(z[..., 0]*np.sin(wd*tmin) + z[..., 1]*np.cos(wd*tmin)))
359 return np.moveaxis(response_peaks, 0, -1)
362def octspace(low, high, points_per_octave):
363 """
364 Constructs octave spacing between low and high values
366 Parameters
367 ----------
368 low : float
369 Starting value for the spacing
370 high : float
371 Upper value for the spacing
372 points_per_octave : int
373 Number of points per octave
375 Returns
376 -------
377 octave_points : np.ndarray
378 Octave-spaced points
379 """
380 num_octaves = np.log2(high/low)
381 num_steps = np.ceil(num_octaves*points_per_octave)
382 point_indices = np.arange(num_steps+1)
383 log_points = np.log2(low) + num_octaves/num_steps*point_indices
384 points = 2**log_points
385 return points
388def sum_decayed_sines(sample_rate, block_size,
389 sine_frequencies=None, sine_tone_range=None, sine_tone_per_octave=None,
390 sine_amplitudes=None, sine_decays=None, sine_delays=None,
391 required_srs=None, srs_breakpoints=None,
392 srs_damping=0.05, srs_type=9,
393 compensation_frequency=None, compensation_decay=0.95,
394 # Paramters for the iteration
395 number_of_iterations=3, convergence=0.8,
396 error_tolerance=0.05,
397 tau=None, num_time_constants=None, decay_resolution=None,
398 scale_factor=1.02,
399 acceleration_factor=1.0,
400 plot_results=False, srs_frequencies=None,
401 ignore_compensation_pulse = False,
402 verbose=False
403 ):
404 """
405 Generate a Sum of Decayed Sines signal given an SRS.
407 Note that there are many approaches to do this, with many optional arguments
408 so please read the documentation carefully to understand which arguments
409 must be passed to the function.
411 Parameters
412 ----------
413 sample_rate : float
414 The sample rate of the generated signal.
415 block_size : int
416 The number of samples in the generated signal.
417 sine_frequencies : np.ndarray, optional
418 The frequencies of the sine tones. If this argument is not specified,
419 then the `sine_tone_range` argument must be specified.
420 sine_tone_range : np.ndarray, optional
421 A length-2 array containing the minimum and maximum sine tone to
422 generate. If this is argument is not specified, then the
423 `sine_frequencies` argument must be specified instead.
424 sine_tone_per_octave : int, optional
425 The number of sine tones per octave. If not specified along with
426 `sine_tone_range`, then a default value of 4 will be used if the
427 `srs_damping` is >= 0.05. Otherwise, the formula of
428 `sine_tone_per_octave = 9 - srs_damping*100` will be used.
429 sine_amplitudes : np.ndarray, optional
430 The initial amplitude of the sine tones used in the optimization. If
431 not specified, they will be set to the value of the SRS at each frequency
432 divided by the quality factor of the SRS.
433 sine_decays : np.ndarray, optional
434 An array of decay value time constants (often represented by variable
435 tau). Tau is the time for the amplitude of motion to decay 63% defined
436 by the equation `1/(2*np.pi*freq*zeta)` where `freq` is the frequency
437 of the sine tone and `zeta` is the fraction of critical damping.
438 If not specified, then either the `tau` or `num_time_constants`
439 arguments must be specified instead.
440 sine_delays : np.ndarray, optional
441 An array of delay values for the sine components. If not specified,
442 all tones will have zero delay.
443 required_srs : np.ndarray, optional
444 The SRS to match defined at each of the `sine_frequencies`. If this
445 argument is not passed, then the `srs_breakpoints` argument must be
446 passed instead. The default is None.
447 srs_breakpoints : np.ndarray, optional
448 A numpy array with shape `(n,2)` where the first column is the
449 frequencies at which each breakpoint occurs, and the second column is
450 the value of the SRS at that breakpoint. SRS values at the
451 `sine_frequencies` will be interpolated in a log-log sense from this
452 breakpoint array. If this argument is not specified, then the
453 `required_srs` must be passed instead.
454 srs_damping : float, optional
455 Fraction of critical damping to use in the SRS calculation (e.g. you
456 should specify 0.03 to represent 3%, not 3). The default is 0.03.
457 srs_type : int
458 The type of spectrum desired:
459 If `srs_type` > 0 (pos) then the SRS will be a base
460 acceleration-absolute acceleration model
461 If `srs_type` < 0 (neg) then the SRS will be a base acceleration-relative
462 displacement model (expressed in equivalent static acceleration units).
463 If abs(`srs_type`) is:
464 1--positive primary, 2--negative primary, 3--absolute maximum primary
465 4--positive residual, 5--negative residual, 6--absolute maximum residual
466 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
467 9 -- maximax, the largest absolute value of 1-8
468 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
469 compensation_frequency : float
470 The frequency of the compensation pulse. If not specified, it will be
471 set to 1/3 of the lowest sine tone
472 compensation_decay : float
473 The decay value for the compensation pulse. If not specified, it will
474 be set to 0.95.
475 number_of_iterations : int, optional
476 The number of iterations to perform. At least two iterations should be
477 performed. 3 iterations is preferred, and will be used if this argument
478 is not specified.
479 convergence : float, optional
480 The fraction of the error corrected each iteration. The default is 0.8.
481 error_tolerance : float, optional
482 Allowable relative error in the SRS. The default is 0.05.
483 tau : float, optional
484 If a floating point number is passed, then this will be used for the
485 `sine_decay` values. Alternatively, a dictionary can be passed with
486 the keys containing a length-2 tuple specifying the minimum and maximum
487 frequency range, and the value specifying the value of `tau` within that
488 frequency range. If this latter approach is used, all `sine_frequencies`
489 must be contained within a frequency range. If this argument is not
490 specified, then either `sine_decays` or `num_time_constants` must be
491 specified instead.
492 num_time_constants : int, optional
493 If an integer is passed, then this will be used to set the `sine_decay`
494 values by ensuring the specified number of time constants occur in the
495 `block_size`. Alternatively, a dictionary can be passed with the keys
496 containing a length-2 tuple specifying the minimum and maximum
497 frequency range, and the value specifying the value of
498 `num_time_constants` over that frequency range. If this latter approach
499 is used, all `sine_frequencies` must be contained within a frequency
500 range. If this argument is not specified, then either `sine_decays` or
501 `tau` must be specified instead.
502 decay_resolution : float, optional
503 A scalar identifying the resolution of the fractional decay rate
504 (often known by the variable `zeta`). The decay parameters will be
505 rounded to this value. The default is to not round.
506 scale_factor : float, optional
507 A scaling applied to the sine tone amplitudes so the achieved SRS better
508 fits the specified SRS, rather than just touching it. The default is 1.02.
509 acceleration_factor : float, optional
510 Optional scale factor to convert acceleration into velocity and
511 displacement. For example, if sine amplitudes are in G and displacement
512 is desired in inches, the acceleration factor should be set to 386.089.
513 If sine amplitudes are in G and displacement is desired in meters, the
514 acceleration factor should be set to 9.80665. The default is 1, which
515 assumes consistent units (e.g. acceleration in m/s^2, velocity in m/s,
516 displacement in m).
517 plot_results : bool, optional
518 If True, a figure will be plotted showing the acceleration, velocity,
519 and displacement signals, as well as the desired and achieved SRS.
520 srs_frequencies : np.ndarray, optional
521 If specified, these frequencies will be used to compute the SRS that
522 will be plotted when the `plot_results` value is `True`.
523 ignore_compensation_pulse : bool, optional
524 If True, the compensation pulse will be ignored.
525 verbose : True, optional
526 If True, additional diagnostics will be printed to the console.
528 Returns
529 -------
530 acceleration_signal : ndarray
531 The acceleration signal that satisfies the SRS.
532 velocity_signal : ndarray
533 The velocity of the acceleration signal.
534 displacement_signal : ndarray
535 The displacement of the acceleration signal.
536 sine_frequencies : ndarray
537 An array of frequencies for each sine tone, including the compensation
538 pulse
539 sine_amplitudes : ndarray
540 An array of amplitudes for each sine tone, including the compensation
541 pulse
542 sine_decays : ndarray
543 An array of decay values for each sine tone, including the compensation
544 pulse
545 sine_delays : ndarray
546 An array of delay values for each sine tone, including the compensation
547 pulse
548 fig : matplotlib.Figure
549 A reference to the plotted figure. Only returned if `plot_results` is
550 `True`.
551 ax : matplotlib.Axes
552 A reference to the plotted axes. Only returned if `plot_results` is
553 `True`.
555 """
556 # Handle the sine tone frequencies
557 if sine_frequencies is None and sine_tone_range is None:
558 raise ValueError('Either `sine_frequencies` or `sine_tone_range` must be specified')
559 if sine_frequencies is not None and sine_tone_range is not None:
560 raise ValueError('`sine_frequencies` can not be specified simultaneously with `sine_tone_range`')
561 if sine_frequencies is None:
562 # Create sine tones
563 if sine_tone_per_octave is None:
564 sine_tone_per_octave = int(np.floor(9-srs_damping*100))
565 sine_frequencies = octspace(sine_tone_range[0], sine_tone_range[1],
566 sine_tone_per_octave)
567 # Now set up the SRS
568 if required_srs is None and srs_breakpoints is None:
569 raise ValueError('Either `required_srs` or `srs_breakpoints` must be specified')
570 if required_srs is not None and srs_breakpoints is not None:
571 raise ValueError('`required_srs` can not be specified simultaneously with `srs_breakpoints`')
572 if required_srs is None:
573 required_srs = loginterp(sine_frequencies,
574 srs_breakpoints[:, 0],
575 srs_breakpoints[:, 1])
576 if sine_amplitudes is None:
577 srs_amplitudes = required_srs.copy()
578 srs_amplitudes[np.arange(srs_amplitudes.size) % 2 == 0] *= -1
579 quality = 1/(2*srs_damping)
580 srs_amplitudes /= quality
581 sine_amplitudes = srs_amplitudes
583 if sine_delays is None:
584 sine_delays = np.zeros(sine_frequencies.size)
586 if compensation_frequency is None:
587 compensation_frequency = np.min(sine_frequencies)/3
589 # Set up decay terms
590 decay_terms_specified = 0
591 if sine_decays is not None:
592 decay_terms_specified += 1
593 if tau is not None:
594 decay_terms_specified += 1
595 if num_time_constants is not None:
596 decay_terms_specified += 1
597 if decay_terms_specified == 0:
598 raise ValueError('One of `sine_decays`, `tau`, or `num_time_constants` must be specified')
599 if decay_terms_specified > 1:
600 raise ValueError('Only one of `sine_decays`, `tau`, or `num_time_constants` can be specified')
602 # Now check and see which is defined
603 if num_time_constants is not None:
604 period = block_size / sample_rate
605 if isinstance(num_time_constants, dict):
606 tau = {}
607 for freq_range, num_time_constant in num_time_constants.items():
608 tau[freq_range] = period / num_time_constant
609 else:
610 tau = period/num_time_constants
611 if tau is not None:
612 sine_decays = []
613 for freq in sine_frequencies:
614 if isinstance(tau, dict):
615 this_decay = None
616 for freq_range, this_tau in tau.items():
617 if freq_range[0] <= freq <= freq_range[1]:
618 this_decay = 1/(2*np.pi*freq*this_tau)
619 break
620 if this_decay is None:
621 raise ValueError('No frequency range matching frequency {:} was found in the specified decay parameters.'.format(freq))
622 sine_decays.append(this_decay)
623 else:
624 sine_decays.append(1/(2*np.pi*freq*tau))
625 sine_decays = np.array(sine_decays)
626 # Otherwise we just keep the specified sine_decays
628 # Now handle the minimum resolution on the decay
629 if decay_resolution is not None:
630 sine_decays = decay_resolution*np.round(sine_decays/decay_resolution)
632 if compensation_frequency is None:
633 compensation_frequency = np.min(sine_frequencies)/3
635 # Now we can actually run the generation process
636 (acceleration_signal, used_frequencies, used_amplitudes, used_decays,
637 used_delays, used_comp_frequency, used_comp_amplitude, used_comp_decay,
638 used_comp_delay) = _sum_decayed_sines(
639 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
640 scale_factor*required_srs, compensation_frequency, compensation_decay,
641 sample_rate, block_size, srs_damping, srs_type, number_of_iterations,
642 convergence, error_tolerance, ignore_compensation_pulse, verbose)
644 all_frequencies = np.concatenate((used_frequencies,
645 [used_comp_frequency]))
646 all_amplitudes = np.concatenate((used_amplitudes,
647 [used_comp_amplitude]))
648 all_decays = np.concatenate((used_decays,
649 [used_comp_decay]))
650 all_delays = np.concatenate((used_delays,
651 [used_comp_delay]))
652 # Now compute displacement and velocity
653 velocity_signal, displacement_signal = sum_decayed_sines_displacement_velocity(
654 all_frequencies, all_amplitudes, all_decays, all_delays, sample_rate,
655 block_size, acceleration_factor)
657 return_vals = (acceleration_signal, velocity_signal, displacement_signal,
658 all_frequencies, all_amplitudes, all_decays, all_delays)
660 # Plot the results
661 if plot_results:
662 fig, ax = plt.subplots(2, 2, figsize=(8, 6))
663 times = np.arange(block_size)/sample_rate
664 ax[0, 0].plot(times, acceleration_signal)
665 ax[0, 0].set_ylabel('Acceleration')
666 ax[0, 0].set_xlabel('Time (s)')
667 ax[0, 1].plot(times, velocity_signal)
668 ax[0, 1].set_ylabel('Velocity')
669 ax[0, 1].set_xlabel('Time (s)')
670 ax[1, 0].plot(times, displacement_signal)
671 ax[1, 0].set_ylabel('Displacement')
672 ax[1, 0].set_xlabel('Time (s)')
673 # Compute SRS
674 if srs_frequencies is None:
675 srs_frequencies = sine_frequencies
676 this_srs, this_frequencies = srs(
677 acceleration_signal, 1/sample_rate, srs_frequencies,
678 srs_damping, srs_type)
679 if srs_breakpoints is None:
680 srs_abscissa = sine_frequencies
681 srs_ordinate = required_srs
682 else:
683 srs_abscissa, srs_ordinate = srs_breakpoints.T
684 ax[1, 1].plot(srs_abscissa, srs_ordinate, 'k--')
685 ax[1, 1].plot(this_frequencies, this_srs)
686 ax[1, 1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
687 ax[1, 1].set_xlabel('Frequency (Hz)')
688 ax[1, 1].set_yscale('log')
689 ax[1, 1].set_xscale('log')
690 fig.tight_layout()
691 return_vals += (fig, ax)
692 ax[1, 1].legend(('Reference', 'Decayed Sine'))
693 return return_vals
696def _sum_decayed_sines(sine_frequencies, sine_amplitudes,
697 sine_decays, sine_delays,
698 required_srs,
699 compensation_frequency, compensation_decay,
700 sample_rate, block_size,
701 srs_damping, srs_type,
702 number_of_iterations=3, convergence=0.8,
703 error_tolerance=0.05, ignore_compensation_pulse = False,
704 verbose=False):
705 """
706 Optimizes the amplitudes of sums of decayed sines
708 Parameters
709 ----------
710 sine_frequencies : ndarray
711 An array of frequencies for each sine tone
712 sine_amplitudes : ndarray
713 An array of amplitudes for each sine tone
714 sine_decays : ndarray
715 An array of decay values for each sine tone
716 sine_delays : ndarray
717 An array of delay values for each sine tone
718 required_srs : ndarray
719 An array of SRS values for each frequency in `sine_frequencies`
720 compensation_frequency : float
721 The frequency of the compensation pulse
722 compensation_decay : float
723 The decay value for the compensation pulse
724 sample_rate : float
725 The sample rate of the signal
726 block_size : int
727 The number of samples in the signal
728 srs_damping : float
729 Fraction of critical damping to use in the SRS calculation (e.g. you
730 should specify 0.03 to represent 3%, not 3). The default is 0.03.
731 srs_type : int
732 The type of spectrum desired:
733 If `spectrum_type` > 0 (pos) then the SRS will be a base
734 acceleration-absolute acceleration model
735 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
736 displacement model (expressed in equivalent static acceleration units).
737 If abs(`spectrum_type`) is:
738 1--positive primary, 2--negative primary, 3--absolute maximum primary
739 4--positive residual, 5--negative residual, 6--absolute maximum residual
740 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
741 9 -- maximax, the largest absolute value of 1-8
742 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
743 number_of_iterations : int
744 The number of iterations that will be performed. The default is 3.
745 convergence : float, optional
746 The fraction of the error corrected each iteration. The default is 0.8.
747 error_tolerance : float, optional
748 Allowable relative error in the SRS. The default is 0.05.
749 ignore_compensation_pulse : bool, optional
750 If True, ignores the compensation pulse. Default is false.
751 verbose : bool, optional
752 If True, information on the interations will be provided. The default
753 is False.
755 Returns
756 -------
757 this_signal : np.ndarray
758 A numpy array containing the generated signal
759 sine_frequencies : ndarray
760 An array of frequencies for each sine tone
761 sine_amplitudes : ndarray
762 An array of amplitudes for each sine tone
763 sine_decays : ndarray
764 An array of decay values for each sine tone
765 sine_delays : ndarray
766 An array of delay values for each sine tone
767 this_compensation_frequency : float
768 The frequency of the compensation term.
769 this_compensation_amplitude : float
770 The amplitude of the compensation term.
771 this_compensation_decay : float
772 The decay constant of the compensation term
773 this_compensation_delay : float
774 The delay constant of the compenstation term
775 """
776 for i in range(number_of_iterations):
777 if verbose:
778 print('Starting Iteration {:}'.format(i+1))
779 # Compute updated sine amplitudes and compensation terms
780 this_sine_amplitudes, this_compensation_amplitude, this_compensation_delay = (
781 _sum_decayed_sines_single_iteration(
782 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
783 required_srs, compensation_frequency, compensation_decay,
784 sample_rate, block_size, srs_damping, srs_type,
785 number_of_iterations=10, convergence=convergence,
786 error_tolerance=error_tolerance, ignore_compensation_pulse=ignore_compensation_pulse,
787 verbose=verbose))
788 # get the SRS error at each frequency term by first computing the signal
789 (this_signal, this_compensation_frequency, this_compensation_amplitude,
790 this_compensation_decay, this_compensation_delay) = sum_decayed_sines_reconstruction_with_compensation(
791 sine_frequencies, this_sine_amplitudes, sine_decays, sine_delays,
792 compensation_frequency, compensation_decay, sample_rate, block_size, ignore_compensation_pulse=ignore_compensation_pulse)
793 # Then computing the SRS of the signal
794 this_srs = srs(this_signal, 1/sample_rate, sine_frequencies, srs_damping,
795 srs_type)[0]
796 srs_error = (this_srs - required_srs)/required_srs
797 sine_amplitudes = this_sine_amplitudes
798 if verbose:
799 print('Sine Table after Iterations:')
800 print('{:>8s}, {:>10s}, {:>10s}, {:>10s}, {:>10s}, {:>10s}'.format(
801 'Sine', 'Frequency', 'Amplitude', 'SRS Val', 'SRS Req', 'Error'))
802 for i, (freq, amp, srs_val, srs_req, srs_err) in enumerate(zip(
803 sine_frequencies, sine_amplitudes, this_srs, required_srs,
804 srs_error)):
805 print('{:>8s}, {:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}'.format(
806 str(i), freq, amp, srs_val, srs_req, srs_err))
807 print('Compensating Pulse:')
808 print('{:>10s}, {:>10s}, {:>10s}, {:>10s}'.format(
809 'Frequency', 'Amplitude', 'Decay', 'Delay'))
810 print('{:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}'.format(
811 this_compensation_frequency, this_compensation_amplitude,
812 this_compensation_decay, this_compensation_delay))
813 return (this_signal, sine_frequencies, sine_amplitudes, sine_decays,
814 sine_delays, this_compensation_frequency, this_compensation_amplitude,
815 this_compensation_decay, this_compensation_delay)
818def _sum_decayed_sines_single_iteration(sine_frequencies, sine_amplitudes,
819 sine_decays, sine_delays,
820 required_srs, compensation_frequency,
821 compensation_decay,
822 sample_rate, block_size,
823 damping_srs, srs_type,
824 number_of_iterations=10, convergence=0.8,
825 error_tolerance=0.05,
826 ignore_compensation_pulse=False,
827 verbose=False):
828 """
829 Iterates on amplitudes of decayed sine waves to match a prescribed SRS
831 Parameters
832 ----------
833 sine_frequencies : ndarray
834 An array of frequencies for each sine tone
835 sine_amplitudes : ndarray
836 An array of amplitudes for each sine tone
837 sine_decays : ndarray
838 An array of decay values for each sine tone
839 sine_delays : ndarray
840 An array of delay values for each sine tone
841 required_srs : ndarray
842 An array of SRS values for each frequency in `sine_frequencies`
843 compensation_frequency : float
844 The frequency of the compensation pulse
845 compensation_decay : float
846 The decay value for the compensation pulse
847 sample_rate : float
848 The sample rate of the signal
849 block_size : int
850 The number of samples in the signal
851 damping_srs : float
852 Fraction of critical damping to use in the SRS calculation (e.g. you
853 should specify 0.03 to represent 3%, not 3). The default is 0.03.
854 srs_type : int
855 The type of spectrum desired:
856 If `spectrum_type` > 0 (pos) then the SRS will be a base
857 acceleration-absolute acceleration model
858 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
859 displacement model (expressed in equivalent static acceleration units).
860 If abs(`spectrum_type`) is:
861 1--positive primary, 2--negative primary, 3--absolute maximum primary
862 4--positive residual, 5--negative residual, 6--absolute maximum residual
863 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
864 9 -- maximax, the largest absolute value of 1-8
865 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
866 number_of_iterations : int
867 Maximum number of iterations that can be performed at each frequency line.
868 Default is 10.
869 convergence : float, optional
870 The fraction of the error corrected each iteration. The default is 0.8.
871 error_tolerance : float, optional
872 Allowable relative error in the SRS. The default is 0.05.
873 ignore_compensation_pulse : bool, optional
874 If True, do not use a compensation pulse. Default is false.
875 verbose : bool, optional
876 If True, information on the interations will be provided. The default
877 is False.
879 Raises
880 ------
881 ValueError
882 If compensation delay is required to be too long.
884 Returns
885 -------
886 sine_amplitudes : np.ndarray
887 An array the same size as the input `sine_amplitudes` array with updated
888 amplitudes.
889 compensation_amplitude : float
890 Amplitude of the compensation pulse
891 compensation_delay : float
892 Delay of the compensation pulse.
893 """
894 # Define some helper variables
895 Amax = np.max(np.abs(sine_amplitudes))
896 # Reduce the error tolerance so there's a bit of room for round-off
897 error_tolerance = error_tolerance*0.9
898 # Get filter weights for SRS calculations
899 b, a = sdof_ramp_invariant_filter_weights(sine_frequencies, sample_rate,
900 damping_srs, srs_type)
901 # Copy the arrays so we don't overwrite anything
902 sine_frequencies = np.array(sine_frequencies).copy()
903 sine_amplitudes = np.array(sine_amplitudes).copy()
904 sine_decays = np.array(sine_decays).copy()
905 sine_delays = np.array(sine_delays).copy()
906 # Iterate one frequency at a time
907 for i, (frequency, amplitude, decay, delay, desired_srs) in enumerate(zip(
908 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
909 required_srs)):
910 # if i == 10:
911 # break
912 srs_error = float('inf')
913 iteration_count = 1
914 increment = 0.1
915 # Get the pulse without this frequency line
916 other_sine_frequencies = np.concatenate((sine_frequencies[:i], sine_frequencies[i+1:]))
917 other_sine_amplitudes = np.concatenate((sine_amplitudes[:i], sine_amplitudes[i+1:]))
918 other_sine_decays = np.concatenate((sine_decays[:i], sine_decays[i+1:]))
919 other_sine_delays = np.concatenate((sine_delays[:i], sine_delays[i+1:]))
920 other_pulse = sum_decayed_sines_reconstruction(
921 other_sine_frequencies, other_sine_amplitudes, other_sine_decays,
922 other_sine_delays, sample_rate, block_size)
923 # Now iterate at the single frequency until convergence
924 while abs(srs_error) > error_tolerance:
925 sine_amplitudes[i] = amplitude
926 # Find the pulse at this frequency line
927 this_pulse = sum_decayed_sines_reconstruction(
928 frequency, amplitude, decay, delay, sample_rate, block_size)
929 # Get the compensating pulse
930 if not ignore_compensation_pulse:
931 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
932 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
933 else:
934 compensation_amplitude = 0
935 compensation_delay = 0
936 # Find the number of samples to shift the waveform
937 num_shift = -int(np.floor(compensation_delay*sample_rate))
938 if abs(num_shift) >= block_size:
939 raise ValueError('The number of samples in the compensation delay ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
940 + ' The entire pulse will consist of part of the compensation pulse.'
941 + ' Please increase the block_size or compensation frequency.')
942 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
943 # Find compensating time history
944 compensation_pulse = sum_decayed_sines_reconstruction(
945 compensation_frequency, compensation_amplitude, compensation_decay,
946 compensation_delay_corrected, sample_rate, block_size)
947 # Build the composite waveform, need to shift the signals to align
948 if num_shift >= 0:
949 composite_pulse = (
950 compensation_pulse
951 + np.concatenate((
952 np.zeros(num_shift),
953 (other_pulse + this_pulse)[:block_size-num_shift])))
954 else:
955 composite_pulse = (
956 other_pulse + this_pulse +
957 + np.concatenate((
958 np.zeros(-num_shift),
959 (compensation_pulse)[:block_size+num_shift])))
960 # Find the SRS at the current frequency line
961 srs_prediction = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
962 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
963 srs_error = (srs_prediction - desired_srs)/desired_srs
964 if verbose:
965 print('Iteration at frequency {:}, {:0.4f}\n Iteration Count: {:}, SRS Error: {:0.4f}'.format(i, frequency, iteration_count, srs_error))
966 # Now we're going to compute the same thing again with a perturbed
967 # amplitude, this will allow us to compute the slope change at the
968 # current amplitude
969 amplitude_change = np.sign(srs_error)*increment*np.sign(amplitude)*Amax
970 # Now check and see if we need to modify the amplitude
971 if abs(srs_error) > error_tolerance:
972 if amplitude_change == 0: # perturb it a bit
973 amplitude_change = np.sign(srs_error)*np.sign(sine_amplitudes[i])*Amax*increment/10
974 new_amplitude = amplitude + amplitude_change
975 # Can't allow the sign of the amplitude to change
976 if np.sign(amplitude) != np.sign(new_amplitude):
977 new_amplitude *= -1
978 if amplitude == new_amplitude:
979 new_amplitude = amplitude * 1.01
980 sine_amplitudes[i] = new_amplitude
981 # Find new component at this frequency
982 this_pulse = sum_decayed_sines_reconstruction(
983 frequency, new_amplitude, decay, delay, sample_rate, block_size)
984 # Get the compensating pulse
985 if not ignore_compensation_pulse:
986 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
987 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
988 else:
989 compensation_amplitude = 0
990 compensation_delay = 0
991 # Find the number of samples to shift the waveform
992 num_shift = -int(np.floor(compensation_delay*sample_rate))
993 if abs(num_shift) >= block_size:
994 raise ValueError('The number of samples in the compensation delay ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
995 + ' The entire pulse will consist of part of the compensation pulse. '
996 + 'Please increase the block_size or compensation frequency.')
997 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
998 # Find compensating time history
999 compensation_pulse = sum_decayed_sines_reconstruction(
1000 compensation_frequency, compensation_amplitude, compensation_decay,
1001 compensation_delay_corrected, sample_rate, block_size)
1002 # Build the composite waveform, need to shift the signals to align
1003 if num_shift >= 0:
1004 composite_pulse = (
1005 compensation_pulse
1006 + np.concatenate((
1007 np.zeros(num_shift),
1008 (other_pulse + this_pulse)[:block_size-num_shift])))
1009 else:
1010 composite_pulse = (
1011 other_pulse + this_pulse +
1012 + np.concatenate((
1013 np.zeros(-num_shift),
1014 (compensation_pulse)[:block_size+num_shift])))
1015 # Find the SRS at the current frequency line
1016 srs_perturbed = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
1017 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
1018 # Get slope of correction
1019 correction_slope = (abs(amplitude)-abs(new_amplitude))/(srs_prediction-srs_perturbed)
1020 if correction_slope > 0:
1021 # this is what we want, it means we are changing in the right
1022 # direction
1023 new_amplitude = convergence*correction_slope*(desired_srs-srs_prediction) + abs(amplitude)
1024 # Never let the change be too large
1025 amplitude_check = max([abs(amplitude), Amax/10])
1026 if new_amplitude > 2*amplitude_check:
1027 new_amplitude = 2*amplitude_check
1028 # Never let it be less than zero
1029 if new_amplitude < 0:
1030 new_amplitude = np.abs(amplitude)/10
1031 # Move the previous amplitude into a different variable
1032 old_amplitude = amplitude
1033 # Never allow new amplitude to be zero
1034 if new_amplitude == 0:
1035 amplitude = 2*np.finfo(float).eps*np.sign(amplitude)
1036 else:
1037 amplitude = new_amplitude*np.sign(amplitude)
1038 # Store it for the next iteration
1039 sine_amplitudes[i] = amplitude
1040 # Find new component at this frequency
1041 this_pulse = sum_decayed_sines_reconstruction(
1042 frequency, amplitude, decay, delay, sample_rate, block_size)
1043 # Get the compensating pulse
1044 if not ignore_compensation_pulse:
1045 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1046 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
1047 else:
1048 compensation_amplitude = 0
1049 compensation_delay = 0
1050 # Find the number of samples to shift the waveform
1051 num_shift = -int(np.floor(compensation_delay*sample_rate))
1052 if abs(num_shift) >= block_size:
1053 raise ValueError('The number of samples in the compensation delay'
1054 + ' ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
1055 + ' The entire pulse will consist of part of the compensation pulse. '
1056 + 'Please increase the block_size or compensation frequency.')
1057 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
1058 # Find compensating time history
1059 compensation_pulse = sum_decayed_sines_reconstruction(
1060 compensation_frequency, compensation_amplitude, compensation_decay,
1061 compensation_delay_corrected, sample_rate, block_size)
1062 # Build the composite waveform, need to shift the signals to align
1063 if num_shift >= 0:
1064 composite_pulse = (
1065 compensation_pulse
1066 + np.concatenate((
1067 np.zeros(num_shift),
1068 (other_pulse + this_pulse)[:block_size-num_shift])))
1069 else:
1070 composite_pulse = (
1071 other_pulse + this_pulse +
1072 + np.concatenate((
1073 np.zeros(-num_shift),
1074 (compensation_pulse)[:block_size+num_shift])))
1075 # Find the SRS at the current frequency line
1076 srs_corrected = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
1077 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
1078 srs_error = (srs_corrected - desired_srs)/desired_srs
1079 # If the SRS error is positive and the amplitude has already been reduced to zero, stop trying
1080 if srs_error > 0 and amplitude == 0:
1081 iteration_count = number_of_iterations + 1
1082 if verbose:
1083 print(' Old Amplitude: {:0.4f}, New Amplitude: {:0.4f}, Error: {:0.4f}'.format(old_amplitude, amplitude, srs_error))
1084 # Don't allow amplitude to change more than a factor of 10 in an iteration
1085 if np.abs(10*old_amplitude) < np.abs(amplitude) or np.abs(0.1*old_amplitude) > np.abs(amplitude):
1086 iteration_count = number_of_iterations + 1
1087 else:
1088 # The slope is negative, so we try a bigger increment
1089 increment *= 2
1090 amplitude_change = np.sign(srs_error)*increment*np.sign(amplitude)*Amax
1091 if verbose:
1092 print('Slope of correction was negative, trying a bigger increment')
1093 iteration_count += 1
1094 if iteration_count > number_of_iterations:
1095 print(' Warning: SRS did not converge for frequency {:}: {:0.4f}'.format(i, frequency))
1096 break
1097 if not ignore_compensation_pulse:
1098 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1099 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
1100 else:
1101 compensation_amplitude = 0
1102 compensation_delay = 0
1103 return sine_amplitudes, compensation_amplitude, compensation_delay
1106def sum_decayed_sines_compensating_pulse_parameters(sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1107 compensation_frequency, compensation_decay):
1108 omegas = sine_frequencies*2*np.pi
1109 omega_comp = compensation_frequency*2*np.pi
1110 var = -np.sum(sine_amplitudes/(sine_frequencies*(sine_decays**2+1)))
1111 compensation_amplitude = compensation_frequency*(compensation_decay**2+1)*var
1112 var0 = (np.sum(sine_amplitudes*sine_delays/(omegas*(sine_decays**2+1))))/compensation_amplitude
1113 # TODO Verify if this is a bug or not, I think this should be matlab code but graflab code had it as a matrix division.
1114 # var1 = (np.sum(2*sine_decays*sine_amplitudes/(omegas**2*(sine_decays**2 + 1)**2)))/compensation_amplitude
1115 var1 = (np.sum((2*sine_decays*sine_amplitudes)[np.newaxis, :]@np.linalg.pinv((omegas**2*(sine_decays**2 + 1)**2)[np.newaxis, :])))/compensation_amplitude
1116 var2 = 2*compensation_decay/(omega_comp*omega_comp*(compensation_decay**2+1)**2)
1117 var3 = omega_comp*(compensation_decay**2+1)
1118 compensation_delay = -var3*(var2 + var1 + var0)
1119 return compensation_amplitude, compensation_delay
1122def sum_decayed_sines_reconstruction(sine_frequencies, sine_amplitudes,
1123 sine_decays, sine_delays, sample_rate,
1124 block_size):
1125 """
1126 Computes a sum of decayed sines signal from a set of frequencies, amplitudes,
1127 decays, and delays.
1129 Parameters
1130 ----------
1131 sine_frequencies : ndarray
1132 An array of frequencies for each sine tone
1133 sine_amplitudes : ndarray
1134 An array of amplitudes for each sine tone
1135 sine_decays : ndarray
1136 An array of decay values for each sine tone
1137 sine_delays : ndarray
1138 An array of delay values for each sine tone
1139 sample_rate : float
1140 The sample rate of the signal
1141 block_size : int
1142 The number of samples in the signal
1144 Returns
1145 -------
1146 ndarray
1147 A signal containing the sum of decayed sinusoids.
1149 """
1150 times = np.arange(block_size)/sample_rate
1151 # See if any delays are below zero and if so then shift all delays forward
1152 if np.any(sine_delays < 0):
1153 sine_delays = sine_delays - np.min(sine_delays)
1154 # Now go through and compute the sine tones
1155 omegas = 2*np.pi*sine_frequencies
1156 this_times = times[:, np.newaxis] - sine_delays
1157 response = sine_amplitudes*np.exp(-sine_decays*omegas*this_times)*np.sin(omegas*this_times)
1158 response[this_times < 0] = 0
1159 return np.sum(response, axis=-1)
1162def sum_decayed_sines_reconstruction_with_compensation(
1163 sine_frequencies, sine_amplitudes,
1164 sine_decays, sine_delays, compensation_frequency, compensation_decay,
1165 sample_rate, block_size, ignore_compensation_pulse = False):
1166 """
1169 Parameters
1170 ----------
1171 sine_frequencies : ndarray
1172 An array of frequencies for each sine tone
1173 sine_amplitudes : ndarray
1174 An array of amplitudes for each sine tone
1175 sine_decays : ndarray
1176 An array of decay values for each sine tone
1177 sine_delays : ndarray
1178 An array of delay values for each sine tone
1179 compensation_frequency : float
1180 The frequency of the compensation pulse
1181 compensation_decay : float
1182 The decay value for the compensation pulse
1183 sample_rate : float
1184 The sample rate of the signal
1185 block_size : int
1186 The number of samples in the signal
1187 ignore_compensation_pulse : bool, optional
1188 If True, ignores the compensation pulse
1190 Returns
1191 -------
1192 signal : ndarray
1193 A signal containing the sum of decayed sinusoids.
1194 compensation_frequency : float
1195 The frequency of the compensation pulse
1196 compensation_amplitude : float
1197 The amplitude value for the compensation pulse
1198 compensation_decay : float
1199 The decay value for the compensation pulse
1200 compensation_delay : float
1201 The delay value for the compensation pulse
1202 """
1203 sine_frequencies = np.array(sine_frequencies).flatten()
1204 sine_amplitudes = np.array(sine_amplitudes).flatten()
1205 sine_delays = np.array(sine_delays).flatten()
1206 sine_decays = np.array(sine_decays).flatten()
1207 if not ignore_compensation_pulse:
1208 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1209 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1210 compensation_frequency, compensation_decay)
1211 else:
1212 compensation_amplitude = 0
1213 compensation_delay = 0
1214 sine_frequencies = np.concatenate((sine_frequencies, [compensation_frequency]))
1215 sine_amplitudes = np.concatenate((sine_amplitudes, [compensation_amplitude]))
1216 sine_delays = np.concatenate((sine_delays, [compensation_delay]))
1217 sine_decays = np.concatenate((sine_decays, [compensation_decay]))
1218 signal = sum_decayed_sines_reconstruction(
1219 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1220 sample_rate, block_size)
1221 return (signal, compensation_frequency, compensation_amplitude,
1222 compensation_decay, compensation_delay)
1225def sum_decayed_sines_displacement_velocity(
1226 sine_frequencies, sine_amplitudes,
1227 sine_decays, sine_delays, sample_rate, block_size,
1228 acceleration_factor=1):
1229 """
1230 Creates velocity and displacement signals from acceleration sinusoids.
1232 Parameters
1233 ----------
1234 sine_frequencies : ndarray
1235 An array of frequencies for each sine tone
1236 sine_amplitudes : ndarray
1237 An array of amplitudes for each sine tone
1238 sine_decays : ndarray
1239 An array of decay values for each sine tone
1240 sine_delays : ndarray
1241 An array of delay values for each sine tone
1242 sample_rate : float
1243 The sample rate of the signal
1244 block_size : int
1245 The number of samples in the signal
1246 acceleration_factor : float, optional
1247 Optional scale factor to convert acceleration into velocity and
1248 displacement. For example, if sine amplitudes are in G and displacement
1249 is desired in inches, the acceleration factor should be set to 386.089.
1250 If sine amplitudes are in G and displacement is desired in meters, the
1251 acceleration factor should be set to 9.80665. The default is 1, which
1252 assumes consistent units (e.g. acceleration in m/s^2, velocity in m/s,
1253 displacement in m).
1255 Returns
1256 -------
1257 v : ndarray
1258 The velocity of the signal.
1259 d : ndarray
1260 The displacement of the signal.
1262 """
1263 # Make sure everything is a numpy array
1264 sine_frequencies = np.array(sine_frequencies).flatten()
1265 sine_amplitudes = np.array(sine_amplitudes).flatten()
1266 sine_delays = np.array(sine_delays).flatten()
1267 sine_decays = np.array(sine_decays).flatten()
1268 # Transform units
1269 sine_amplitudes = sine_amplitudes * acceleration_factor
1270 # See if any delays are below zero and if so then shift all delays forward
1271 if np.any(sine_delays < 0):
1272 sine_delays = sine_delays - np.min(sine_delays)
1273 f = sine_frequencies
1274 A = sine_amplitudes
1275 z = sine_decays
1276 tau = sine_delays
1277 x = np.zeros(block_size)
1278 tmp1 = x.copy()
1279 tmp2 = x.copy()
1280 tmp3 = x.copy()
1281 x1 = np.ones(x.shape)
1282 v = x.copy()
1283 d = x.copy()
1284 w = 2*np.pi*f
1285 w2 = w*w
1286 zw = z*w
1287 zp1 = z*z + 1
1288 zm1 = z*z - 1
1289 t = np.arange(block_size)/sample_rate
1290 for k in range(len(sine_frequencies)):
1291 indices = t-tau[k] > 0 # index's where vel and disp are evaluated
1292 Awz1 = A[k]/(w[k]*zp1[k])
1293 tmp1[indices] = Awz1*np.exp(-zw[k]*(t[indices]-tau[k]))
1294 tmp2[indices] = z[k]*np.sin(w[k]*(t[indices]-tau[k])) + np.cos(w[k]*(t[indices]-tau[k]))
1295 v[indices] = v[indices] - tmp1[indices]*tmp2[indices] + Awz1*x1[indices]
1296 Awz2 = A[k]/(w2[k]*zp1[k]*zp1[k])
1297 tmp1[indices] = Awz2*np.exp(-zw[k]*(t[indices]-tau[k]))
1298 tmp2[indices] = zm1[k]*np.sin(w[k]*(t[indices]-tau[k])) + 2*z[k]*np.cos(w[k]*(t[indices]-tau[k]))
1299 tmp3[indices] = Awz1*(t[indices]-tau[k])
1300 tmp4 = 2*z[k]*Awz2
1301 d[indices] = d[indices] + tmp1[indices]*tmp2[indices] + tmp3[indices] - tmp4
1302 return v, d
1305def loginterp(x, xp, fp):
1306 return 10**np.interp(np.log10(x), np.log10(xp), np.log10(fp))
1309def optimization_error_function(
1310 amplitude_lin_scales, sample_rate, block_size,
1311 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1312 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1313 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1314 b=None, a=None, frequency_index=None):
1315 if frequency_index is None:
1316 frequency_index = slice(None)
1317 # Apply the scale factors
1318 new_sine_amplitudes = sine_amplitudes.copy()
1319 new_sine_amplitudes[frequency_index] *= amplitude_lin_scales
1320 # Compute a new time history signal
1321 (input_signal, new_compensation_frequency, new_compensation_amplitude,
1322 new_compensation_decay, new_compensation_delay) = sum_decayed_sines_reconstruction_with_compensation(
1323 sine_frequencies, new_sine_amplitudes, sine_decays, sine_delays,
1324 compensation_frequency, compensation_decay, sample_rate, block_size)
1325 # Transform the signal to the control degrees of freedom
1326 control_responses = np.array([oaconvolve(control_irf, input_signal) for control_irf in control_irfs])[..., :block_size]
1327 # Compute SRSs at all frequencies
1328 predicted_control_srs, frequencies = srs(control_responses, 1/sample_rate, sine_frequencies[frequency_index],
1329 srs_damping, srs_type,
1330 None if b is None else b[frequency_index],
1331 None if a is None else a[frequency_index])
1332 # Compute error
1333 mean_control_error = np.mean(
1334 (control_weights[:, np.newaxis]
1335 * ((predicted_control_srs - control_srs[..., frequency_index])
1336 / control_srs[..., frequency_index])))
1337 if limit_srs is not None:
1338 limit_responses = np.array([oaconvolve(limit_irf, input_signal) for limit_irf in limit_irfs])[..., :block_size]
1339 predicted_limit_srs, frequencies = srs(limit_responses, 1/sample_rate, sine_frequencies[frequency_index],
1340 srs_damping, srs_type,
1341 None if b is None else b[frequency_index],
1342 None if a is None else a[frequency_index])
1343 max_limit_error = np.max(
1344 ((predicted_limit_srs - limit_srs[..., frequency_index])
1345 / limit_srs[..., frequency_index]))
1346 if max_limit_error >= 0 and mean_control_error >= 0:
1347 srs_error = np.max((mean_control_error, max_limit_error))
1348 elif max_limit_error <= 0 and mean_control_error <= 0:
1349 srs_error = np.max((mean_control_error, max_limit_error))
1350 elif max_limit_error >= 0 and mean_control_error <= 0:
1351 srs_error = max_limit_error
1352 elif max_limit_error <= 0 and mean_control_error >= 0:
1353 srs_error = mean_control_error
1354 else:
1355 limit_responses = None
1356 predicted_limit_srs = None
1357 srs_error = mean_control_error
1359 return (np.abs(srs_error), input_signal, control_responses,
1360 predicted_control_srs, limit_responses, predicted_limit_srs,
1361 new_sine_amplitudes, new_compensation_frequency,
1362 new_compensation_amplitude, new_compensation_decay,
1363 new_compensation_delay)
1366def optimization_callback(intermediate_result, rms_error_threshold=0.02,
1367 verbose=True):
1368 if verbose:
1369 print('Amplitude Scale: {:}, error is {:}'.format(intermediate_result.x.squeeze(),
1370 intermediate_result.fun))
1371 if rms_error_threshold is not None:
1372 if intermediate_result.fun < rms_error_threshold:
1373 raise StopIteration
1376def sum_decayed_sines_minimize(sample_rate, block_size,
1377 sine_frequencies=None, sine_tone_range=None, sine_tone_per_octave=None,
1378 sine_amplitudes=None, sine_decays=None, sine_delays=None,
1379 control_srs=None, control_breakpoints=None,
1380 srs_damping=0.05, srs_type=9,
1381 compensation_frequency=None, compensation_decay=0.95,
1382 # Parameters for defining decays
1383 tau=None, num_time_constants=None, decay_resolution=None,
1384 scale_factor=1.02,
1385 acceleration_factor=1.0,
1386 # Parameters for imposing limits
1387 limit_breakpoints=None, limit_transfer_functions=None,
1388 control_transfer_functions=None, control_weights=None,
1389 # Parameters for the optimizer
1390 minimize_iterations=1, rms_error_threshold=None,
1391 optimization_passes=3,
1392 plot_results=False, verbose=False
1393 ):
1394 # Handle the sine tone frequencies
1395 if sine_frequencies is None and sine_tone_range is None:
1396 raise ValueError('Either `sine_frequencies` or `sine_tone_range` must be specified')
1397 if sine_frequencies is not None and sine_tone_range is not None:
1398 raise ValueError('`sine_frequencies` can not be specified simultaneously with `sine_tone_range`')
1399 if sine_frequencies is None:
1400 # Create sine tones
1401 if sine_tone_per_octave is None:
1402 sine_tone_per_octave = int(np.floor(9-srs_damping*100))
1403 sine_frequencies = octspace(sine_tone_range[0], sine_tone_range[1],
1404 sine_tone_per_octave)
1405 # Now set up the SRS
1406 if control_srs is None and control_breakpoints is None:
1407 raise ValueError('Either `control_srs` or `control_breakpoints` must be specified')
1408 if control_srs is not None and control_breakpoints is not None:
1409 raise ValueError('`control_srs` can not be specified simultaneously with `control_breakpoints`')
1410 if control_srs is None:
1411 frequencies = control_breakpoints[:, 0]
1412 breakpoint_curves = control_breakpoints[:, 1:].T
1413 control_srs = np.array([loginterp(sine_frequencies,
1414 frequencies,
1415 breakpoint_curve) for breakpoint_curve in breakpoint_curves])
1416 else:
1417 control_srs = np.atleast_2d(control_srs)
1418 if control_transfer_functions is None:
1419 control_transfer_functions = np.ones(((control_srs.shape[0]-1)*2, block_size//2+1))
1420 tf_frequencies = np.fft.rfftfreq(control_transfer_functions.shape[-1]*2-1, 1/sample_rate)
1421 if sine_amplitudes is None:
1422 tf_at_frequencies = np.array([
1423 np.interp(sine_frequencies, tf_frequencies, np.abs(control_transfer_function))
1424 for control_transfer_function in control_transfer_functions])
1425 srs_amplitudes = np.array([nnls(ai[:, np.newaxis], bi)[0][0] for ai, bi in zip(tf_at_frequencies.T, control_srs.T)])
1426 srs_amplitudes[np.arange(srs_amplitudes.size) % 2 == 0] *= -1
1427 quality = 1/(2*srs_damping)
1428 srs_amplitudes /= quality
1429 sine_amplitudes = srs_amplitudes
1431 if sine_delays is None:
1432 sine_delays = np.zeros(sine_frequencies.size)
1434 if compensation_frequency is None:
1435 compensation_frequency = np.min(sine_frequencies)/3
1437 # Set up decay terms
1438 decay_terms_specified = 0
1439 if sine_decays is not None:
1440 decay_terms_specified += 1
1441 if tau is not None:
1442 decay_terms_specified += 1
1443 if num_time_constants is not None:
1444 decay_terms_specified += 1
1445 if decay_terms_specified == 0:
1446 raise ValueError('One of `sine_decays`, `tau`, or `num_time_constants` must be specified')
1447 if decay_terms_specified > 1:
1448 raise ValueError('Only one of `sine_decays`, `tau`, or `num_time_constants` can be specified')
1450 # Now check and see which is defined
1451 if num_time_constants is not None:
1452 period = block_size / sample_rate
1453 if isinstance(num_time_constants, dict):
1454 tau = {}
1455 for freq_range, num_time_constant in num_time_constants.items():
1456 tau[freq_range] = period / num_time_constant
1457 else:
1458 tau = period/num_time_constants
1459 if tau is not None:
1460 sine_decays = []
1461 for freq in sine_frequencies:
1462 if isinstance(tau, dict):
1463 this_decay = None
1464 for freq_range, this_tau in tau.items():
1465 if freq_range[0] <= freq <= freq_range[1]:
1466 this_decay = 1/(2*np.pi*freq*this_tau)
1467 break
1468 if this_decay is None:
1469 raise ValueError('No frequency range matching frequency {:} was found in the specified decay parameters.'.format(freq))
1470 sine_decays.append(this_decay)
1471 else:
1472 sine_decays.append(1/(2*np.pi*freq*tau))
1473 sine_decays = np.array(sine_decays)
1474 # Otherwise we just keep the specified sine_decays
1476 # Now handle the minimum resolution on the decay
1477 if decay_resolution is not None:
1478 sine_decays = decay_resolution*np.round(sine_decays/decay_resolution)
1480 if compensation_frequency is None:
1481 compensation_frequency = np.min(sine_frequencies)/3
1483 # Now set up limits and transfer functions
1484 if control_weights is None:
1485 control_weights = np.ones((control_srs.shape[0]))
1486 if limit_breakpoints is not None:
1487 frequencies = limit_breakpoints[:, 0]
1488 breakpoint_curves = limit_breakpoints[:, 1:].T
1489 limit_srs = np.array([loginterp(sine_frequencies,
1490 frequencies,
1491 breakpoint_curve) for breakpoint_curve in breakpoint_curves])
1492 else:
1493 limit_srs = None
1495 # Compute impulse responses
1496 control_irfs = np.fft.irfft(control_transfer_functions, axis=-1)
1497 if limit_transfer_functions is not None:
1498 limit_irfs = np.fft.irfft(limit_transfer_functions, axis=-1)
1499 else:
1500 limit_irfs = None
1501 b, a = sdof_ramp_invariant_filter_weights(sine_frequencies, sample_rate, srs_damping, srs_type)
1503 # Normalize control weights
1504 control_weights = control_weights/np.linalg.norm(control_weights)
1506 # Copy things so we don't overwrite
1507 sine_amplitudes = sine_amplitudes.copy()
1509 # Now we will iterate over all of the frequencies and compute the new
1510 # amplitudes
1511 for j in range(optimization_passes):
1512 for i in range(len(sine_frequencies)):
1513 if verbose:
1514 print('Pass {:}, Analyzing Frequency {:}: {:0.2f}'.format(j+1, i, sine_frequencies[i]))
1516 # We will now go through and optimize the frequency line
1517 def error_function(x):
1518 return optimization_error_function(
1519 x, sample_rate, block_size,
1520 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1521 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1522 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1523 b, a, frequency_index=[i])[0]
1525 def callback(intermediate_result):
1526 return optimization_callback(
1527 intermediate_result, rms_error_threshold, verbose)
1529 optimization_result = minimize(error_function,
1530 np.ones(1),
1531 method='Powell',
1532 bounds=[(0, np.inf)],
1533 callback=callback,
1534 options={'maxiter': minimize_iterations})
1536 # Populate the amplitudes with the updated result
1537 if verbose:
1538 print('Initial Amplitude: {:}, Updated Amplitude: {:}\n'.format(sine_amplitudes[i], sine_amplitudes[i]*optimization_result.x.squeeze()))
1539 sine_amplitudes[i] *= optimization_result.x.squeeze()
1541 (error, input_signal, control_responses, predicted_control_srs,
1542 limit_responses, predicted_limit_srs,
1543 sine_amplitudes, compensation_frequency,
1544 compensation_amplitude, compensation_decay,
1545 compensation_delay) = optimization_error_function(
1546 np.ones(sine_frequencies.size),
1547 sample_rate, block_size,
1548 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1549 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1550 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1551 b, a)
1553 all_frequencies = np.concatenate((sine_frequencies,
1554 [compensation_frequency]))
1555 all_amplitudes = np.concatenate((sine_amplitudes,
1556 [compensation_amplitude]))
1557 all_decays = np.concatenate((sine_decays,
1558 [compensation_decay]))
1559 all_delays = np.concatenate((sine_delays,
1560 [compensation_delay]))
1562 # Now compute displacement and velocity
1563 velocity_signal, displacement_signal = sum_decayed_sines_displacement_velocity(
1564 all_frequencies, all_amplitudes, all_decays, all_delays, sample_rate,
1565 block_size, acceleration_factor)
1567 return_vals = (input_signal, velocity_signal, displacement_signal,
1568 control_responses, predicted_control_srs,
1569 limit_responses, predicted_limit_srs,
1570 all_frequencies, all_amplitudes, all_decays, all_delays,
1571 )
1573 if plot_results:
1574 fig, ax = plt.subplots(2, 2, figsize=(8, 6))
1575 times = np.arange(block_size)/sample_rate
1576 ax[0, 0].plot(times, input_signal)
1577 ax[0, 0].set_ylabel('Acceleration')
1578 ax[0, 0].set_xlabel('Time (s)')
1579 ax[0, 1].plot(times, velocity_signal)
1580 ax[0, 1].set_ylabel('Velocity')
1581 ax[0, 1].set_xlabel('Time (s)')
1582 ax[1, 0].plot(times, displacement_signal)
1583 ax[1, 0].set_ylabel('Displacement')
1584 ax[1, 0].set_xlabel('Time (s)')
1585 # Compute SRS
1586 this_srs, this_frequencies = srs(
1587 input_signal, 1/sample_rate, sine_frequencies,
1588 srs_damping, srs_type)
1589 if control_breakpoints is None:
1590 srs_abscissa = sine_frequencies
1591 srs_ordinate = control_srs.T
1592 else:
1593 srs_abscissa = control_breakpoints[:, 0]
1594 srs_ordinate = control_breakpoints[:, 1:]
1595 ax[1, 1].plot(srs_abscissa, srs_ordinate, 'k--')
1596 ax[1, 1].plot(this_frequencies, this_srs)
1597 ax[1, 1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1598 ax[1, 1].set_xlabel('Frequency (Hz)')
1599 ax[1, 1].set_yscale('log')
1600 ax[1, 1].set_xscale('log')
1601 ax[1, 1].legend(('Reference', 'Decayed Sine'))
1602 fig.tight_layout()
1603 return_vals += (fig, ax)
1604 # Now plot all control channels
1605 for i, (predicted, control, response) in enumerate(zip(predicted_control_srs, control_srs, control_responses)):
1606 fig, ax = plt.subplots(1, 2, figsize=(8, 3))
1607 ax[0].set_title('Control Signal {:}'.format(i))
1608 ax[0].plot(times, response)
1609 ax[0].set_label('Acceleration')
1610 ax[0].set_xlabel('Time (s)')
1611 ax[1].set_title('Control SRS {:}'.format(i))
1612 ax[1].plot(sine_frequencies, control, 'k--')
1613 ax[1].plot(sine_frequencies, predicted)
1614 ax[1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1615 ax[1].set_xlabel('Frequency (Hz)')
1616 ax[1].set_yscale('log')
1617 ax[1].set_xscale('log')
1618 ax[1].legend(('Desired', 'Achieved'))
1619 fig.tight_layout()
1620 return_vals += (fig, ax)
1621 # Now plot all limit channels
1622 if limit_srs is not None:
1623 for i, (predicted, limit, response) in enumerate(zip(predicted_limit_srs, limit_srs, limit_responses)):
1624 fig, ax = plt.subplots(1, 2, figsize=(8, 3))
1625 ax[0].set_title('Limit Signal {:}'.format(i))
1626 ax[0].plot(times, response)
1627 ax[0].set_label('Acceleration')
1628 ax[0].set_xlabel('Time (s)')
1629 ax[1].set_title('Limit SRS {:}'.format(i))
1630 ax[1].plot(sine_frequencies, limit, 'k--')
1631 ax[1].plot(sine_frequencies, predicted)
1632 ax[1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1633 ax[1].set_xlabel('Frequency (Hz)')
1634 ax[1].set_yscale('log')
1635 ax[1].set_xscale('log')
1636 ax[1].legend(('Desired', 'Achieved'))
1637 fig.tight_layout()
1638 return_vals += (fig, ax)
1640 return return_vals