简介一套基于MATLAB的扫频法开环传递函数求解程序面向控制工程专业学生、科研人员及系统调试工程师用于通过频率响应实验确定线性时不变系统的开环传递函数模型。压缩包内仅含1个m脚本文件包体大小约2KB代码结构紧凑但功能完整涵盖扫频输入信号生成、系统激励、响应采集、幅频与相频数据处理、传递函数估计等步骤并预留了使用bode、freqs等函数绘制Bode图或奈奎斯特图的接口便于直观查看系统频率特性。扫频法原理清晰易懂代码注释完整支持针对不同被控对象修改扫频范围和采样参数后直接复用。目前该资源已有5697人浏览学习程序经过实际运行验证。下载后即可获得可运行的MATLAB源码既能深入理解扫频法求取开环传递函数的完整思路也能学习利用实验数据进行系统辨识的实操方法对控制系统建模、动态特性分析以及后续控制器设计均有较强的实用价值。1. 扫频法求开环传递函数的现场价值调试伺服系统或电源环路时手头没有网络分析仪是常态。你有一个被控对象、一块数据采集卡和一套 MATLAB却要在今天下班前给出开环传递函数的 Bode 图判断穿越频率和相位裕度是否符合指标。扫频法就是干这个的给系统输入端施加频率连续变化的正弦激励同步采集输入和输出信号在频域里逐点算出增益和相位最后拼成完整的频率响应曲线再用 MATLAB 的辨识工具把它变成G(s) num/den形式的传递函数。整个过程只需要一个激励通道、一个反馈采集通道不给系统附加任何硬件误差主要取决于激励设计和后期提取算法。这套方法适合三类人做运动控制或电源环路补偿的工程师需要确认理论建模与实际对象偏差的算法工程师以及刚接触系统辨识、想从实验数据里得到传递函数模型的学生。标题里的“MATLAB程序”不是指某个现成脚本而是指从激励生成、数据采集到频响提取、模型拟合的一条完整链路。接下来按这个链路走一遍重点说清楚哪些参数决定成败哪些坑在相位提取和直流偏置上等着你。2. 扫频法求传递函数的基本原理与激励信号制作2.1 为什么用扫频而不是直接给阶跃开环传递函数的一个重要表达形式是频率响应把正弦信号x(t)A·sin(2πft)送入线性时不变系统稳态输出仍然是同频正弦只是幅度变成A·|G(j2πf)|相位移动了∠G(j2πf)。遍历频率 f就能得到G(jω)在整个频带上的幅频和相频特性。理论上给一个阶跃信号也能通过拉普拉斯变换求传递函数但阶跃的能量分布在高频段迅速衰减信噪比差而且对积分环节和振荡环节的辨识精度很低。扫频信号的优点是把能量均匀或按对数规律分配到每个关心的频点上每个频点的激励时间又足够长让系统达到稳态从而用窄带提取的方式获得干净的幅值和相位。这也是为什么实际工程里扫频法求开环传递函数成为首选。2.2 对数扫频与步进正弦两种生成方案MATLAB 程序里生成激励信号有两种常见做法对应不同的求传递函数策略。方案 A 是线性或对数 chirp 信号用chirp函数一次性生成整个时域序列激励时间短适合在线快速测试。缺点是每个频点的驻留时间不均匀低频段可能还没到稳态就扫过去了相位提取误差偏大。方案 B 是步进正弦stepped sine把关心的频率范围按对数间隔分成若干频点每个频点单独生成正弦波持续若干个周期后再切到下一个频点。它虽然耗时但每个频点都能保证稳态抗噪声能力强更适合需要精确求开环传递函数的场景。实际做系统辨识时我一般用方案 B因为它和后面的峰值提取算法配合最好。下面是一个生成多频点步进正弦激励的 MATLAB 程序骨架fs 100e3; % 采样率 100 kHz f_start 20; % 起始频率 20 Hz f_end 20e3; % 终止频率 20 kHz pts_per_dec 20; % 每十倍频程点数 cycles 10; % 每个频点激励周期数 freqs logspace(log10(f_start), log10(f_end), ... round(pts_per_dec * log10(f_end/f_start))); sig []; t_total []; for k 1:length(freqs) fk freqs(k); t (0:ceil(cycles*fs/fk)-1) / fs; sig [sig, 0.5 * sin(2*pi*fk*t)]; t_total [t_total, t sum(1:0)]; % 仅示意记录全局时间 end这段代码里fs决定了每个频点的波形分辨率必须大于最高频率的 10 倍以上否则谐波混叠会污染提取结果。cycles取 10 到 20 比较稳妥太少达不到稳态太多浪费时间。pts_per_dec决定频响曲线的点数20 点/十倍频程已经能覆盖常见的二阶振荡峰。2.3 激励信号输出前的环节设计生成的正弦序列不能直接送进被控对象。先检查幅值是否在系统线性区内——开环增益高的系统幅值过大容易让输出饱和导致实测增益偏低幅值过小则信噪比不足。通常先给定一个标称幅值的 20%观察输出波形无明显畸变再逐步加大。还要考虑功率放大器或驱动器的输出阻抗如果激励通路存在二次低通其影响必须从实测结果里扣掉否则求出来的传递函数会额外包含前向通路特性。把激励序列写入数据采集卡的模拟输出前建议先做一次零填充。头部加一段静音用于同步触发采集尾部加静音保证最后一个频点的稳态响应被完整记录。这个细节在后面的相位提取阶段会省很多事因为 MATLAB 程序需要对输入信号和输出信号做精确对齐两端的延时代价会直接反映在相位曲线上的线性倾斜里。3. 从扫频数据中提取频率响应的 MATLAB 实现3.1 用 FFT 还是用 Goertzel 提取各频点幅值相位拿到输入和输出两路时域信号后求开环传递函数的关键是从每个频点附近提取出该频率分量的幅度和相位。FFT 是整个频段一次性算完频率分辨率受采样点数和窗函数限制扫频各频点频率通常不在 FFT 的频率网格上会引入泄漏误差需要加窗和插值修正。Goertzel 算法更适合单频提取场景。它每次只计算某一个指定频率的 DFT 系数算法简单占用资源低而且可以直接给出该频点的实部和虚部换算幅度和相位非常方便。MATLAB 从 R2018b 起提供了goertzel函数输入一段时域序列和目标频率索引返回对应频点的 DFT 值。实际使用场景里采集到的输入和输出信号可以用同一个goertzel调用分别处理保证两者的参考相位基准一致。3.2 逐频点计算增益和相位的完整程序下面这段 MATLAB 程序展示从原始采集数据u_in激励和y_out响应中逐频点计算幅频和相频的完整流程freqs logspace(log10(20), log10(20000), 60); n_freq length(freqs); gain_db zeros(1, n_freq); phase_deg zeros(1, n_freq); N length(u_in); for k 1:n_freq fk freqs(k); % 只取该频点附近且包含整数周期的数据段 seg_len round(10 * fs / fk); start_idx round(2 * fs / fk) 1; if start_idx seg_len - 1 N seg_len N - start_idx 1; end u_seg u_in(start_idx : start_idx seg_len - 1); y_seg y_out(start_idx : start_idx seg_len - 1); % Goertzel 提取指定频率的复数分量 dft_u goertzel(u_seg, fk, fs); dft_y goertzel(y_seg, fk, fs); amp_u abs(dft_u) * 2 / seg_len; amp_y abs(dft_y) * 2 / seg_len; phase_u angle(dft_u); phase_y angle(dft_y); gain_db(k) 20 * log10(amp_y / max(amp_u, 1e-12)); phase_deg(k) (phase_y - phase_u) * 180 / pi; end这里goertzel的调用方式是自定义写法实际使用时需要根据 MATLAB 版本确认参数格式。核心思想是取 10 个整周期的数据段舍弃前 2 个周期的瞬态响应再对输入和输出在同一频点上做窄带提取用输出复数除以输入复数得到该频点的频率响应。取整周期段的做法能显著降低频谱泄漏比直接整段 FFT 加窗的精度高一个量级。3.3 相位连续化处理angle函数返回的相位在 -π 到 π 之间跨越 ±180° 时会出现跳变直接绘制 Bode 图会看到锯齿状折线。处理方法是对相频曲线做 unwrapphase_deg_cont unwrap(phase_deg * pi / 180) * 180 / pi;unwrap会自动检测相邻点之间超过 180° 的跳变并补上 2π 的整数倍。这里要注意如果扫频点数太少两个相邻频点间的真实相位变化超过 180°unwrap 会补错方向。所以每十倍频程至少要有 10 个以上的测点带宽较宽或系统阶次较高时建议加到 20 点以上。另一个常见问题是起始相位参考不一致造成整体偏置可以在低频段先和理论模型对比一次把固定偏差去掉。3.4 数据组织与保存格式扫频法得到的结果数据建议保存为三列文本频率、增益、相位为后续拟合做准备。打开方式不限但要注意编码问题。常见做法是用writematrix直接导出或者用save存成 mat 文件。如果现场只有 CSV 格式的示波器数据也可以用 MATLAB 的readmatrix读进来再按列切分出输入输出信号。数据组织上最重要的一个原则是确保采样率信息随数据一起保存因为求传递函数时所有频点计算都依赖 fs丢了 fs 后面只能靠猜。这里有一个工程细节值得注意采集到的信号如果带有直流偏置Goertzel 提取结果会包含直流泄漏尤其是在低频段幅度误差会被放大。处理办法是先把均值减掉u_in u_in - mean(u_in)y_out y_out - mean(y_out)。对开环传递函数辨识来说直流增益通常无穷大减去偏置不会丢信息反而能消除 ADC 失调带来的相位污染。4. 由 Bode 数据拟合开环传递函数模型4.1 拟合前的频域加权策略拿到频响数据后求传递函数的下一步是把它拟合成有理传递函数G(s) b_m s^m ... / (s^n a_{n-1} s^{n-1} ...)。MATLAB 的系统辨识工具箱提供了tfest函数直接输入频率响应数据就能估计零极点。但直接调用往往会得到差强人意的结果因为 Bode 图低频段增益很高高频段接近零普通最小二乘会被低频大数值支配导致高频段拟合失真。拟合前要对频域数据做加权处理。常见做法是用对数频率间隔重新插值使每个十倍频程内的采样点数一致或者给tfest传入频率响应的标准差向量在低频段给更小的权重。另一种做法是先把增益转成线性幅度再做加权最小二乘避免对数变换带来的噪声放大。4.2 用 tfest 估计传递函数的程序骨架下面这段代码演示从实测频率响应数据出发用tfest估计开环传递函数的过程% freq_hz: 频点向量单位 Hz % resp: 复数频率响应resp(k) gain(k) * exp(1i * phase_rad(k)) data idfrd(resp, freq_hz, 0, FrequencyUnit, Hz); np 3; % 极点个数 nz 1; % 零点个数 options tfestOptions(WeightingFilter, inv, EnforceStability, true); sys_est tfest(data, np, nz, options); % 与实测数据对比 [mag, phase, wout] bode(sys_est, freq_hz * 2 * pi);tfest的第一个参数是频域数据对象idfrd函数把实测频率响应包装成辨识工具能识别的格式。EnforceStability强制极点落在左半平面这对开环传递函数很有用——如果被控对象本身是积分环节加惯性环节拟合时很容易跑出右半平面极点来补偿相位得到物理上不存在的模型。WeightingFilter设为inv时按幅值倒加权相当于在对数坐标上进行拟合更贴合 Bode 图的视觉直觉。4.3 开环传递函数增益按首一还是尾一写热词里提到的“开环传递函数增益是首一还是尾一”是一个值得说清楚的问题。传递函数的增益写法不是纯粹的数学格式问题它直接影响辨识参数的物理意义和数值稳定性。首一形式是把分子多项式写成最高次项系数为 1G(s) (s z1) / (s^2 a1·s a2)尾一形式是把常数项归为 1G(s) K · (1 s/z1) / (1 s/p1)(1 s/p2)对开环传递函数而言尾一形式更适合描述实际系统的增益预算。理由很直接开环 Bode 图的低频段增益由K/(s^n)决定K是环路增益设计的核心参数写成尾一后 K 直接代表了直流增益或积分增益工程上可以直接和理论计算的环路增益对照。首一形式下的分子最高次系数为 1低频增益隐含在其他参数里无法一眼看出开环增益的大小。MATLAB 的tfest输出默认是首一形式对含积分环节的系统首一形式会让低频增益表征不直观数值也容易偏大。如果你更关心开环增益是否达标辨识后可以用zpk转换成零极点增益形式再手工整理成尾一形式。下面这行代码把估计结果转换并显示增益 Ksys_zpk zpk(sys_est); % 查看零极点与增益 K用zpk形式观察时如果 K 的数值和理论计算的环路增益差一个量级以上先检查是不是辨识时频率范围没覆盖到穿越频率附近。开环增益在穿越频率处的信息量最大扫频范围应该至少覆盖穿越频率典型 0.1 倍到 10 倍否则拟合出的 K 不可信。4.4 模型阶次怎么定tfest的np和nz参数需要事先估计。给一个参考做法看实测相位曲线在扫频范围内最终下降了多少。一个极点贡献 -90°一个零点贡献 90°。如果相位从低频到高频下降了大约 270°那系统的极点个数大约是 3 到 4 个积分环节也算一个极点。零点个数看相位是否在中频段有回升回升意味着存在零点。阶次宁可少一个也不要多。多一个零极点对会让拟合程序用零极点对消来补偿噪声得到的模型阶次虚高后续做补偿器设计时反而碍事。拟合完成后对比模型的 Bode 曲线和实测数据点如果偏差在 2 dB 和 10° 以内就可以接受。5. 扫频法求开环传递函数的几个高频坑与验证技巧5.1 采样率、频率分配与相位校准fs的选择直接决定扫频上限。奈奎斯特限制下采样率至少是最高激励频率的 10 到 20 倍主要给谐波留空间。被控对象如果存在较强的非线性输出中会有 2 次、3 次谐波这些谐波如果不被采样率滤除会混叠到基频附近的提取结果里增益和相位都会失真。频率点分配要用对数间隔低频段频点间距远小于高频段。这是由传递函数本身的特性决定的系统的转折频率在对数坐标上是均匀分布的用logspace生成频率点就能保证每个转折频率附近都有足够的测点。另一个相关坑是最低频率和最高频率的选取最低频率应该比系统最慢的转折频率低 5 倍以上否则低频段斜率和增益拟合不出来最高频率比穿越频率高 5 到 10 倍即可过高只会增加测试时间。相位提取时经常遇到相位曲线整体偏移的情况这通常是信号采集通道之间的延时差造成的。每个通道的模拟滤波器和调理电路会引入不同的群延时表现为相频曲线叠加了一个线性项-2πf·τ。校准方法是把采集卡的输入和输出短接直接测一条直通频响得到的相位偏差曲线从实测相位里扣除。5.2 最小时间成本验证法每次正式测试前花 30 秒做一个快速验证只扫 10 个频点覆盖目标穿越频率的 0.3 倍到 3 倍把实测相位曲线和模型对比。如果这 10 个点相位趋势正确、穿越频率附近的增益斜率合理再开始全频段扫频。这个做法可以在采集系统出错、接线反相、幅值饱和等问题上快速止损比跑完整个扫频再发现数据不可用高效得多。5.3 加一个在线闭环交叉验证开环传递函数的最终验证方式是闭环阶跃响应。把辨识出的开环模型补上控制器构成单位负反馈仿真闭环阶跃和实际系统的阶跃响应对比超调量和上升时间。幅度偏差在 10% 以内说明扫频法求开环传递函数的数据链路基本可靠。如果阶跃响应实测超调明显大于模型预测优先怀疑相位提取偏乐观实际相位裕度小于模型值。这时把扫频激励的幅值降低三分之一重新测一次如果相位曲线比原来更差说明之前的数据已经触碰到系统非线性区。线性度检查就这么简单减小激励幅值模型不变才是合格数据。本文还有配套的精品资源点击获取 MATLAB扫频法求开环传递函数全流程解析
简介一套基于MATLAB的扫频法开环传递函数求解程序面向控制工程专业学生、科研人员及系统调试工程师用于通过频率响应实验确定线性时不变系统的开环传递函数模型。压缩包内仅含1个m脚本文件包体大小约2KB代码结构紧凑但功能完整涵盖扫频输入信号生成、系统激励、响应采集、幅频与相频数据处理、传递函数估计等步骤并预留了使用bode、freqs等函数绘制Bode图或奈奎斯特图的接口便于直观查看系统频率特性。扫频法原理清晰易懂代码注释完整支持针对不同被控对象修改扫频范围和采样参数后直接复用。目前该资源已有5697人浏览学习程序经过实际运行验证。下载后即可获得可运行的MATLAB源码既能深入理解扫频法求取开环传递函数的完整思路也能学习利用实验数据进行系统辨识的实操方法对控制系统建模、动态特性分析以及后续控制器设计均有较强的实用价值。1. 扫频法求开环传递函数的现场价值调试伺服系统或电源环路时手头没有网络分析仪是常态。你有一个被控对象、一块数据采集卡和一套 MATLAB却要在今天下班前给出开环传递函数的 Bode 图判断穿越频率和相位裕度是否符合指标。扫频法就是干这个的给系统输入端施加频率连续变化的正弦激励同步采集输入和输出信号在频域里逐点算出增益和相位最后拼成完整的频率响应曲线再用 MATLAB 的辨识工具把它变成G(s) num/den形式的传递函数。整个过程只需要一个激励通道、一个反馈采集通道不给系统附加任何硬件误差主要取决于激励设计和后期提取算法。这套方法适合三类人做运动控制或电源环路补偿的工程师需要确认理论建模与实际对象偏差的算法工程师以及刚接触系统辨识、想从实验数据里得到传递函数模型的学生。标题里的“MATLAB程序”不是指某个现成脚本而是指从激励生成、数据采集到频响提取、模型拟合的一条完整链路。接下来按这个链路走一遍重点说清楚哪些参数决定成败哪些坑在相位提取和直流偏置上等着你。2. 扫频法求传递函数的基本原理与激励信号制作2.1 为什么用扫频而不是直接给阶跃开环传递函数的一个重要表达形式是频率响应把正弦信号x(t)A·sin(2πft)送入线性时不变系统稳态输出仍然是同频正弦只是幅度变成A·|G(j2πf)|相位移动了∠G(j2πf)。遍历频率 f就能得到G(jω)在整个频带上的幅频和相频特性。理论上给一个阶跃信号也能通过拉普拉斯变换求传递函数但阶跃的能量分布在高频段迅速衰减信噪比差而且对积分环节和振荡环节的辨识精度很低。扫频信号的优点是把能量均匀或按对数规律分配到每个关心的频点上每个频点的激励时间又足够长让系统达到稳态从而用窄带提取的方式获得干净的幅值和相位。这也是为什么实际工程里扫频法求开环传递函数成为首选。2.2 对数扫频与步进正弦两种生成方案MATLAB 程序里生成激励信号有两种常见做法对应不同的求传递函数策略。方案 A 是线性或对数 chirp 信号用chirp函数一次性生成整个时域序列激励时间短适合在线快速测试。缺点是每个频点的驻留时间不均匀低频段可能还没到稳态就扫过去了相位提取误差偏大。方案 B 是步进正弦stepped sine把关心的频率范围按对数间隔分成若干频点每个频点单独生成正弦波持续若干个周期后再切到下一个频点。它虽然耗时但每个频点都能保证稳态抗噪声能力强更适合需要精确求开环传递函数的场景。实际做系统辨识时我一般用方案 B因为它和后面的峰值提取算法配合最好。下面是一个生成多频点步进正弦激励的 MATLAB 程序骨架fs 100e3; % 采样率 100 kHz f_start 20; % 起始频率 20 Hz f_end 20e3; % 终止频率 20 kHz pts_per_dec 20; % 每十倍频程点数 cycles 10; % 每个频点激励周期数 freqs logspace(log10(f_start), log10(f_end), ... round(pts_per_dec * log10(f_end/f_start))); sig []; t_total []; for k 1:length(freqs) fk freqs(k); t (0:ceil(cycles*fs/fk)-1) / fs; sig [sig, 0.5 * sin(2*pi*fk*t)]; t_total [t_total, t sum(1:0)]; % 仅示意记录全局时间 end这段代码里fs决定了每个频点的波形分辨率必须大于最高频率的 10 倍以上否则谐波混叠会污染提取结果。cycles取 10 到 20 比较稳妥太少达不到稳态太多浪费时间。pts_per_dec决定频响曲线的点数20 点/十倍频程已经能覆盖常见的二阶振荡峰。2.3 激励信号输出前的环节设计生成的正弦序列不能直接送进被控对象。先检查幅值是否在系统线性区内——开环增益高的系统幅值过大容易让输出饱和导致实测增益偏低幅值过小则信噪比不足。通常先给定一个标称幅值的 20%观察输出波形无明显畸变再逐步加大。还要考虑功率放大器或驱动器的输出阻抗如果激励通路存在二次低通其影响必须从实测结果里扣掉否则求出来的传递函数会额外包含前向通路特性。把激励序列写入数据采集卡的模拟输出前建议先做一次零填充。头部加一段静音用于同步触发采集尾部加静音保证最后一个频点的稳态响应被完整记录。这个细节在后面的相位提取阶段会省很多事因为 MATLAB 程序需要对输入信号和输出信号做精确对齐两端的延时代价会直接反映在相位曲线上的线性倾斜里。3. 从扫频数据中提取频率响应的 MATLAB 实现3.1 用 FFT 还是用 Goertzel 提取各频点幅值相位拿到输入和输出两路时域信号后求开环传递函数的关键是从每个频点附近提取出该频率分量的幅度和相位。FFT 是整个频段一次性算完频率分辨率受采样点数和窗函数限制扫频各频点频率通常不在 FFT 的频率网格上会引入泄漏误差需要加窗和插值修正。Goertzel 算法更适合单频提取场景。它每次只计算某一个指定频率的 DFT 系数算法简单占用资源低而且可以直接给出该频点的实部和虚部换算幅度和相位非常方便。MATLAB 从 R2018b 起提供了goertzel函数输入一段时域序列和目标频率索引返回对应频点的 DFT 值。实际使用场景里采集到的输入和输出信号可以用同一个goertzel调用分别处理保证两者的参考相位基准一致。3.2 逐频点计算增益和相位的完整程序下面这段 MATLAB 程序展示从原始采集数据u_in激励和y_out响应中逐频点计算幅频和相频的完整流程freqs logspace(log10(20), log10(20000), 60); n_freq length(freqs); gain_db zeros(1, n_freq); phase_deg zeros(1, n_freq); N length(u_in); for k 1:n_freq fk freqs(k); % 只取该频点附近且包含整数周期的数据段 seg_len round(10 * fs / fk); start_idx round(2 * fs / fk) 1; if start_idx seg_len - 1 N seg_len N - start_idx 1; end u_seg u_in(start_idx : start_idx seg_len - 1); y_seg y_out(start_idx : start_idx seg_len - 1); % Goertzel 提取指定频率的复数分量 dft_u goertzel(u_seg, fk, fs); dft_y goertzel(y_seg, fk, fs); amp_u abs(dft_u) * 2 / seg_len; amp_y abs(dft_y) * 2 / seg_len; phase_u angle(dft_u); phase_y angle(dft_y); gain_db(k) 20 * log10(amp_y / max(amp_u, 1e-12)); phase_deg(k) (phase_y - phase_u) * 180 / pi; end这里goertzel的调用方式是自定义写法实际使用时需要根据 MATLAB 版本确认参数格式。核心思想是取 10 个整周期的数据段舍弃前 2 个周期的瞬态响应再对输入和输出在同一频点上做窄带提取用输出复数除以输入复数得到该频点的频率响应。取整周期段的做法能显著降低频谱泄漏比直接整段 FFT 加窗的精度高一个量级。3.3 相位连续化处理angle函数返回的相位在 -π 到 π 之间跨越 ±180° 时会出现跳变直接绘制 Bode 图会看到锯齿状折线。处理方法是对相频曲线做 unwrapphase_deg_cont unwrap(phase_deg * pi / 180) * 180 / pi;unwrap会自动检测相邻点之间超过 180° 的跳变并补上 2π 的整数倍。这里要注意如果扫频点数太少两个相邻频点间的真实相位变化超过 180°unwrap 会补错方向。所以每十倍频程至少要有 10 个以上的测点带宽较宽或系统阶次较高时建议加到 20 点以上。另一个常见问题是起始相位参考不一致造成整体偏置可以在低频段先和理论模型对比一次把固定偏差去掉。3.4 数据组织与保存格式扫频法得到的结果数据建议保存为三列文本频率、增益、相位为后续拟合做准备。打开方式不限但要注意编码问题。常见做法是用writematrix直接导出或者用save存成 mat 文件。如果现场只有 CSV 格式的示波器数据也可以用 MATLAB 的readmatrix读进来再按列切分出输入输出信号。数据组织上最重要的一个原则是确保采样率信息随数据一起保存因为求传递函数时所有频点计算都依赖 fs丢了 fs 后面只能靠猜。这里有一个工程细节值得注意采集到的信号如果带有直流偏置Goertzel 提取结果会包含直流泄漏尤其是在低频段幅度误差会被放大。处理办法是先把均值减掉u_in u_in - mean(u_in)y_out y_out - mean(y_out)。对开环传递函数辨识来说直流增益通常无穷大减去偏置不会丢信息反而能消除 ADC 失调带来的相位污染。4. 由 Bode 数据拟合开环传递函数模型4.1 拟合前的频域加权策略拿到频响数据后求传递函数的下一步是把它拟合成有理传递函数G(s) b_m s^m ... / (s^n a_{n-1} s^{n-1} ...)。MATLAB 的系统辨识工具箱提供了tfest函数直接输入频率响应数据就能估计零极点。但直接调用往往会得到差强人意的结果因为 Bode 图低频段增益很高高频段接近零普通最小二乘会被低频大数值支配导致高频段拟合失真。拟合前要对频域数据做加权处理。常见做法是用对数频率间隔重新插值使每个十倍频程内的采样点数一致或者给tfest传入频率响应的标准差向量在低频段给更小的权重。另一种做法是先把增益转成线性幅度再做加权最小二乘避免对数变换带来的噪声放大。4.2 用 tfest 估计传递函数的程序骨架下面这段代码演示从实测频率响应数据出发用tfest估计开环传递函数的过程% freq_hz: 频点向量单位 Hz % resp: 复数频率响应resp(k) gain(k) * exp(1i * phase_rad(k)) data idfrd(resp, freq_hz, 0, FrequencyUnit, Hz); np 3; % 极点个数 nz 1; % 零点个数 options tfestOptions(WeightingFilter, inv, EnforceStability, true); sys_est tfest(data, np, nz, options); % 与实测数据对比 [mag, phase, wout] bode(sys_est, freq_hz * 2 * pi);tfest的第一个参数是频域数据对象idfrd函数把实测频率响应包装成辨识工具能识别的格式。EnforceStability强制极点落在左半平面这对开环传递函数很有用——如果被控对象本身是积分环节加惯性环节拟合时很容易跑出右半平面极点来补偿相位得到物理上不存在的模型。WeightingFilter设为inv时按幅值倒加权相当于在对数坐标上进行拟合更贴合 Bode 图的视觉直觉。4.3 开环传递函数增益按首一还是尾一写热词里提到的“开环传递函数增益是首一还是尾一”是一个值得说清楚的问题。传递函数的增益写法不是纯粹的数学格式问题它直接影响辨识参数的物理意义和数值稳定性。首一形式是把分子多项式写成最高次项系数为 1G(s) (s z1) / (s^2 a1·s a2)尾一形式是把常数项归为 1G(s) K · (1 s/z1) / (1 s/p1)(1 s/p2)对开环传递函数而言尾一形式更适合描述实际系统的增益预算。理由很直接开环 Bode 图的低频段增益由K/(s^n)决定K是环路增益设计的核心参数写成尾一后 K 直接代表了直流增益或积分增益工程上可以直接和理论计算的环路增益对照。首一形式下的分子最高次系数为 1低频增益隐含在其他参数里无法一眼看出开环增益的大小。MATLAB 的tfest输出默认是首一形式对含积分环节的系统首一形式会让低频增益表征不直观数值也容易偏大。如果你更关心开环增益是否达标辨识后可以用zpk转换成零极点增益形式再手工整理成尾一形式。下面这行代码把估计结果转换并显示增益 Ksys_zpk zpk(sys_est); % 查看零极点与增益 K用zpk形式观察时如果 K 的数值和理论计算的环路增益差一个量级以上先检查是不是辨识时频率范围没覆盖到穿越频率附近。开环增益在穿越频率处的信息量最大扫频范围应该至少覆盖穿越频率典型 0.1 倍到 10 倍否则拟合出的 K 不可信。4.4 模型阶次怎么定tfest的np和nz参数需要事先估计。给一个参考做法看实测相位曲线在扫频范围内最终下降了多少。一个极点贡献 -90°一个零点贡献 90°。如果相位从低频到高频下降了大约 270°那系统的极点个数大约是 3 到 4 个积分环节也算一个极点。零点个数看相位是否在中频段有回升回升意味着存在零点。阶次宁可少一个也不要多。多一个零极点对会让拟合程序用零极点对消来补偿噪声得到的模型阶次虚高后续做补偿器设计时反而碍事。拟合完成后对比模型的 Bode 曲线和实测数据点如果偏差在 2 dB 和 10° 以内就可以接受。5. 扫频法求开环传递函数的几个高频坑与验证技巧5.1 采样率、频率分配与相位校准fs的选择直接决定扫频上限。奈奎斯特限制下采样率至少是最高激励频率的 10 到 20 倍主要给谐波留空间。被控对象如果存在较强的非线性输出中会有 2 次、3 次谐波这些谐波如果不被采样率滤除会混叠到基频附近的提取结果里增益和相位都会失真。频率点分配要用对数间隔低频段频点间距远小于高频段。这是由传递函数本身的特性决定的系统的转折频率在对数坐标上是均匀分布的用logspace生成频率点就能保证每个转折频率附近都有足够的测点。另一个相关坑是最低频率和最高频率的选取最低频率应该比系统最慢的转折频率低 5 倍以上否则低频段斜率和增益拟合不出来最高频率比穿越频率高 5 到 10 倍即可过高只会增加测试时间。相位提取时经常遇到相位曲线整体偏移的情况这通常是信号采集通道之间的延时差造成的。每个通道的模拟滤波器和调理电路会引入不同的群延时表现为相频曲线叠加了一个线性项-2πf·τ。校准方法是把采集卡的输入和输出短接直接测一条直通频响得到的相位偏差曲线从实测相位里扣除。5.2 最小时间成本验证法每次正式测试前花 30 秒做一个快速验证只扫 10 个频点覆盖目标穿越频率的 0.3 倍到 3 倍把实测相位曲线和模型对比。如果这 10 个点相位趋势正确、穿越频率附近的增益斜率合理再开始全频段扫频。这个做法可以在采集系统出错、接线反相、幅值饱和等问题上快速止损比跑完整个扫频再发现数据不可用高效得多。5.3 加一个在线闭环交叉验证开环传递函数的最终验证方式是闭环阶跃响应。把辨识出的开环模型补上控制器构成单位负反馈仿真闭环阶跃和实际系统的阶跃响应对比超调量和上升时间。幅度偏差在 10% 以内说明扫频法求开环传递函数的数据链路基本可靠。如果阶跃响应实测超调明显大于模型预测优先怀疑相位提取偏乐观实际相位裕度小于模型值。这时把扫频激励的幅值降低三分之一重新测一次如果相位曲线比原来更差说明之前的数据已经触碰到系统非线性区。线性度检查就这么简单减小激励幅值模型不变才是合格数据。本文还有配套的精品资源点击获取