Coverage for src / sdynpy / modal / sdynpy_modeshape.py: 6%

1# -*- coding: utf-8 -*- 

2""" 

3Functions for computing shapes after the poles are known 

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. 

9 

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. 

14 

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. 

19 

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""" 

23 

24from enum import Enum 

25import numpy as np 

26from ..core.sdynpy_data import TransferFunctionArray 

27from ..core.sdynpy_shape import shape_array 

28from ..core.sdynpy_coordinate import CoordinateArray 

29from ..signal_processing.sdynpy_complex import collapse_complex_to_real 

30import warnings 

31 

32 

33class ShapeSelection(Enum): 

34 ALL = 0 

35 DRIVE_POINT_COEFFICIENT = 1 

36 PARTICIPATION_FACTOR = 2 

37 

38 

39def compute_residues(experimental_frf: TransferFunctionArray, 

40 natural_frequencies: np.ndarray, 

41 damping_ratios: np.ndarray, 

42 real_modes: bool = False, 

43 residuals: bool = True, 

44 min_frequency: float = None, 

45 max_frequency: float = None, 

46 weighting: np.ndarray = 'uniform', 

47 displacement_derivative: int = 0, 

48 frequency_lines_at_resonance: int = None, 

49 frequency_lines_for_residuals: int = None): 

50 """ 

51 Fit residues to FRF data given frequency and damping values 

52 

53 Parameters 

54 ---------- 

55 experimental_frf : TransferFunctionArray 

56 Experimental FRF data to which modes will be fit 

57 natural_frequencies : np.ndarray 

58 Natural Frequencies (in Hz) at which modes will be fit 

59 damping_ratios : np.ndarray 

60 Damping Ratios at which modes will be fit 

61 real_modes : bool, optional 

62 If true, fit residues will be real-valued. False allows complex modes. 

63 The default is False. 

64 residuals : bool, optional 

65 Use residuals in the FRF fit. The default is True. 

66 min_frequency : float, optional 

67 Minimum frequency to use in the shape fit. The default is the lowest 

68 frequency in the experimental FRF. 

69 max_frequency : float, optional 

70 Maximum frequency to use in the shape fit. The default is the highest 

71 frequency in the experimental FRF. 

72 weighting : np.ndarray or string, optional 

73 A weighting array to use to fit shapes better at specific frequencies. 

74 The default is weighted by the log magnitude of the FRF matrix. Can be 

75 defined as 'magnitude','uniform', or an ndarray with shape identical 

76 to the ordinate of the experimental frfs 

77 displacement_derivative : int, optional 

78 Defines the type of data in the FRF based on the number of derivatives 

79 from displacement (0 - displacement, 1 - velocity, 2 - acceleration). 

80 The default is 0 (displacement). 

81 frequency_lines_at_resonance : int, optional 

82 Defines the number of frequency lines to look at around the specified 

83 natural frequencies for computing residues. If not specified, all 

84 frequency lines are used for computing shapes. 

85 frequency_lines_for_residuals : int, optional 

86 Defines the number of frequency lines at the low and high frequency to 

87 use in computing shapes. Only used if frequency_lines_at_resonance is 

88 specified. If not specified, the lower 10% and upper 10% of frequency 

89 lines will be kept. 

90 

91 Returns 

92 ------- 

93 shape_residues : np.ndarray 

94 A (..., n_modes) shaped np.ndarray where ... is the shape of the input 

95 experimental_frf array. There will be one residue for each experimental 

96 frf (reference and response) for each mode. 

97 synthesized_frf : TransferFunctionArray 

98 Transfer function array containing the analytical fits using the 

99 residues. 

100 residual_frf : TransferFunctionArray 

101 Transfer function array containing the residual data from the analytical 

102 fits. 

103 """ 

104 flat_frf = experimental_frf.flatten() 

105 frequencies = flat_frf[0].abscissa.copy() 

106 if min_frequency is None: 

107 min_frequency = np.min(frequencies) 

108 if max_frequency is None: 

109 max_frequency = np.max(frequencies) 

110 abscissa_indices = np.ones(frequencies.shape, dtype=bool) 

111 abscissa_indices &= (frequencies >= min_frequency) 

112 abscissa_indices &= (frequencies <= max_frequency) 

113 frequencies = frequencies[abscissa_indices] 

114 frf_matrix = flat_frf.ordinate[:, abscissa_indices].T.copy() 

115 angular_frequencies = 2 * np.pi * frequencies[:, np.newaxis] 

116 angular_natural_frequencies = 2 * np.pi * np.array(natural_frequencies).flatten() 

117 damping_ratios = np.array(damping_ratios).flatten() 

118 

119 # Reduce to the kept frequency lines 

120 if frequency_lines_at_resonance is not None: 

121 solve_indices = np.argmin(np.abs(angular_natural_frequencies - angular_frequencies), axis=0) 

122 # print(solve_indices) 

123 solve_indices = np.unique( 

124 solve_indices[:, np.newaxis] + np.arange(frequency_lines_at_resonance) - frequency_lines_at_resonance // 2) 

125 solve_indices = solve_indices[(solve_indices >= 0) & ( 

126 solve_indices < angular_frequencies.size)] 

127 # Add the residual indices 

128 if residuals: 

129 if frequency_lines_for_residuals is None: 

130 low_freq_indices = np.arange(angular_frequencies.size // 10) 

131 high_freq_indices = angular_frequencies.size - \ 

132 np.arange(angular_frequencies.size // 10) - 1 

133 else: 

134 low_freq_indices = np.arange(frequency_lines_for_residuals) 

135 high_freq_indices = angular_frequencies.size - \ 

136 np.arange(frequency_lines_for_residuals) - 1 

137 solve_indices = np.unique(np.concatenate( 

138 (solve_indices, low_freq_indices, high_freq_indices))) 

139 kernel_indices = np.concatenate((solve_indices, solve_indices + angular_frequencies.size)) 

140 # print(kernel_indices) 

141 

142 # Set up the kernel to solve the least squares residue problem 

143 denominator = ((angular_natural_frequencies**2 - angular_frequencies**2)**2 

144 + (2 * damping_ratios * angular_natural_frequencies * angular_frequencies)**2) 

145 kernel_rr = (angular_natural_frequencies**2 - angular_frequencies**2) / denominator 

146 kernel_ri = 2 * damping_ratios * angular_natural_frequencies * angular_frequencies / denominator 

147 kernel_ir = -2 * damping_ratios * angular_natural_frequencies * angular_frequencies / denominator 

148 kernel_ii = (angular_natural_frequencies**2 - angular_frequencies**2) / denominator 

149 low_frequency_residual = (-1 / angular_frequencies**2) 

150 high_frequency_residual = (np.ones(angular_frequencies.shape)) 

151 zeros = np.zeros(angular_frequencies.shape) 

152 

153 kernel = np.concatenate((np.concatenate((kernel_rr, kernel_ri, low_frequency_residual, zeros, high_frequency_residual, zeros), axis=-1), 

154 np.concatenate((kernel_ir, kernel_ii, zeros, low_frequency_residual, zeros, high_frequency_residual), axis=-1)), axis=0) 

155 

156 # print(kernel.shape) 

157 

158 # Reduce kernel depending on whether or not we want real modes or complex 

159 # modes 

160 if real_modes: 

161 kernel[..., angular_natural_frequencies.size:2 * angular_natural_frequencies.size] = 0 

162 

163 # Reduce kernel depending on whether or not we want residuals 

164 if not residuals: 

165 kernel[..., -4:] = 0 

166 

167 # Perform the solution in acceleration space 

168 kernel *= np.tile(-angular_frequencies**2, (2, 1)) 

169 frf_matrix *= (1j * angular_frequencies)**(2 - displacement_derivative) 

170 

171 # print(kernel.shape) 

172 

173 # Weighting matrix 

174 if isinstance(weighting, str): 

175 if weighting.lower() == 'uniform': 

176 weighting = np.ones(frf_matrix.shape) 

177 elif weighting.lower() == 'magnitude': 

178 max_frf = np.max(np.abs(frf_matrix), axis=0, keepdims=True) 

179 max_frf[max_frf == 0.0] = 1.0 

180 weighting = np.abs(frf_matrix) / max_frf 

181 min_weighting = np.min(weighting, axis=0, keepdims=True) 

182 min_weighting[min_weighting == 0] = 1 

183 weighting_ratio = weighting / min_weighting 

184 weighting_ratio[weighting_ratio == 0] = 1 

185 weighting = np.log10(weighting_ratio) 

186 else: 

187 raise ValueError( 

188 'If weighting is specified as a string, must be "uniform" or "magnitude"') 

189 

190 # Now assemble the FRF matrix 

191 frf_matrix_to_fit = np.concatenate((frf_matrix.real, 

192 frf_matrix.imag), axis=0) 

193 

194 # Now perform the weighting 

195 weighting = np.tile(weighting, (2, 1)) 

196 weighted_kernel = weighting.T[..., np.newaxis] * kernel 

197 weighted_frf_matrix_to_fit = (weighting * frf_matrix_to_fit).T[..., np.newaxis] 

198 

199 # print(weighted_kernel.shape) 

200 # print(weighted_frf_matrix_to_fit.shape) 

201 

202 if frequency_lines_at_resonance is not None: 

203 weighted_kernel = weighted_kernel[:, kernel_indices] 

204 weighted_frf_matrix_to_fit = weighted_frf_matrix_to_fit[:, kernel_indices] 

205 

206 # now solve 

207 residues = (np.linalg.pinv(weighted_kernel) @ weighted_frf_matrix_to_fit).squeeze().T 

208 

209 frf_fit_matrix = kernel @ residues 

210 frf_fit_matrix = (frf_fit_matrix[:frf_fit_matrix.shape[0] // 2, :] 

211 + 1j * frf_fit_matrix[frf_fit_matrix.shape[0] // 2:, :]) 

212 output_frf = flat_frf.extract_elements(abscissa_indices) 

213 output_frf.ordinate = (frf_fit_matrix / (1j * angular_frequencies) 

214 ** (2 - displacement_derivative)).T 

215 output_frf = output_frf.reshape(experimental_frf.shape) 

216 

217 # Extract shape residues and residuals 

218 shape_residues = (residues[:angular_natural_frequencies.size] 

219 + 1j * residues[angular_natural_frequencies.size:2 * angular_natural_frequencies.size]) 

220 shape_residues = np.moveaxis(shape_residues.reshape(-1, *experimental_frf.shape), 0, -1) 

221 if real_modes: 

222 shape_residues = np.real(shape_residues) 

223 

224 # Extract residuals 

225 residual_frf = flat_frf.extract_elements(abscissa_indices) 

226 residual_matrix = kernel[:, -4:] @ residues[-4:] 

227 residual_matrix = (residual_matrix[:residual_matrix.shape[0] // 2, :] 

228 + 1j * residual_matrix[residual_matrix.shape[0] // 2:, :]) 

229 residual_frf.ordinate = (residual_matrix / (1j * angular_frequencies) 

230 ** (2 - displacement_derivative)).T 

231 residual_frf = residual_frf.reshape(experimental_frf.shape) 

232 

233 return shape_residues, output_frf, residual_frf 

234 

235 

236def compute_shapes(natural_frequencies: np.ndarray, 

237 damping_ratios: np.ndarray, 

238 coordinates: CoordinateArray, 

239 residue_matrix: np.ndarray, 

240 shape_selection=ShapeSelection.ALL, 

241 participation_factors: np.ndarray = None): 

242 abs_coordinates = abs(coordinates) 

243 sign_matrix = np.prod(coordinates.sign(), axis=-1) 

244 equality_matrix = np.all(abs_coordinates == abs_coordinates[..., 0, np.newaxis], axis=-1) 

245 drive_point_indices = np.argmax(equality_matrix, axis=0) 

246 residue_scales = np.array([sign_matrix[response_index, reference_index] 

247 for reference_index, response_index in enumerate(drive_point_indices)]) 

248 signed_residue_matrix = residue_matrix * residue_scales[:, np.newaxis] 

249 drive_point_residues = np.array([signed_residue_matrix[response_index, reference_index, :] 

250 for reference_index, response_index in enumerate(drive_point_indices)]) 

251 negative_drive_points = np.where(drive_point_residues < 0) 

252 if np.isrealobj(residue_matrix) and np.any(drive_point_residues < 0): 

253 print('Negative Drive Point Residues Found!') 

254 for mode_index, reference_index in np.sort(np.array(negative_drive_points).T[..., ::-1], axis=0): 

255 print(' Mode {:} Reference {:}'.format( 

256 mode_index + 1, str(coordinates[0, reference_index, -1]))) 

257 mode_shape_scaling = np.sqrt(np.abs(drive_point_residues) if np.isrealobj( 

258 residue_matrix) else drive_point_residues) 

259 mode_shape_matrix = signed_residue_matrix / mode_shape_scaling 

260 if isinstance(shape_selection, str): 

261 if shape_selection.lower() in ['all']: 

262 shape_selection = ShapeSelection.ALL 

263 elif shape_selection.lower() in ['drive', 'drive point', 'dp']: 

264 shape_selection = ShapeSelection.DRIVE_POINT_COEFFICIENT 

265 elif shape_selection.lower() in ['part', 'participation', 'participation factor']: 

266 shape_selection = ShapeSelection.PARTICIPATION_FACTOR 

267 if not shape_selection == ShapeSelection.ALL: 

268 if shape_selection == ShapeSelection.DRIVE_POINT_COEFFICIENT: 

269 shape_selection_indices = np.argmax(drive_point_residues, axis=0) 

270 elif shape_selection == ShapeSelection.PARTICIPATION_FACTOR: 

271 shape_selection_indices = np.argmax(np.abs(participation_factors), axis=-1) 

272 else: 

273 raise ValueError('Invalid Shape Selection Technique') 

274 mode_shape_matrix = np.array([mode_shape_matrix[:, reference_index, mode_index] 

275 for mode_index, reference_index in enumerate(shape_selection_indices)]).T 

276 else: 

277 shape_selection_indices = np.arange(mode_shape_matrix.shape[1])[ 

278 :, np.newaxis] * np.ones(mode_shape_matrix.shape[2], dtype=int) 

279 ref_array = np.ndarray 

280 shapes = shape_array(coordinate=coordinates[..., 0, 0], 

281 shape_matrix=np.moveaxis(mode_shape_matrix, 0, -1), 

282 frequency=natural_frequencies, damping=damping_ratios, 

283 comment1=coordinates[0, :, -1].string_array()[shape_selection_indices]) 

284 return shapes, negative_drive_points 

285 

286 

287def compute_shapes_multireference(experimental_frf: TransferFunctionArray, 

288 natural_frequencies: np.ndarray, 

289 damping_ratios: np.ndarray, 

290 participation_factors: np.ndarray, 

291 real_modes: bool = False, 

292 lower_residuals: bool = True, 

293 upper_residuals: bool = True, 

294 min_frequency: float = None, 

295 max_frequency: float = None, 

296 displacement_derivative: int = 0, 

297 frequency_lines_at_resonance: int = None, 

298 frequency_lines_for_residuals: int = None): 

299 """ 

300 Computes mode shapes from multireference datasets. 

301 

302 Uses the modal participation factor as a constraint on the mode shapes to 

303 solve for the shapes in one pass, rather than solving for residues and 

304 subsequently solving for shapes. 

305 

306 Parameters 

307 ---------- 

308 experimental_frf : TransferFunctionArray 

309 Experimental FRF data to which modes will be fit 

310 natural_frequencies : np.ndarray 

311 Natural Frequencies (in Hz) at which modes will be fit 

312 damping_ratios : np.ndarray 

313 Damping Ratios at which modes will be fit 

314 participation_factors : np.ndarray 

315 Mode participation factors from which the shapes can be computed. 

316 Should have shape (n_modes x n_inputs) 

317 real_modes : bool, optional 

318 Specifies whether to solve for real modes or complex modes (default). 

319 lower_residuals : bool, optional 

320 Use lower residuals in the FRF fit. The default is True. 

321 upper_residuals : bool, optional 

322 Use upper residuals in the FRF fit. The default is True. 

323 min_frequency : float, optional 

324 Minimum frequency to use in the shape fit. The default is the lowest 

325 frequency in the experimental FRF. 

326 max_frequency : float, optional 

327 Maximum frequency to use in the shape fit. The default is the highest 

328 frequency in the experimental FRF. 

329 displacement_derivative : int, optional 

330 Defines the type of data in the FRF based on the number of derivatives 

331 from displacement (0 - displacement, 1 - velocity, 2 - acceleration). 

332 The default is 0 (displacement). 

333 frequency_lines_at_resonance : int, optional 

334 Defines the number of frequency lines to look at around the specified 

335 natural frequencies for computing residues. If not specified, all 

336 frequency lines are used for computing shapes. 

337 frequency_lines_for_residuals : int, optional 

338 Defines the number of frequency lines at the low and high frequency to 

339 use in computing shapes. Only used if frequency_lines_at_resonance is 

340 specified. If not specified, the lower 10% and upper 10% of frequency 

341 lines will be kept for computing residuals. 

342 

343 Raises 

344 ------ 

345 ValueError 

346 If the FRF is not 2-dimensional with references on the columns and 

347 responses on the rows. 

348 

349 Returns 

350 ------- 

351 output_shape : ShapeArray 

352 ShapeArray containing the mode shapes of the system 

353 frfs_resynthesized : TransferFunctionArray 

354 FRFs resynthesized from the fit shapes and residuals 

355 residual_frfs : TransferFunctionArray 

356 FRFs resynthesized only from the residuals used in the calculation 

357 kernel_frfs: TransferFunctionArray 

358 FRFs synthesized by the kernel solution without taking into account 

359 construction of mode shapes and reciprocity with participation factors 

360 """ 

361 original_coordinates = experimental_frf.coordinate 

362 if not experimental_frf.ndim == 2: 

363 raise ValueError('FRF must be shaped n_outputs x n_inputs') 

364 abs_coordinate = abs(original_coordinates) 

365 experimental_frf = experimental_frf[abs_coordinate] 

366 # Also need to adjust the participation factors 

367 participation_factors = participation_factors*np.sign(original_coordinates[0, :, 1].direction) 

368 if real_modes: 

369 participation_factors = collapse_complex_to_real(participation_factors) 

370 frequencies = experimental_frf[0, 0].abscissa.copy() 

371 if min_frequency is None: 

372 min_frequency = np.min(frequencies) 

373 if max_frequency is None: 

374 max_frequency = np.max(frequencies) 

375 abscissa_indices = np.ones(frequencies.shape, dtype=bool) 

376 abscissa_indices &= (frequencies >= min_frequency) 

377 abscissa_indices &= (frequencies <= max_frequency) 

378 frequencies = frequencies[abscissa_indices] 

379 frf_matrix = experimental_frf.ordinate[..., abscissa_indices].copy() 

380 angular_frequencies = 2 * np.pi * frequencies 

381 reconstruction_angular_frequencies = frequencies*2*np.pi 

382 angular_natural_frequencies = 2 * np.pi * np.array(natural_frequencies).flatten() 

383 damping_ratios = np.array(damping_ratios).flatten() 

384 

385 # Reduce to the kept frequency lines 

386 if frequency_lines_at_resonance is not None: 

387 solve_indices = np.argmin(np.abs(angular_natural_frequencies - angular_frequencies[:, np.newaxis]), axis=0) 

388 # print(solve_indices) 

389 solve_indices = np.unique( 

390 solve_indices[:, np.newaxis] + np.arange(frequency_lines_at_resonance) - frequency_lines_at_resonance // 2) 

391 solve_indices = solve_indices[(solve_indices >= 0) & ( 

392 solve_indices < angular_frequencies.size)] 

393 # Add the residual indices 

394 if lower_residuals: 

395 if frequency_lines_for_residuals is None: 

396 low_freq_indices = np.arange(angular_frequencies.size // 10) 

397 else: 

398 low_freq_indices = np.arange(frequency_lines_for_residuals) 

399 else: 

400 low_freq_indices = np.array([],dtype=int) 

401 if upper_residuals: 

402 if frequency_lines_for_residuals is None: 

403 high_freq_indices = angular_frequencies.size - \ 

404 np.arange(angular_frequencies.size // 10) - 1 

405 else: 

406 high_freq_indices = angular_frequencies.size - \ 

407 np.arange(frequency_lines_for_residuals) - 1 

408 else: 

409 high_freq_indices = np.array([],dtype=int) 

410 solve_indices = np.unique(np.concatenate( 

411 (solve_indices, low_freq_indices, high_freq_indices))) 

412 print(solve_indices) 

413 frf_matrix = frf_matrix[..., solve_indices] 

414 angular_frequencies = angular_frequencies[solve_indices] 

415 

416 poles = -damping_ratios*angular_natural_frequencies + 1j*np.sqrt(1-damping_ratios**2)*angular_natural_frequencies 

417 

418 if real_modes: 

419 kernel = generate_kernel_real( 

420 angular_frequencies, poles, 

421 participation_factors, lower_residuals, 

422 upper_residuals, displacement_derivative) 

423 if angular_frequencies.size != reconstruction_angular_frequencies.size: 

424 reconstruction_kernel = generate_kernel_real( 

425 reconstruction_angular_frequencies, poles, 

426 participation_factors, lower_residuals, 

427 upper_residuals, displacement_derivative) 

428 else: 

429 reconstruction_kernel = kernel 

430 else: 

431 kernel = generate_kernel_complex( 

432 angular_frequencies, poles, participation_factors, lower_residuals, 

433 upper_residuals, displacement_derivative) 

434 if angular_frequencies.size != reconstruction_angular_frequencies.size: 

435 reconstruction_kernel = generate_kernel_complex( 

436 reconstruction_angular_frequencies, poles, participation_factors, lower_residuals, 

437 upper_residuals, displacement_derivative) 

438 else: 

439 reconstruction_kernel = kernel 

440 

441 coefficients, (full_reconstruction, residual_reconstruction) = stack_and_lstsq( 

442 frf_matrix, kernel, 

443 return_reconstruction=True, 

444 reconstruction_partitions = [slice(None), 

445 slice(poles.size*(1 if real_modes else 2),None)], 

446 reconstruction_kernel = reconstruction_kernel) 

447 

448 residual_frfs = experimental_frf.copy().extract_elements(abscissa_indices) 

449 residual_frfs.ordinate = residual_reconstruction 

450 kernel_frfs = experimental_frf.copy().extract_elements(abscissa_indices) 

451 kernel_frfs.ordinate = full_reconstruction 

452 

453 # Extract the shapes and residues 

454 if real_modes: 

455 shapes = (coefficients[:poles.size]).T 

456 else: 

457 shapes = (coefficients[:poles.size] + 1j*coefficients[poles.size:2*poles.size]).T 

458 

459 # Now we have to go in and find the scale factor to scale the shapes correctly 

460 drive_points = np.where(experimental_frf.response_coordinate == experimental_frf.reference_coordinate) 

461 if drive_points[0].size == 0: 

462 print('Warning: Drive Points Not Found in Dataset, Shapes are Unscaled.') 

463 output_shape = shape_array(abs_coordinate[:, 0, 0], shapes.T, 

464 angular_natural_frequencies/(2*np.pi), 

465 damping_ratios, 

466 1) 

467 frfs_resynthesized = kernel_frfs 

468 else: 

469 shapes_normalized = [] 

470 drive_residues = shapes[drive_points[0]].T*participation_factors[:,drive_points[1]] 

471 for k, drive_residue in enumerate(drive_residues): 

472 if real_modes: 

473 # Throw away negative drive points 

474 bad_indices = drive_residue < 0 

475 else: 

476 # Throw away non-negative imaginary parts 

477 bad_indices_positive = drive_residue.imag > 0 

478 # Throw away small values compared to the average value (only considering negative imaginary parts) 

479 bad_indices_small = np.abs(drive_residue) < 0.01*np.mean(np.abs(drive_residue[~bad_indices_positive])) 

480 # Combine into a single criteria 

481 bad_indices = bad_indices_positive | bad_indices_small 

482 # Get the good indices that are remaining 

483 remaining_indices = np.where(~bad_indices)[0] 

484 if len(remaining_indices) == 0: 

485 print('Warning: Mode {:} had no valid drive points and is therefore unscaled.'.format(k+1)) 

486 shapes_normalized.append(shapes[:,k]) 

487 continue 

488 # We will then construct the least squares solution 

489 shape_coefficients = shapes[drive_points[0][remaining_indices],k][:,np.newaxis] 

490 residue_coefficients = np.sqrt(drive_residue[remaining_indices])[:,np.newaxis] 

491 # Before we compute the scale, we need to make sure that we have all of the signs the same way. 

492 # This is because the square root can give you +/- root where root**2 = complex number 

493 # This will mess up the least squares at it will try to find something between the 

494 # two vectors. 

495 scale_vector = (residue_coefficients/shape_coefficients).flatten() 

496 sign_vector = np.array((scale_vector.real,scale_vector.imag)) 

497 # Get the signs 

498 signs = np.sign(np.dot(sign_vector[:,0],sign_vector)) 

499 residue_coefficients = residue_coefficients*signs[:,np.newaxis] 

500 # Now compute the least-squares solution 

501 scale = np.linalg.lstsq(shape_coefficients,residue_coefficients)[0].squeeze() 

502 print('Scale for mode {:}: {:}'.format(k+1,scale)) 

503 shapes_normalized.append(shapes[:,k]*scale) 

504 shapes_normalized = np.array(shapes_normalized).T 

505 

506 output_shape = shape_array(abs_coordinate[:, 0, 0], shapes_normalized.T, 

507 angular_natural_frequencies/(2*np.pi), 

508 damping_ratios, 

509 1) 

510 

511 frfs_resynthesized = (output_shape.compute_frf( 

512 frequencies, 

513 responses=experimental_frf.response_coordinate[:, 0], 

514 references=experimental_frf.reference_coordinate[0, :], 

515 displacement_derivative=displacement_derivative) 

516 + residual_frfs) 

517 

518 frfs_resynthesized = frfs_resynthesized[original_coordinates] 

519 residual_frfs = residual_frfs[original_coordinates] 

520 kernel_frfs = kernel_frfs[original_coordinates] 

521 

522 return output_shape, frfs_resynthesized, residual_frfs, kernel_frfs 

523 

524def generate_kernel_complex(omegas, poles, participation_factors, 

525 lower_residuals = False, upper_residuals = False, 

526 displacement_derivative = 0): 

527 """ 

528  

529 

530 Parameters 

531 ---------- 

532 omegas : np.ndarray 

533 The angular frequencies (in radians/s) at which the kernel matrix will 

534 be computed. This should be a 1D array with length num_freqs. 

535 poles : np.ndarray 

536 An array of poles corresponding to the modes of the structure. This 

537 should be a 1D array with length num_modes. 

538 participation_factors : np.ndarray 

539 A 2D array of participation factors corresponding to the reference 

540 degrees of freedom. This should have shape (num_modes, num_inputs). 

541 lower_residuals : bool, optional 

542 If True, construct the kernel matrix such that lower residuals will be 

543 computed in the least-squares operation. The default is False. 

544 upper_residuals : bool, optional 

545 If True, construct the kernel matrix such that upper residuals will be 

546 computed in the least-squares operation. The default is False. 

547 displacement_derivative : int, optional 

548 The derivative of displacement used to construct the frequency response 

549 functions. Should be 0 for a receptance (displacement/force) frf, 1 

550 for a mobility (velocity/force) frf, or 2 for an accelerance 

551 (acceleration/force) frf. 

552 

553 Returns 

554 ------- 

555 kernel_matrix : np.ndarray 

556 A 3D matrix that represents the kernel that can be inverted to solve 

557 for mode shapes. The size of the output array will be num_inputs x 

558 num_freq*2 x num_modes*2. The top partition of the num_freq dimension 

559 corresponds to the real part of the frf, and the 

560 bottom portion corresponds to the imaginary part of the frf. The top 

561 partition of the num_modes dimension corresponds to the real part of 

562 the mode shape matrix, and the bottom partition corresponds to the 

563 imaginary part. If residuals are included, there will be an extra two 

564 entries along the num_modes dimension corresponding to real and 

565 imaginary parts of the residual matrix for each of the residuals 

566 included. 

567 """ 

568 omegas = np.array(omegas) 

569 poles = np.array(poles) 

570 participation_factors = np.array(participation_factors) 

571 if omegas.ndim != 1: 

572 raise ValueError('`omegas` should be 1-dimensional') 

573 if poles.ndim != 1: 

574 raise ValueError('`poles` should be 1-dimensional') 

575 if participation_factors.ndim != 2: 

576 raise ValueError('`participation_factors` should be 2-dimensional (num_modes x num_inputs)') 

577 n_modes,n_inputs = participation_factors.shape 

578 if poles.size != n_modes: 

579 raise ValueError('number of poles ({:}) is not equal to the number of rows in the participation_factors array ({:})'.format( 

580 poles.size,n_modes)) 

581 # Set up broadcasting 

582 # We want the output array to be n_input x n_freq*2 x n_modes*2 

583 # So let's adjust the terms so they have the right shapes 

584 # We want frequency lines to be the middle dimension 

585 omega = omegas[np.newaxis,:,np.newaxis] 

586 # We want inputs to be the first dimension and modes the 

587 # last dimension 

588 participation_factors = participation_factors.T[:,np.newaxis,:] 

589 # Split up terms into real and imaginary parts 

590 p_r = poles.real 

591 p_i = poles.imag 

592 l_r = participation_factors.real 

593 l_i = participation_factors.imag 

594 

595 # Now go through the different derivatives and compute the kernel plus any 

596 # residuals 

597 if displacement_derivative == 0: 

598 kernel = np.array([ 

599 [(-l_i*omega - l_i*p_i - l_r*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

600 + (l_i*omega - l_i*p_i - l_r*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2), 

601 (l_i*p_r - l_r*omega - l_r*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

602 + (l_i*p_r + l_r*omega - l_r*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2)], 

603 [(-l_i*p_r - l_r*omega + l_r*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2) 

604 + (l_i*p_r - l_r*omega - l_r*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2), 

605 (l_i*omega - l_i*p_i - l_r*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2) 

606 + (l_i*omega + l_i*p_i + l_r*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2)]]) 

607 if lower_residuals: 

608 kernel_RL = np.concatenate([ 

609 np.kron(-1/omegas[:,np.newaxis]**2,np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0])), 

610 np.kron(-1/omegas[:,np.newaxis]**2,np.kron(np.eye(n_inputs)[:,np.newaxis],[0,1]))],axis=1) 

611 if upper_residuals: 

612 kernel_RU = np.concatenate([ 

613 np.kron(np.ones((omegas.size,1)),np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0])), 

614 np.kron(np.ones((omegas.size,1)),np.kron(np.eye(n_inputs)[:,np.newaxis],[0,1]))],axis=1) 

615 elif displacement_derivative == 1: 

616 kernel = np.array([ 

617 [(-l_i*omega*p_r + l_r*omega**2 + l_r*omega*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

618 + (l_i*omega*p_r + l_r*omega**2 - l_r*omega*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2), 

619 (-l_i*omega**2 - l_i*omega*p_i - l_r*omega*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

620 + (-l_i*omega**2 + l_i*omega*p_i + l_r*omega*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2)], 

621 [(-l_i*omega**2 - l_i*omega*p_i - l_r*omega*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

622 + (l_i*omega**2 - l_i*omega*p_i - l_r*omega*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2), 

623 (l_i*omega*p_r - l_r*omega**2 - l_r*omega*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

624 + (l_i*omega*p_r + l_r*omega**2 - l_r*omega*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2)]]) 

625 if lower_residuals: 

626 kernel_RL = np.concatenate([ 

627 np.kron( 1/omegas[:,np.newaxis],np.kron(np.eye(n_inputs)[:,np.newaxis],[0,1])), 

628 np.kron(-1/omegas[:,np.newaxis],np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0]))],axis=1) 

629 if upper_residuals: 

630 kernel_RU = np.concatenate([ 

631 np.kron(-omegas[:,np.newaxis],np.kron(np.eye(n_inputs)[:,np.newaxis],[0,1])), 

632 np.kron( omegas[:,np.newaxis],np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0]))],axis=1) 

633 elif displacement_derivative == 2: 

634 kernel = np.array([ 

635 [(-l_i*omega**3 + l_i*omega**2*p_i + l_r*omega**2*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2) 

636 + (l_i*omega**3 + l_i*omega**2*p_i + l_r*omega**2*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2), 

637 (-l_i*omega**2*p_r - l_r*omega**3 + l_r*omega**2*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2) 

638 + (-l_i*omega**2*p_r + l_r*omega**3 + l_r*omega**2*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2)], 

639 [(-l_i*omega**2*p_r + l_r*omega**3 + l_r*omega**2*p_i)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

640 + (l_i*omega**2*p_r + l_r*omega**3 - l_r*omega**2*p_i)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2), 

641 (-l_i*omega**3 - l_i*omega**2*p_i - l_r*omega**2*p_r)/(omega**2 + 2*omega*p_i + p_i**2 + p_r**2) 

642 + (-l_i*omega**3 + l_i*omega**2*p_i + l_r*omega**2*p_r)/(omega**2 - 2*omega*p_i + p_i**2 + p_r**2)]]) 

643 if lower_residuals: 

644 kernel_RL = np.concatenate([ 

645 np.kron(np.ones((omegas.size,1)),np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0])), 

646 np.kron(np.ones((omegas.size,1)),np.kron(np.eye(n_inputs)[:,np.newaxis],[0,1]))],axis=1) 

647 if upper_residuals: 

648 kernel_RU = np.concatenate([ 

649 np.kron(-omegas[:,np.newaxis]**2,np.kron(np.eye(n_inputs)[:,np.newaxis],[1,0])),