企业官网建设流程全解析

简介：面向信号处理与控制系统方向的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决定谱线间隔，fs=1000、nfft=4096时理论谱线间隔约为 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 窗是默认选择，但不是所有场景都合适。窗函数决定了主瓣宽度和旁瓣衰减之间的权衡，主瓣越宽，频率分辨能力越差；旁瓣越低，泄漏到大频率区间的能量越少。

其中N是窗长度。实际操作中，如果两个频率成分只差几赫兹，主瓣要尽量窄，Hann 和 Hamming 都行；如果关心的是比主信号低 60 dB 的小幅值分量，Blackman 能减少大信号泄漏对小信号的掩盖，代价是峰值被展宽。我一般默认用 Hann，只有当分析结果里出现明显的“裙边”时才换成 Kaiser 并增大 β。

2.4 数据预处理不能省

pwelch虽然能处理绝大部分信号，但对异常段很敏感。一个数据丢失片段或一个尖峰，会抬升整个频段的噪声底。常规流程是先在Powerspectrum外面对数据做三件事：

1. 用detrend(x)去掉线性趋势，尤其是采集卡零漂明显时。
2. 截断掉首尾异常段，比如由继电器吸合产生的尖峰。
3. 若关心的是几百赫兹以上的成分，先做一次低通滤波，避免高频混叠分量折叠回低频带。

需要注意，滤波会改变信号频谱形状，不能想当然地“滤干净”。如果目标就是看原始信号的频率构成，预处理越少越好；如果目标是给系统辨识铺路，滤波尽量只在已知干扰频带附近使用。

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 而不是 linspace

Bode 图横轴是对数刻度，频率点也应当在对数空间均匀分布。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 的选择主要看你想关注的频率间隔。两个频率相隔 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);

这里FreqUnits='Hz'和前面手动除以2*pi的效果一样，但不再需要改数据。MagUnits='abs'适合观察增益接近 1 或接近 0 的窄带系统，如果看宽频响应还是默认的 dB 更直观。PhaseWrapping='on'会把相位限制在 ±180 度以内，和 unwrap 后的曲线是两种展示风格，通常在发布给团队的报告里用 unwrap 版本，在控制设计讨论里用缠绕版本。

5.2 导出矢量图避免放大失真

频谱图和 Bode 图经常要放到 PDF 报告甚至论文里，位图放大后线条和文字会糊。MATLAB 里最稳妥的方式是用exportgraphics，它支持矢量输出且自动裁剪空白边。

exportgraphics(gcf, 'bode_response.pdf', 'ContentType', 'vector');

ContentType='vector'确保坐标轴、曲线和网格都保存为矢量，缩放不损失细节。老版本 MATLAB 没有exportgraphics时，用print('bode_response','-dpdf','-vector')也可以达到同等效果。最后单独说一句：如果你的脚本要换电脑跑，bodeoptions里的对象结构在不同版本 Toolbox 下略有差异，报错时优先检查'cstprefs'这个默认预置组是否存在，不存在时直接去掉第一个参数即可。

本文还有配套的精品资源，点击获取