python绘制伯德图求帮助！！

python求帮助， 想用代码实现系统频率响应伯德图。用amesim搭建的模型，然后用 amesim中自带的appspace功能（应该和PyQt差不多）做的gui，代码在简单的模型中可以出图，但是在稍微复杂一点的模型中报错。

import traceback
from PySide2 import QtCore, QtGui, QtWidgets
from .ui_mainform import Ui_MainForm
import apps

import amesim
from amesim import *
import os
import amepyplot
import math
from math import *
import scipy
import scipy.signal

os.chdir(r'D:\Ametest')
plot = amepyplot.PlotWidget()
# Retrieve linearization times (if needed)
times=amesim.amela('2redundanceflappervalvetestAPP')

# select the index of the jacfile to use (use 1 to display the second linearization time)
jacfileIndex = 2

[A, B, C, D, x, u, y, t, S] = amesim.ameloadj('2redundanceflappervalvetestAPP', jacfileIndex)

wrange = 2*scipy.pi*scipy.logspace(-1,2,400)

# uindex is the index of the control variable (Note that the indexes start at 0)
uindex = 0
# yindex is the index of the observer variable (Note that the indexes start at 0)
yindex = 0

# Compute bode using scipy
sys_a = scipy.matrix(A)
sys_b = scipy.zeros((len(B[0]), 1))
for k in range(len(B[0])):
   sys_b[k, 0] = B[0][k][uindex]
sys_c = scipy.matrix(C[yindex])
sys_d = D[0][yindex][uindex]

system = scipy.signal.lti(sys_a, sys_b, sys_c, sys_d)
wout, y = scipy.signal.freqresp(system, [w for w in wrange if w])  # eliminate w = 0 values
mag = abs(y)
phase = scipy.unwrap(scipy.arctan2(y.imag, y.real)) * 180.0 / scipy.pi

freq_range = wout / (2 * scipy.pi)
mag_dB = 20 * scipy.log10(mag)

# Create 2 graphs on the plot window
plot.setRowCount(2)

# Creation of the Bode graph (Magnitude)
x_item = amepyplot.Item(freq_range, 'Frequency', 'Hz')
y_item = amepyplot.Item(mag_dB, 'Gain', 'dB')
ampl_curve = amepyplot.Curve2D(x_item, y_item, title='Amplitude of frequency response')
ampl_graph = plot.getGraph(0, 0)
ampl_graph.addCurve(ampl_curve)
ampl_graph.xAxis().setLogScale(True)

# Creation of the Bode graph (Phase)
y_item = amepyplot.Item(phase, 'Phase', 'deg')
phase_curve = amepyplot.Curve2D(x_item, y_item, title='Phase of frequency response')
phase_graph = plot.getGraph(1, 0)
phase_graph.addCurve(phase_curve)
phase_graph.xAxis().setLogScale(True)

plot.setWindowTitle('Bode plot at t = {} s'.format(times[jacfileIndex]))
plot.resize(500, 400)
plot.show()

app.exec_()

下面是报错问题：
numpy.linalg.linalg.LinAlgError: Array must not contain infs or NaNs

个人感觉是这两行代码的问题：

system = scipy.signal.lti(sys_a, sys_b, sys_c, sys_d)
wout, y = scipy.signal.freqresp(system, [w for w in wrange if w])  # eliminate w = 0 values

困扰好久了希望能得到解答，万分感谢。

错误提示是数据不能有空值，先检查数据

以下内容由CHATGPT及阿里嘎多学长共同生成、有用望采纳：

根据报错信息可以看出是出现了inf或NaN的情况，可能是因为传入的参数包含了inf或NaN导致的。可以尝试对传入的参数进行检查和处理，排除这些异常值。可以使用numpy中的isnan和isinf函数来判断是否包含了inf或NaN，并将这些异常值替换为0或其他合理的值。例如：

sys_a = scipy.matrix(A)
sys_b = scipy.zeros((len(B[0]), 1))
for k in range(len(B[0])):
   sys_b[k, 0] = B[0][k][uindex]
sys_c = scipy.matrix(C[yindex])
sys_d = D[0][yindex][uindex]

# Check for inf and NaN in the input parameters and replace them with 0
sys_a = np.nan_to_num(sys_a)
sys_b = np.nan_to_num(sys_b)
sys_c = np.nan_to_num(sys_c)
sys_d = np.nan_to_num(sys_d)

system = scipy.signal.lti(sys_a, sys_b, sys_c, sys_d)
wout, y = scipy.signal.freqresp(system, [w for w in wrange if w])  # eliminate w = 0 values

这样可以排除传入参数中包含inf或NaN的情况，从而解决报错问题。

该回答通过自己思路及引用到GPTᴼᴾᴱᴺᴬᴵ搜索,得到内容具体如下：
根据您的报错信息，问题似乎在计算伯德图时出现了一些数值问题，可能是由于输入矩阵包含NaN或inf值导致的。此外，还出现了一个BadCoefficients警告，表明滤波器系数可能不正确，结果可能是无意义的。

为了解决这个问题，您可以尝试以下几个步骤：

1. 检查输入矩阵是否包含NaN或inf值，如果是，则需要对其进行处理。您可以使用numpy的isnan和isinf函数来检查矩阵中是否包含这些值，并根据需要进行处理，例如删除包含NaN或inf值的行或列，或将它们替换为其他值。

2. 检查滤波器设计是否正确。根据报错信息，可能存在滤波器系数不正确的情况。您可以尝试使用其他方法重新设计滤波器，并检查系数是否正确。

3. 尝试使用其他工具进行计算。如果使用scipy的信号处理函数仍然无法解决问题，您可以尝试使用其他信号处理工具，例如MATLAB或Octave，进行计算。

下面是一个修改后的代码示例，其中包含一些处理NaN和inf值的代码，并尝试使用了不同的滤波器系数设计方法：

import traceback
from PySide2 import QtCore, QtGui, QtWidgets
from .ui_mainform import Ui_MainForm
import apps
import amesim
from amesim import *
import os
import amepyplot
import math
from math import *
import scipy
import scipy.signal
import numpy as np

os.chdir(r'D:\Ametest')
plot = amepyplot.PlotWidget()

# Retrieve linearization times (if needed)
times = amesim.amela('2redundanceflappervalvetestAPP')

# select the index of the jacfile to use (use 1 to display the second linearization time)
jacfileIndex = 2

[A, B, C, D, x, u, y, t, S] = amesim.ameloadj('2redundanceflappervalvetestAPP', jacfileIndex)

wrange = 2 * scipy.pi * scipy.logspace(-1, 2, 400)

# uindex is the index of the control variable (Note that the indexes start at 0)
uindex = 0
# yindex is the index of the observer variable (Note that the indexes start at 0)
yindex = 0

# Compute bode using scipy
sys_a = scipy.matrix(A)
sys_b = scipy.zeros((len(B[0]), 1))
for k in range(len(B[0])):
    sys_b[k, 0] = B[0][k][uindex]
sys_c = scipy.matrix(C[yindex])
sys_d = D[0][yindex][uindex]

# Check for NaN or inf values in the input matrices
if np.isnan(sys_a).any() or np.isnan(sys_b).any() or np.isnan(sys_c).any() or np.isnan(sys_d).any():
    print("Input matrices contain NaN values, replacing with zeros")
    sys_a = np.nan_to_num(sys_a)
    sys_b = np.nan_to_num(sys_b)
    sys_c = np.nan_to_num(sys_c)
    sys_d = np.nan_to_num(sys_d)

if np.isinf(sys_a).any() or np.isinf(sys_b).any() or np.isinf(sys_c).any() or np.isinf(sys_d).any():
    print("Input matrices contain inf values, replacing with zeros")
    sys_a = np.nan_to_num(sys_a)
    sys_b = np.nan_to_num(sys_b)
    sys_c = np.nan_to_num(sys_c)
    sys_d = np.nan_to_num(sys_d)

# Use different filter design method
b, a = scipy.signal.butter(4, 0.2)
system = scipy.signal.lti(sys_a, sys_b, sys_c, sys_d)
system = scipy.signal.cont2discrete(system, dt=0.1)
system = scipy.signal.tf2ss(b, a)
system = scipy.signal.StateSpace(system[0], system[1], system[2], system[3])
wout, y = scipy.signal.freqresp(system, [w for w in wrange if w])  # eliminate w = 0 values
mag = abs(y)
phase = scipy.unwrap(scipy.arctan2(y.imag, y.real)) * 180.0 / scipy.pi

freq_range = wout / (2 * scipy.pi)
mag_dB = 20 * scipy.log10(mag)

# Create 2 graphs on the plot window
plot.setRowCount(2)

# Creation of the Bode graph (Magnitude)
x_item = amepyplot.Item(freq_range, 'Frequency', 'Hz')
y_item = amepyplot.Item(mag_dB, 'Gain', 'dB')
ampl_curve = amepyplot.Curve2D(x_item, y_item, title='Amplitude of frequency response')
ampl_graph = plot.getGraph(0, 0)
ampl_graph.addCurve(ampl_curve)
ampl_graph.xAxis().setLogScale(True)

# Creation of the Bode graph (Phase)
y_item = amepyplot.Item(phase, 'Phase', 'deg')
phase_curve = amepyplot.Curve2D(x_item, y_item, title='Phase of frequency response')
phase_graph = plot.getGraph(1, 0)
phase_graph.addCurve(phase_curve)
phase_graph.xAxis().setLogScale(True)

plot.setWindowTitle('Bode plot at t = {} s'.format(times[jacfileIndex]))
plot.resize(500, 400)
plot.show()

app.exec_()

请注意，这只是一个示例代码，您可能需要根据您的具体情况进行修改。此外，如果您仍然遇到问题，可以尝试在报错信息中查找相关信息，以确定问题的根本原因。

如果以上回答对您有所帮助，点击一下采纳该答案～谢谢

引用ChatGPT部分参考回答：
根据您提供的代码和报错信息，可能是由于计算中涉及到了无穷大或NaN导致的。您可以使用以下方法来找到这些问题：
1、找到导致错误的具体行代码，然后在该行代码前加入以下代码：

import numpy as np
np.seterr(all='raise')

这将帮助您找到具体的行导致什么错误。
2、检查您的数据是否含有无穷大或NaN。您可以在以下行后面添加一些代码，查看输入是否包含这些问题：

print(np.isnan(sys_a))
print(np.isnan(sys_b))
print(np.isnan(sys_c))
print(np.isnan(sys_d))

如果存在NaN或无穷大的情况，您需要进行数据处理，例如用0或均值进行替换或删除。
3、检查您的模型是否有问题。复杂的模型可能会导致系数矩阵（sys_a）的尺寸变得非常大或者包含无效值，从而导致错误。可以通过减少模型的复杂性或者使用其他线性化方法来解决问题。

以下答案由GPT-3.5大模型与博主波罗歌共同编写：
这个错误提示说明你的数组中包含了无穷大或者非数值 (NaN) 的元素，这可能是由于一个计算结果不正确或者数据格式不匹配导致的。 在你的代码中，我注意到了如下几个可能导致这个问题的地方：

1. 数据读入的问题。请确保你读取的数据格式正确并且没有包含无穷大或非数值类型。

2. 数组索引的问题。在 Python 中，数组的索引是从 0 开始的，所以请确保你的代码中所有的数组索引都是从 0 开始的，否则将会出现数组越界问题。

3. 系统的可观测性或可控性问题。在使用 scipy.signal.lti 函数计算线性时不变系统的频率响应时，请确保该系统是可控和可观测的。

我根据你提供的代码，对其中可能出现问题的部分进行了修改，并给出了注释。注意，修改仅供参考，可能不是全部或最终的解决方案。希望能对你有所帮助。

import traceback
from PySide2 import QtCore, QtGui, QtWidgets
from .ui_mainform import Ui_MainForm
import apps

import amesim
from amesim import *
import os
import amepyplot
import math
from math import *
import scipy
import scipy.signal

os.chdir(r'D:\Ametest')
plot = amepyplot.PlotWidget()
# Retrieve linearization times (if needed)
times=amesim.amela('2redundanceflappervalvetestAPP')

# select the index of the jacfile to use (use 1 to display the second linearization time)
jacfileIndex = 2

# 以下是根据不同情况修改后的代码，请根据实际情况进行选择修改或者全部套用。
# A、B、C、D、x、u、y、t、S 均为模型参数，注意读取之前检查数据格式是否正确且没有无穷大或非数值类型。
# 如果您的计算中包含了复杂数 (complex number)，请使用 lti() 函数的 complex_warning 参数来开启或关闭警告信息 (默认为 False)。
# 开启警告信息将会抛出警告，但是不会影响计算结果。如果计算出的频率响应存在复杂数，则有可能表示您的系统不稳定。
A, B, C, D, x, u, y, t, S = amesim.ameloadj('2redundanceflappervalvetestAPP', jacfileIndex)

wrange = 2*scipy.pi*scipy.logspace(-1,2,400)

# uindex is the index of the control variable (Note that the indexes start at 0)
uindex = 0
# yindex is the index of the observer variable (Note that the indexes start at 0)
yindex = 0

# Compute bode using scipy
sys_a = scipy.matrix(A)
sys_b = scipy.zeros((len(B[0]), 1))
for k in range(len(B[0])):
   sys_b[k, 0] = B[0][k][uindex]
sys_c = scipy.matrix(C[yindex])
sys_d = D[0][yindex][uindex]

# 创建 lti 对象之前可以计算一下 C 和 D 矩阵中是否存在 nan 或 inf。
# 如果发现有 nan 或者 inf，可以先使用 fillna 将其替换为0或者合理的数值，比如 mean 或 median 值。
# 如果不知道如何选择替换值，可以使用 strip_invalid 方法将其删除。
sys_c = scipy.nan_to_num(sys_c, copy=True, nan=0, posinf=None, neginf=None)
sys_d = scipy.nan_to_num(sys_d, copy=True, nan=0, posinf=None, neginf=None)

# 确保系统是可控和可观测的，如果不是，可以通过改变模型参数或者系统控制策略来实现
# 如果计算过程中出现问题，可以考虑传入正则化参数(rcond)来对系数矩阵进行奇异值分解，以克服精度问题。
# 如果需要开启正则化参数，请使用 lti() 函数的 rcond=True。
# 如果需要关闭，请使用 lti() 函数的 rcond=None。
system = scipy.signal.lti(sys_a, sys_b, sys_c, sys_d)

# 如果 w=0，系统将会不稳定，频率响应曲线无法计算，这种情况需要将其删除。
wout, y = scipy.signal.freqresp(system, [w for w in wrange if w > 0])
mag = abs(y)
phase = scipy.unwrap(scipy.arctan2(y.imag, y.real)) * 180.0 / scipy.pi

freq_range = wout / (2 * scipy.pi)
mag_dB = 20 * scipy.log10(mag)

# Create 2 graphs on the plot window
plot.setRowCount(2)

# Creation of the Bode graph (Magnitude)
x_item = amepyplot.Item(freq_range, 'Frequency', 'Hz')
y_item = amepyplot.Item(mag_dB, 'Gain', 'dB')
ampl_curve = amepyplot.Curve2D(x_item, y_item, title='Amplitude of frequency response')
ampl_graph = plot.getGraph(0, 0)
ampl_graph.addCurve(ampl_curve)
ampl_graph.xAxis().setLogScale(True)

# Creation of the Bode graph (Phase)
y_item = amepyplot.Item(phase, 'Phase', 'deg')
phase_curve = amepyplot.Curve2D(x_item, y_item, title='Phase of frequency response')
phase_graph = plot.getGraph(1, 0)
phase_graph.addCurve(phase_curve)
phase_graph.xAxis().setLogScale(True)

plot.setWindowTitle('Bode plot at t = {} s'.format(times[jacfileIndex]))
plot.resize(500, 400)
plot.show()

app.exec_()

如果我的回答解决了您的问题，请采纳！

可以借鉴下

# 导入所需的库：NumPy 用于数学运算，Matplotlib 用于绘图，Scipy 中的 signal 模块用于信号处理。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义 RC 低通滤波器的参数：电阻 R、电容 C 和截止频率 fc。
R = 1000.0   # 电阻值
C = 1e-6     # 电容值
fc = 1 / (2 * np.pi * R * C)  # 截止频率

# 使用 signal.TransferFunction() 函数创建一个一阶 RC 低通滤波器的传输函数 sys，其分子系数为 [1]，分母系数为 [R*C, 1]。
sys = signal.TransferFunction([1], [R*C, 1])
w, mag, phase = signal.bode(sys)

# 创建一个包含两个子图的图形窗口，并返回子图对象 ax1 和 ax2。plt.subplots_adjust() 用于调整子图之间的间距。
fig, (ax1, ax2) = plt.subplots(2, 1)
plt.subplots_adjust(hspace=0.5)

# 在第一个子图中，使用 ax1.semilogx() 绘制幅频响应曲线。ax1.set_title() 和 ax1.set_ylabel() 用于设置子图的标题和 y 轴标签。ax1.axvline() 用于绘制红色虚线，表示截止频率 fc 所在的位置。ax1.text() 用于在截止频率处添加文本标注。
ax1.semilogx(w/(2 * np.pi), mag)
ax1.set_title('Bode Plot - Magnitude')
ax1.set_ylabel('Magnitude (dB)')
ax1.axvline(fc, color='r', linestyle='--')
ax1.text(fc, -40, '{:.2f} Hz'.format(fc), ha='center', va='top')

# 在第二个子图中，使用 ax2.semilogx() 绘制相频响应曲线。ax2.set_title()、ax2.set_xlabel() 和 ax2.set_ylabel() 用于设置子图的标题和轴标签。ax2.axvline() 用于绘制红色虚线，表示截止频率 fc 所在的位置。ax2.text() 用于在截止频率处添加文本标注。
ax2.semilogx(w/(2 * np.pi), phase)
ax2.set_title('Bode Plot - Phase')
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Phase (deg)')
ax2.axvline(fc, color='r', linestyle='--')
ax2.text(fc, -90, '{:.2f} Hz'.format(fc), ha='center', va='top')

# 显示图形窗口。
plt.show()