ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

MATLAB频谱分析与Bode图绘制:pwelch与频率响应实践

MATLAB频谱分析与Bode图绘制:pwelch与频率响应实践 简介面向信号处理与控制系统方向的MATLAB学习者这份资源聚焦频谱绘制与Bode图绘制两大实用技能适合需要快速上手频域分析的本科生、研究生或工程师。压缩包内共2个文件均为.m脚本分别是基于pwelch函数实现功率谱估计与绘图的频谱分析脚本以及调用bode/bodeplot展示系统频率响应的Bode图脚本整体仅1KB精简而直接。目前已有1648人学习下载对该类型入门任务有较高参考价值。通过学习这两个脚本可掌握数据预处理、窗函数选择、频率分辨率设定以及传递函数建模、频率范围调整与幅频/相频呈现的完整链路并能迁移到噪声分析、滤波器设计及控制系统稳定性评估等实际场景。1. 为什么这两个脚本总是一起出现实际项目里只要压缩包里同时出现Powerspectrum.m和plotmakebode.m说明后半段工作不是“画个波形”这么简单而是要回答两个问题采集到的信号里有哪些频率成分系统在这些频率上的增益和相位分别是什么。前者把时域采样用pwelch转成功率谱密度后者把传递函数或状态空间模型画成标准Bode图。伺服系统调试、振动噪声分析、滤波器设计验证都会用到这组工具。直接调用fft或freqz虽然快但坐标单位、窗函数增益、模型形式不一致拿到真实数据时很难和理论曲线对上。这篇把两个脚本从函数选型到参数调整、再到联用时的对比逻辑讲清楚重点是每一步哪个参数在起作用改它会发生什么。2. 频谱绘制背后的pwelch周期图、窗函数和分辨率2.1 从周期图到 Welch 功率谱估计直接用plot(abs(fft(x)))画频谱在演示中能看工程上不行。原因有两层一是 FFT 离散化带来的栅栏效应让峰值落在两条谱线之间时幅度偏低二是单次周期图的方差很大相邻频率点起伏明显把真实谱形状淹没掉。Welch 把数据分成若干段每段加窗后做 FFT再对功率谱取平均段数越多估计方差越小但每个段越短频率分辨率越差。这一对矛盾就是pwelch参数设计的核心。Powerspectrum.m的角色就是把这段逻辑封装成固定模板。不同项目里数据来源不同但核心调用方式是一致的准备好列向量x和采样率fs然后交给pwelch。2.2 Powerspectrum.m 的骨架与 pwelch 参数拆解一个典型脚本长这样function [pxx, f] Powerspectrum(x, fs) % 输入: x 时域信号向量, fs 采样率(Hz) % 输出: pxx 单边功率谱密度, f 频率轴(Hz) x x(:); % 统一为列向量 x x - mean(x); % 去直流避免 0 Hz 处尖峰 nfft 4096; % FFT 点数 win hann(nfft, periodic); % 周期 Hann 窗频谱泄漏较小 nol round(nfft * 0.75); % 75% 重叠数据量不够时降为 50% [pxx, f] pwelch(x, win, nol, nfft, fs); pxx 10 * log10(pxx); % 转 dB 便于观察动态范围 plot(f, pxx, LineWidth, 1); xlabel(频率 (Hz)); ylabel(功率谱密度 (dB/Hz)); grid on; end代码里值得说的有几点。x(:)是为了防止行向量和列向量混用导致后续计算维度错误这是数据采集脚本最常见的低级问题。mean(x)去直流非常关键如果不去掉0 Hz 处会有一个幅度极大的谱峰它会掩盖低频段的真实细节。nfft决定谱线间隔fs1000、nfft4096时理论谱线间隔约为 0.244 Hz但这不代表分辨率一定有 0.244 Hz实际分辨率受窗长度限制不能靠单纯增大nfft来“看清”更近的频率。nol是重叠采样点加窗后每段两端被压小重叠能让数据被重复利用一般取 50% 到 75%。最后转成 dB 是因为功率谱动态范围可能达到几个数量级线性坐标只能看出最大值附近的情况。pwelch默认输出单边功率谱密度单位是幅值平方/Hz。如果采样率fs用 Hz频率轴f就是 Hz如果数据来自旋转机械有人喜欢用fs*60转成 CPM单位会跟着变。这个脚本把结果存进pxx和f而不是直接打印是为了方便后面跟 Bode 图做对比。2.3 窗函数怎么选Powerspectrum.m里用的 Hann 窗是默认选择但不是所有场景都合适。窗函数决定了主瓣宽度和旁瓣衰减之间的权衡主瓣越宽频率分辨能力越差旁瓣越低泄漏到大频率区间的能量越少。窗函数主瓣宽度近似旁瓣衰减适用场景Hann4π/N-31 dB多数信号分析泄漏与分辨率折中Hamming4π/N-43 dB与 Hann 接近旁瓣更低但第一旁瓣衰减快Kaiser(β6)约 4.5π/N-50 dB 以下需要手动调节主瓣/旁瓣权衡Blackman6π/N-58 dB强干扰附近需要压低旁瓣分辨率让步其中N是窗长度。实际操作中如果两个频率成分只差几赫兹主瓣要尽量窄Hann 和 Hamming 都行如果关心的是比主信号低 60 dB 的小幅值分量Blackman 能减少大信号泄漏对小信号的掩盖代价是峰值被展宽。我一般默认用 Hann只有当分析结果里出现明显的“裙边”时才换成 Kaiser 并增大 β。2.4 数据预处理不能省pwelch虽然能处理绝大部分信号但对异常段很敏感。一个数据丢失片段或一个尖峰会抬升整个频段的噪声底。常规流程是先在Powerspectrum外面对数据做三件事用detrend(x)去掉线性趋势尤其是采集卡零漂明显时。截断掉首尾异常段比如由继电器吸合产生的尖峰。若关心的是几百赫兹以上的成分先做一次低通滤波避免高频混叠分量折叠回低频带。需要注意滤波会改变信号频谱形状不能想当然地“滤干净”。如果目标就是看原始信号的频率构成预处理越少越好如果目标是给系统辨识铺路滤波尽量只在已知干扰频带附近使用。3. Bode 图绘制plotmakebode.m里的模型与频率向量3.1 Bode 图读什么以及模型怎么进Bode 图是频域里最直观的系统描述方式横轴是对数频率上图画幅值单位常用 dB下图画相位单位是度。闭环稳定性讨论里的增益裕度和相位裕度就是从这两条曲线上读出来的。MATLAB 的bode函数接受tf、ss、zpk三种模型对象离散系统也可以直接用z tf(z, ts)构造后传入。这里有个容易忽略的点bode计算的是频率响应采样值它需要知道沿虚轴取哪些点。如果不给频率向量MATLAB 会按照系统动态特性自动选点结果是曲线很平滑但多次对比时可能对不上。plotmakebode.m存在的意义就是把频率向量、线型、坐标范围固定下来让每次画出的 Bode 图风格一致并且能和Powerspectrum.m的输出放在同一张图里。3.2 频率向量用 logspace 而不是 linspaceBode 图横轴是对数刻度频率点也应当在对数空间均匀分布。linspace(0, 1000, 300)会把大量点集中在高频段低频段每十倍频程只有几个点画出来的渐近线是折线。正确做法是用logspace。function plotmakebode(G, w) % 把系统 G 的 Bode 图按自定义样式画出来 % w 为对数均匀分布的角频率向量单位 rad/s if nargin 2 || isempty(w) w logspace(-2, 4, 300); % 默认 0.01 ~ 10000 rad/s end [mag, phase] bode(G, w); % 返回幅度与相位数组 mag squeeze(mag); % SISO 时变成一维向量 phase squeeze(phase); subplot(2,1,1); semilogx(w, 20*log10(mag), b-, LineWidth, 1.4); grid on; ylabel(幅值 (dB)); title(Bode 图 - 幅度); subplot(2,1,2); semilogx(w, phase, r-, LineWidth, 1.4); grid on; ylabel(相位 (deg)); xlabel(频率 (rad/s)); title(Bode 图 - 相位); end代码里squeeze很关键。bode(G, w)在单输入单输出系统上返回的是多维数组维度是(输出数, 输入数, 频率点数)不压缩的话semilogx会报维度错误。20*log10(mag)是把线性幅值转成 dB这一变换让增益从 0.001 到 1000 的变化范围能画在同一张图上。logspace(-2, 4, 300)表示从 0.01 到 10000 rad/s 取 300 个对数间隔点每十年约 50 个点足够画出平滑曲线。如果你已经知道转折频率w0更聪明的方式是只取转折频率附近的范围比如w logspace(log10(w0*0.2), log10(w0*5), 400)。这样既能看到低频渐近线又能看到高频滚降图的宽度不会被无关频段占满。3.3 相位跳变与 unwrap连续系统的相位通常从 0 度渐变到负几百度的某个值但 MATLAB 默认把相位限制在 ±180 度之间。结果是从 -180 度再往下画会跳到 180 度图上是锯齿形容易误判相位裕度。比较理论 Bode 和实测频谱之前先对相位做一次展开。phase squeeze(phase); phase unwrap(phase * pi / 180) * 180 / pi;unwrap的本质是检测相邻采样点之间的跳变如果跳变超过 180 度就加减 360 度把它接成连续曲线。这个操作对按频率顺序采样的w有效如果w是乱序的必须先排序再 unwrap。plotmakebode.m没有默认做这一步因为展开后的曲线在某些情况下反而不直观写脚本时建议把两版都试一下。4. 把频谱和 Bode 图对上号数据换算、参数调节与典型坑4.1 实测频谱和理论 Bode 图如何对齐频谱分析的对象是“信号”Bode 图描述的是“系统”。两者要在同一张图上对比必须满足一个前提已知输入激励且系统近似线性。常见做法是给系统一个白噪声或扫频激励采集输入输出信号然后用tfestimate直接估计频率响应。但在脚本较简单的场景里用户往往只有输出数据此时可以用理论关系输出功率谱约等于输入功率谱乘以系统幅频响应的平方。最直接的对比方式是把 Bode 图的横轴从 rad/s 改成 Hz因为Powerspectrum.m里pwelch输出的是 Hz。代码这样写% 把角频率转换成 Hz方便与 pwelch 的结果叠加 w logspace(1, 3, 200); f_hz w / (2*pi); [mag, ~] bode(G, w); mag_db 20 * log10(squeeze(mag)); plot(f_hz, mag_db, LineWidth, 1.5); hold on;这里w/(2*pi)是把角频率转换成循环频率。转换之后Bode 图横轴的单位和pwelch返回的f一致低频段和高频段的启停位置也能对齐。需要注意的是转换横轴不会改变系统特性但会影响你观察转折频率的直觉1 Hz 对应 6.28 rad/s不要在两条曲线上错位比较。4.2 参数调整速查表两个脚本联用时的参数选择有规律可循整理成一张表便于在真机上快速试参数作用典型值NFFT决定 FFT 插值密度谱线间隔 fs/NFFT1024 ~ 8192noverlap重叠采样数提升平均次数50% ~ 75% 的 NFFTnfft 分段长度频域插值不增加真实分辨率不推荐超出太多w 下限Bode 图最低频点0.1 倍转折频率w 上限Bode 图最高频点10 倍转折频率每十年点数曲线平滑度50 ~ 200NFFT 的选择主要看你想关注的频率间隔。两个频率相隔 10 Hz采样率 1000 Hz 时NFFT 至少要 256实际因为加窗会再模糊一些建议直接取 1024 以上。重叠率不是越高越好75% 对随机噪声够用再高计算量翻倍但方差改善有限。Bode 图的频率范围如果太宽转折频率附近的分辨率会被压缩所以先通过一次粗略扫描找到转折区再缩小范围细化。4.3 实际系统里常见的三个坑第一个坑是相位图锯齿。上一节已经提到unwrap的必要性实际工程里还会遇到另一种情况系统本身有纯延迟环节相位随频率持续下降即使 unwrap 后也可能降到上千度。这时候不能简单比较表观相位要看减去延迟项w*Td之后的残余相位。第二个坑是直流分量被去掉后低频段对不上。Powerspectrum.m里做了x - mean(x)这等于把直流到极低频的功率全部扔掉。如果理论 Bode 图在低频段增益很高而实测频谱在 0.1 Hz 以下几乎没有能量两者无法直接对比。解决方法是保留直流或只去掉缓慢趋势并把对比的最低频率设到转折频率以下 10 倍即可。第三个坑是窗函数增益导致绝对幅值偏移。pwelch对每段数据加窗后窗形状会改变段总能量。如果两个脚本分别采用不同窗比较时会发现整体平移了几个 dB。修正方法是把谱密度除以窗均值让能量恢复为接近原始数据win hann(4096, periodic); [pxx, f] pwelch(x, win, [], 4096, fs); pxx pxx / mean(win); % 能量修正补偿窗函数对幅值的加权这个修正只影响绝对电平不影响峰值的相对位置。做滤波器通带增益对比时除不除差很多只做故障特征频率识别时不影响结论。我的习惯是两种都保留修正后的版本用于和 Bode 图幅值对齐未修正的版本用于噪声底监测。5. 用 bodeoptions 和 exportgraphics 把两张图沉淀成可交付文件5.1 bodeoptions 一次设定全部样式plotmakebode.m手动写semilogx灵活度高但如果要嵌入 Simulink 或 Control System Toolbox 环境直接使用bodeplot加bodeoptions更省事。它可以一次性设置坐标、单位、相位缠绕和标题避免每张图都重复写绘图属性。opts bodeoptions(cstprefs); opts.FreqUnits Hz; % 坐标系用 Hz便于与频谱图拼接 opts.MagUnits abs; % 不转 dB直接看线性增益 opts.PhaseWrapping on; % 相位限制在 ±180 度 opts.XLim {[0.5 500]}; opts.Title.String 闭环系统 Bode 图; opts.Grid on; opts.XLabel.String 频率 (Hz); h bodeplot(G, w, opts);这里FreqUnitsHz和前面手动除以2*pi的效果一样但不再需要改数据。MagUnitsabs适合观察增益接近 1 或接近 0 的窄带系统如果看宽频响应还是默认的 dB 更直观。PhaseWrappingon会把相位限制在 ±180 度以内和 unwrap 后的曲线是两种展示风格通常在发布给团队的报告里用 unwrap 版本在控制设计讨论里用缠绕版本。5.2 导出矢量图避免放大失真频谱图和 Bode 图经常要放到 PDF 报告甚至论文里位图放大后线条和文字会糊。MATLAB 里最稳妥的方式是用exportgraphics它支持矢量输出且自动裁剪空白边。exportgraphics(gcf, bode_response.pdf, ContentType, vector);ContentTypevector确保坐标轴、曲线和网格都保存为矢量缩放不损失细节。老版本 MATLAB 没有exportgraphics时用print(bode_response,-dpdf,-vector)也可以达到同等效果。最后单独说一句如果你的脚本要换电脑跑bodeoptions里的对象结构在不同版本 Toolbox 下略有差异报错时优先检查cstprefs这个默认预置组是否存在不存在时直接去掉第一个参数即可。本文还有配套的精品资源点击获取
返回列表