Coverage for src / sdynpy / signal_processing / sdynpy_frf_inverse.py: 9%
89 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 inverses of 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
25from scipy.interpolate import interp1d
26import warnings
29def frf_inverse(frf_matrix,
30 method='standard',
31 response_weighting_matrix=None,
32 reference_weighting_matrix=None,
33 regularization_weighting_matrix=None,
34 regularization_parameter=None,
35 cond_num_threshold=None,
36 num_retained_values=None):
37 """
38 Computes the inverse of an FRF matrix for source estimation problems.
40 Parameters
41 ----------
42 frf_matrix : NDArray
43 Transfer function as an np.ndarray, should be organized such that the
44 size is [number of lines, number of responses, number of references]
45 method : str, optional
46 The method to be used for the FRF matrix inversions. The available
47 methods are:
48 - standard - basic pseudo-inverse via numpy.linalg.pinv with the
49 default rcond parameter, this is the default method
50 - threshold - pseudo-inverse via numpy.linalg.pinv with a specified
51 condition number threshold
52 - tikhonov - pseudo-inverse using the Tikhonov regularization method
53 - truncation - pseudo-inverse where a fixed number of singular values
54 are retained for the inverse
55 response_weighting_matrix : sdpy.Matrix or np.ndarray, optional
56 Diagonal matrix used to weight response degrees of freedom (to solve the
57 problem as a weight least squares) by multiplying the rows of the FRF
58 matrix by a scalar weights. This matrix can also be a 3D matrix such that
59 the the weights are different for each frequency line. The matrix should
60 be sized [number of lines, number of references, number of references],
61 where the number of lines either be one (the same weights at all frequencies)
62 or the length of the abscissa (for the case where a 3D matrix is supplied).
63 reference_weighting_matrix : sdpy.Matrix or np.ndarray, optional
64 Diagonal matrix used to weight reference degrees of freedom (generally for
65 normalization) by multiplying the columns of the FRF matrix by a scalar weights.
66 This matrix can also be a 3D matrix such that the the weights are different
67 for each frequency line. The matrix should be sized
68 [number of lines, number of references, number of references], where the number
69 of lines either be one (the same weights at all frequencies) or the length
70 of the abscissa (for the case where a 3D matrix is supplied).
71 regularization_weighting_matrix : sdpy.Matrix or np.ndarray, optional
72 Matrix used to weight input degrees of freedom via Tikhonov regularization.
73 This matrix can also be a 3D matrix such that the the weights are different
74 for each frequency line. The matrix should be sized
75 [number of lines, number of references, number of references], where the number
76 of lines either be one (the same weights at all frequencies) or the length
77 of the abscissa (for the case where a 3D matrix is supplied).
78 regularization_parameter : float or np.ndarray, optional
79 Scaling parameter used on the regularization weighting matrix when the tikhonov
80 method is chosen. A vector of regularization parameters can be provided so the
81 regularization is different at each frequency line. The vector must match the
82 length of the abscissa in this case (either be size [num_lines,] or [num_lines, 1]).
83 cond_num_threshold : float or np.ndarray, optional
84 Condition number used for SVD truncation when the threshold method is chosen.
85 A vector of condition numbers can be provided so it varies as a function of
86 frequency. The vector must match the length of the abscissa in this case.
87 num_retained_values : float or np.ndarray, optional
88 Number of singular values to retain in the pseudo-inverse when the truncation
89 method is chosen. A vector of can be provided so the number of retained values
90 can change as a function of frequency. The vector must match the length of the
91 abscissa in this case.
93 Returns
94 -------
95 np.ndarray
96 Inverse of the supplied FRF matrix
98 Raises
99 ------
100 Warning
101 If regularization_weighting_matrix is supplied without selecting the tikhonov method
102 Exception
103 If the threshold method is chosen but a condition number threshold isn't supplied
104 Exception
105 If the declared method is not one of the available options
107 Notes
108 -----
109 This function solves the inverse problem for the supplied FRF matrix. All of the inverse
110 methods use the SVD (or modified SVD) to compute the pseudo-inverse.
112 References
113 ----------
114 .. [1] Wikipedia, "Moore-Penrose inverse".
115 https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse
116 .. [2] A.N. Tithe, D.J. Thompson, The quantification of structure-borne transmission pathsby inverse methods. Part 2: Use of regularization techniques,
117 Journal of Sound and Vibration, Volume 264, Issue 2, 2003, Pages 433-451, ISSN 0022-460X,
118 https://doi.org/10.1016/S0022-460X(02)01203-8.
119 .. [3] Wikipedia, "Ridge regression".
120 https://en.wikipedia.org/wiki/Ridge_regression
121 """
122 if regularization_weighting_matrix is not None and method != 'tikhonov':
123 warnings.warn("Regularization weighting matrix is only used during Tikhonov regularization and is not being used")
125 if response_weighting_matrix is not None:
126 frf_matrix = response_weighting_matrix@frf_matrix
127 if reference_weighting_matrix is not None:
128 frf_matrix = frf_matrix@reference_weighting_matrix
130 if method == 'standard':
131 frf_pinv = np.linalg.pinv(frf_matrix)
132 elif method == 'threshold':
133 if cond_num_threshold is None:
134 raise ValueError('A condition number threshold must be supplied when using the threshold method')
135 rcond_thresh = 1/cond_num_threshold
136 frf_pinv = np.linalg.pinv(frf_matrix, rcond=rcond_thresh)
137 elif method == 'tikhonov':
138 frf_pinv = pinv_by_tikhonov(frf_matrix,
139 regularization_weighting_matrix=regularization_weighting_matrix,
140 regularization_parameter=regularization_parameter)
141 elif method == 'truncation':
142 frf_pinv = pinv_by_truncation(frf_matrix,
143 num_retained_values=num_retained_values)
144 else:
145 error_string_start = 'Declared method "'
146 error_string_end = '" is not an available option'
147 raise NameError(error_string_start+method+error_string_end)
149 return frf_pinv
152def pinv_by_tikhonov(frf_matrix,
153 regularization_weighting_matrix=None,
154 regularization_parameter=None):
155 """
156 Computes the pseudo-inverse of an FRF matrix via Tikhonov regularization
158 Parameters
159 ----------
160 frf_matrix : NDArray
161 Transfer function as an np.ndarray, should be organized such that the
162 frequency is the first axis, the responses are the second axis, and the
163 references are the third axis
164 regularization_weighting_matrix : sdpy.Matrix or np.ndarray, optional
165 Matrix used to weight input degrees of freedom in the regularization, the
166 default is identity. This can be a 3D matrix such that the the weights are
167 different for each frequency line. The matrix should be sized
168 [number of lines, number of references, number of references], where the number
169 of lines either be one (the same weights at all frequencies) or the length
170 of the abscissa (for the case where a 3D matrix is supplied).
171 regularization_parameter : float or np.ndarray, optional
172 Scaling parameter used on the regularization weighting matrix. A vector of
173 regularization parameters can be provided so the regularization is different at
174 each frequency line. The vector must match the length of the abscissa in this
175 case (either be size [num_lines,] or [num_lines, 1]).
177 Returns
178 -------
179 np.ndarray
180 Inverse of the supplied FRF matrix
182 Raises
183 ------
184 Exception
185 If a regularization parameter isn't supplied
186 Warning
187 If a 2D regularization weighting matrix is supplied when a vector of regularization parameters is supplied
189 Notes
190 -----
191 Tikhonov regularization is being done via the SVD method (details on Wikipedia or
192 in the paper by Tithe and Thompson).
194 References
195 ----------
196 .. [1] A.N. Tithe, D.J. Thompson, The quantification of structure-borne transmission pathsby inverse methods. Part 2: Use of regularization techniques,
197 Journal of Sound and Vibration, Volume 264, Issue 2, 2003, Pages 433-451, ISSN 0022-460X,
198 https://doi.org/10.1016/S0022-460X(02)01203-8.
199 .. [2] Wikipedia, "Ridge regression".
200 https://en.wikipedia.org/wiki/Ridge_regression
201 """
203 if regularization_parameter is None:
204 raise ValueError('A regularization parameter must be supplied when using Tikhonov regularization')
206 regularization_parameter = np.asarray(regularization_parameter, dtype=np.float64)
208 u, s, vh = np.linalg.svd(frf_matrix, full_matrices=False)
209 vh_conjT = np.swapaxes(vh.conj(), -2, -1)
210 u_conjT = np.swapaxes(u.conj(), -2, -1)
212 # Need to add a second axis to the regularization parameter if it's a vector
213 # so it broadcasts appropriately.
214 if regularization_parameter.size > 1 and regularization_parameter.ndim == 1:
215 regularization_parameter = regularization_parameter[..., np.newaxis]
217 if regularization_parameter.size > 1 and regularization_parameter.shape[1] > 1:
218 raise ValueError('The regularization parameter is not a scalar or vector and cannot be interpreted')
220 if regularization_parameter.size > 1 and regularization_parameter.size != frf_matrix.shape[0]:
221 raise ValueError('The regularization parameter is a vector but the size'
222 + ' does not match the number of frequency lines and cannot be broadcasted to the FRFs')
224 # This is adding a third dimension to the regularization weighting matrix (if
225 # it doesn't already exist) in the case that the regularization parameter is a
226 # vector, this is required so the weighting matrix and parameter will broadcast
227 # together appropriately.
228 if regularization_weighting_matrix is not None and regularization_parameter.size > 1 and regularization_weighting_matrix.ndim == 2:
229 warnings.warn('The regularization weighting matrix is 2D while the regularization'
230 + ' parameter is a vector, the regularization weighting matrix is '
231 + 'being expanded to 3D. Consider adding a third dimension to the regularization weighting matrix as appropriate')
232 desired_shape = [regularization_parameter.shape[0], regularization_weighting_matrix.shape[0], regularization_weighting_matrix.shape[1]]
233 regularization_weighting_matrix = np.broadcast_to(regularization_weighting_matrix, desired_shape)
235 # This is adding a third dimension to the regularization parameter vector when a
236 # regularization weighting matrix is being used. This is required so the broadcasting
237 # works correctly. I.E., if the regularization weighting matrix is sized
238 # [num_lines, num_refs, num_refs] (which it always is), the regularization parameter
239 # vector has to be sized [num_lines, 1, 1]. This is not required when the regularization
240 # parameter is a scaler.
241 if (regularization_weighting_matrix is not None
242 and regularization_weighting_matrix.ndim == 3
243 and regularization_parameter.size > 1 and regularization_parameter.ndim != 3):
244 regularization_parameter = regularization_parameter[..., np.newaxis]
246 if regularization_weighting_matrix is None:
247 regularized_sInv = s/(s*s+regularization_parameter)
248 if frf_matrix.ndim == 3:
249 # This turns s into a 3D matrix with the singular values on the diagonal
250 regularized_sInv = np.apply_along_axis(np.diag, 1, regularized_sInv)
251 else:
252 # This is for the case where a single frequency line is supplied for the FRF
253 regularized_sInv = np.diag(regularized_sInv)
254 else:
255 if frf_matrix.ndim == 3:
256 # This turns s into a 3D matrix with the singular values on the diagonal
257 s = np.apply_along_axis(np.diag, 1, s)
258 regularized_sInv = np.linalg.pinv(np.moveaxis(s, -1, -2)@s+regularization_weighting_matrix*regularization_parameter)@np.moveaxis(s, -1, -2)
259 else:
260 # This is for the case where a single frequency line is supplied for the FRF
261 regularized_sInv = s/(s*s+regularization_parameter*regularization_weighting_matrix)
263 frf_pinv = vh_conjT@regularized_sInv@u_conjT
265 return frf_pinv
268def pinv_by_truncation(frf_matrix,
269 num_retained_values):
270 """
271 Computes the pseudo-inverse of an FRF matrix where a fixed number of singular
272 values are retained for the inverse
274 Parameters
275 ----------
276 frf_matrix : NDArray
277 Transfer function as an np.ndarray, should be organized such that the
278 frequency is the first axis, the responses are the second axis, and the
279 references are the third axis
280 num_retained_values : float or np.ndarray, optional
281 Number of singular values to retain in the pseudo-inverse when the truncation
282 method is chosen. A vector of can be provided so the number of retained values
283 can change as a function of frequency. The vector must match the length of the
284 abscissa in this case.
286 Returns
287 -------
288 np.ndarray
289 Inverse of the supplied FRF matrix
291 Raises
292 ------
293 Exception
294 If a number of retained values isn't supplied
295 Warning
296 if the number of retained values is greater than the actual number of singular
297 values
298 """
299 if num_retained_values is None:
300 raise ValueError('The number of retained values must be supplied when using the truncation method')
302 num_retained_values = np.asarray(num_retained_values, dtype=np.intc)
304 # Need to add a second axis to the number of retained values if it's a vector
305 # so it broadcasts appropriately.
306 if num_retained_values.size > 1 and num_retained_values.ndim == 1:
307 num_retained_values = num_retained_values[..., np.newaxis]
309 if num_retained_values.size > 1 and num_retained_values.shape[1] > 1:
310 raise ValueError('The number of retained values parameter is not a scalar or vector and cannot be interpreted')
312 if num_retained_values.size > 1 and num_retained_values.size != frf_matrix.shape[0]:
313 raise ValueError('The number of retained values parameter is a vector but '
314 + 'the size does not match the number of frequency lines and cannot be broadcasted to the FRFs')
316 u, s, vh = np.linalg.svd(frf_matrix, full_matrices=False)
317 vh_conjT = np.swapaxes(vh.conj(), -2, -1)
318 u_conjT = np.swapaxes(u.conj(), -2, -1)
320 # The weird looking logic is to handle if num_retained_values is either a string
321 # or a scaler. It also changes the number of retained values in the case that
322 # it is greater than the actual number of singular values (it changes to retain
323 # all the singular values).
324 if num_retained_values.size > 1 and any(num_retained_values > s.shape[1]) or num_retained_values.size == 1 and num_retained_values > s.shape[1]:
325 warnings.warn('The requested number of retained singular values is greater than the '
326 + 'actual number of singular values (for at least one frequency line). All singular values are being retained.')
327 num_retained_values[num_retained_values > s.shape[1]] = s.shape[1]
329 s_truncated_inv = np.zeros((s.shape[0], s.shape[1], s.shape[1]))
330 num_lines = s.shape[0]
332 # Making num_retained_values a vector if an integer was supplied
333 if not num_retained_values.size > 1:
334 num_retained_values = (np.ones(num_lines) * num_retained_values)[..., np.newaxis]
336 num_retained_values = num_retained_values.astype(int)
338 for ii in range(num_lines):
339 if np.any(s[ii, :num_retained_values[ii, 0]] == 0):
340 continue
341 else:
342 s_truncated_inv[ii, :num_retained_values[ii, 0], :num_retained_values[ii, 0]] = np.diag(1/s[ii, :num_retained_values[ii, 0]])
344 frf_pinv = vh_conjT@s_truncated_inv@u_conjT
346 return frf_pinv
349def compute_tikhonov_modified_singular_values(original_singular_values, regularization_parameter):
350 """
351 Computes the modified singular values that would be seen in Tikhonov
352 regularization. This is only intended for visualization purposes
353 only.
355 Parameters
356 ----------
357 original_singular_values : np.ndarray
358 The singular values from the oringinal FRF matrix. Should be organized
359 [frequency lines, singular values.]
360 regularization_parameter : float or np.ndarray
361 The regularization parameter that would be used in an inverse
362 via Tikhonov regularization.
364 Returns
365 -------
366 modified_singular_values : np.ndarray
367 The modified singular values from the Tikhonov regularization.
369 Notes
370 -----
371 This is not a mathematically rigorous version of the modified singular
372 values and should only be used for subjective evaluation of the effect of
373 the regularization parameter.
374 """
375 regularization_parameter = np.asarray(regularization_parameter, dtype=np.float64)
377 if regularization_parameter.size != 1 and regularization_parameter.ndim == 1:
378 regularization_parameter = regularization_parameter[..., np.newaxis]
380 modified_singular_values = (original_singular_values**2 + regularization_parameter)/original_singular_values
382 return modified_singular_values