如何使用Python获取数值黑盒函数的频率响应

随机噪声激励方案可行性 #

可以使用随机噪声作为激励计算频率响应，相比脉冲法不需要对齐黑盒内部更新时序，结果更稳定。直接使用fft(out)/fft(in)的方法存在频谱泄露、低能量频点误差放大的问题，推荐使用功率谱比值法计算，这是线性系统频率响应测量的标准工程方案。

原理是通过计算输入信号的自功率谱Pxx和输入输出的互功率谱Pxy，得到系统频率响应H = Pxy / Pxx，该方法天然抗干扰，结果平滑度远高于直接FFT相除。

实现步骤与代码 #

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

class blackBox:
    tick = 0
    result = 0
    er1 = er2 = er3 = er4 = er5 = er6 = er7 =  0

    def __call__(s,input):
        s.tick+=1        
        if (s.tick==5):        
            s.tick = 0
            s.result = input + 5.0*s.er1 + 9.0*s.er2 + 5.0*s.er3 -5.0*s.er4 -9.0*s.er5 -5.0*s.er6 - s.er7
            s.er7 = s.er6
            s.er6 = s.er5
            s.er5 = s.er4
            s.er4 = s.er3
            s.er3 = s.er2
            s.er2 = s.er1
            s.er1 = input
        return s.result        

# 参数配置
sample_rate = 10000  # 对应100微秒调用一次，可按需修改
input_len = 200000   # 输入序列越长，频率分辨率越高、结果越平滑
nperseg = 1024       # welch法分段长度，可平衡分辨率和平滑度

# 生成白噪声输入，喂入黑盒得到输出
x = np.random.randn(input_len)
p = blackBox()
y = np.array([p(val) for val in x])

# 计算功率谱与频率响应
freq, Pxx = signal.welch(x, fs=sample_rate, nperseg=nperseg)
freq, Pxy = signal.csd(x, y, fs=sample_rate, nperseg=nperseg)
H = Pxy / Pxx

# 绘制幅频、相频响应
plt.figure(figsize=(10,6))
plt.subplot(211)
plt.plot(freq, 20*np.log10(np.abs(H)))
plt.ylabel('幅值 [dB]')
plt.grid(True)

plt.subplot(212)
plt.plot(freq, np.angle(H, deg=True))
plt.ylabel('相位 [度]')
plt.xlabel('频率 [Hz]')
plt.grid(True)

plt.tight_layout()
plt.show()

原有脉冲法的问题修正 #

你原有代码存在频率轴错位问题：使用rfft计算单边频谱，但调用fftfreq生成双边频率轴，替换为rfftfreq即可修正该问题：

# 原有代码修正片段
f = np.fft.rfft(out)
f1 = abs(f)
freq = np.fft.rfftfreq(n, d=1./sample_rate)  # 替换fftfreq为rfftfreq
plt.plot(freq, f1)

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