Coverage for src / sdynpy / signal_processing / sdynpy_frf.py: 10%
169 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Tools for computing frequency response functions
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import numpy as np
25import scipy.signal as sig
26import matplotlib.pyplot as plt
27from .sdynpy_lrm import frf_local_model
30def sysmat2frf(frequencies, M, C, K, frf_type='disp'):
31 '''Compute Frequency Response Functions given system matrices M, C, and K
33 This function computes frequency response functions (FRFs) given the system
34 mass matrix, M, damping matrix, C, and stiffness matrix, K, in the equation
35 of motion
37 M x_dd + C x_d + K x = F
39 This will return the frequency response function
41 H = x/F
43 as a function of frequency.
45 Parameters
46 ----------
47 frequencies : ndarray
48 frequencies should be a 1D np.ndarray consisting of the frequency lines
49 at which to evaluate the frequency response function. They should be
50 in units of cycles per second or hertz (Hz), rather than in angular
51 frequency of radians/s.
52 M : ndarray
53 M should be a 2D np.ndarray consisting of the system mass matrix.
54 C : ndarray
55 C should be a 2D np.ndarray consisting of the system damping matrix.
56 K : ndarray
57 K should be a 2D np.ndarray consisting of the system stiffness matrix.
58 frf_type : str
59 frf_type should be one of ['disp','vel','accel'] or ['displacement',
60 'velocity','acceleration'] to specify which "type" of frequency
61 response function to compute. By default it computes a displacement or
62 "receptance" FRF. However, if an acceleration or "Accelerance" FRF is
63 desired, specify 'accel' instead. The displacement, velocity, and
64 acceleration FRFs differ by a factor of 1j*omega where omega is the
65 angular frequency at a given frequency line.
67 Returns
68 -------
69 H : ndarray
70 A 3D np array with shape (nf,no,ni), where nf is the number of
71 frequency lines, no is the number of outputs, and ni is the number of
72 inputs. Since M, C, and K should be square, ni should equal no.
73 Values in H are complex.
75 Notes
76 -----
77 This performs a direct inversion of the system matrix, therefore it is not
78 advisable to compute FRFs of large systems using this method. An
79 alternative would be to compute modes first then compute the FRFs from
80 the modes.
81 '''
82 if not frf_type in ['displacement', 'velocity', 'acceleration', 'disp', 'vel',
83 'accel']:
84 raise ValueError('frf_type must be one of {:}'.format(['displacement',
85 'velocity', 'acceleration', 'disp', 'vel', 'accel']))
86 Z = (-(2 * np.pi * frequencies[:, np.newaxis, np.newaxis])**2 * M
87 + 1j * (2 * np.pi * frequencies[:, np.newaxis, np.newaxis]) * C
88 + K)
89 H = np.linalg.inv(Z)
90 if frf_type in ['vel', 'velocity']:
91 H = 1j * (2 * np.pi * frequencies[:, np.newaxis, np.newaxis]) * H
92 elif frf_type in ['accel', 'acceleration']:
93 H = -(2 * np.pi * frequencies[:, np.newaxis, np.newaxis])**2 * H
95 return H
98def modes2frf(frequencies, natural_frequencies, damping_ratios, mode_shapes=None,
99 input_mode_shapes=None, frf_type='disp'):
100 '''Compute Frequency Response Functions given modal properties
102 This function computes frequency responses given modal parameters.
104 Parameters
105 ----------
106 frequencies : ndarray
107 frequencies should be a 1D np.ndarray consisting of the frequency lines
108 at which to evaluate the frequency response function. They should be
109 in units of cycles per second or hertz (Hz), rather than in angular
110 frequency of radians/s.
111 natural_frequencies : ndarray
112 Natural frequencies of the structure in cycles per second or hertz (Hz)
113 rather than in angular frequency of radians/s.
114 damping_ratios : ndarray
115 Critical damping ratios of the structure in ratio form rather than
116 percentange, e.g. 2% damping would be specified as 0.02 rather than 2.
117 mode_shapes : ndarray, optional
118 A 2D mode shape matrix with shape (no,nm) where no is the number of
119 output responses and nm is the number of modes. If the optional
120 argument input_mode_shapes is not specified, this mode shape matrix
121 will also be used for the inputs, resulting in a square FRF matrix.
122 If mode_shapes is not specified, the values of the modal FRF matrix
123 are returned as 2D array.
124 input_mode_shapes : ndarray, optional
125 A 2D mode shape matrix with shape (ni,nm) where ni is the number of
126 input forces and nm is the number of modes. If the optional argument
127 input_mode_shapes is specified, it be used for the inputs, resulting in
128 a potentially nonsquare FRF matrix.
129 frf_type : str, optional
130 frf_type should be one of ['disp','vel','accel'] or ['displacement',
131 'velocity','acceleration'] to specify which "type" of frequency
132 response function to compute. By default it computes a displacement or
133 "receptance" FRF. However, if an acceleration or "Accelerance" FRF is
134 desired, specify 'accel' instead. The displacement, velocity, and
135 acceleration FRFs differ by a factor of 1j*omega where omega is the
136 angular frequency at a given frequency line.
138 Returns
139 -------
140 H : ndarray
141 A 3D (or 2D) np array depending on whether or not mode_shapes is
142 defined with shape (nf,no,ni) or (nf,nm), where nf is the number of
143 frequency lines, no is the number of outputs, ni is the number of
144 inputs, and nm is the number of modes. Values in H are complex.
146 Notes
147 -----
148 This function assumes mass normalized mode shapes.
150 '''
151 if not frf_type in ['displacement', 'velocity', 'acceleration', 'disp', 'vel',
152 'accel']:
153 raise ValueError('frf_type must be one of {:}'.format(['displacement',
154 'velocity', 'acceleration', 'disp', 'vel', 'accel']))
156 Z = (-(2 * np.pi * frequencies[:, np.newaxis])**2
157 + 1j * (2 * np.pi * frequencies[:, np.newaxis]) * 2 *
158 damping_ratios * (2 * np.pi * natural_frequencies)
159 + (2 * np.pi * natural_frequencies)**2)
160 if mode_shapes is not None:
161 if input_mode_shapes is None:
162 input_mode_shapes = mode_shapes
163 H = np.einsum('ij,fj,lj->fil', mode_shapes, 1 / Z, input_mode_shapes)
164 else:
165 H = 1/Z
167 if frf_type in ['vel', 'velocity']:
168 if np.ndim(H) == 3:
169 H = 1j * (2 * np.pi * frequencies[:, np.newaxis, np.newaxis]) * H
170 elif np.ndim(H) == 2:
171 H = 1j * (2 * np.pi * frequencies[:, np.newaxis]) * H
172 # num = [1,0]
173 elif frf_type in ['accel', 'acceleration']:
174 if np.ndim(H) == 3:
175 H = -(2 * np.pi * frequencies[:, np.newaxis, np.newaxis])**2 * H
176 elif np.ndim(H) == 2:
177 H = -(2 * np.pi * frequencies[:, np.newaxis])**2 * H
178 return H
180def timedata2frf(references, responses, dt=1, samples_per_average=None,
181 overlap=0.0, method='H1', window=np.array((1.0,)),
182 response_fft=lambda array: np.fft.rfft(array, axis=-1),
183 reference_fft=lambda array: np.fft.rfft(array, axis=-1),
184 response_fft_array=None, reference_fft_array=None,
185 return_model_data=False, **lrm_kwargs):
186 # TODO Update DOCSTRING with changes after they are all made.
187 '''Creates an FRF matrix given time histories of responses and references
189 This function creates a nf x no x ni FRF matrix from the time histories
190 provided.
192 Parameters
193 ----------
194 references : ndarray
195 A ni x nt or nt array where ni is the number of references and nt is
196 the number of time steps in the signal. If averaging is specified, nt
197 should be divisible by the number of averages.
198 responses : ndarray
199 A no x nt or nt array where no is the number of responses and nt is the
200 number of time steps in the signal. If averaging is specified, nt
201 should be divisible by the number of averages.
202 dt : float
203 The time between samples
204 samples_per_average : int
205 The number of time samples per average. If not specified, it is set
206 to the number of samples in the time signal, and no averaging is
207 performed
208 overlap : float
209 The overlap as a fraction of the frame (e.g. 0.5 specifies 50% overlap).
210 If not specified, no overlap is used.
211 method : str in ['H1','H2','Hv','Hs','LRM']
212 The method for creating the frequency response function. 'H1' is
213 default if not specified.
214 window : ndarray or str
215 A 1D ndarray with length samples_per_average that specifies the
216 coefficients of the window. No window is applied if not specified.
217 If a string is specified, then the window will be obtained from scipy
218 fft : function
219 FFT Function that should be used. FFT must take the fft over axis -1.
220 response_fft_array : np.ndarray
221 Array to store the data into before taking the FFT. Should be
222 size number_of_responses, n_averages, samples_per_average
223 reference_fft_array : np.ndarray
224 Array to store the data into before taking the FFT. Should be
225 size number_of_references, n_averages, samples_per_average
226 return_model_data : boolean
227 Wheter to return selected model orders used in the 'LRM' method.
228 default is False
229 **lrm_kwargs
230 Additional keyword arguments specify parameters if method == 'LRM'.
231 Possible arguments are: f_out (ndarray), bandwidth (float), transient
232 (boolean), modelset (iterable of (3,1) ndarray), export_ratio (float),
233 max_parallel (int), print_time (bool). See sdynpy_lrm.frf_local_model.
234 samples_per_average, overlap, and window are void if method=='LRM'.
236 Returns
237 -------
238 frequencies : ndarray
239 A nf array of the frequency values associated with H
240 H : ndarray
241 A nf x no x ni array where nf is the number of frequency lines, no is
242 the number of outputs, and ni is the number of inputs.
243 model_data: dict
244 Contains None if method != 'LRM'. See sdynpy_lrm.frf_local_model.
246 Notes
247 -----
248 There are requirements for the shapes of references and responses for some
249 FRF computations.
251 '''
252 if references.shape[-1] != responses.shape[-1]:
253 raise ValueError(
254 'reference and responses must have the same number of time steps (last dimension)!')
255 if not method in ['H1', 'H2', 'H3', 'Hs', 'Hv', 'Hcd','LRM']:
256 raise ValueError('method parameter must be one of ["H1","H2","H3","Hv","Hcd"]')
257 if references.ndim == 1:
258 references = references[np.newaxis, :]
259 if responses.ndim == 1:
260 responses = responses[np.newaxis, :]
261 # Set up time indices
262 ntimes_total = references.shape[-1]
263 if samples_per_average is None or method == 'LRM':
264 samples_per_average = ntimes_total
265 overlap_samples = int(samples_per_average * overlap)
266 time_starts = np.arange(0, ntimes_total - samples_per_average +
267 1, samples_per_average - overlap_samples)
268 time_indices_for_averages = time_starts[:, np.newaxis] + np.arange(samples_per_average)
269 if method == 'LRM':
270 window=np.array((1.0,)) # local modeling does not use windowing
271 if isinstance(window, str):
272 window = sig.windows.get_window(window.lower(), samples_per_average, fftbins=True)
273 # Sort the references and responses into their time arrays
274 if response_fft_array is None:
275 response_fft_array = responses[..., time_indices_for_averages] * window
276 else:
277 response_fft_array[...] = responses[:, time_indices_for_averages] * window
278 if reference_fft_array is None:
279 reference_fft_array = references[:, time_indices_for_averages] * window
280 else:
281 reference_fft_array[...] = references[:, time_indices_for_averages] * window
282 # Compute FFTs
283 frequencies = np.fft.rfftfreq(samples_per_average, dt)
284 response_ffts = response_fft(response_fft_array)
285 reference_ffts = reference_fft(reference_fft_array)
286 # Now start computing FRFs
287 if method == 'H1':
288 # We want to compute X*F^H = [X1;X2;X3][F1^H F2^H F3^H]
289 Gxf = np.einsum('...iaf,...jaf->...fij', response_ffts, np.conj(reference_ffts))
290 Gff = np.einsum('...iaf,...jaf->...fij', reference_ffts, np.conj(reference_ffts))
291 # Add small values to any matrices that are singular
292 singular_matrices = np.abs(np.linalg.det(Gff)) < 2 * np.finfo(Gff.dtype).eps
293 Gff[singular_matrices] += np.eye(Gff.shape[-1]) * np.finfo(Gff.dtype).eps
294 H = np.moveaxis(np.linalg.solve(np.moveaxis(Gff, -2, -1), np.moveaxis(Gxf, -2, -1)), -2, -1)
295 elif method == 'H2':
296 if (response_ffts.shape != reference_ffts.shape):
297 raise ValueError('For H2, Number of inputs must equal number of outputs')
298 Gxx = np.einsum('...iaf,...jaf->...fij', response_ffts, np.conj(response_ffts))
299 Gfx = np.einsum('...iaf,...jaf->...fij', reference_ffts, np.conj(response_ffts))
300 singular_matrices = np.abs(np.linalg.det(Gfx)) < 2 * np.finfo(Gfx.dtype).eps
301 Gfx[singular_matrices] += np.eye(Gfx.shape[-1]) * np.finfo(Gfx.dtype).eps
302 H = np.moveaxis(np.linalg.solve(np.moveaxis(Gfx, -2, -1), np.moveaxis(Gxx, -2, -1)), -2, -1)
303 elif method == 'H3':
304 if (response_ffts.shape != reference_ffts.shape):
305 raise ValueError('For H3, Number of inputs must equal number of outputs')
306 Gxf = np.einsum('...iaf,...jaf->...fij', response_ffts, np.conj(reference_ffts))
307 Gff = np.einsum('...iaf,...jaf->...fij', reference_ffts, np.conj(reference_ffts))
308 # Add small values to any matrices that are singular
309 singular_matrices = np.abs(np.linalg.det(Gff)) < 2 * np.finfo(Gff.dtype).eps
310 Gff[singular_matrices] += np.eye(Gff.shape[-1]) * np.finfo(Gff.dtype).eps
311 Gxx = np.einsum('...iaf,...jaf->...fij', response_ffts, np.conj(response_ffts))
312 Gfx = np.einsum('...iaf,...jaf->...fij', reference_ffts, np.conj(response_ffts))
313 singular_matrices = np.abs(np.linalg.det(Gfx)) < 2 * np.finfo(Gfx.dtype).eps
314 Gfx[singular_matrices] += np.eye(Gfx.shape[-1]) * np.finfo(Gfx.dtype).eps
315 H = (np.moveaxis(np.linalg.solve(np.moveaxis(Gfx, -2, -1), np.moveaxis(Gxx, -2, -1)), -2, -1) +
316 np.moveaxis(np.linalg.solve(np.moveaxis(Gff, -2, -1), np.moveaxis(Gxf, -2, -1)), -2, -1)) / 2
317 elif method == 'Hcd':
318 Gxf = np.einsum('...iaf,...jaf->...fij', response_ffts, np.conj(reference_ffts))
319 Gff = np.einsum('...iaf,...jaf->...fij', reference_ffts, np.conj(reference_ffts))
320 # Add small values to any matrices that are singular
321 singular_matrices = np.abs(np.linalg.det(Gff)) < 2 * np.finfo(Gff.dtype).eps
322 Gff[singular_matrices] += np.eye(Gff.shape[-1]) * np.finfo(Gff.dtype).eps
323 Lfz = np.linalg.cholesky(Gff)
324 Lzf = np.conj(np.moveaxis(Lfz, -2, -1))
325 Gxz = np.moveaxis(np.linalg.solve(np.moveaxis(
326 Lzf, -2, -1), np.moveaxis(Gxf, -2, -1)), -2, -1)
327 H = np.moveaxis(np.linalg.solve(np.moveaxis(Lfz, -2, -1), np.moveaxis(Gxz, -2, -1)), -2, -1)
328 elif method == 'Hv':
329 Gxx = np.einsum('...iaf,...iaf->...if', response_ffts,
330 np.conj(response_ffts))[..., np.newaxis, np.newaxis]
331 Gxf = np.einsum('...iaf,...jaf->...ifj', response_ffts,
332 np.conj(reference_ffts))[..., np.newaxis, :]
333 Gff = np.einsum('...iaf,...jaf->...fij', reference_ffts,
334 np.conj(reference_ffts))[..., np.newaxis, :, :, :]
335 # Broadcast over all responses
336 Gff = np.broadcast_to(Gff,Gxx.shape[:-2]+Gff.shape[-2:])
337 Gffx = np.block([[Gff, np.conj(np.moveaxis(Gxf, -2, -1))],
338 [Gxf, Gxx]])
339 # Compute eigenvalues
340 lam, evect = np.linalg.eigh(np.moveaxis(Gffx, -2, -1))
341 # Get the evect corresponding to the minimum eigenvalue
342 evect = evect[..., 0] # Assumes evals are sorted ascending
343 H = np.moveaxis(-evect[..., :-1] / evect[..., -1:], # Scale so last value is -1
344 -3, -2)
345 elif method == 'LRM':
346 frequencies, H, model_data = frf_local_model(reference_ffts,response_ffts,
347 frequencies,**lrm_kwargs)
348 else:
349 raise NotImplementedError('Method {:} has not been implemented yet!'.format(method))
351 if return_model_data:
352 if method != 'LRM':
353 model_data = {'info': None,'modelset': None,'model_selected': None}
354 return frequencies, H, model_data
355 else:
356 return frequencies, H
359def fft2frf(references, responses, method='H1',freqs_in=None,**kwargs):
360 '''Creates an FRF matrix given ffts of responses and references
362 This function creates a nf x no x ni FRF matrix from the ffts
363 provided.
365 Parameters
366 ----------
367 references : ndarray
368 A ni x nf or nf array where ni is the number of references and nt is
369 the number of frequencies in the fft.
370 responses : ndarray
371 A no x nf or nf array where no is the number of responses and nf is the
372 number of time frequencies in the fft.
373 method : str in ['H1','H2','Hv','Hs','LRM']
374 The method for creating the frequency response function.
375 freqs_in : ndarray
376 Only used if method == 'LRM'. A nf or nf x 1 array of frquency values.
377 **lrm_kwargs
378 Additional keyword arguments specify parameters if method == 'LRM'.
379 Possible arguments are: f_out (ndarray), bandwidth (float), transient
380 (boolean), modelset (iterable of (3,1) ndarray), max_parallel (int),
381 export_ratio (float). See sdynpy_lrm.frf_local_model.
383 Returns
384 -------
385 H : ndarray
386 A nf x no x ni array where nf is the number of frequency lines, no is
387 the number of outputs, and ni is the number of inputs. The output
388 frequency lines will correspond to the frequency lines in the ffts.
389 freqs_out : None or ndarray
390 None unless method == 'LRM'. See sdynpy_lrm.frf_local_model.
391 model_data: None dict
392 None unless method == 'LRM'. See sdynpy_lrm.frf_local_model.
394 Notes
395 -----
396 There are requirements for the shapes of references and responses for some
397 FRF computation methods. No averaging is performed by this function.
399 '''
400 if references.ndim == 1:
401 references = references[np.newaxis, :]
402 elif references.ndim > 2:
403 raise ValueError('references should be at maximum a 2 dimensional array')
404 if responses.ndim == 1:
405 responses = responses[np.newaxis, :]
406 elif responses.ndim > 2:
407 raise ValueError('responses should be at maximum a 2 dimensional array')
408 if not method in ['H1', 'H2', 'Hs', 'Hv', 'LRM']:
409 raise ValueError('method parameter must be one of ["H1","H2","Hs","Hv","LRM"]')
410# H = np.zeros((references.shape[1],responses.shape[0],references.shape[0]),dtype=complex)
411# for i in range(H.shape[1]):
412# for j in range(H.shape[2]):
413# if method == 'H1':
414# H[:,i,j] = (responses[i,:]*references[j,:].conj())/(references[j,:]*references[j,:].conj())
415# else:
416# raise NotImplementedError('Method {:} has not been implemented'.format(method))
417 if method == 'H1':
418 Gxf = np.einsum('if,jf->fij', responses, references.conj())
419 Gff = np.einsum('if,jf->fij', references, references.conj())
420 H = np.linalg.solve(Gff.transpose(0, 2, 1), Gxf.transpose((0, 2, 1))).transpose((0, 2, 1))
421 elif method == 'LRM':
422 if freqs_in is None:
423 raise ValueError('if method == \'LRM\', freqs_in must be specified')
424 import sdynpy.signal_processing.sdynpy_lrm as frflm
425 freqs_out, H, model_data = frflm.frf_local_model(references, responses,
426 freqs_in, **lrm_kwargs)
427 else:
428 raise NotImplementedError('Method {:} has not been implemented'.format(method))
429 if method != 'LRM':
430 freqs_out = None
431 model_data = None
432 return H, freqs_out, model_data
435def plot(H, f, responses=None, references=None, real_imag=False):
436 fig = plt.figure()
437 phase_axis = fig.add_subplot(2, 1, 1)
438 mag_axis = fig.add_subplot(2, 1, 2)
440 if responses is None:
441 responses = np.arange(H.shape[1])[:, np.newaxis]
442 if references is None:
443 references = np.arange(H.shape[2])
445 H_to_plot = H[:, responses, references]
446 H_to_plot = H_to_plot.transpose(*np.arange(1, H_to_plot.ndim), 0)
448 for index in np.ndindex(H_to_plot.shape[:-1]):
449 if real_imag:
450 phase_axis.plot(f, np.imag(H_to_plot[index]))
451 mag_axis.plot(f, np.real(H_to_plot[index]))
452 else:
453 phase_axis.plot(f, np.angle(H_to_plot[index]) * 180 / np.pi)
454 mag_axis.plot(f, np.abs(H_to_plot[index]))
456 if real_imag:
457 phase_axis.set_ylabel('Imag(H)')
458 mag_axis.set_ylabel('Real(H)')
459 else:
460 phase_axis.set_ylabel('Angle(H)')
461 mag_axis.set_ylabel('Abs(H)')
462 mag_axis.set_yscale('log')
463 mag_axis.set_xlabel('Frequency')
464 fig.tight_layout()
467def delay_signal(times, signal, dt):
468 '''Delay a time signal by the specified amount
470 This function takes a signal and delays it by a specified amount of time
471 that need not be an integer number of samples. It does this by adjusting
472 the phaes of the signal's FFT.
474 Parameters
475 ----------
476 times : np.ndarray
477 A signal specifying the time values at which the samples in signal
478 occur.
479 signal : np.ndarray
480 A signal that is to be delayed (n_signals x len(times))
481 dt : float
482 The amount of time to delay the signal
484 Returns
485 -------
486 updated_signal : np.ndarray
487 The time-shifted signal.
488 '''
489 fft_omegas = np.fft.fftfreq(len(times), np.mean(np.diff(times))) * 2 * np.pi
490 signal_fft = np.fft.fft(signal, axis=-1)
491 signal_fft *= np.exp(-1j * fft_omegas * dt)
492 return np.fft.ifft(signal_fft, axis=-1).real.astype(signal.dtype)