尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

MATLAB数字信号处理实战:从FFT频谱分析到包络谱故障诊断

MATLAB数字信号处理实战:从FFT频谱分析到包络谱故障诊断 简介《数字信号处理的MATLAB实现》万永革著高清扫描版配齐了目录适合数理基础薄弱、重视动手实践的本科生、研究生和工程技术人员自学信号处理。压缩包共476个文件、约54.24MB其中364个m脚本对应各章可运行示例62个06数据文件与7个dat、5个mat文件提供实验数据24个txt和3个tex辅助说明与源码另有pdf主文档和doc说明便于对照学习。资源已有527人浏览学习内容结合大量实例阐述数字信号处理的实际应用。读者可通过配套程序直接运行或稍作修改快速实现滤波、谱分析、时频分析等典型算法如tfrpage.m涉及的时频处理既能验证教材理论也能迁移到自己的课题或工程问题中是一份动手型学习资料。1. 从万永革这本书说开去MATLAB做数字信号处理的基本盘如果你手头也有一本《数字信号处理的MATLAB实现》大概率你和我一样刚翻开时信心满满翻到第三章开始怀疑人生到了滤波器那一章干脆直接复制代码跑一下就算“学完了”。万永革这本书的好处是例程密集、覆盖全面从离散时间系统、傅里叶变换到自适应滤波几乎都有涉及而且很多章节可以直接抄到一段能运行的MATLAB代码。但问题是抄能跑不代表你理解它为什么这么写等你自己拿一段实测信号做频谱分析、滤波、包络谱的时候各种隐藏的坑就全冒出来了。这篇文章不是什么读书笔记而是在我折腾了几年信号处理项目之后回头看这本书里那些例程时觉得真正值钱、也最容易卡住人的几个技术点。适合刚补完信号与系统基础、准备用MATLAB处理实测数据的人也适合已经跑过书里代码但总觉得“差点意思”的朋友。我不按目录讲就挑实际工作中高频出现的四个问题展开频谱分析怎么画才对、IIR滤波器怎么用才稳、希尔伯特变换怎么提取包络、以及怎么把书写式的例程改造成自己能用的工具。1.1 这本书教了什么没教什么万永革这本书的核心价值是把数字信号处理的抽象公式和MATLAB函数一一对应起来。比如conv对应卷积fft对应DFTfreqz对应系统频率响应butter对应模拟/数字滤波器设计。对初学者来说这种“公式-函数”的映射关系比死磕教科书高效得多。但它的局限也很明显书上每个例子都是孤立的信号是干净的、参数是预设好的。真实数据往往带噪声、有趋势项、长度不是理想的2的幂这些问题书上很少展开。所以我的建议是把这本书当成“函数用法词典”和“算法概念入门”而不要指望抄几段代码就能直接处理现场数据。下面要讲的几个点正是从“抄例程”到“做项目”之间最容易摔跟头的地方。2. 频谱分析最容易踩的坑采样率、频率轴和点乘2.1 FFT之后到底该画什么很多人仿照书上的代码写完X fft(x); plot(abs(X));就认为频谱出来了。这个图往往是一堆对称的尖峰横轴是“点数”而不是“赫兹”而且峰值高度也和信号的真实幅值对不上。问题出在三个地方没取模、没做单边谱、没标定频率轴和幅值。正确的做法是先把频率轴建立起来。假设采样率fs信号长度N那么FFT结果的第k个点对应的频率是k*fs/N范围从0到fs。对于实信号频谱关于fs/2对称所以一般只看前半段。fs 1000; % 采样率 1000 Hz t (0:1023)/fs; % 1.024 秒1024 个点 x 2*sin(2*pi*50*t) 1.5*sin(2*pi*120*t) 0.8*randn(size(t)); N length(x); X fft(x); f (0:N-1)*fs/N; % 完整频率轴 % 取单边频谱 X_half X(1:N/2); f_half f(1:N/2); mag abs(X_half)/N; % 先除以 N mag(2:end) mag(2:end)*2; % 除直流外单边谱幅值乘 2 plot(f_half, mag); xlabel(频率 (Hz)); ylabel(幅值);这段代码执行完50Hz处峰值接近2120Hz处接近1.5噪声底大约在0.05以下。为什么交流分量要乘2因为单边谱把原来对称分布在正负频率上的能量合并到了正频率一侧所以幅值加倍直流分量本身不翻倍。如果不做这一步你画出来的幅值永远是真实值的一半。2.2 FFT结果的复数到底怎么理解FFT输出的每个点都是一个复数实部对应余弦分量的幅度虚部对应正弦分量的幅度。取模就是该频率分量的合成幅值angle(X)得到的是初始相位。实际分析中幅值谱最重要相位谱往往只在需要重构信号时才用。另外要注意FFT点数不一定是信号长度。如果信号很短可以做补零FFT来“插值”频谱X fft(x, 4096)这样频率轴更密曲线更平滑。补零不改变频谱分辨率分辨率只取决于真实信号时长N/fs但会让谱线看起来更细腻。书里很少讲这个区别我在项目里经常用它来定位窄带峰值。2.3 点乘和直接乘的分水岭MATLAB里*是矩阵乘法.*才是逐元素乘法。信号处理中绝大部分操作都是逐元素的比如加窗、混频、信号相乘。很多人第一次写加窗代码时会写成x * window结果要么报维度不匹配的错要么意外得到矩阵外积。window hann(N); % 行向量窗函数 x_windowed x .* window; % 正确逐元素乘我见过最隐蔽的坑是列向量和行向量混用。x是列向量window是行向量用.*也会触发广播机制生成一个N×N矩阵而不是N×1向量。代码不报错但结果数据量瞬间膨胀几十倍内存小一点直接卡死。规矩很简单做逐元素运算前先确认两个向量的维度一致统一用x(:)把它变成列向量或者用transpose统一方向。3. 滤波器设计实战IIR函数背后那些书上不写的门道3.1 归一化频率和滤波器阶数书里介绍IIR滤波器通常从butter开始[b,a] butter(2, Wn, low)。这个Wn不是赫兹而是相对于奈奎斯特频率fs/2的归一化频率。比如采样率1000Hz想截止到100HzWn 100/(1000/2) 0.2。忘记归一化是最常见的错有人直接填Wn 100因为频率轴超界filter返回的结果完全是乱的。阶数选择也需要概念。butter(2,...)是二阶滤波器每倍频程衰减约12dB如果信号里的带外干扰很强二阶往往不够。但阶数不是越高越好高阶IIR滤波器的数值稳定性会变差尤其是截止频率很低的时候。3.2 filter与filtfilt零相位不是免费的同样是滤波filter和filtfilt结果差异肉眼可见。filter是因果滤波器输出相对输入存在相位延迟阶数越高延迟越大。filtfilt先正向滤波再反向滤波两次滤波的相位延迟相互抵消实现零相位偏移波形上的特征点位置不会漂移。离线数据分析场景我几乎总是用filtfilt。比如处理振动信号要定位故障冲击发生的时刻相位延迟会直接导致特征位置错位。代价是filtfilt不能用于实时处理因为它需要整段数据才能做反向滤波。[b, a] butter(2, 0.2, low); y_causal filter(b, a, x); y_zero filtfilt(b, a, x);对比一下滤波后的上升沿y_causal明显滞后于原始信号y_zero的边缘和原始信号对齐。另外filtfilt的等效滤波器阶数是原来的两倍边界效应也更明显数据两端会出现小幅振铃处理短序列时要警惕。3.3 高阶滤波器的数值稳定性高阶IIR直接写成[b,a]再调用filter在极端截止频率下可能因为极点位置偏差导致输出爆炸。更稳妥的做法是用二阶节SOS形式[b, a] butter(10, 0.2, low); sos tf2sos(b, a); y sosfilt(sos, x);或者直接用designfiltd designfilt(lowpassiir, FilterOrder, 10, ... HalfPowerFrequency, 100, SampleRate, fs); y filtfilt(d, x);designfilt返回滤波器对象内部会自动处理数值优化代码也更可读。我自己的习惯是简单场景用butterSOS复杂设计带通、带阻、多频段直接用designfilt少碰裸[b,a]能省掉大量调试时间。4. 希尔伯特变换与包络谱故障诊断里最常用的解调手段4.1 解析信号与包络提取旋转机械的故障信号通常表现为“高频载波低频调制”比如滚动轴承内圈故障会产生以固有频率为载波、以故障特征频率为调制波的非平稳振动。直接对原始信号做FFT载波频率附近会出现一堆边带特征频率很难分辨。这时候就要先提取包络再对包络做频谱分析得到包络谱。MATLAB里hilbert(x)返回的是解析信号不是传统意义上的希尔伯特变换结果。解析信号的实部是原信号虚部是原信号的90度相移版本所以包络就是解析信号的模fs 2000; t (0:1999)/fs; fc 500; % 载波频率 fm 10; % 调制频率 x (1 0.5*cos(2*pi*fm*t)) .* cos(2*pi*fc*t); x x 0.05*randn(size(t)); % 加一点噪声 analytical hilbert(x); envelope abs(analytical);到这里得到的是包含直流分量的包络。做包络谱前通常要去均值envelope envelope - mean(envelope)否则零频处会有一个巨大的峰值。然后用和第二章一样的FFT流程处理包络就能在10Hz处看到一个清晰的谱峰这就是调制频率。4.2 完整包络谱处理的建议流程实际项目里我不会直接对原始信号做hilbert而是先做带通滤波把关心的高频载波频带选出来再做包络解调。原因很简单hilbert对宽带噪声很敏感如果信号里除了目标频带还有其他强干扰包络里会混入大量无用成分包络谱变成一锅粥。推荐的流程是带通滤波 → 截取稳定段 → 去均值 →hilbert取包络 → 包络去均值 → FFT得到包络谱。例如先带通到450~550Hz再提取包络得到的谱线会干净得多。d designfilt(bandpassiir, FilterOrder, 8, ... HalfPowerFrequency1, 450, HalfPowerFrequency2, 550, ... SampleRate, fs); x filtered filtfilt(d, x); analytic hilbert(x_filtered); env abs(analytic); env env - mean(env); % 接下来对 env 做 FFT画包络谱写代码请注意hilbert默认沿第一个非单一维度运算。如果你的信号是行向量输出解析信号也是行向量方向要一致。我遇到过因为维度不一致导致后续abs结果长度不对的情况排查了很久才发现是x和x_filtered一个行一个列。统一用x x(:)转成列向量能省掉这个烦恼。5. 把脚本变成工具调试、绘图和工程化的经验5.1 绘图里不起眼但很影响效率的几个细节信号处理离不开可视化但很多人在图上浪费的时间比调算法还多。第一个是局部放大别用鼠标缩放完再截屏直接在代码里控制xlim和ylim比如看0~50Hz范围内的频谱就写xlim([0 50])。如果要对比多个通道可以用subplot统一坐标范围或者用linkaxes同步缩放。第二个是颜色映射。MATLAB默认的parula色系比老版jet更适合显示连续数据但很多人习惯用colormap(jet)喷气色在高亮区域会产生伪轮廓给人视觉上造成虚假的边界。做热力图或时频图时建议用colormap(parula)或colormap(turbo)。这点在书里基本没提但对判断频谱图细节影响很大。第三个是坐标轴层级问题。用polarplot画方向谱时刻度文字经常被数据曲线盖住一般可以通过ax gca; ax.Layer top;把坐标轴放到最上层或者单独设置刻度字体大小。这个冷门技巧我是在某次画阵列方向图时被折腾了很久才找到的。5.2 MATLAB环境问题排查除了算法环境问题也常卡人。安装MATLAB后双击图标一闪而过或直接打不开通常不是软件坏了而是许可证未激活、临时目录权限不够或者显卡驱动和OpenGL不兼容。我试过最快的定位方法是从命令行启动MATLABmatlab -desktop -logfile /tmp/matlab.log启动失败时日志里一般会给出具体原因。常见解决方案包括更新显卡驱动、设置环境变量MATLAB_USE_OPENGL0、把临时目录改到有权限的位置。做完这些90%的启动问题都能解决。5.3 函数封装和调试习惯书里的例程大多是一段脚本变量全堆在基础工作区。数据量小还没什么项目一复杂就乱套。我的做法是把每个功能封装成函数并用inputParser做参数解析function [f, mag] myFFT(x, fs) % 输入 x: 信号, fs: 采样率 % 输出 f: 频率轴, mag: 单边幅值谱 N length(x); X fft(x); f (0:N-1)*fs/N; X_half X(1:floor(N/2)); f_half f(1:floor(N/2)); mag abs(X_half)/N; mag(2:end) mag(2:end)*2; end调试代码时养成两个习惯一是dbstop if error程序报错时自动停在错误行不用在终端里翻栈二是关键中间变量用assert检查维度比如assert(length(x)length(window))能在错误发生后第一时间暴露来源而不是等到画图发现数据不对再去猜。6. 怎样把万永革这本教材真正变成自己的知识体系6.1 从抄代码变成改代码我推荐的用法不是从头到尾读而是把它当工具书带着问题去查。比如你手头有一段振动信号想看有没有轴承故障直接翻到希尔伯特变换那一章把例程里的正弦替换成自己的数据然后观察输出。跑通之后再思考三个“为什么”为什么例程里要先滤波为什么包络谱的横轴单位是Hz为什么0Hz处有个大峰值把这三个问题想明白你对这个工具的理解就超过只抄代码的人了。6.2 用一个端到端的小项目收尾我建议你动手做一个完整的流程生成一段调幅信号加噪声带通滤波提取包络画包络谱识别调制频率。这个流程走完基本涵盖了数字信号处理最常用的几个模块而且每一步都能和万永革书里的章节对应上信号生成对应第三章FFT对应第四章滤波对应第五章希尔伯特变换对应第七八章。做完之后再回到书里看那些推导你会突然发现原来公式和代码是同一件事。我自己在做这个流程时最大的体会是MATLAB只是表达工具真正的门槛在于你知道每个函数返回什么、每个参数物理上代表什么。书不会替你踩这些坑但你可以拿它当一个可靠的起点。遇到问题先跑一遍例程再动手改参数再解释变化原因这个循环坚持下来数字信号处理的能力提升会非常快。本文还有配套的精品资源点击获取
返回列表