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 图会看到锯齿状折线。处理方法是对相频曲线做 unwrap：

phase_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 开环传递函数增益按首一还是尾一写

热词里提到的“开环传递函数增益是首一还是尾一”是一个值得说清楚的问题。传递函数的增益写法不是纯粹的数学格式问题，它直接影响辨识参数的物理意义和数值稳定性。

首一形式是把分子多项式写成最高次项系数为 1：

G(s) = (s + z1) / (s^2 + a1·s + a2)

尾一形式是把常数项归为 1：

G(s) = K · (1 + s/z1) / (1 + s/p1)(1 + s/p2)

对开环传递函数而言，尾一形式更适合描述实际系统的增益预算。理由很直接：开环 Bode 图的低频段增益由K/(s^n)决定，K是环路增益设计的核心参数，写成尾一后 K 直接代表了直流增益或积分增益，工程上可以直接和理论计算的环路增益对照。首一形式下的分子最高次系数为 1，低频增益隐含在其他参数里，无法一眼看出开环增益的大小。

MATLAB 的tfest输出默认是首一形式，对含积分环节的系统，首一形式会让低频增益表征不直观，数值也容易偏大。如果你更关心开环增益是否达标，辨识后可以用zpk转换成零极点增益形式，再手工整理成尾一形式。下面这行代码把估计结果转换并显示增益 K：

sys_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 个点相位趋势正确、穿越频率附近的增益斜率合理，再开始全频段扫频。这个做法可以在采集系统出错、接线反相、幅值饱和等问题上快速止损，比跑完整个扫频再发现数据不可用高效得多。