
简介面向控制系统分析与设计人员这份MATLAB源码演示了如何利用扫频法实验求解开环传递函数适合自动化、电子工程及相关专业学生、工程师参考。扫频法通过按预设频率范围生成正弦激励信号记录系统在不同频点下的输出响应依据幅度比与相位差估算频率特性进而推算开环传递函数G(s)。程序在2KB的m文件中完整覆盖信号生成、系统激励、响应采集、数据处理与频域估计流程预留了幅频特性/相频特性分析接口并可借助Bode图或Nyquist图进行结果校验。压缩包仅含1个m文件体量小巧、结构清晰便于学习理解与二次修改。已有5697人学习下载资源聚焦频率响应法的核心实现能够帮助读者掌握基于MATLAB的系统辨识与频域建模方法加深对开环传递函数物理意义的工程认知。1. 扫频法求开环传递函数阶跃测不出的对象用它反而容易给含积分环节的开环对象加阶跃输出根本不收敛会一直匀速爬升。想从这段时域波形里拟合传递函数结果对截取长度和噪声都极其敏感往往拟合出一个“看起来能跟阶跃对得上、但 Bode 图完全不对”的模型。扫频法没有这个毛病它把频率连续变化的正弦信号送进对象逐点观察稳态输出的幅值比和相位差直接得到频响再反推 s 域传递函数。对伺服电机、电源环路这类需要交越频率和相位裕度的对象用扫频法实测几乎是默认路径也是我平时写 MATLAB 辨识程序最频繁的场景。下面按“先立原理、再写程序、后设参数、最后验证”展开程序可以直接改着用。2. 从传递函数扫频到频响动手前先理清三件事2.1 频响函数与传递函数的关系扫频直接测到的是什么对一个线性时不变对象输入u A*sin(2*pi*f*t)稳态输出一定是同频率的正弦只有幅值和相位发生变化y A*|G(j*2*pi*f)|*sin(2*pi*f*t φ)。扫频把所有关心的频率都走一遍就得到一组复数点G(jω)这就是对象在激励带宽内的频响函数。这组点本质上是G(s)在虚轴上的采样G(jω) G(s)|_{sjω}。有了足够多的复数点再用有理多项式去拟合就能还原出开环传递函数。扫频法的优势在于它是窄带激励每个时刻主要激发一个频率非线性干扰和宽带噪声不像阶跃法那样全部混在一起。对含积分环节的对象阶跃响应没有稳态但正弦响应在每个频率点都是确定的积分环节只表现为相位滞后 90 度所以照样能测。有一点必须提醒如果对象本身开环不稳定不要在开环状态下直接扫。先用机械限位、内部反馈或闭环把工作点稳住再把扫频信号注入到参考点或误差点原理不变但数据处理时要把已知的闭环关系扣除这属于闭环辨识的范畴。2.2 为什么不能直接拿 FFT 相除从 DFT 到 H1 估计见过不少人直接把Y./U当频响这在理论上成立但工程上基本不能用。原因有二第一输入信号在某些频点能量很低除法会把噪声放大到不可用第二输出端总有测量噪声直接比值没有对噪声做任何平均算出来是单次样本的有偏估计。正确的做法是用互谱和自谱。常见公式是Suy conj(U) * Y / N互功率谱Suu conj(U) * U / N输入自功率谱频响估计用H1 Suy / Suu。互谱计算里输出噪声与输入不相关在多次分段平均后趋近于零所以 H1 对输出噪声是稳健的。在 MATLAB 里不需要手写循环分段直接用cpsd和pwelch组合即可同一段数据可以得到频响和相干函数。2.3 相干函数判断扫频数据可不可信的关键指标相干函数定义为gamma2 |Suy|^2 / (Suu * Syy)取值在 0 到 1 之间。等于 1 表示输出完全由输入线性决定明显低于 1 说明该频点被噪声、非线性或泄漏污染了。这个指标比任何误差棒都实用因为扫频拟合的权重可以直接按相干值来定相干高的频点多给权重相干低的频点剔除。后面第 5 章会用到它做验证但建议在动手写程序前就把“每个频点相干值”作为结果的一部分一起输出。相频曲线也值得单独检查。幅频曲线在某些频点掉下去只是增益问题相频如果有不规则的毛刺或跳变往往意味着该频段信噪比差或存在未对齐的延迟。扫频数据处理完先看相干再看相频最后才看拟合。3. MATLAB程序从生成扫频信号到拟合开环传递函数3.1 生成对数扫频激励比 chirp 函数更可控的写法对数扫频比线性扫频更适合系统辨识因为它在每个十倍频程分配相同时间低频段不会一闪而过。下面这段代码用相位积分的方式生成信号比直接调用chirp更直观也方便改成任意扫频规律。fs 10000; % 采样率与采集卡一致Hz f0 0.5; % 起始频率低于首个转折频率的 1/10Hz f1 200; % 终止频率高于对象带宽 3~5 倍Hz T 60; % 总扫频时长s按 20/f0 起调 A 0.1; % 激励幅值V以输出不进入饱和为准 t (0:1/fs:T-1/fs).; inst_f f0 * (f1/f0).^(t/T); % 对数扫频瞬时频率 phase 2*pi*cumsum(inst_f)/fs; % 相位积分避免相位跳变 u A * sin(phase);cumsum是对瞬时频率做时间积分得到连续变化的相位这样信号在任意时刻都不存在相位突变。瞬时频率表达式f0*(f1/f0)^(t/T)让log(f)随时间线性增长保证每个十倍频程的时间占比一致。如果对象低频响应很慢把T加大即可不需要改结构。生成的u就是给对象的激励。实际测量时用 DAC 播放这个信号同时用 ADC 同步录制输入u和输出y存成两列 csv。输入输出必须用同一时钟源采样否则相位差会随扫描时间漂移后面对齐会很麻烦。3.2 用 readmatrix 导入 csv 并估计频响csv 数据导入 MATLAB 后做 FFT 分析推荐用readmatrix一次性读入两列波形然后分段估计功率谱。分段是必要的如果不分段、直接对整段做 FFTH1 估计退化成单次除法噪声收敛不到零。data readmatrix(sweep_data.csv); % 第 1 列输入 u第 2 列输出 y u data(:,1); y data(:,2); u u - mean(u); % 只去均值不要 detrend 默认的线性趋势 y y - mean(y); % 低频扫频段容易被线性拟合吃掉 seg round(fs / f0 * 5); % 每段包含约 5 个最低频周期 nfft 2^nextpow2(seg); noverlap round(seg * 0.5); win hann(seg, periodic); % 周期 Hann 窗适合功率谱估计 [Suy, f] cpsd(u, y, win, noverlap, nfft, fs); [Suu, ~] pwelch(u, win, noverlap, nfft, fs); [Syy, ~] pwelch(y, win, noverlap, nfft, fs); G Suy ./ Suu; % H1 频响估计 gamma2 abs(Suy).^2 ./ (Suu .* Syy); % 相干函数cpsd计算输入输出互功率谱与pwelch使用相同的窗和分段方式因此Suy ./ Suu的频率轴对齐。段长seg取“最低频 5 个周期”是一个兼顾频率分辨率和分段数量的折中段太长平均次数少段太短低频分辨率不足。如果低频相干仍然不好把seg加倍再试。通道延迟这里不建议用xcorr自动对齐。互相关峰值会把对象本身的群延迟也算进去结果等于把对象的相位特性削掉一块。采集卡通道延迟是固定值直接按采样点数平移即可对象自身的延迟要保留在传递函数里。3.3 用 invfreqs 把频响拟合为传递函数首一和尾一先对齐频响曲线拿到后用invfreqs拟合连续时间传递函数。频率轴要换算成角频率并且只取扫频范围内的频点边界处的泄漏点不要进拟合。sel f f0 * 0.9 f f1 * 0.9; w 2 * pi * f(sel); Gsel G(sel); nb 4; % 分子阶数 na 5; % 分母阶数有积分环节时 na nb 1 [B, A] invfreqs(Gsel, w, nb, na, [], [], 50); Gfit tf(B, A);invfreqs使用迭代算法最后一个参数是迭代次数。默认初始值在多数情况下能收敛但如果拟合结果明显偏离原始频响可以传入权重向量把拟合重点压到交越频率附近例如wt 1 ./ abs(Gsel)可以让低频大增益不再支配误差函数。拟合完要做一次规范化这涉及“开环传递函数增益是首一还是尾一”这个问题的根源。MATLAB 返回的分子分母都是按 s 降幂排列分母最高次项系数不一定为 1。两种常见写法B B / A(1); A A / A(1); % 尾一形式分母最高次项为 1 % 此时直流增益 K B(end) / A(end) K A(end); B B / K; A A / K; % 常数项归一分母写成时间常数形式 % 此时直流增益直接读 B(end)两种写法的传递函数完全相同但增益数值不一样。别人给你的开环增益如果对不上先确认对方用的是哪种分母归一化。含积分环节时没有直流增益K 应理解为误差系数读法同样随分母写法变化。如果对象存在纯延迟观察高频段相位是否持续下滑而幅频平坦。用高频段polyfit(w, unwrap(angle(Gsel)), 1)求斜率tau -dφ/dω然后Gfit tf(B, A, ioDelay, tau)把延迟补进去。4. 扫频参数怎么设先查表再按结果微调4.1 四个关键参数的取值表扫频参数之间相互影响实际调试时按下面这张表起调再根据相干曲线微调。参数起始取值判断依据与微调方向激励幅值 A输出达到额定值的 5%~10%相干低且无谐波就加大输出出现限幅就减小起始频率 f0低于首个转折频率的 1/10低频相干差就继续降 f0或加大 T终止频率 f1高于对象带宽 3~5 倍高频相位平坦但相干骤降说明扫超了收 f1扫频时长 T按 20/f0 起调低频相干低于 0.9 就加倍分段长度 seg约 5 个最低频周期低频频率分辨率不足时加倍幅值是最容易出问题的参数。激励太小高频段输出淹没在噪声里激励太大对象的饱和与非线性会把谐波带到输出端。一个稳妥做法是先用小幅值扫一遍观察输出峰值和相干曲线再逐步加大 A直到相干整体抬升而输出不出现平顶。分段长度很容易被忽略。seg取 5 个最低频周期时0.5 Hz 对应 10 秒一段如果 f0 是 0.05 Hz每段就要 100 秒。整段数据只有几分钟时分段数会很少相干曲线噪声很大这时要接受较低的分段数量不要强行缩短段长。4.2 三个必避的坑起始瞬态、折返频点与谐波污染第一个坑是起始瞬态。扫频信号从零幅值开始虽然sin(0)0没有幅值跳变但信号从无到有本身就是一个宽频激励对象在低频段的响应需要时间建立。处理方法是正式扫频前先以固定频率f0激励10/f0秒让对象的低频响应进入稳态然后再开始扫频。分析数据时从扫频起点截取前面的定频激励段直接丢弃。第二个坑是扫频折返点。扫到f1后信号突然变化这个时刻会产生宽带泄漏污染附近频点的频响估计。不要把分析频率范围取满到f1一般截到f1*0.9就停。如果采用循环扫频每圈结束点的跳变依然存在Welch 分段平均能在一定程度上摊平泄漏但边界频点仍然建议剔除。第三个坑是非线性谐波。大信号激励下饱和、死区、摩擦都会产生 2 次、3 次谐波。H1 估计会把谐波当作与输入不相关的分量摊到噪声里表现为相干函数在某些频段周期性凹陷。如果相干曲线的凹陷间隔恰好等于某个基频的整数倍基本就是非线性在起作用。优先降低 A如果必须在大信号下辨识应改用对数扫频加解调的 Farina 方法分离各次谐波而不是用平均法硬扛。5. 验证扫频结果相干性核查、回代与闭环对照5.1 用相干函数给每个频点打分数据处理好后把幅频、相频、相干三条曲线一起画出来这是扫频法求开环传递函数的第一步验收sel f f0 f f1 * 0.9; figure; subplot(3,1,1); semilogx(f(sel), 20*log10(abs(G(sel)))); grid on; ylabel(幅值 dB); subplot(3,1,2); semilogx(f(sel), rad2deg(unwrap(angle(G(sel))))); grid on; ylabel(相位 deg); subplot(3,1,3); semilogx(f(sel), gamma2(sel)); ylim([0 1]); grid on; ylabel(相干); xlabel(Hz);相干高于 0.9 的频点可以放心使用0.8 到 0.9 之间要看相位是否平滑低于 0.8 的频点建议剔除后再拟合。如果低频段相干大面积偏低优先加长 T 或加大 A高频段相干低多数是输出接近噪声底把 f1 下调或者做分段扫频更实际。相位曲线记得先unwrap否则在 -180 度附近会来回跳人眼很难判断是否平滑。5.2 回代仿真用拟合模型反跑一次扫频把拟合出的Gfit放回 Simulink 或者直接用lsim重新跑一遍同样的扫频输入时域重合度能反映模型在整个频带的整体质量t_rec (0:length(y)-1)./fs; y_sim lsim(Gfit, u, t_rec); plot(t_rec, y, t_rec, y_sim);如果 Bode 图对得上但时域波形尾部偏差大多半是起始瞬态没有处理干净如果整体相位跟不上检查是否遗漏了 ioDelay。相位延迟可以从高频段unwrap(angle(G))对 ω 的直线斜率粗估用tf(B, A, ioDelay, tau)补上再拟合一次。5.3 用 margin 和 pzmap 做最后一道检查开环传递函数拟合出来是为了算交越频率和相位裕度直接看结果[Gm, Pm, Wcg, Wcp] margin(Gfit);对照扫频原始数据检查Wcp处的相位是否与实测一致。不一致时用invfreqs的权重向量把拟合重点压到交越频率附近重新拟合。最后跑一次pzmap(Gfit)如果出现离原点很远、且零点极点几乎重合的“凑对”结构说明模型阶数给高了把 na 降 1 再拟合通常比硬扛高阶模型更干净。本文还有配套的精品资源点击获取