使用pyFRF计算锤击试验FRF时结果异常及报错的排查求助

我正在开展移动锤击试验的实验模态诊断数据分析，已将时域加速度与外力数据通过FFT转换至频域，所得频域图形符合预期。但使用pyFRF包计算并绘制FRF（频响函数）时，生成的图形不符合预期，还出现了运行警告。

所用代码 #

from pyFRF import FRF
import pandas as pd
import numpy as np
import scipy
import matplotlib.pyplot as plt
import pyExSi

# data import
skiprow_list = [0,2,3,4,5,6,7,8,9,10]
column_list = np.arange(1,15)

data = pd.read_csv("2024_01_03_WTB_drgania/a12.CSV",
                    delimiter=';',
                    decimal=',',
                    skiprows=skiprow_list,
                    usecols=column_list)

t = data[['C1:  Time[s]']]                  # time
a = data[['C1:  [m/s^2]']]                  # acceleration (acc sensor)
p = data[['C7:  [N]']]                      # force (hammer)


# date frame to 1D arrays
t = t.values.flatten()
a = a.values.flatten()
p = p.values.flatten()


N = len(t)
dt = (t[N-1] - t[0])/N                      # time step
fs = 1/dt                                   # sampling frequency
f0 = fs/N 

freq_syn = np.arange(0.0, 2000, f0)
omega = 2 * np.pi * freq_syn


# create an object
frf_object = FRF(sampling_freq=int(fs), exc=p, resp=a, resp_type='a', frf_type='H1')

# get frf
H1_pyFRF = frf_object.get_FRF()  # returns the requested FRF estimator matrix of shape (response DOF, excitation DOF, frequency points)
freq_pyFRF= frf_object.get_f_axis()  # returns frequency series
H1_pyFRF.shape  # single input, single output (SISO) system


plt.semilogy(freq_pyFRF[1:], np.abs(H1_pyFRF)[0,0,1:], color="b", label='pyFRF')
plt.xlim(left=0, right=2000)
plt.xlabel('Frequency [Hz]')
plt.ylabel('FRF H1')
plt.legend();

运行警告 #

RuntimeWarning: invalid value encountered in divide
  self.frf_conversion = np.power(-1.j / self.get_w_axis(), 2)

预期效果 #

预期FRF图形为带有多个明显共振峰的半对数曲线，是典型的模态FRF特征曲线。

问题排查与修正方案 #

1. 警告原因

警告源于0Hz处的除法运算：当resp_type='a'（加速度响应）时，pyFRF会进行响应类型转换（加速度转位移需要除以角频率ω²），而0Hz对应的ω=0，导致除以0的无效运算。

2. 数据预处理优化

• 信号同步检查：锤击试验中力信号与加速度信号必须严格同步，可先绘制时域力脉冲与加速度响应曲线，确认两者起始点对齐。
• 去除直流分量：直流分量会干扰0Hz附近的FRF结果，对信号做去直流处理：

  p = p - np.mean(p)
  a = a - np.mean(a)
  

3. pyFRF参数调整

• 采样频率避免强制转int：代码中int(fs)会丢失采样频率的精度，直接使用浮点型fs：

  frf_object = FRF(sampling_freq=fs, exc=p, resp=a, resp_type='a', frf_type='H1')
  

• 指定频率范围：通过freq_range参数限定计算到2000Hz，减少不必要的计算点：

  frf_object = FRF(
      sampling_freq=fs, 
      exc=p, 
      resp=a, 
      resp_type='a', 
      frf_type='H1',
      freq_range=(0, 2000)
  )
  

4. 修正后示例代码

from pyFRF import FRF
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 数据导入
skiprow_list = [0,2,3,4,5,6,7,8,9,10]
column_list = np.arange(1,15)

data = pd.read_csv("2024_01_03_WTB_drgania/a12.CSV",
                    delimiter=';',
                    decimal=',',
                    skiprows=skiprow_list,
                    usecols=column_list)

t = data[['C1:  Time[s]']].values.flatten()
a = data[['C1:  [m/s^2]']].values.flatten()
p = data[['C7:  [N]']].values.flatten()

# 去除直流分量
p = p - np.mean(p)
a = a - np.mean(a)

N = len(t)
dt = (t[-1] - t[0])/N
fs = 1/dt

# 创建FRF对象，指定频率范围
frf_object = FRF(
    sampling_freq=fs, 
    exc=p, 
    resp=a, 
    resp_type='a', 
    frf_type='H1',
    freq_range=(0, 2000)
)

H1_pyFRF = frf_object.get_FRF()
freq_pyFRF = frf_object.get_f_axis()

# 绘图跳过0Hz点
plt.semilogy(freq_pyFRF[1:], np.abs(H1_pyFRF)[0,0,1:], color="b", label='pyFRF')
plt.xlim(left=0, right=2000)
plt.xlabel('频率 [Hz]')
plt.ylabel('FRF H1')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Fjotolf Hansen