Coverage for src/sdynpy/signal_processing/sdynpy_srs.py: 8%
880 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-18 21:32 +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, windows
12from scipy.optimize import minimize, NonlinearConstraint, nnls
13from scipy.interpolate import interp1d
14from scipy.integrate import cumulative_trapezoid
15from dataclasses import dataclass
16from typing import Optional, List, Dict, Any
19def srs(signal, dt, frequencies=None, damping=0.05, spectrum_type=9,
20 b_filter_weights=None, a_filter_weights=None):
21 """
22 Computes shock response spectrum of a signal
24 Computes the shock response spectrum using a ramp invariant digital filter
25 simulation of the single degree-of-freedom system.
27 Parameters
28 ----------
29 signal : np.ndarray
30 A shape (..., num_samples) ndarray containing the sampled signals to
31 compute SRSs from.
32 dt : float
33 The time between samples (1/sample rate)
34 frequencies : np.ndarray
35 An iterable of frequency lines at which the shock response spectrum
36 will be calculated. If not specified, it will be computed from the
37 sample rate/10000 to the sample rate/4 over 50 lines.
38 damping : float
39 Fraction of critical damping to use in the SRS calculation (e.g. you
40 should specify 0.03 to represent 3%, not 3). The default is 0.03.
41 spectrum_type : int
42 The type of spectrum desired:
43 If `spectrum_type` > 0 (pos) then the SRS will be a base
44 acceleration-absolute acceleration model
45 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
46 displacement model (expressed in equivalent static acceleration units).
47 If abs(`spectrum_type`) is:
48 1--positive primary, 2--negative primary, 3--absolute maximum primary
49 4--positive residual, 5--negative residual, 6--absolute maximum residual
50 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
51 9 -- maximax, the largest absolute value of 1-8
52 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
53 b_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`.
57 a_filter_weights : np.ndarray, optional
58 Optional filter weights with shape (frequencies.size,3).
59 The default is to automatically compute filter weights using
60 `sdof_ramp_invariant_filter_weights`.
62 Returns
63 -------
64 srs : np.ndarray
65 The shock response spectrum at the frequencies specified. If
66 `spectrum_type` == 10 or -10, then this will have shape
67 (... x 9 x len(frequencies)) where each row is a different type of SRS
68 frequencies : np.ndarray
69 The frequencies at which the SRSs were computed.
70 """
71 # Compute default parameters
72 sample_rate = 1/dt
73 if frequencies is None:
74 frequencies = np.logspace(np.log10(sample_rate/1e4),
75 np.log10(sample_rate/4),
76 50)
77 else:
78 frequencies = np.array(frequencies).flatten()
79 if frequencies.size == 0:
80 raise ValueError('Frequencies must have nonzero size')
81 signal = np.array(signal)
82 num_samples = signal.shape[-1]
84 # Compute bad parameters
85 abs_type = abs(spectrum_type)
86 if abs_type > 10:
87 raise ValueError("`spectrum_type` must be in the range [-10,10]")
88 if abs_type == 0:
89 raise ValueError('`spectrum_type` must not be 0')
90 if sample_rate <= 0:
91 raise ValueError('`dt` must be positive')
92 if frequencies.min() < 0:
93 raise ValueError('`frequencies` must all be positive')
94 if damping < 0:
95 raise ValueError('`damping` must be positive')
96 if damping == 0:
97 print('Warning: Damping is Zero in SRS calculation!')
99 # Get the filter coefficients if not specified
100 if b_filter_weights is None or a_filter_weights is None:
101 b_filter_weights, a_filter_weights = sdof_ramp_invariant_filter_weights(
102 frequencies, sample_rate, damping, spectrum_type)
103 else:
104 b_filter_weights = np.atleast_2d(np.array(b_filter_weights))
105 a_filter_weights = np.atleast_2d(np.array(a_filter_weights))
107 if b_filter_weights.shape != (frequencies.size, 3):
108 raise ValueError('`b_filter_weights` and `a_filter_weights` must have shape (frequencies.size,3)')
110 srs_output = np.zeros(signal.shape[:-1]
111 + ((9,) if np.abs(spectrum_type) == 10 else ())
112 + (frequencies.size,))
113 for i_freq, (frequency, b, a) in enumerate(zip(
114 frequencies, b_filter_weights, a_filter_weights)):
115 # Append enough zeros to get 1/100 of the residual response
116 zeros_to_append = int(np.max([2, np.ceil(0.01*sample_rate/frequency)]))
117 signals_extended = np.concatenate((
118 signal, np.zeros(signal.shape[:-1] + (zeros_to_append,))), axis=-1)
120 # Filter the signal
121 response = sdof_filter(b, a, signals_extended)
123 # Now compute the SRS parameters
124 # Primary response
125 primary_maximum = np.array(np.max(response[..., :num_samples], axis=-1))
126 primary_maximum[primary_maximum < 0] = 0
127 primary_minimum = np.array(np.abs(np.min(response[..., :num_samples], axis=-1)))
128 primary_minimum[primary_minimum < 0] = 0
129 primary_abs = np.array(np.max([primary_maximum, primary_minimum], axis=0))
131 # Now compute the residual response. we need time steps and
132 # responses at two points
133 residual_times = [0, (zeros_to_append-1)/sample_rate]
134 residual_responses = np.concatenate([response[..., num_samples, np.newaxis],
135 response[..., num_samples+zeros_to_append-1, np.newaxis]],
136 axis=-1)
138 # Compute the peak response after the structure has stopped responding
139 peak_residual_responses = sdof_free_decay_peak_response(
140 residual_responses, residual_times, frequency, damping)
142 # Pull out the residual peaks
143 residual_maximum = np.array(np.max(peak_residual_responses, axis=-1))
144 residual_maximum[residual_maximum < 0] = 0
145 residual_minimum = np.array(np.abs(np.min(peak_residual_responses, axis=-1)))
146 residual_minimum[residual_minimum < 0] = 0
147 residual_abs = np.array(np.max([residual_maximum, residual_minimum], axis=0))
149 all_srs_this_frequency_line = np.array([
150 primary_maximum, primary_minimum, primary_abs,
151 residual_maximum, residual_minimum, residual_abs,
152 np.max([primary_maximum, residual_maximum], axis=0),
153 np.max([primary_minimum, residual_minimum], axis=0),
154 np.max([primary_abs, residual_abs], axis=0)])
156 if abs_type < 10:
157 srs_output[..., i_freq] = all_srs_this_frequency_line[abs_type-1]
158 else:
159 srs_output[..., i_freq] = np.moveaxis(all_srs_this_frequency_line, 0, -1)
160 return srs_output, frequencies
163def sdof_ramp_invariant_filter_weights(frequencies, sample_rate, damping,
164 spectrum_type):
165 """
166 Computes filter weights for SDOF resonators using a ramp-invariant filter.
168 The weights are used in conjunction with the function
169 `sdof_filter` and `srs` to calculate the shock response
170 spectrum of an acceleration time history using a ramp
171 invarient simulation.
173 Parameters
174 ----------
175 frequencies : np.ndarray
176 The frequencies at which filters should be computed
177 sample_rate : float
178 The sample rate of the measurement
179 damping : float
180 The fraction of critical damping for the system (e.g. 0.03, not 3 for 3%)
181 spectrum_type : int
182 See `spectrum_type` argument for compute_srs
184 Returns
185 -------
186 b : np.ndarray
187 Filter coefficients with shape (frequencies.shape,3)
188 a : np.ndarray
189 Filter coefficients with shape (frequencies.shape,3)
190 """
191 # Make sure it's a numpy array
192 frequencies = np.array(frequencies)
193 # Preallocate arrays
194 a = np.ones(frequencies.shape+(3,))
195 b = np.ones(frequencies.shape+(3,))
197 normalized_frequencies = 2*np.pi*frequencies/sample_rate
199 small_freq_indices = normalized_frequencies < 0.001
201 # Small frequencies
202 small_freqs = normalized_frequencies[small_freq_indices]
204 # 2*z*w + w*w*(1-2*z*z) where z is damping and w is normalized frequencies
205 a[small_freq_indices, 1] = (
206 2*damping*small_freqs
207 + small_freqs**2*(1-2*damping**2)
208 )
209 # -2*z*w + 2*z*z*w*w where z is damping and w is normalized frequencies
210 a[small_freq_indices, 2] = (
211 -2*damping*small_freqs
212 + 2*damping**2*small_freqs**2
213 )
215 if spectrum_type > 0:
216 # Absolute Acceleration Model
217 # z*w + (w*w)*( (1/6) - 2*z*z/3 );
218 b[small_freq_indices, 0] = (
219 damping*small_freqs
220 + small_freqs**2
221 * (1/6 - 2*damping**2/3)
222 )
223 # 2*w*w*(1-z*z)/3
224 b[small_freq_indices, 1] = (
225 2*small_freqs**2*(1-damping**2)/3
226 )
227 # -z*w + w*w*( (1/6) - 4*z*z/3 );
228 b[small_freq_indices, 2] = (
229 - damping*small_freqs
230 + small_freqs**2
231 * (1/6 - 4*damping**2/3)
232 )
233 else:
234 # Relative Displacement Model
235 # -w*w/6;
236 b[small_freq_indices, 0] = -small_freqs**2/6
237 # -2*w*w/3;
238 b[small_freq_indices, 1] = -2*small_freqs**2/3
239 # -w*w/6;
240 b[small_freq_indices, 2] = -small_freqs**2/6
242 # Large frequencies
243 large_freqs = normalized_frequencies[~small_freq_indices]
245 # Define some helper variables
246 sq = np.sqrt(1-damping**2)
247 e = np.exp(-damping*large_freqs)
248 wd = large_freqs*sq
249 sp = e*np.sin(wd)
250 fact = (2*damping**2 - 1)*sp/sq
251 c = e*np.cos(wd)
253 a[~small_freq_indices, 1] = 2*(1-c)
254 a[~small_freq_indices, 2] = -1 + e**2
256 if spectrum_type > 0:
257 # Absolute Acceleration Model
258 spwd = sp/wd
259 b[~small_freq_indices, 0] = 1-spwd
260 b[~small_freq_indices, 1] = 2*(spwd - c)
261 b[~small_freq_indices, 2] = e**2 - spwd
262 else:
263 # Relative Displacement Model
264 b[~small_freq_indices, 0] = -(2*damping*(c-1) + fact + large_freqs)/large_freqs
265 b[~small_freq_indices, 1] = -(-2*c*large_freqs + 2*damping*(1-e**2) - 2*fact)/large_freqs
266 b[~small_freq_indices, 2] = -(e**2*(large_freqs+2*damping) - 2*damping*c + fact)/large_freqs
268 return b, a
271def sdof_filter(b, a, signal, zi=None):
272 """
273 Applies a filter to simulate a single degree of freedom system
275 Parameters
276 ----------
277 b : np.ndarray
278 Size 3 array representing filter coefficients used by scipy.signal.lfilter
279 a : np.array
280 Size 3 array representing filter coefficients used by scipy.signal.lfilter
281 signal : np.ndarray
282 The signals that are to be filtered
283 zi : np.ndarray, optional
284 Optional initial state for the filters, having length
285 `max(len(a), len(b)) - 1`. If not specified, zero initial conditions
286 are assumed.
288 Returns
289 -------
290 filtered_signal : np.ndarray
291 The filtered signal
293 """
294 aa = np.array([a[0], a[1]-2, a[2] + 1])
295 return lfilter(b, aa, signal, axis=-1, zi=zi)
298def sdof_free_decay_peak_response(responses, times_at_responses, frequency, damping):
299 """
300 Calculates peak response of a freely-decaying sdof system.
302 The residual response is the peak response of the sdof system
303 as it decays after the input has ended, i.e., the input is zero.
304 The first two peaks in the decaying response will be the largest.
305 One peak will be negative and one positive. We don't know before
306 the calculation if the first and largest peak in amplitude will be
307 positive or negative.
309 Parameters
310 ----------
311 responses : np.array
312 A (...,2) shape array giving the response at two time steps
313 times_at_responses : np.array
314 A length-2 array giving the time values that the responses occur at
315 frequency : float
316 The frequency at which the computation is performed
317 damping : float
318 The fraction of critical damping for the system (e.g. 0.03, not 3 for 3%)
320 Returns
321 -------
322 response_peaks : np.ndarray
323 A shape (...,2) array giving the amplitude at the two peaks after the
324 response has decayed
326 Notes
327 -----
328 The response has the form
329 a(t) = exp(-zeta wn t)[z[0]sin(wd t) + z[1]cos(wd t)].
330 If I know a(t) at two values of time, t, I can calculate the constants
331 z[0] and z[1]. The general form of the response can then be solved for
332 the maximum by finding the time of the first maximum by setting the
333 derivative to zero. Then substituting the time of the maximum response
334 back into the general equation to find the maximum response.
335 The second peak will occur half a cycle later.
337 """
338 responses = np.array(responses)
339 times_at_responses = np.array(times_at_responses)
340 # Set up some initial helper variables
341 fd = frequency*np.sqrt(1-damping*damping)
342 wd = 2*np.pi*fd
343 wn = 2*np.pi*frequency
344 wdt = wd*times_at_responses
345 e = np.exp(-wn*damping*times_at_responses)
346 sd = np.sin(wdt)
347 cd = np.cos(wdt)
348 c = e[:, np.newaxis]*np.array((sd, cd)).T
349 z = np.linalg.solve(c, responses[..., np.newaxis])[..., 0]
351 # the response has the form
352 # a(t) = exp(-zeta wn t)[z[0]sin(wd t) + z[1]cos(wd t)]
353 # the first maximum response will occur at the time tmax
354 tmax = np.array((1/wd) * np.arctan((wd*z[..., 0] - damping*wn*z[..., 1])/(wd*z[..., 1] - damping*wn*z[..., 0])))
355 tmax[tmax < 0] = tmax[tmax < 0] + np.pi/wd
356 tmax[tmax > np.pi/wd] = tmax[tmax > np.pi/wd] - np.pi/wd
358 # We can now plug it into the equation to find the maximum
359 response_peaks = []
360 response_peaks.append(np.exp(-damping*wn*tmax)*(z[..., 0]*np.sin(wd*tmax) + z[..., 1]*np.cos(wd*tmax)))
361 tmin = tmax + np.pi/wd
362 response_peaks.append(np.exp(-damping*wn*tmin)*(z[..., 0]*np.sin(wd*tmin) + z[..., 1]*np.cos(wd*tmin)))
363 return np.moveaxis(response_peaks, 0, -1)
366def octspace(low, high, points_per_octave):
367 """
368 Constructs octave spacing between low and high values
370 Parameters
371 ----------
372 low : float
373 Starting value for the spacing
374 high : float
375 Upper value for the spacing
376 points_per_octave : int
377 Number of points per octave
379 Returns
380 -------
381 octave_points : np.ndarray
382 Octave-spaced points
383 """
384 num_octaves = np.log2(high/low)
385 num_steps = np.ceil(num_octaves*points_per_octave)
386 point_indices = np.arange(num_steps+1)
387 log_points = np.log2(low) + num_octaves/num_steps*point_indices
388 points = 2**log_points
389 return points
392def sum_decayed_sines(sample_rate, block_size,
393 sine_frequencies=None, sine_tone_range=None, sine_tone_per_octave=None,
394 sine_amplitudes=None, sine_decays=None, sine_delays=None,
395 required_srs=None, srs_breakpoints=None,
396 srs_damping=0.05, srs_type=9,
397 compensation_frequency=None, compensation_decay=0.95,
398 # Paramters for the iteration
399 number_of_iterations=3, convergence=0.8,
400 error_tolerance=0.05,
401 tau=None, num_time_constants=None, decay_resolution=None,
402 scale_factor=1.02,
403 acceleration_factor=1.0,
404 plot_results=False, srs_frequencies=None,
405 ignore_compensation_pulse = False,
406 verbose=False
407 ):
408 """
409 Generate a Sum of Decayed Sines signal given an SRS.
411 Note that there are many approaches to do this, with many optional arguments
412 so please read the documentation carefully to understand which arguments
413 must be passed to the function.
415 Parameters
416 ----------
417 sample_rate : float
418 The sample rate of the generated signal.
419 block_size : int
420 The number of samples in the generated signal.
421 sine_frequencies : np.ndarray, optional
422 The frequencies of the sine tones. If this argument is not specified,
423 then the `sine_tone_range` argument must be specified.
424 sine_tone_range : np.ndarray, optional
425 A length-2 array containing the minimum and maximum sine tone to
426 generate. If this is argument is not specified, then the
427 `sine_frequencies` argument must be specified instead.
428 sine_tone_per_octave : int, optional
429 The number of sine tones per octave. If not specified along with
430 `sine_tone_range`, then a default value of 4 will be used if the
431 `srs_damping` is >= 0.05. Otherwise, the formula of
432 `sine_tone_per_octave = 9 - srs_damping*100` will be used.
433 sine_amplitudes : np.ndarray, optional
434 The initial amplitude of the sine tones used in the optimization. If
435 not specified, they will be set to the value of the SRS at each frequency
436 divided by the quality factor of the SRS.
437 sine_decays : np.ndarray, optional
438 An array of decay value time constants (often represented by variable
439 tau). Tau is the time for the amplitude of motion to decay 63% defined
440 by the equation `1/(2*np.pi*freq*zeta)` where `freq` is the frequency
441 of the sine tone and `zeta` is the fraction of critical damping.
442 If not specified, then either the `tau` or `num_time_constants`
443 arguments must be specified instead.
444 sine_delays : np.ndarray, optional
445 An array of delay values for the sine components. If not specified,
446 all tones will have zero delay.
447 required_srs : np.ndarray, optional
448 The SRS to match defined at each of the `sine_frequencies`. If this
449 argument is not passed, then the `srs_breakpoints` argument must be
450 passed instead. The default is None.
451 srs_breakpoints : np.ndarray, optional
452 A numpy array with shape `(n,2)` where the first column is the
453 frequencies at which each breakpoint occurs, and the second column is
454 the value of the SRS at that breakpoint. SRS values at the
455 `sine_frequencies` will be interpolated in a log-log sense from this
456 breakpoint array. If this argument is not specified, then the
457 `required_srs` must be passed instead.
458 srs_damping : float, optional
459 Fraction of critical damping to use in the SRS calculation (e.g. you
460 should specify 0.03 to represent 3%, not 3). The default is 0.03.
461 srs_type : int
462 The type of spectrum desired:
463 If `srs_type` > 0 (pos) then the SRS will be a base
464 acceleration-absolute acceleration model
465 If `srs_type` < 0 (neg) then the SRS will be a base acceleration-relative
466 displacement model (expressed in equivalent static acceleration units).
467 If abs(`srs_type`) is:
468 1--positive primary, 2--negative primary, 3--absolute maximum primary
469 4--positive residual, 5--negative residual, 6--absolute maximum residual
470 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
471 9 -- maximax, the largest absolute value of 1-8
472 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
473 compensation_frequency : float
474 The frequency of the compensation pulse. If not specified, it will be
475 set to 1/3 of the lowest sine tone
476 compensation_decay : float
477 The decay value for the compensation pulse. If not specified, it will
478 be set to 0.95.
479 number_of_iterations : int, optional
480 The number of iterations to perform. At least two iterations should be
481 performed. 3 iterations is preferred, and will be used if this argument
482 is not specified.
483 convergence : float, optional
484 The fraction of the error corrected each iteration. The default is 0.8.
485 error_tolerance : float, optional
486 Allowable relative error in the SRS. The default is 0.05.
487 tau : float, optional
488 If a floating point number is passed, then this will be used for the
489 `sine_decay` values. Alternatively, a dictionary can be passed with
490 the keys containing a length-2 tuple specifying the minimum and maximum
491 frequency range, and the value specifying the value of `tau` within that
492 frequency range. If this latter approach is used, all `sine_frequencies`
493 must be contained within a frequency range. If this argument is not
494 specified, then either `sine_decays` or `num_time_constants` must be
495 specified instead.
496 num_time_constants : int, optional
497 If an integer is passed, then this will be used to set the `sine_decay`
498 values by ensuring the specified number of time constants occur in the
499 `block_size`. Alternatively, a dictionary can be passed with the keys
500 containing a length-2 tuple specifying the minimum and maximum
501 frequency range, and the value specifying the value of
502 `num_time_constants` over that frequency range. If this latter approach
503 is used, all `sine_frequencies` must be contained within a frequency
504 range. If this argument is not specified, then either `sine_decays` or
505 `tau` must be specified instead.
506 decay_resolution : float, optional
507 A scalar identifying the resolution of the fractional decay rate
508 (often known by the variable `zeta`). The decay parameters will be
509 rounded to this value. The default is to not round.
510 scale_factor : float, optional
511 A scaling applied to the sine tone amplitudes so the achieved SRS better
512 fits the specified SRS, rather than just touching it. The default is 1.02.
513 acceleration_factor : float, optional
514 Optional scale factor to convert acceleration into velocity and
515 displacement. For example, if sine amplitudes are in G and displacement
516 is desired in inches, the acceleration factor should be set to 386.089.
517 If sine amplitudes are in G and displacement is desired in meters, the
518 acceleration factor should be set to 9.80665. The default is 1, which
519 assumes consistent units (e.g. acceleration in m/s^2, velocity in m/s,
520 displacement in m).
521 plot_results : bool, optional
522 If True, a figure will be plotted showing the acceleration, velocity,
523 and displacement signals, as well as the desired and achieved SRS.
524 srs_frequencies : np.ndarray, optional
525 If specified, these frequencies will be used to compute the SRS that
526 will be plotted when the `plot_results` value is `True`.
527 ignore_compensation_pulse : bool, optional
528 If True, the compensation pulse will be ignored.
529 verbose : True, optional
530 If True, additional diagnostics will be printed to the console.
532 Returns
533 -------
534 acceleration_signal : ndarray
535 The acceleration signal that satisfies the SRS.
536 velocity_signal : ndarray
537 The velocity of the acceleration signal.
538 displacement_signal : ndarray
539 The displacement of the acceleration signal.
540 sine_frequencies : ndarray
541 An array of frequencies for each sine tone, including the compensation
542 pulse
543 sine_amplitudes : ndarray
544 An array of amplitudes for each sine tone, including the compensation
545 pulse
546 sine_decays : ndarray
547 An array of decay values for each sine tone, including the compensation
548 pulse
549 sine_delays : ndarray
550 An array of delay values for each sine tone, including the compensation
551 pulse
552 fig : matplotlib.Figure
553 A reference to the plotted figure. Only returned if `plot_results` is
554 `True`.
555 ax : matplotlib.Axes
556 A reference to the plotted axes. Only returned if `plot_results` is
557 `True`.
559 """
560 # Handle the sine tone frequencies
561 if sine_frequencies is None and sine_tone_range is None:
562 raise ValueError('Either `sine_frequencies` or `sine_tone_range` must be specified')
563 if sine_frequencies is not None and sine_tone_range is not None:
564 raise ValueError('`sine_frequencies` can not be specified simultaneously with `sine_tone_range`')
565 if sine_frequencies is None:
566 # Create sine tones
567 if sine_tone_per_octave is None:
568 sine_tone_per_octave = int(np.floor(9-srs_damping*100))
569 sine_frequencies = octspace(sine_tone_range[0], sine_tone_range[1],
570 sine_tone_per_octave)
571 # Now set up the SRS
572 if required_srs is None and srs_breakpoints is None:
573 raise ValueError('Either `required_srs` or `srs_breakpoints` must be specified')
574 if required_srs is not None and srs_breakpoints is not None:
575 raise ValueError('`required_srs` can not be specified simultaneously with `srs_breakpoints`')
576 if required_srs is None:
577 required_srs = loginterp(sine_frequencies,
578 srs_breakpoints[:, 0],
579 srs_breakpoints[:, 1])
580 if sine_amplitudes is None:
581 srs_amplitudes = required_srs.copy()
582 srs_amplitudes[np.arange(srs_amplitudes.size) % 2 == 0] *= -1
583 quality = 1/(2*srs_damping)
584 srs_amplitudes /= quality
585 sine_amplitudes = srs_amplitudes
587 if sine_delays is None:
588 sine_delays = np.zeros(sine_frequencies.size)
590 if compensation_frequency is None:
591 compensation_frequency = np.min(sine_frequencies)/3
593 # Set up decay terms
594 decay_terms_specified = 0
595 if sine_decays is not None:
596 decay_terms_specified += 1
597 if tau is not None:
598 decay_terms_specified += 1
599 if num_time_constants is not None:
600 decay_terms_specified += 1
601 if decay_terms_specified == 0:
602 raise ValueError('One of `sine_decays`, `tau`, or `num_time_constants` must be specified')
603 if decay_terms_specified > 1:
604 raise ValueError('Only one of `sine_decays`, `tau`, or `num_time_constants` can be specified')
606 # Now check and see which is defined
607 if num_time_constants is not None:
608 period = block_size / sample_rate
609 if isinstance(num_time_constants, dict):
610 tau = {}
611 for freq_range, num_time_constant in num_time_constants.items():
612 tau[freq_range] = period / num_time_constant
613 else:
614 tau = period/num_time_constants
615 if tau is not None:
616 sine_decays = []
617 for freq in sine_frequencies:
618 if isinstance(tau, dict):
619 this_decay = None
620 for freq_range, this_tau in tau.items():
621 if freq_range[0] <= freq <= freq_range[1]:
622 this_decay = 1/(2*np.pi*freq*this_tau)
623 break
624 if this_decay is None:
625 raise ValueError('No frequency range matching frequency {:} was found in the specified decay parameters.'.format(freq))
626 sine_decays.append(this_decay)
627 else:
628 sine_decays.append(1/(2*np.pi*freq*tau))
629 sine_decays = np.array(sine_decays)
630 # Otherwise we just keep the specified sine_decays
632 # Now handle the minimum resolution on the decay
633 if decay_resolution is not None:
634 sine_decays = decay_resolution*np.round(sine_decays/decay_resolution)
636 if compensation_frequency is None:
637 compensation_frequency = np.min(sine_frequencies)/3
639 # Now we can actually run the generation process
640 (acceleration_signal, used_frequencies, used_amplitudes, used_decays,
641 used_delays, used_comp_frequency, used_comp_amplitude, used_comp_decay,
642 used_comp_delay) = _sum_decayed_sines(
643 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
644 scale_factor*required_srs, compensation_frequency, compensation_decay,
645 sample_rate, block_size, srs_damping, srs_type, number_of_iterations,
646 convergence, error_tolerance, ignore_compensation_pulse, verbose)
648 all_frequencies = np.concatenate((used_frequencies,
649 [used_comp_frequency]))
650 all_amplitudes = np.concatenate((used_amplitudes,
651 [used_comp_amplitude]))
652 all_decays = np.concatenate((used_decays,
653 [used_comp_decay]))
654 all_delays = np.concatenate((used_delays,
655 [used_comp_delay]))
656 # Now compute displacement and velocity
657 velocity_signal, displacement_signal = sum_decayed_sines_displacement_velocity(
658 all_frequencies, all_amplitudes, all_decays, all_delays, sample_rate,
659 block_size, acceleration_factor)
661 return_vals = (acceleration_signal, velocity_signal, displacement_signal,
662 all_frequencies, all_amplitudes, all_decays, all_delays)
664 # Plot the results
665 if plot_results:
666 fig, ax = plt.subplots(2, 2, figsize=(8, 6))
667 times = np.arange(block_size)/sample_rate
668 ax[0, 0].plot(times, acceleration_signal)
669 ax[0, 0].set_ylabel('Acceleration')
670 ax[0, 0].set_xlabel('Time (s)')
671 ax[0, 1].plot(times, velocity_signal)
672 ax[0, 1].set_ylabel('Velocity')
673 ax[0, 1].set_xlabel('Time (s)')
674 ax[1, 0].plot(times, displacement_signal)
675 ax[1, 0].set_ylabel('Displacement')
676 ax[1, 0].set_xlabel('Time (s)')
677 # Compute SRS
678 if srs_frequencies is None:
679 srs_frequencies = sine_frequencies
680 this_srs, this_frequencies = srs(
681 acceleration_signal, 1/sample_rate, srs_frequencies,
682 srs_damping, srs_type)
683 if srs_breakpoints is None:
684 srs_abscissa = sine_frequencies
685 srs_ordinate = required_srs
686 else:
687 srs_abscissa, srs_ordinate = srs_breakpoints.T
688 ax[1, 1].plot(srs_abscissa, srs_ordinate, 'k--')
689 ax[1, 1].plot(this_frequencies, this_srs)
690 ax[1, 1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
691 ax[1, 1].set_xlabel('Frequency (Hz)')
692 ax[1, 1].set_yscale('log')
693 ax[1, 1].set_xscale('log')
694 fig.tight_layout()
695 return_vals += (fig, ax)
696 ax[1, 1].legend(('Reference', 'Decayed Sine'))
697 return return_vals
700def _sum_decayed_sines(sine_frequencies, sine_amplitudes,
701 sine_decays, sine_delays,
702 required_srs,
703 compensation_frequency, compensation_decay,
704 sample_rate, block_size,
705 srs_damping, srs_type,
706 number_of_iterations=3, convergence=0.8,
707 error_tolerance=0.05, ignore_compensation_pulse = False,
708 verbose=False):
709 """
710 Optimizes the amplitudes of sums of decayed sines
712 Parameters
713 ----------
714 sine_frequencies : ndarray
715 An array of frequencies for each sine tone
716 sine_amplitudes : ndarray
717 An array of amplitudes for each sine tone
718 sine_decays : ndarray
719 An array of decay values for each sine tone
720 sine_delays : ndarray
721 An array of delay values for each sine tone
722 required_srs : ndarray
723 An array of SRS values for each frequency in `sine_frequencies`
724 compensation_frequency : float
725 The frequency of the compensation pulse
726 compensation_decay : float
727 The decay value for the compensation pulse
728 sample_rate : float
729 The sample rate of the signal
730 block_size : int
731 The number of samples in the signal
732 srs_damping : float
733 Fraction of critical damping to use in the SRS calculation (e.g. you
734 should specify 0.03 to represent 3%, not 3). The default is 0.03.
735 srs_type : int
736 The type of spectrum desired:
737 If `spectrum_type` > 0 (pos) then the SRS will be a base
738 acceleration-absolute acceleration model
739 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
740 displacement model (expressed in equivalent static acceleration units).
741 If abs(`spectrum_type`) is:
742 1--positive primary, 2--negative primary, 3--absolute maximum primary
743 4--positive residual, 5--negative residual, 6--absolute maximum residual
744 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
745 9 -- maximax, the largest absolute value of 1-8
746 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
747 number_of_iterations : int
748 The number of iterations that will be performed. The default is 3.
749 convergence : float, optional
750 The fraction of the error corrected each iteration. The default is 0.8.
751 error_tolerance : float, optional
752 Allowable relative error in the SRS. The default is 0.05.
753 ignore_compensation_pulse : bool, optional
754 If True, ignores the compensation pulse. Default is false.
755 verbose : bool, optional
756 If True, information on the interations will be provided. The default
757 is False.
759 Returns
760 -------
761 this_signal : np.ndarray
762 A numpy array containing the generated signal
763 sine_frequencies : ndarray
764 An array of frequencies for each sine tone
765 sine_amplitudes : ndarray
766 An array of amplitudes for each sine tone
767 sine_decays : ndarray
768 An array of decay values for each sine tone
769 sine_delays : ndarray
770 An array of delay values for each sine tone
771 this_compensation_frequency : float
772 The frequency of the compensation term.
773 this_compensation_amplitude : float
774 The amplitude of the compensation term.
775 this_compensation_decay : float
776 The decay constant of the compensation term
777 this_compensation_delay : float
778 The delay constant of the compenstation term
779 """
780 for i in range(number_of_iterations):
781 if verbose:
782 print('Starting Iteration {:}'.format(i+1))
783 # Compute updated sine amplitudes and compensation terms
784 this_sine_amplitudes, this_compensation_amplitude, this_compensation_delay = (
785 _sum_decayed_sines_single_iteration(
786 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
787 required_srs, compensation_frequency, compensation_decay,
788 sample_rate, block_size, srs_damping, srs_type,
789 number_of_iterations=10, convergence=convergence,
790 error_tolerance=error_tolerance, ignore_compensation_pulse=ignore_compensation_pulse,
791 verbose=verbose))
792 # get the SRS error at each frequency term by first computing the signal
793 (this_signal, this_compensation_frequency, this_compensation_amplitude,
794 this_compensation_decay, this_compensation_delay) = sum_decayed_sines_reconstruction_with_compensation(
795 sine_frequencies, this_sine_amplitudes, sine_decays, sine_delays,
796 compensation_frequency, compensation_decay, sample_rate, block_size, ignore_compensation_pulse=ignore_compensation_pulse)
797 # Then computing the SRS of the signal
798 this_srs = srs(this_signal, 1/sample_rate, sine_frequencies, srs_damping,
799 srs_type)[0]
800 srs_error = (this_srs - required_srs)/required_srs
801 sine_amplitudes = this_sine_amplitudes
802 if verbose:
803 print('Sine Table after Iterations:')
804 print('{:>8s}, {:>10s}, {:>10s}, {:>10s}, {:>10s}, {:>10s}'.format(
805 'Sine', 'Frequency', 'Amplitude', 'SRS Val', 'SRS Req', 'Error'))
806 for i, (freq, amp, srs_val, srs_req, srs_err) in enumerate(zip(
807 sine_frequencies, sine_amplitudes, this_srs, required_srs,
808 srs_error)):
809 print('{:>8s}, {:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}'.format(
810 str(i), freq, amp, srs_val, srs_req, srs_err))
811 print('Compensating Pulse:')
812 print('{:>10s}, {:>10s}, {:>10s}, {:>10s}'.format(
813 'Frequency', 'Amplitude', 'Decay', 'Delay'))
814 print('{:>10.4f}, {:>10.4f}, {:>10.4f}, {:>10.4f}'.format(
815 this_compensation_frequency, this_compensation_amplitude,
816 this_compensation_decay, this_compensation_delay))
817 return (this_signal, sine_frequencies, sine_amplitudes, sine_decays,
818 sine_delays, this_compensation_frequency, this_compensation_amplitude,
819 this_compensation_decay, this_compensation_delay)
822def _sum_decayed_sines_single_iteration(sine_frequencies, sine_amplitudes,
823 sine_decays, sine_delays,
824 required_srs, compensation_frequency,
825 compensation_decay,
826 sample_rate, block_size,
827 damping_srs, srs_type,
828 number_of_iterations=10, convergence=0.8,
829 error_tolerance=0.05,
830 ignore_compensation_pulse=False,
831 verbose=False):
832 """
833 Iterates on amplitudes of decayed sine waves to match a prescribed SRS
835 Parameters
836 ----------
837 sine_frequencies : ndarray
838 An array of frequencies for each sine tone
839 sine_amplitudes : ndarray
840 An array of amplitudes for each sine tone
841 sine_decays : ndarray
842 An array of decay values for each sine tone
843 sine_delays : ndarray
844 An array of delay values for each sine tone
845 required_srs : ndarray
846 An array of SRS values for each frequency in `sine_frequencies`
847 compensation_frequency : float
848 The frequency of the compensation pulse
849 compensation_decay : float
850 The decay value for the compensation pulse
851 sample_rate : float
852 The sample rate of the signal
853 block_size : int
854 The number of samples in the signal
855 damping_srs : float
856 Fraction of critical damping to use in the SRS calculation (e.g. you
857 should specify 0.03 to represent 3%, not 3). The default is 0.03.
858 srs_type : int
859 The type of spectrum desired:
860 If `spectrum_type` > 0 (pos) then the SRS will be a base
861 acceleration-absolute acceleration model
862 If `spectrum_type` < 0 (neg) then the SRS will be a base acceleration-relative
863 displacement model (expressed in equivalent static acceleration units).
864 If abs(`spectrum_type`) is:
865 1--positive primary, 2--negative primary, 3--absolute maximum primary
866 4--positive residual, 5--negative residual, 6--absolute maximum residual
867 7--largest of 1&4, maximum positive, 8--largest of 2&5, maximum negative
868 9 -- maximax, the largest absolute value of 1-8
869 10 -- returns a matrix s(9,length(fn)) with all the types 1-9.
870 number_of_iterations : int
871 Maximum number of iterations that can be performed at each frequency line.
872 Default is 10.
873 convergence : float, optional
874 The fraction of the error corrected each iteration. The default is 0.8.
875 error_tolerance : float, optional
876 Allowable relative error in the SRS. The default is 0.05.
877 ignore_compensation_pulse : bool, optional
878 If True, do not use a compensation pulse. Default is false.
879 verbose : bool, optional
880 If True, information on the interations will be provided. The default
881 is False.
883 Raises
884 ------
885 ValueError
886 If compensation delay is required to be too long.
888 Returns
889 -------
890 sine_amplitudes : np.ndarray
891 An array the same size as the input `sine_amplitudes` array with updated
892 amplitudes.
893 compensation_amplitude : float
894 Amplitude of the compensation pulse
895 compensation_delay : float
896 Delay of the compensation pulse.
897 """
898 # Define some helper variables
899 Amax = np.max(np.abs(sine_amplitudes))
900 # Reduce the error tolerance so there's a bit of room for round-off
901 error_tolerance = error_tolerance*0.9
902 # Get filter weights for SRS calculations
903 b, a = sdof_ramp_invariant_filter_weights(sine_frequencies, sample_rate,
904 damping_srs, srs_type)
905 # Copy the arrays so we don't overwrite anything
906 sine_frequencies = np.array(sine_frequencies).copy()
907 sine_amplitudes = np.array(sine_amplitudes).copy()
908 sine_decays = np.array(sine_decays).copy()
909 sine_delays = np.array(sine_delays).copy()
910 # Iterate one frequency at a time
911 for i, (frequency, amplitude, decay, delay, desired_srs) in enumerate(zip(
912 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
913 required_srs)):
914 # if i == 10:
915 # break
916 srs_error = float('inf')
917 iteration_count = 1
918 increment = 0.1
919 # Get the pulse without this frequency line
920 other_sine_frequencies = np.concatenate((sine_frequencies[:i], sine_frequencies[i+1:]))
921 other_sine_amplitudes = np.concatenate((sine_amplitudes[:i], sine_amplitudes[i+1:]))
922 other_sine_decays = np.concatenate((sine_decays[:i], sine_decays[i+1:]))
923 other_sine_delays = np.concatenate((sine_delays[:i], sine_delays[i+1:]))
924 other_pulse = sum_decayed_sines_reconstruction(
925 other_sine_frequencies, other_sine_amplitudes, other_sine_decays,
926 other_sine_delays, sample_rate, block_size)
927 # Now iterate at the single frequency until convergence
928 while abs(srs_error) > error_tolerance:
929 sine_amplitudes[i] = amplitude
930 # Find the pulse at this frequency line
931 this_pulse = sum_decayed_sines_reconstruction(
932 frequency, amplitude, decay, delay, sample_rate, block_size)
933 # Get the compensating pulse
934 if not ignore_compensation_pulse:
935 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
936 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
937 else:
938 compensation_amplitude = 0
939 compensation_delay = 0
940 # Find the number of samples to shift the waveform
941 num_shift = -int(np.floor(compensation_delay*sample_rate))
942 if abs(num_shift) >= block_size:
943 raise ValueError('The number of samples in the compensation delay ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
944 + ' The entire pulse will consist of part of the compensation pulse.'
945 + ' Please increase the block_size or compensation frequency.')
946 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
947 # Find compensating time history
948 compensation_pulse = sum_decayed_sines_reconstruction(
949 compensation_frequency, compensation_amplitude, compensation_decay,
950 compensation_delay_corrected, sample_rate, block_size)
951 # Build the composite waveform, need to shift the signals to align
952 if num_shift >= 0:
953 composite_pulse = (
954 compensation_pulse
955 + np.concatenate((
956 np.zeros(num_shift),
957 (other_pulse + this_pulse)[:block_size-num_shift])))
958 else:
959 composite_pulse = (
960 other_pulse + this_pulse +
961 + np.concatenate((
962 np.zeros(-num_shift),
963 (compensation_pulse)[:block_size+num_shift])))
964 # Find the SRS at the current frequency line
965 srs_prediction = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
966 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
967 srs_error = (srs_prediction - desired_srs)/desired_srs
968 if verbose:
969 print('Iteration at frequency {:}, {:0.4f}\n Iteration Count: {:}, SRS Error: {:0.4f}'.format(i, frequency, iteration_count, srs_error))
970 # Now we're going to compute the same thing again with a perturbed
971 # amplitude, this will allow us to compute the slope change at the
972 # current amplitude
973 amplitude_change = np.sign(srs_error)*increment*np.sign(amplitude)*Amax
974 # Now check and see if we need to modify the amplitude
975 if abs(srs_error) > error_tolerance:
976 if amplitude_change == 0: # perturb it a bit
977 amplitude_change = np.sign(srs_error)*np.sign(sine_amplitudes[i])*Amax*increment/10
978 new_amplitude = amplitude + amplitude_change
979 # Can't allow the sign of the amplitude to change
980 if np.sign(amplitude) != np.sign(new_amplitude):
981 new_amplitude *= -1
982 if amplitude == new_amplitude:
983 new_amplitude = amplitude * 1.01
984 sine_amplitudes[i] = new_amplitude
985 # Find new component at this frequency
986 this_pulse = sum_decayed_sines_reconstruction(
987 frequency, new_amplitude, decay, delay, sample_rate, block_size)
988 # Get the compensating pulse
989 if not ignore_compensation_pulse:
990 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
991 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
992 else:
993 compensation_amplitude = 0
994 compensation_delay = 0
995 # Find the number of samples to shift the waveform
996 num_shift = -int(np.floor(compensation_delay*sample_rate))
997 if abs(num_shift) >= block_size:
998 raise ValueError('The number of samples in the compensation delay ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
999 + ' The entire pulse will consist of part of the compensation pulse. '
1000 + 'Please increase the block_size or compensation frequency.')
1001 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
1002 # Find compensating time history
1003 compensation_pulse = sum_decayed_sines_reconstruction(
1004 compensation_frequency, compensation_amplitude, compensation_decay,
1005 compensation_delay_corrected, sample_rate, block_size)
1006 # Build the composite waveform, need to shift the signals to align
1007 if num_shift >= 0:
1008 composite_pulse = (
1009 compensation_pulse
1010 + np.concatenate((
1011 np.zeros(num_shift),
1012 (other_pulse + this_pulse)[:block_size-num_shift])))
1013 else:
1014 composite_pulse = (
1015 other_pulse + this_pulse +
1016 + np.concatenate((
1017 np.zeros(-num_shift),
1018 (compensation_pulse)[:block_size+num_shift])))
1019 # Find the SRS at the current frequency line
1020 srs_perturbed = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
1021 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
1022 # Get slope of correction
1023 correction_slope = (abs(amplitude)-abs(new_amplitude))/(srs_prediction-srs_perturbed)
1024 if correction_slope > 0:
1025 # this is what we want, it means we are changing in the right
1026 # direction
1027 new_amplitude = convergence*correction_slope*(desired_srs-srs_prediction) + abs(amplitude)
1028 # Never let the change be too large
1029 amplitude_check = max([abs(amplitude), Amax/10])
1030 if new_amplitude > 2*amplitude_check:
1031 new_amplitude = 2*amplitude_check
1032 # Never let it be less than zero
1033 if new_amplitude < 0:
1034 new_amplitude = np.abs(amplitude)/10
1035 # Move the previous amplitude into a different variable
1036 old_amplitude = amplitude
1037 # Never allow new amplitude to be zero
1038 if new_amplitude == 0:
1039 amplitude = 2*np.finfo(float).eps*np.sign(amplitude)
1040 else:
1041 amplitude = new_amplitude*np.sign(amplitude)
1042 # Store it for the next iteration
1043 sine_amplitudes[i] = amplitude
1044 # Find new component at this frequency
1045 this_pulse = sum_decayed_sines_reconstruction(
1046 frequency, amplitude, decay, delay, sample_rate, block_size)
1047 # Get the compensating pulse
1048 if not ignore_compensation_pulse:
1049 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1050 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
1051 else:
1052 compensation_amplitude = 0
1053 compensation_delay = 0
1054 # Find the number of samples to shift the waveform
1055 num_shift = -int(np.floor(compensation_delay*sample_rate))
1056 if abs(num_shift) >= block_size:
1057 raise ValueError('The number of samples in the compensation delay'
1058 + ' ({:}) is larger than the block_size ({:}).'.format(num_shift, block_size)
1059 + ' The entire pulse will consist of part of the compensation pulse. '
1060 + 'Please increase the block_size or compensation frequency.')
1061 compensation_delay_corrected = compensation_delay + num_shift/sample_rate
1062 # Find compensating time history
1063 compensation_pulse = sum_decayed_sines_reconstruction(
1064 compensation_frequency, compensation_amplitude, compensation_decay,
1065 compensation_delay_corrected, sample_rate, block_size)
1066 # Build the composite waveform, need to shift the signals to align
1067 if num_shift >= 0:
1068 composite_pulse = (
1069 compensation_pulse
1070 + np.concatenate((
1071 np.zeros(num_shift),
1072 (other_pulse + this_pulse)[:block_size-num_shift])))
1073 else:
1074 composite_pulse = (
1075 other_pulse + this_pulse +
1076 + np.concatenate((
1077 np.zeros(-num_shift),
1078 (compensation_pulse)[:block_size+num_shift])))
1079 # Find the SRS at the current frequency line
1080 srs_corrected = srs(composite_pulse, 1/sample_rate, frequency, damping_srs,
1081 srs_type, b[i], a[i])[0][0] # only get the SRS and there should only be one value due to one frequency line
1082 srs_error = (srs_corrected - desired_srs)/desired_srs
1083 # If the SRS error is positive and the amplitude has already been reduced to zero, stop trying
1084 if srs_error > 0 and amplitude == 0:
1085 iteration_count = number_of_iterations + 1
1086 if verbose:
1087 print(' Old Amplitude: {:0.4f}, New Amplitude: {:0.4f}, Error: {:0.4f}'.format(old_amplitude, amplitude, srs_error))
1088 # Don't allow amplitude to change more than a factor of 10 in an iteration
1089 if np.abs(10*old_amplitude) < np.abs(amplitude) or np.abs(0.1*old_amplitude) > np.abs(amplitude):
1090 iteration_count = number_of_iterations + 1
1091 else:
1092 # The slope is negative, so we try a bigger increment
1093 increment *= 2
1094 amplitude_change = np.sign(srs_error)*increment*np.sign(amplitude)*Amax
1095 if verbose:
1096 print('Slope of correction was negative, trying a bigger increment')
1097 iteration_count += 1
1098 if iteration_count > number_of_iterations:
1099 print(' Warning: SRS did not converge for frequency {:}: {:0.4f}'.format(i, frequency))
1100 break
1101 if not ignore_compensation_pulse:
1102 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1103 sine_frequencies, sine_amplitudes, sine_decays, sine_delays, compensation_frequency, compensation_decay)
1104 else:
1105 compensation_amplitude = 0
1106 compensation_delay = 0
1107 return sine_amplitudes, compensation_amplitude, compensation_delay
1110def sum_decayed_sines_compensating_pulse_parameters(sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1111 compensation_frequency, compensation_decay):
1112 omegas = sine_frequencies*2*np.pi
1113 omega_comp = compensation_frequency*2*np.pi
1114 var = -np.sum(sine_amplitudes/(sine_frequencies*(sine_decays**2+1)))
1115 compensation_amplitude = compensation_frequency*(compensation_decay**2+1)*var
1116 var0 = (np.sum(sine_amplitudes*sine_delays/(omegas*(sine_decays**2+1))))/compensation_amplitude
1117 # 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.
1118 # var1 = (np.sum(2*sine_decays*sine_amplitudes/(omegas**2*(sine_decays**2 + 1)**2)))/compensation_amplitude
1119 var1 = (np.sum((2*sine_decays*sine_amplitudes)[np.newaxis, :]@np.linalg.pinv((omegas**2*(sine_decays**2 + 1)**2)[np.newaxis, :])))/compensation_amplitude
1120 var2 = 2*compensation_decay/(omega_comp*omega_comp*(compensation_decay**2+1)**2)
1121 var3 = omega_comp*(compensation_decay**2+1)
1122 compensation_delay = -var3*(var2 + var1 + var0)
1123 return compensation_amplitude, compensation_delay
1126def sum_decayed_sines_reconstruction(sine_frequencies, sine_amplitudes,
1127 sine_decays, sine_delays, sample_rate,
1128 block_size):
1129 """
1130 Computes a sum of decayed sines signal from a set of frequencies, amplitudes,
1131 decays, and delays.
1133 Parameters
1134 ----------
1135 sine_frequencies : ndarray
1136 An array of frequencies for each sine tone
1137 sine_amplitudes : ndarray
1138 An array of amplitudes for each sine tone
1139 sine_decays : ndarray
1140 An array of decay values for each sine tone
1141 sine_delays : ndarray
1142 An array of delay values for each sine tone
1143 sample_rate : float
1144 The sample rate of the signal
1145 block_size : int
1146 The number of samples in the signal
1148 Returns
1149 -------
1150 ndarray
1151 A signal containing the sum of decayed sinusoids.
1153 """
1154 times = np.arange(block_size)/sample_rate
1155 # See if any delays are below zero and if so then shift all delays forward
1156 if np.any(sine_delays < 0):
1157 sine_delays = sine_delays - np.min(sine_delays)
1158 # Now go through and compute the sine tones
1159 omegas = 2*np.pi*sine_frequencies
1160 this_times = times[:, np.newaxis] - sine_delays
1161 response = sine_amplitudes*np.exp(-sine_decays*omegas*this_times)*np.sin(omegas*this_times)
1162 response[this_times < 0] = 0
1163 return np.sum(response, axis=-1)
1166def sum_decayed_sines_reconstruction_with_compensation(
1167 sine_frequencies, sine_amplitudes,
1168 sine_decays, sine_delays, compensation_frequency, compensation_decay,
1169 sample_rate, block_size, ignore_compensation_pulse = False):
1170 """
1173 Parameters
1174 ----------
1175 sine_frequencies : ndarray
1176 An array of frequencies for each sine tone
1177 sine_amplitudes : ndarray
1178 An array of amplitudes for each sine tone
1179 sine_decays : ndarray
1180 An array of decay values for each sine tone
1181 sine_delays : ndarray
1182 An array of delay values for each sine tone
1183 compensation_frequency : float
1184 The frequency of the compensation pulse
1185 compensation_decay : float
1186 The decay value for the compensation pulse
1187 sample_rate : float
1188 The sample rate of the signal
1189 block_size : int
1190 The number of samples in the signal
1191 ignore_compensation_pulse : bool, optional
1192 If True, ignores the compensation pulse
1194 Returns
1195 -------
1196 signal : ndarray
1197 A signal containing the sum of decayed sinusoids.
1198 compensation_frequency : float
1199 The frequency of the compensation pulse
1200 compensation_amplitude : float
1201 The amplitude value for the compensation pulse
1202 compensation_decay : float
1203 The decay value for the compensation pulse
1204 compensation_delay : float
1205 The delay value for the compensation pulse
1206 """
1207 sine_frequencies = np.array(sine_frequencies).flatten()
1208 sine_amplitudes = np.array(sine_amplitudes).flatten()
1209 sine_delays = np.array(sine_delays).flatten()
1210 sine_decays = np.array(sine_decays).flatten()
1211 if not ignore_compensation_pulse:
1212 compensation_amplitude, compensation_delay = sum_decayed_sines_compensating_pulse_parameters(
1213 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1214 compensation_frequency, compensation_decay)
1215 else:
1216 compensation_amplitude = 0
1217 compensation_delay = 0
1218 sine_frequencies = np.concatenate((sine_frequencies, [compensation_frequency]))
1219 sine_amplitudes = np.concatenate((sine_amplitudes, [compensation_amplitude]))
1220 sine_delays = np.concatenate((sine_delays, [compensation_delay]))
1221 sine_decays = np.concatenate((sine_decays, [compensation_decay]))
1222 signal = sum_decayed_sines_reconstruction(
1223 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1224 sample_rate, block_size)
1225 return (signal, compensation_frequency, compensation_amplitude,
1226 compensation_decay, compensation_delay)
1229def sum_decayed_sines_displacement_velocity(
1230 sine_frequencies, sine_amplitudes,
1231 sine_decays, sine_delays, sample_rate, block_size,
1232 acceleration_factor=1):
1233 """
1234 Creates velocity and displacement signals from acceleration sinusoids.
1236 Parameters
1237 ----------
1238 sine_frequencies : ndarray
1239 An array of frequencies for each sine tone
1240 sine_amplitudes : ndarray
1241 An array of amplitudes for each sine tone
1242 sine_decays : ndarray
1243 An array of decay values for each sine tone
1244 sine_delays : ndarray
1245 An array of delay values for each sine tone
1246 sample_rate : float
1247 The sample rate of the signal
1248 block_size : int
1249 The number of samples in the signal
1250 acceleration_factor : float, optional
1251 Optional scale factor to convert acceleration into velocity and
1252 displacement. For example, if sine amplitudes are in G and displacement
1253 is desired in inches, the acceleration factor should be set to 386.089.
1254 If sine amplitudes are in G and displacement is desired in meters, the
1255 acceleration factor should be set to 9.80665. The default is 1, which
1256 assumes consistent units (e.g. acceleration in m/s^2, velocity in m/s,
1257 displacement in m).
1259 Returns
1260 -------
1261 v : ndarray
1262 The velocity of the signal.
1263 d : ndarray
1264 The displacement of the signal.
1266 """
1267 # Make sure everything is a numpy array
1268 sine_frequencies = np.array(sine_frequencies).flatten()
1269 sine_amplitudes = np.array(sine_amplitudes).flatten()
1270 sine_delays = np.array(sine_delays).flatten()
1271 sine_decays = np.array(sine_decays).flatten()
1272 # Transform units
1273 sine_amplitudes = sine_amplitudes * acceleration_factor
1274 # See if any delays are below zero and if so then shift all delays forward
1275 if np.any(sine_delays < 0):
1276 sine_delays = sine_delays - np.min(sine_delays)
1277 f = sine_frequencies
1278 A = sine_amplitudes
1279 z = sine_decays
1280 tau = sine_delays
1281 x = np.zeros(block_size)
1282 tmp1 = x.copy()
1283 tmp2 = x.copy()
1284 tmp3 = x.copy()
1285 x1 = np.ones(x.shape)
1286 v = x.copy()
1287 d = x.copy()
1288 w = 2*np.pi*f
1289 w2 = w*w
1290 zw = z*w
1291 zp1 = z*z + 1
1292 zm1 = z*z - 1
1293 t = np.arange(block_size)/sample_rate
1294 for k in range(len(sine_frequencies)):
1295 indices = t-tau[k] > 0 # index's where vel and disp are evaluated
1296 Awz1 = A[k]/(w[k]*zp1[k])
1297 tmp1[indices] = Awz1*np.exp(-zw[k]*(t[indices]-tau[k]))
1298 tmp2[indices] = z[k]*np.sin(w[k]*(t[indices]-tau[k])) + np.cos(w[k]*(t[indices]-tau[k]))
1299 v[indices] = v[indices] - tmp1[indices]*tmp2[indices] + Awz1*x1[indices]
1300 Awz2 = A[k]/(w2[k]*zp1[k]*zp1[k])
1301 tmp1[indices] = Awz2*np.exp(-zw[k]*(t[indices]-tau[k]))
1302 tmp2[indices] = zm1[k]*np.sin(w[k]*(t[indices]-tau[k])) + 2*z[k]*np.cos(w[k]*(t[indices]-tau[k]))
1303 tmp3[indices] = Awz1*(t[indices]-tau[k])
1304 tmp4 = 2*z[k]*Awz2
1305 d[indices] = d[indices] + tmp1[indices]*tmp2[indices] + tmp3[indices] - tmp4
1306 return v, d
1309def loginterp(x, xp, fp):
1310 return 10**np.interp(np.log10(x), np.log10(xp), np.log10(fp))
1313def optimization_error_function(
1314 amplitude_lin_scales, sample_rate, block_size,
1315 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1316 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1317 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1318 b=None, a=None, frequency_index=None):
1319 if frequency_index is None:
1320 frequency_index = slice(None)
1321 # Apply the scale factors
1322 new_sine_amplitudes = sine_amplitudes.copy()
1323 new_sine_amplitudes[frequency_index] *= amplitude_lin_scales
1324 # Compute a new time history signal
1325 (input_signal, new_compensation_frequency, new_compensation_amplitude,
1326 new_compensation_decay, new_compensation_delay) = sum_decayed_sines_reconstruction_with_compensation(
1327 sine_frequencies, new_sine_amplitudes, sine_decays, sine_delays,
1328 compensation_frequency, compensation_decay, sample_rate, block_size)
1329 # Transform the signal to the control degrees of freedom
1330 control_responses = np.array([oaconvolve(control_irf, input_signal) for control_irf in control_irfs])[..., :block_size]
1331 # Compute SRSs at all frequencies
1332 predicted_control_srs, frequencies = srs(control_responses, 1/sample_rate, sine_frequencies[frequency_index],
1333 srs_damping, srs_type,
1334 None if b is None else b[frequency_index],
1335 None if a is None else a[frequency_index])
1336 # Compute error
1337 mean_control_error = np.mean(
1338 (control_weights[:, np.newaxis]
1339 * ((predicted_control_srs - control_srs[..., frequency_index])
1340 / control_srs[..., frequency_index])))
1341 if limit_srs is not None:
1342 limit_responses = np.array([oaconvolve(limit_irf, input_signal) for limit_irf in limit_irfs])[..., :block_size]
1343 predicted_limit_srs, frequencies = srs(limit_responses, 1/sample_rate, sine_frequencies[frequency_index],
1344 srs_damping, srs_type,
1345 None if b is None else b[frequency_index],
1346 None if a is None else a[frequency_index])
1347 max_limit_error = np.max(
1348 ((predicted_limit_srs - limit_srs[..., frequency_index])
1349 / limit_srs[..., frequency_index]))
1350 if max_limit_error >= 0 and mean_control_error >= 0:
1351 srs_error = np.max((mean_control_error, max_limit_error))
1352 elif max_limit_error <= 0 and mean_control_error <= 0:
1353 srs_error = np.max((mean_control_error, max_limit_error))
1354 elif max_limit_error >= 0 and mean_control_error <= 0:
1355 srs_error = max_limit_error
1356 elif max_limit_error <= 0 and mean_control_error >= 0:
1357 srs_error = mean_control_error
1358 else:
1359 limit_responses = None
1360 predicted_limit_srs = None
1361 srs_error = mean_control_error
1363 return (np.abs(srs_error), input_signal, control_responses,
1364 predicted_control_srs, limit_responses, predicted_limit_srs,
1365 new_sine_amplitudes, new_compensation_frequency,
1366 new_compensation_amplitude, new_compensation_decay,
1367 new_compensation_delay)
1370def optimization_callback(intermediate_result, rms_error_threshold=0.02,
1371 verbose=True):
1372 if verbose:
1373 print('Amplitude Scale: {:}, error is {:}'.format(intermediate_result.x.squeeze(),
1374 intermediate_result.fun))
1375 if rms_error_threshold is not None:
1376 if intermediate_result.fun < rms_error_threshold:
1377 raise StopIteration
1380def sum_decayed_sines_minimize(sample_rate, block_size,
1381 sine_frequencies=None, sine_tone_range=None, sine_tone_per_octave=None,
1382 sine_amplitudes=None, sine_decays=None, sine_delays=None,
1383 control_srs=None, control_breakpoints=None,
1384 srs_damping=0.05, srs_type=9,
1385 compensation_frequency=None, compensation_decay=0.95,
1386 # Parameters for defining decays
1387 tau=None, num_time_constants=None, decay_resolution=None,
1388 scale_factor=1.02,
1389 acceleration_factor=1.0,
1390 # Parameters for imposing limits
1391 limit_breakpoints=None, limit_transfer_functions=None,
1392 control_transfer_functions=None, control_weights=None,
1393 # Parameters for the optimizer
1394 minimize_iterations=1, rms_error_threshold=None,
1395 optimization_passes=3,
1396 plot_results=False, verbose=False
1397 ):
1398 # Handle the sine tone frequencies
1399 if sine_frequencies is None and sine_tone_range is None:
1400 raise ValueError('Either `sine_frequencies` or `sine_tone_range` must be specified')
1401 if sine_frequencies is not None and sine_tone_range is not None:
1402 raise ValueError('`sine_frequencies` can not be specified simultaneously with `sine_tone_range`')
1403 if sine_frequencies is None:
1404 # Create sine tones
1405 if sine_tone_per_octave is None:
1406 sine_tone_per_octave = int(np.floor(9-srs_damping*100))
1407 sine_frequencies = octspace(sine_tone_range[0], sine_tone_range[1],
1408 sine_tone_per_octave)
1409 # Now set up the SRS
1410 if control_srs is None and control_breakpoints is None:
1411 raise ValueError('Either `control_srs` or `control_breakpoints` must be specified')
1412 if control_srs is not None and control_breakpoints is not None:
1413 raise ValueError('`control_srs` can not be specified simultaneously with `control_breakpoints`')
1414 if control_srs is None:
1415 frequencies = control_breakpoints[:, 0]
1416 breakpoint_curves = control_breakpoints[:, 1:].T
1417 control_srs = np.array([loginterp(sine_frequencies,
1418 frequencies,
1419 breakpoint_curve) for breakpoint_curve in breakpoint_curves])
1420 else:
1421 control_srs = np.atleast_2d(control_srs)
1422 if control_transfer_functions is None:
1423 control_transfer_functions = np.ones(((control_srs.shape[0]-1)*2, block_size//2+1))
1424 tf_frequencies = np.fft.rfftfreq(control_transfer_functions.shape[-1]*2-1, 1/sample_rate)
1425 if sine_amplitudes is None:
1426 tf_at_frequencies = np.array([
1427 np.interp(sine_frequencies, tf_frequencies, np.abs(control_transfer_function))
1428 for control_transfer_function in control_transfer_functions])
1429 srs_amplitudes = np.array([nnls(ai[:, np.newaxis], bi)[0][0] for ai, bi in zip(tf_at_frequencies.T, control_srs.T)])
1430 srs_amplitudes[np.arange(srs_amplitudes.size) % 2 == 0] *= -1
1431 quality = 1/(2*srs_damping)
1432 srs_amplitudes /= quality
1433 sine_amplitudes = srs_amplitudes
1435 if sine_delays is None:
1436 sine_delays = np.zeros(sine_frequencies.size)
1438 if compensation_frequency is None:
1439 compensation_frequency = np.min(sine_frequencies)/3
1441 # Set up decay terms
1442 decay_terms_specified = 0
1443 if sine_decays is not None:
1444 decay_terms_specified += 1
1445 if tau is not None:
1446 decay_terms_specified += 1
1447 if num_time_constants is not None:
1448 decay_terms_specified += 1
1449 if decay_terms_specified == 0:
1450 raise ValueError('One of `sine_decays`, `tau`, or `num_time_constants` must be specified')
1451 if decay_terms_specified > 1:
1452 raise ValueError('Only one of `sine_decays`, `tau`, or `num_time_constants` can be specified')
1454 # Now check and see which is defined
1455 if num_time_constants is not None:
1456 period = block_size / sample_rate
1457 if isinstance(num_time_constants, dict):
1458 tau = {}
1459 for freq_range, num_time_constant in num_time_constants.items():
1460 tau[freq_range] = period / num_time_constant
1461 else:
1462 tau = period/num_time_constants
1463 if tau is not None:
1464 sine_decays = []
1465 for freq in sine_frequencies:
1466 if isinstance(tau, dict):
1467 this_decay = None
1468 for freq_range, this_tau in tau.items():
1469 if freq_range[0] <= freq <= freq_range[1]:
1470 this_decay = 1/(2*np.pi*freq*this_tau)
1471 break
1472 if this_decay is None:
1473 raise ValueError('No frequency range matching frequency {:} was found in the specified decay parameters.'.format(freq))
1474 sine_decays.append(this_decay)
1475 else:
1476 sine_decays.append(1/(2*np.pi*freq*tau))
1477 sine_decays = np.array(sine_decays)
1478 # Otherwise we just keep the specified sine_decays
1480 # Now handle the minimum resolution on the decay
1481 if decay_resolution is not None:
1482 sine_decays = decay_resolution*np.round(sine_decays/decay_resolution)
1484 if compensation_frequency is None:
1485 compensation_frequency = np.min(sine_frequencies)/3
1487 # Now set up limits and transfer functions
1488 if control_weights is None:
1489 control_weights = np.ones((control_srs.shape[0]))
1490 if limit_breakpoints is not None:
1491 frequencies = limit_breakpoints[:, 0]
1492 breakpoint_curves = limit_breakpoints[:, 1:].T
1493 limit_srs = np.array([loginterp(sine_frequencies,
1494 frequencies,
1495 breakpoint_curve) for breakpoint_curve in breakpoint_curves])
1496 else:
1497 limit_srs = None
1499 # Compute impulse responses
1500 control_irfs = np.fft.irfft(control_transfer_functions, axis=-1)
1501 if limit_transfer_functions is not None:
1502 limit_irfs = np.fft.irfft(limit_transfer_functions, axis=-1)
1503 else:
1504 limit_irfs = None
1505 b, a = sdof_ramp_invariant_filter_weights(sine_frequencies, sample_rate, srs_damping, srs_type)
1507 # Normalize control weights
1508 control_weights = control_weights/np.linalg.norm(control_weights)
1510 # Copy things so we don't overwrite
1511 sine_amplitudes = sine_amplitudes.copy()
1513 # Now we will iterate over all of the frequencies and compute the new
1514 # amplitudes
1515 for j in range(optimization_passes):
1516 for i in range(len(sine_frequencies)):
1517 if verbose:
1518 print('Pass {:}, Analyzing Frequency {:}: {:0.2f}'.format(j+1, i, sine_frequencies[i]))
1520 # We will now go through and optimize the frequency line
1521 def error_function(x):
1522 return optimization_error_function(
1523 x, sample_rate, block_size,
1524 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1525 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1526 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1527 b, a, frequency_index=[i])[0]
1529 def callback(intermediate_result):
1530 return optimization_callback(
1531 intermediate_result, rms_error_threshold, verbose)
1533 optimization_result = minimize(error_function,
1534 np.ones(1),
1535 method='Powell',
1536 bounds=[(0, np.inf)],
1537 callback=callback,
1538 options={'maxiter': minimize_iterations})
1540 # Populate the amplitudes with the updated result
1541 if verbose:
1542 print('Initial Amplitude: {:}, Updated Amplitude: {:}\n'.format(sine_amplitudes[i], sine_amplitudes[i]*optimization_result.x.squeeze()))
1543 sine_amplitudes[i] *= optimization_result.x.squeeze()
1545 (error, input_signal, control_responses, predicted_control_srs,
1546 limit_responses, predicted_limit_srs,
1547 sine_amplitudes, compensation_frequency,
1548 compensation_amplitude, compensation_decay,
1549 compensation_delay) = optimization_error_function(
1550 np.ones(sine_frequencies.size),
1551 sample_rate, block_size,
1552 sine_frequencies, sine_amplitudes, sine_decays, sine_delays,
1553 compensation_frequency, compensation_decay, control_irfs, limit_irfs,
1554 srs_damping, srs_type, control_srs, control_weights, limit_srs,
1555 b, a)
1557 all_frequencies = np.concatenate((sine_frequencies,
1558 [compensation_frequency]))
1559 all_amplitudes = np.concatenate((sine_amplitudes,
1560 [compensation_amplitude]))
1561 all_decays = np.concatenate((sine_decays,
1562 [compensation_decay]))
1563 all_delays = np.concatenate((sine_delays,
1564 [compensation_delay]))
1566 # Now compute displacement and velocity
1567 velocity_signal, displacement_signal = sum_decayed_sines_displacement_velocity(
1568 all_frequencies, all_amplitudes, all_decays, all_delays, sample_rate,
1569 block_size, acceleration_factor)
1571 return_vals = (input_signal, velocity_signal, displacement_signal,
1572 control_responses, predicted_control_srs,
1573 limit_responses, predicted_limit_srs,
1574 all_frequencies, all_amplitudes, all_decays, all_delays,
1575 )
1577 if plot_results:
1578 fig, ax = plt.subplots(2, 2, figsize=(8, 6))
1579 times = np.arange(block_size)/sample_rate
1580 ax[0, 0].plot(times, input_signal)
1581 ax[0, 0].set_ylabel('Acceleration')
1582 ax[0, 0].set_xlabel('Time (s)')
1583 ax[0, 1].plot(times, velocity_signal)
1584 ax[0, 1].set_ylabel('Velocity')
1585 ax[0, 1].set_xlabel('Time (s)')
1586 ax[1, 0].plot(times, displacement_signal)
1587 ax[1, 0].set_ylabel('Displacement')
1588 ax[1, 0].set_xlabel('Time (s)')
1589 # Compute SRS
1590 this_srs, this_frequencies = srs(
1591 input_signal, 1/sample_rate, sine_frequencies,
1592 srs_damping, srs_type)
1593 if control_breakpoints is None:
1594 srs_abscissa = sine_frequencies
1595 srs_ordinate = control_srs.T
1596 else:
1597 srs_abscissa = control_breakpoints[:, 0]
1598 srs_ordinate = control_breakpoints[:, 1:]
1599 ax[1, 1].plot(srs_abscissa, srs_ordinate, 'k--')
1600 ax[1, 1].plot(this_frequencies, this_srs)
1601 ax[1, 1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1602 ax[1, 1].set_xlabel('Frequency (Hz)')
1603 ax[1, 1].set_yscale('log')
1604 ax[1, 1].set_xscale('log')
1605 ax[1, 1].legend(('Reference', 'Decayed Sine'))
1606 fig.tight_layout()
1607 return_vals += (fig, ax)
1608 # Now plot all control channels
1609 for i, (predicted, control, response) in enumerate(zip(predicted_control_srs, control_srs, control_responses)):
1610 fig, ax = plt.subplots(1, 2, figsize=(8, 3))
1611 ax[0].set_title('Control Signal {:}'.format(i))
1612 ax[0].plot(times, response)
1613 ax[0].set_label('Acceleration')
1614 ax[0].set_xlabel('Time (s)')
1615 ax[1].set_title('Control SRS {:}'.format(i))
1616 ax[1].plot(sine_frequencies, control, 'k--')
1617 ax[1].plot(sine_frequencies, predicted)
1618 ax[1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1619 ax[1].set_xlabel('Frequency (Hz)')
1620 ax[1].set_yscale('log')
1621 ax[1].set_xscale('log')
1622 ax[1].legend(('Desired', 'Achieved'))
1623 fig.tight_layout()
1624 return_vals += (fig, ax)
1625 # Now plot all limit channels
1626 if limit_srs is not None:
1627 for i, (predicted, limit, response) in enumerate(zip(predicted_limit_srs, limit_srs, limit_responses)):
1628 fig, ax = plt.subplots(1, 2, figsize=(8, 3))
1629 ax[0].set_title('Limit Signal {:}'.format(i))
1630 ax[0].plot(times, response)
1631 ax[0].set_label('Acceleration')
1632 ax[0].set_xlabel('Time (s)')
1633 ax[1].set_title('Limit SRS {:}'.format(i))
1634 ax[1].plot(sine_frequencies, limit, 'k--')
1635 ax[1].plot(sine_frequencies, predicted)
1636 ax[1].set_ylabel('SRS ({:0.2f}% damping)'.format(srs_damping*100))
1637 ax[1].set_xlabel('Frequency (Hz)')
1638 ax[1].set_yscale('log')
1639 ax[1].set_xscale('log')
1640 ax[1].legend(('Desired', 'Achieved'))
1641 fig.tight_layout()
1642 return_vals += (fig, ax)
1644 return return_vals
1647@dataclass
1648class WindowedRandomMetrics:
1649 """
1650 Scalar metrics describing one generated realization.
1652 Parameters
1653 ----------
1654 lanl_rms_db_error : float
1655 LANL-style root-mean-square dB error between the achieved and target
1656 SRS over the matching frequencies.
1657 energy : float
1658 Signal energy computed from the acceleration time history.
1659 rea : float
1660 Root energy amplitude.
1661 peak_accel : float
1662 Maximum absolute acceleration in the compensated windowed history.
1663 peak_vel : float
1664 Maximum absolute velocity derived by integrating the compensated
1665 acceleration history.
1666 peak_disp : float
1667 Maximum absolute displacement derived by integrating the compensated
1668 acceleration history.
1670 Notes
1671 -----
1672 These metrics are derived from the compensated windowed realization and are
1673 intended for ranking or selecting among multiple random realizations.
1674 """
1675 lanl_rms_db_error: float
1676 energy: float
1677 rea: float
1678 peak_accel: float
1679 peak_vel: float
1680 peak_disp: float
1683@dataclass
1684class WindowedRandomResults:
1685 """
1686 Collection of generated windowed-random realizations and associated results.
1688 Parameters
1689 ----------
1690 time : ndarray
1691 Time vector with shape `(n_samples,)`.
1692 sample_rate : float
1693 Sampling rate in samples per second.
1694 gravity : float
1695 Acceleration-of-gravity conversion factor used when integrating
1696 acceleration to obtain velocity and displacement.
1697 window : ndarray
1698 Applied time-domain window with shape `(n_samples,)`.
1699 ref_srs : ndarray
1700 Reference SRS breakpoint array with shape `(n_breakpoints, 2)`. The
1701 first column contains frequency and the second contains spectral
1702 amplitude.
1703 srs_frequency : ndarray
1704 Frequency vector corresponding to `srs_matrix`, with shape `(n_freq,)`.
1705 srs_matrix : ndarray
1706 Computed SRS values for each realization, with shape
1707 `(n_freq, n_realizations)`.
1708 xc_matrix : ndarray
1709 Compensated windowed acceleration histories, with shape
1710 `(n_samples, n_realizations)`.
1711 rr_matrix : ndarray
1712 Stationary random acceleration histories prior to windowing, with shape
1713 `(n_samples, n_realizations)`.
1714 metrics : list of WindowedRandomMetrics
1715 Metrics for each realization.
1716 metadata : dict
1717 Dictionary of generation settings and other bookkeeping information.
1719 Notes
1720 -----
1721 Each column of `xc_matrix`, `rr_matrix`, and `srs_matrix` corresponds to
1722 one realization.
1723 """
1724 time: np.ndarray
1725 sample_rate: float
1726 gravity: float
1727 window: np.ndarray
1728 ref_srs: np.ndarray
1729 srs_frequency: np.ndarray
1730 srs_matrix: np.ndarray
1731 xc_matrix: np.ndarray
1732 rr_matrix: np.ndarray
1733 metrics: List[WindowedRandomMetrics]
1734 metadata: Dict[str, Any]
1736 @property
1737 def metrics_array(self) -> np.ndarray:
1738 """
1739 Return the realization metrics as a dense numeric array.
1741 Returns
1742 -------
1743 ndarray
1744 Array with shape `(n_realizations, 6)` containing columns:
1746 1. LANL RMS dB error
1747 2. energy
1748 3. REA
1749 4. peak acceleration
1750 5. peak velocity
1751 6. peak displacement
1752 """
1753 return np.array([
1754 [
1755 m.lanl_rms_db_error,
1756 m.energy,
1757 m.rea,
1758 m.peak_accel,
1759 m.peak_vel,
1760 m.peak_disp,
1761 ]
1762 for m in self.metrics
1763 ], dtype=float)
1766@dataclass
1767class SelectedRealization:
1768 """
1769 One realization selected from an ensemble.
1771 Parameters
1772 ----------
1773 index : int
1774 Zero-based index of the selected realization.
1775 time : ndarray
1776 Time vector with shape `(n_samples,)`.
1777 xc : ndarray
1778 Compensated windowed acceleration history with shape `(n_samples,)`.
1779 rr : ndarray
1780 Stationary random acceleration history with shape `(n_samples,)`.
1781 srs : ndarray
1782 Shock response spectrum of the selected realization with shape
1783 `(n_freq, 2)`, where the first column is frequency and the second is
1784 SRS amplitude.
1785 metrics : WindowedRandomMetrics
1786 Scalar metrics associated with the selected realization.
1787 """
1788 index: int
1789 time: np.ndarray
1790 xc: np.ndarray
1791 rr: np.ndarray
1792 srs: np.ndarray
1793 metrics: WindowedRandomMetrics
1796def breakpoint_spectrum(
1797 frequencies: np.ndarray,
1798 break_frequencies: np.ndarray,
1799 break_values: np.ndarray,
1800 floor: Optional[float] = None,
1801) -> np.ndarray:
1802 """
1803 Evaluate a breakpoint-defined spectrum using log-log interpolation.
1805 Parameters
1806 ----------
1807 frequencies : ndarray
1808 Frequencies at which the spectrum should be evaluated.
1809 break_frequencies : ndarray
1810 Monotonically increasing breakpoint frequencies.
1811 break_values : ndarray
1812 Spectrum amplitudes at `break_frequencies`.
1813 floor : float, optional
1814 Value used outside the breakpoint range. If not provided, the default
1815 is `min(break_values) / 1000`.
1817 Returns
1818 -------
1819 ndarray
1820 Spectrum amplitudes evaluated at `frequencies`.
1822 Raises
1823 ------
1824 ValueError
1825 If any breakpoint frequencies or breakpoint amplitudes are not
1826 positive.
1828 Notes
1829 -----
1830 This function reproduces the numerical intent of MATLAB ``bpspec``:
1831 straight lines on log-log axes between breakpoint pairs, with a constant
1832 floor outside the breakpoint range.
1833 """
1834 frequencies = np.asarray(frequencies, dtype=float)
1835 break_frequencies = np.asarray(break_frequencies, dtype=float)
1836 break_values = np.asarray(break_values, dtype=float)
1838 if np.any(break_frequencies <= 0):
1839 raise ValueError("break_frequencies must be positive")
1840 if np.any(break_values <= 0):
1841 raise ValueError("break_values must be positive")
1843 if floor is None:
1844 floor = np.min(break_values) / 1000.0
1846 y = np.full_like(frequencies, floor, dtype=float)
1847 inside = (frequencies >= break_frequencies[0]) & (frequencies <= break_frequencies[-1])
1849 if np.any(inside):
1850 f = interp1d(
1851 np.log(break_frequencies),
1852 np.log(break_values),
1853 kind="linear",
1854 bounds_error=False,
1855 fill_value=np.log(floor),
1856 )
1857 y[inside] = np.exp(f(np.log(frequencies[inside])))
1859 return y
1862def temporal_moments(x: np.ndarray, dt: float) -> Dict[str, float]:
1863 """
1864 Compute temporal moments of a sampled signal.
1866 Parameters
1867 ----------
1868 x : ndarray
1869 Input signal amplitudes with shape `(n_samples,)`.
1870 dt : float
1871 Sampling interval.
1873 Returns
1874 -------
1875 dict
1876 Dictionary containing:
1878 - ``"energy"`` : float
1879 Signal energy.
1880 - ``"tau"`` : float
1881 Time centroid.
1882 - ``"duration_rms"`` : float
1883 Root-mean-square duration.
1884 - ``"skewness_metric"`` : float
1885 Cube-root form of the third central temporal moment.
1886 - ``"kurtosis_metric"`` : float
1887 Fourth-root form of the fourth central temporal moment.
1888 - ``"rea"`` : float
1889 Root energy amplitude.
1891 Notes
1892 -----
1893 This function follows the intent of MATLAB ``gtempmom`` and is used to
1894 compute metrics for each generated acceleration history.
1895 """
1896 x = np.asarray(x, dtype=float).reshape(-1)
1897 t = dt * np.arange(len(x))
1898 z = x ** 2
1900 E = abs(dt * np.sum(z))
1901 tau = float(dt * np.dot(t, z) / E)
1903 t_centered = t - tau
1904 D2 = float((dt / E) * np.dot(t_centered ** 2, z))
1905 D = np.sqrt(max(D2, 0.0))
1907 S3 = float((dt / E) * np.dot(t_centered ** 3, z))
1908 K4 = float((dt / E) * np.dot(t_centered ** 4, z))
1910 S = np.sign(S3) * abs(S3) ** (1 / 3) if S3 != 0 else 0.0
1911 K = abs(K4) ** 0.25
1912 R = np.sqrt(E / D) if D > 0 else np.inf
1914 return {
1915 "energy": E,
1916 "tau": tau,
1917 "duration_rms": D,
1918 "skewness_metric": S,
1919 "kurtosis_metric": K,
1920 "rea": R,
1921 }
1924def lanl_rms_db_error(
1925 db_error: np.ndarray,
1926 frequency: np.ndarray,
1927 octave_resolution: float,
1928) -> Dict[str, np.ndarray | float]:
1929 """
1930 Compute LANL-style RMS dB error over octave-spaced frequency lines.
1932 Parameters
1933 ----------
1934 db_error : ndarray
1935 dB error values. May be one-dimensional with shape `(n_freq,)` or
1936 two-dimensional with shape `(n_freq, n_series)`.
1937 frequency : ndarray
1938 Frequency vector with shape `(n_freq,)`.
1939 octave_resolution : float
1940 Number of points per octave.
1942 Returns
1943 -------
1944 dict
1945 Dictionary containing:
1947 - ``"rms"`` : float
1948 Average RMS dB error across series.
1949 - ``"std"`` : float
1950 Average standard deviation of dB error across series.
1951 - ``"avg"`` : ndarray
1952 Mean dB error across series at each frequency.
1954 Raises
1955 ------
1956 ValueError
1957 If the first dimension of `db_error` does not match the length of
1958 `frequency`.
1960 Notes
1961 -----
1962 This follows the octave-band weighting used in the MATLAB
1963 ``lanl_rmsdberror`` routine.
1964 """
1965 db_error = np.asarray(db_error, dtype=float)
1966 frequency = np.asarray(frequency, dtype=float).reshape(-1)
1968 if db_error.ndim == 1:
1969 db_error = db_error[:, None]
1971 if db_error.shape[0] != len(frequency):
1972 raise ValueError("frequency length must match first dimension of db_error")
1974 upper = (2 ** (1 / (2 * octave_resolution))) * frequency
1975 lower = np.concatenate([[frequency[0] / (2 ** (1 / (2 * octave_resolution)))], upper[:-1]])
1976 band = upper - lower
1977 frange = upper[-1] - lower[0]
1979 avg = np.mean(db_error, axis=1)
1981 rms_vals = []
1982 std_vals = []
1983 for j in range(db_error.shape[1]):
1984 rms_vals.append(np.sqrt((1.0 / frange) * np.sum((db_error[:, j] ** 2) * band)))
1985 std_vals.append(np.sqrt((1.0 / frange) * np.sum(((db_error[:, j] - avg) ** 2) * band)))
1987 return {
1988 "rms": float(np.mean(rms_vals)),
1989 "std": float(np.mean(std_vals)),
1990 "avg": avg,
1991 }
1994_NUTTALL_COEFFS = {
1995 0: (1.0, 0.0, 0.0, 0.0),
1996 1: (0.5, 0.5, 0.0, 0.0),
1997 2: (0.42, 0.50, 0.08, 0.0),
1998 3: (7938 / 18608, 9240 / 18608, 1430 / 18608, 0.0),
1999 4: (0.42323, 0.49755, 0.07922, 0.0),
2000 5: (0.44959, 0.49364, 0.05677, 0.0),
2001 6: (0.35875, 0.48829, 0.14128, 0.01168),
2002 7: (0.40217, 0.49703, 0.09892, 0.00188),
2003 8: (0.375, 0.5, 0.125, 0.0),
2004 9: (0.40897, 0.5, 0.09103, 0.0),
2005 10: (10 / 32, 15 / 32, 6 / 32, 1 / 32),
2006 11: (0.338936, 0.481973, 0.161054, 0.018027),
2007 12: (0.355768, 0.487396, 0.144232, 0.012604),
2008 13: (0.53836, 0.46164, 0.0, 0.0),
2009 14: (0.4243801, 0.4973406, 0.0782793, 0.0),
2010 15: (0.3635819, 0.4891775, 0.1365995, 0.0106411),
2011}
2014def nuttall_window(n: int, window_type: int = 12) -> np.ndarray:
2015 """
2016 Generate a Nuttall-family window.
2018 Parameters
2019 ----------
2020 n : int
2021 Number of samples in the window.
2022 window_type : int, default=12
2023 Window type identifier. Only the subset required by the compensation
2024 algorithm is supported here.
2026 Returns
2027 -------
2028 ndarray
2029 Window values with shape `(n,)`.
2031 Raises
2032 ------
2033 ValueError
2034 If `n` is not positive or if `window_type` is unsupported.
2036 Notes
2037 -----
2038 This is a minimal implementation derived from the MATLAB ``nuttall``
2039 routine and is included only to support waveform compensation.
2040 """
2041 if n <= 0:
2042 raise ValueError("n must be positive")
2044 if window_type == 16:
2045 return windows.flattop(n, sym=False)
2046 if window_type == 17:
2047 if n % 2 == 0:
2048 t = np.arange(-n // 2, n // 2)
2049 else:
2050 m = (n - 1) / 2
2051 t = np.arange(-m - 0.5, m + 0.5, 1.0)
2052 return np.cos(np.pi * t / n)
2053 if window_type == 19:
2054 return windows.kaiser(n, beta=9.0, sym=False)
2056 if window_type not in _NUTTALL_COEFFS:
2057 raise ValueError(f"Unsupported Nuttall window type {window_type}")
2059 if n % 2 == 0:
2060 t = np.arange(-n // 2, n // 2)
2061 else:
2062 m = (n - 1) / 2
2063 t = np.arange(-m - 0.5, m + 0.5, 1.0)
2065 a0, a1, a2, a3 = _NUTTALL_COEFFS[window_type]
2066 return (
2067 a0
2068 + a1 * np.cos(2 * np.pi * t / n)
2069 + a2 * np.cos(4 * np.pi * t / n)
2070 + a3 * np.cos(6 * np.pi * t / n)
2071 )
2074def window_compensation(
2075 x: np.ndarray,
2076 window: np.ndarray | float,
2077 sample_rate: float,
2078 delay: float = 0.0,
2079) -> Dict[str, np.ndarray]:
2080 """
2081 Apply a compensation waveform to enforce zero terminal cumulative moments.
2083 Parameters
2084 ----------
2085 x : ndarray
2086 Input waveform with shape `(n_samples,)`.
2087 window : ndarray or float
2088 If an array, it is used directly as the compensation window. If a
2089 scalar, it is interpreted as a compensation frequency and a Nuttall
2090 window is generated internally.
2091 sample_rate : float
2092 Sampling rate in samples per second.
2093 delay : float, default=0.0
2094 Delay applied to the compensating waveform, in seconds.
2096 Returns
2097 -------
2098 dict
2099 Dictionary containing:
2101 - ``"comp"`` : ndarray
2102 The compensating waveform.
2103 - ``"xc"`` : ndarray
2104 The compensated waveform.
2105 - ``"coefficients"`` : ndarray
2106 The three coefficients multiplying the DC, sine, and cosine
2107 compensation basis terms.
2109 Notes
2110 -----
2111 This function follows the numerical approach of MATLAB ``sdwcomp``. The
2112 compensation is constructed from three basis terms: a windowed DC term,
2113 a windowed sine term, and a windowed cosine term. Their amplitudes are
2114 chosen so that the final cumulative sums corresponding to acceleration,
2115 velocity, and displacement-like integrals vanish.
2116 """
2117 x = np.asarray(x, dtype=float).reshape(-1)
2119 if np.max(np.abs(x)) <= np.finfo(float).eps:
2120 return {
2121 "comp": np.zeros_like(x),
2122 "xc": np.zeros_like(x),
2123 "coefficients": np.zeros(3),
2124 }
2126 if np.isscalar(window):
2127 freq = float(window)
2128 lc = int(np.fix(sample_rate / freq))
2129 w = nuttall_window(lc, window_type=12)
2130 else:
2131 w = np.asarray(window, dtype=float).reshape(-1)
2132 lc = len(w)
2133 freq = sample_rate / lc
2135 t = np.arange(lc) / sample_rate
2136 acd = w
2137 acs = w * np.sin(2 * np.pi * freq * t)
2138 acc = w * np.cos(2 * np.pi * freq * t)
2140 ld = int(np.fix(delay * sample_rate))
2141 lx = len(x)
2143 if ld < 0:
2144 if lc + ld <= lx:
2145 l = -ld + lx
2146 xw = np.concatenate([np.zeros(l - lx), x])
2147 acd = np.concatenate([acd, np.zeros(l - lc)])
2148 acs = np.concatenate([acs, np.zeros(l - lc)])
2149 acc = np.concatenate([acc, np.zeros(l - lc)])
2150 else:
2151 l = lc
2152 xw = np.concatenate([np.zeros(-ld), x, np.zeros(l - lx + ld)])
2153 else:
2154 if ld + lc <= lx:
2155 l = lx
2156 xw = x.copy()
2157 acd = np.concatenate([np.zeros(ld), acd, np.zeros(lx - ld - lc)])
2158 acs = np.concatenate([np.zeros(ld), acs, np.zeros(lx - ld - lc)])
2159 acc = np.concatenate([np.zeros(ld), acc, np.zeros(lx - ld - lc)])
2160 else:
2161 l = ld + lc
2162 xw = np.concatenate([x, np.zeros(l - lx)])
2163 acd = np.concatenate([np.zeros(ld), acd])
2164 acs = np.concatenate([np.zeros(ld), acs])
2165 acc = np.concatenate([np.zeros(ld), acc])
2167 xv = np.cumsum(xw)
2168 xd = np.cumsum(xv)
2169 xe = np.cumsum(xd)
2171 vcd = np.cumsum(acd)
2172 vcs = np.cumsum(acs)
2173 vcc = np.cumsum(acc)
2175 dcd = np.cumsum(vcd)
2176 dcs = np.cumsum(vcs)
2177 dcc = np.cumsum(vcc)
2179 ecd = np.cumsum(dcd)
2180 ecs = np.cumsum(dcs)
2181 ecc = np.cumsum(dcc)
2183 rhs = -np.array([xv[-1], xd[-1], xe[-1]], dtype=float)
2184 mat = np.array([
2185 [vcd[-1], vcs[-1], vcc[-1]],
2186 [dcd[-1], dcs[-1], dcc[-1]],
2187 [ecd[-1], ecs[-1], ecc[-1]],
2188 ], dtype=float)
2190 coeff = np.linalg.solve(mat, rhs)
2191 comp = coeff[0] * acd + coeff[1] * acs + coeff[2] * acc
2192 xc = xw + comp
2194 return {
2195 "comp": comp,
2196 "xc": xc,
2197 "coefficients": coeff,
2198 }
2201def _inverse_real_spectrum(one_sided_spectrum: np.ndarray) -> np.ndarray:
2202 """
2203 Reconstruct a real time history from a one-sided complex spectrum.
2205 Parameters
2206 ----------
2207 one_sided_spectrum : ndarray
2208 One-sided complex spectrum of length `floor(N/2) + 1`.
2210 Returns
2211 -------
2212 ndarray
2213 Real-valued time history reconstructed by inverse FFT.
2215 Notes
2216 -----
2217 This matches the convention used by MATLAB ``rffti``: the omitted
2218 negative-frequency half of the spectrum is inferred by Hermitian
2219 symmetry.
2220 """
2221 X = np.asarray(one_sided_spectrum, dtype=complex).reshape(-1)
2222 n = len(X)
2224 if np.abs(np.imag(X[-1])) != 0:
2225 Y = np.concatenate([X, np.conj(X[n - 1:0:-1])])
2226 else:
2227 Y = np.concatenate([X, np.conj(X[n - 2:0:-1])])
2229 return np.fft.ifft(Y).real
2232def generate_single_windowed_random_realization(
2233 ref_srs: np.ndarray,
2234 window: np.ndarray,
2235 sample_rate: float,
2236 gravity: float = 386.0,
2237 damping: float = 0.03,
2238 srs_type: int = 9,
2239 match_points_per_octave: float = 4.0,
2240 output_points_per_octave: float = 12.0,
2241 amp_init: Optional[np.ndarray] = None,
2242 phase_init: Optional[np.ndarray] = None,
2243 iterations: int = 10,
2244 correction: float = 0.8,
2245 error_tolerance: float = 0.05,
2246 randomize_phase: bool = False,
2247 max_frequency: Optional[float] = None,
2248 rng: Optional[np.random.Generator] = None,
2249) -> Dict[str, Any]:
2250 """
2251 Generate a single windowed random realization matched to a target SRS.
2253 Parameters
2254 ----------
2255 ref_srs : ndarray
2256 Target SRS breakpoint array with shape `(n_breakpoints, 2)`. The first
2257 column is frequency and the second is SRS amplitude.
2258 window : ndarray
2259 Time-domain window with shape `(n_samples,)`.
2260 sample_rate : float
2261 Sampling rate in samples per second.
2262 gravity : float, default=386.0
2263 Acceleration-of-gravity conversion factor used when integrating to
2264 velocity and displacement.
2265 damping : float, default=0.03
2266 Fraction of critical damping used for SRS calculations.
2267 srs_type : int, default=9
2268 SRS type passed to :func:`sdynpy_srs.srs`.
2269 match_points_per_octave : float, default=4.0
2270 Frequency resolution used during the iterative SRS matching process.
2271 output_points_per_octave : float, default=12.0
2272 Frequency resolution used for the returned SRS curve.
2273 amp_init : ndarray, optional
2274 Initial spectrum amplitude guess as a two-column array
2275 `[frequency, amplitude]`.
2276 phase_init : ndarray, optional
2277 Initial one-sided phase vector in radians.
2278 iterations : int, default=10
2279 Maximum number of matching iterations.
2280 correction : float, default=0.8
2281 Fraction of the SRS error corrected at each iteration.
2282 error_tolerance : float, default=0.05
2283 Relative SRS error tolerance for convergence.
2284 randomize_phase : bool, default=False
2285 If True, ignore `phase_init` and generate a random phase realization.
2286 max_frequency : float, optional
2287 Maximum frequency for nonzero Fourier amplitude. If not specified, the
2288 default is `min(sample_rate/4, ref_srs[-1, 0])`.
2289 rng : numpy.random.Generator, optional
2290 Random number generator used when random phase is required.
2292 Returns
2293 -------
2294 dict
2295 Dictionary containing:
2297 - ``"time"`` : ndarray
2298 Time vector.
2299 - ``"xc"`` : ndarray
2300 Compensated windowed acceleration history.
2301 - ``"rr"`` : ndarray
2302 Stationary random acceleration history before windowing.
2303 - ``"srs_freq"`` : ndarray
2304 Frequency vector for the returned SRS.
2305 - ``"srs_values"`` : ndarray
2306 SRS values for the returned realization.
2307 - ``"srsc"`` : ndarray
2308 Matching-frequency SRS array with columns `[frequency, amplitude]`.
2309 - ``"metrics"`` : WindowedRandomMetrics
2310 Scalar metrics for the realization.
2311 - ``"amp"`` : ndarray
2312 Final matched spectral amplitudes at the matching frequencies.
2313 - ``"phase"`` : ndarray
2314 Final phase vector used for waveform synthesis.
2316 Raises
2317 ------
2318 ValueError
2319 If the FFT line spacing is too coarse to resolve the lowest matched SRS
2320 frequency, or if the internally defined frequency grid is invalid.
2322 Notes
2323 -----
2324 This function is the numerical core corresponding to one pass through the
2325 MATLAB ``wrcomp_tdh`` / ``gwinrand`` workflow.
2326 """
2327 ref_srs = np.asarray(ref_srs, dtype=float)
2328 window = np.asarray(window, dtype=float).reshape(-1)
2329 rng = np.random.default_rng() if rng is None else rng
2331 n = len(window)
2332 t = np.arange(n) / sample_rate
2334 fmin = ref_srs[0, 0]
2335 fmax = ref_srs[-1, 0]
2337 if max_frequency is None:
2338 max_frequency = min(sample_rate / 4.0, fmax)
2340 freq_match = octspace(fmin, fmax, match_points_per_octave).reshape(-1)
2342 if phase_init is None:
2343 block_size = int(2 ** np.ceil(np.log2(n)))
2344 else:
2345 block_size = 2 * (len(phase_init) - 1)
2347 delta_f = sample_rate / block_size
2348 if delta_f > freq_match[0]:
2349 raise ValueError(
2350 "FFT line spacing is too coarse to match the lowest SRS frequency. "
2351 f"delta_f={delta_f}, fmin={freq_match[0]}"
2352 )
2354 freq_line = np.arange(block_size // 2 + 1) * delta_f
2355 if len(freq_line) < 2:
2356 raise ValueError("Frequency line vector too short")
2357 freq_line[0] = freq_line[1] / 10.0
2359 if amp_init is not None and len(amp_init) > 0:
2360 amp_init = np.asarray(amp_init, dtype=float)
2361 amp_match = breakpoint_spectrum(freq_match, amp_init[:, 0], amp_init[:, 1])
2362 else:
2363 amp_match = None
2365 if amp_match is None:
2366 amp_line = 2 * np.pi * (1 / delta_f) * breakpoint_spectrum(
2367 freq_line, ref_srs[:, 0], ref_srs[:, 1]
2368 )
2369 amp_line[0] = 0.0
2370 amp_match = np.interp(freq_match, freq_line, amp_line, left=0.0, right=0.0) / freq_match
2371 else:
2372 amp_line = breakpoint_spectrum(freq_line, freq_match, amp_match)
2374 if randomize_phase or phase_init is None or len(phase_init) == 0:
2375 phase = 2 * np.pi * rng.random(len(amp_line))
2376 phase[0] = 0.0
2377 phase[-1] = 0.0
2378 else:
2379 phase = np.asarray(phase_init, dtype=float).reshape(-1)
2381 srs_required = breakpoint_spectrum(freq_match, ref_srs[:, 0], ref_srs[:, 1])
2383 spectrum = amp_line * np.exp(1j * phase)
2384 rr_full = _inverse_real_spectrum(spectrum)
2385 rr = rr_full[:n]
2386 xr = window * rr
2388 xc = window_compensation(xr, window, sample_rate, delay=0.0)["xc"][:n]
2390 srsc_values = srs(xc, 1.0 / sample_rate, freq_match, damping, srs_type)[0]
2391 srsc = np.column_stack([freq_match, srsc_values])
2393 maxline = np.where(freq_line > max_frequency)[0]
2394 amp_iter = amp_match.copy()
2396 for _ in range(1, iterations):
2397 amp_match = amp_iter + amp_iter * correction * (srs_required - srsc[:, 1]) / srsc[:, 1]
2398 amp_match = np.maximum(0.0, amp_match)
2399 amp_iter = amp_match.copy()
2401 f1, f2 = freq_match[0], freq_match[1]
2402 f3, f4 = freq_match[-2], freq_match[-1]
2403 freqq = np.concatenate([
2404 [0.0],
2405 [f1 + f1 / 2 - f2 / 2],
2406 freq_match,
2407 [f4 + f4 - f3],
2408 [sample_rate / 2.0],
2409 ])
2410 ampp = np.concatenate([
2411 [amp_match[0] / 1000.0],
2412 [amp_match[0]],
2413 amp_match,
2414 [amp_match[-1]],
2415 [amp_match[-1] / 1000.0],
2416 ])
2418 amp_line = np.interp(freq_line, freqq, ampp, left=0.0, right=0.0)
2419 amp_line[maxline] = 0.0
2421 nz = np.abs(amp_line[np.abs(amp_line) > 0])
2422 if len(nz):
2423 amp_line = np.maximum(0.01 * np.min(nz), amp_line)
2425 spectrum = amp_line * np.exp(1j * phase)
2426 rr_full = _inverse_real_spectrum(spectrum)
2427 rr = rr_full[:n]
2428 xr = window * rr
2430 xc = window_compensation(xr, window, sample_rate, delay=0.0)["xc"][:n]
2431 srsc[:, 1] = srs(xc, 1.0 / sample_rate, freq_match, damping, srs_type)[0]
2433 srs_error = (srsc[:, 1] - srs_required) / srs_required
2434 if np.max(np.abs(srs_error)) < error_tolerance:
2435 break
2437 srs_freq = octspace(fmin, fmax, output_points_per_octave)
2438 srs_values = srs(xc, 1.0 / sample_rate, srs_freq, damping, srs_type)[0]
2440 vc = (gravity / sample_rate) * cumulative_trapezoid(xc, initial=0.0)
2441 dc = (1.0 / sample_rate) * cumulative_trapezoid(vc, initial=0.0)
2443 sref_interp = breakpoint_spectrum(srsc[:, 0], ref_srs[:, 0], ref_srs[:, 1])
2444 db_error = 20.0 * np.log10(srsc[:, 1] / sref_interp)
2445 err = lanl_rms_db_error(db_error, srsc[:, 0], output_points_per_octave)
2446 tm = temporal_moments(xc, 1.0 / sample_rate)
2448 metrics = WindowedRandomMetrics(
2449 lanl_rms_db_error=err["rms"],
2450 energy=tm["energy"],
2451 rea=tm["rea"],
2452 peak_accel=float(np.max(np.abs(xc))),
2453 peak_vel=float(np.max(np.abs(vc))),
2454 peak_disp=float(np.max(np.abs(dc))),
2455 )
2457 return {
2458 "time": t,
2459 "xc": xc,
2460 "rr": rr,
2461 "srs_freq": srs_freq,
2462 "srs_values": srs_values,
2463 "srsc": srsc,
2464 "metrics": metrics,
2465 "amp": amp_match,
2466 "phase": phase,
2467 }
2470def generate_windowed_random(
2471 ref_srs: np.ndarray,
2472 window: np.ndarray,
2473 sample_rate: float,
2474 gravity: float = 386.0,
2475 damping: float = 0.03,
2476 srs_type: int = 9,
2477 match_points_per_octave: float = 4.0,
2478 output_points_per_octave: float = 12.0,
2479 n_realizations: int = 1,
2480 amp_init: Optional[np.ndarray] = None,
2481 phase_init: Optional[np.ndarray] = None,
2482 iterations: int = 10,
2483 correction: float = 0.8,
2484 error_tolerance: float = 0.05,
2485 randomize_phase: bool = False,
2486 max_frequency: Optional[float] = None,
2487 seed: Optional[int] = None,
2488) -> WindowedRandomResults:
2489 """
2490 Generate an ensemble of windowed random time histories matched to a target SRS.
2492 Parameters
2493 ----------
2494 ref_srs : ndarray
2495 Target SRS breakpoint array with shape `(n_breakpoints, 2)`. The first
2496 column is frequency and the second is amplitude.
2497 window : ndarray
2498 Time-domain window with shape `(n_samples,)`.
2499 sample_rate : float
2500 Sampling rate in samples per second.
2501 gravity : float, default=386.0
2502 Acceleration-of-gravity conversion factor used when integrating to
2503 velocity and displacement.
2504 damping : float, default=0.03
2505 Fraction of critical damping used for SRS calculations.
2506 srs_type : int, default=9
2507 SRS type passed to :func:`sdynpy_srs.srs`.
2508 match_points_per_octave : float, default=4.0
2509 Frequency resolution used for SRS matching during iteration.
2510 output_points_per_octave : float, default=12.0
2511 Frequency resolution used for the returned SRS curves.
2512 n_realizations : int, default=1
2513 Number of random realizations to generate.
2514 amp_init : ndarray, optional
2515 Initial spectrum amplitude guess as a two-column array
2516 `[frequency, amplitude]`.
2517 phase_init : ndarray, optional
2518 Initial one-sided phase vector in radians.
2519 iterations : int, default=10
2520 Maximum number of matching iterations per realization.
2521 correction : float, default=0.8
2522 Fraction of the SRS error corrected at each iteration.
2523 error_tolerance : float, default=0.05
2524 Relative SRS error tolerance for convergence.
2525 randomize_phase : bool, default=False
2526 If True, ignore `phase_init` and randomize the phase realization.
2527 max_frequency : float, optional
2528 Maximum frequency for nonzero Fourier amplitude. If not specified, the
2529 default is `min(sample_rate/4, ref_srs[-1, 0])`.
2530 seed : int, optional
2531 Seed for the random number generator.
2533 Returns
2534 -------
2535 WindowedRandomResults
2536 Ensemble result object containing all realizations, their SRS curves,
2537 and selection metrics.
2539 Notes
2540 -----
2541 Each realization is generated independently using the same configuration
2542 except for the random phase realization when randomization is enabled.
2543 """
2544 ref_srs = np.asarray(ref_srs, dtype=float)
2545 window = np.asarray(window, dtype=float).reshape(-1)
2547 rng = np.random.default_rng(seed)
2548 t = np.arange(len(window)) / sample_rate
2550 xc_list = []
2551 rr_list = []
2552 srs_list = []
2553 metrics = []
2554 srs_freq = None
2556 for _ in range(n_realizations):
2557 out = generate_single_windowed_random_realization(
2558 ref_srs=ref_srs,
2559 window=window,
2560 sample_rate=sample_rate,
2561 gravity=gravity,
2562 damping=damping,
2563 srs_type=srs_type,
2564 match_points_per_octave=match_points_per_octave,
2565 output_points_per_octave=output_points_per_octave,
2566 amp_init=amp_init,
2567 phase_init=phase_init,
2568 iterations=iterations,
2569 correction=correction,
2570 error_tolerance=error_tolerance,
2571 randomize_phase=randomize_phase,
2572 max_frequency=max_frequency,
2573 rng=rng,
2574 )
2576 xc_list.append(out["xc"])
2577 rr_list.append(out["rr"])
2578 srs_list.append(out["srs_values"])
2579 metrics.append(out["metrics"])
2581 if srs_freq is None:
2582 srs_freq = out["srs_freq"]
2584 return WindowedRandomResults(
2585 time=t,
2586 sample_rate=sample_rate,
2587 gravity=gravity,
2588 window=window,
2589 ref_srs=ref_srs,
2590 srs_frequency=srs_freq,
2591 srs_matrix=np.column_stack(srs_list),
2592 xc_matrix=np.column_stack(xc_list),
2593 rr_matrix=np.column_stack(rr_list),
2594 metrics=metrics,
2595 metadata={
2596 "damping": damping,
2597 "srs_type": srs_type,
2598 "match_points_per_octave": match_points_per_octave,
2599 "output_points_per_octave": output_points_per_octave,
2600 "iterations": iterations,
2601 "correction": correction,
2602 "error_tolerance": error_tolerance,
2603 "randomize_phase": randomize_phase,
2604 "max_frequency": max_frequency,
2605 "n_realizations": n_realizations,
2606 },
2607 )
2610def select_realization(
2611 results: WindowedRandomResults,
2612 criterion: str = "lanl_rms_db_error",
2613 index: Optional[int] = None,
2614) -> SelectedRealization:
2615 """
2616 Select a realization from an ensemble.
2618 Parameters
2619 ----------
2620 results : WindowedRandomResults
2621 Ensemble results returned by :func:`generate_windowed_random`.
2622 criterion : str, default="lanl_rms_db_error"
2623 Name of the metric field to minimize if `index` is not provided.
2624 index : int, optional
2625 Explicit zero-based realization index to select. If provided, this
2626 takes precedence over `criterion`.
2628 Returns
2629 -------
2630 SelectedRealization
2631 Selected realization containing the time history, SRS, and scalar
2632 metrics.
2634 Raises
2635 ------
2636 AttributeError
2637 If `criterion` is not a valid attribute of
2638 :class:`WindowedRandomMetrics`.
2640 Notes
2641 -----
2642 This function provides the deterministic equivalent of the MATLAB
2643 ``gwinrand_select`` helper without interactive point-picking.
2644 """
2645 if index is None:
2646 values = np.array([getattr(m, criterion) for m in results.metrics], dtype=float)
2647 index = int(np.argmin(values))
2649 return SelectedRealization(
2650 index=index,
2651 time=results.time,
2652 xc=results.xc_matrix[:, index],
2653 rr=results.rr_matrix[:, index],
2654 srs=np.column_stack([results.srs_frequency, results.srs_matrix[:, index]]),
2655 metrics=results.metrics[index],
2656 )