傅里叶变换从原理到工程实践:频谱分析、泄漏与故障诊断
说实话傅里叶变换这四个字我前前后后系统学过不下三遍。本科《信号与系统》一遍考研复习又一遍工作后做振动故障诊断还被迫啃了一遍。前两遍是为了考试第三遍才是真正为了用。直到我把一个轴承故障从频谱图里揪出来才忽然明白当年背得滚瓜烂熟的公式本质上就是一个极其朴素的思想把信号拆成不同频率的“纯音”叠在一起。这篇文章不打算把公式推导从头再来一遍而是想用做项目的方式把“傅里叶变换”这套东西串起来。适合那些对傅里叶变换有点概念、但总觉得隔着一层窗户纸的人也适合做信号处理、音频分析、振动检测、数据科学的朋友想在实际工程里把它用明白、用正确。1. 傅里叶变换到底在干什么换一个坐标系看问题1.1 果汁、合唱团和坐标变换三个直观类比很多人一上来就被 (X(f)\int_{-\infty}^{\infty} x(t)e^{-j2\pi ft}dt) 这个积分吓住了觉得傅里叶变换是数学家的玩具。其实它干的事情特别接地气——做成分分析。你面前有一杯混合果汁橙汁和苹果汁混在一起。你喝一口能感觉出“又橙又苹果”但很难直接说出橙汁占百分之几、苹果汁占百分之几。傅里叶变换就是那个帮你做成分分析的仪器它把一段信号拆解成不同频率成分然后告诉你每个频率的“浓度”有多大也就是幅度以及这个成分在时间上偏移了多少也就是相位。另一个例子是合唱团。几十个人同时唱歌你闭着眼也能听出其中有没有人走音能分辨出高音部和低音部。你的耳朵本质上就在做傅里叶变换——把叠加在一起的声波拆成不同频率的音高。人脑天生有这个能力只是大家没意识到这也叫“频域分析”。我再加一个数学一点的类比坐标系。同一个人站在不同坐标系里坐标值完全不同但他还是他。信号也一样你在示波器上看到的波形是它在“时间坐标”下的表示经过傅里叶变换换到“频率坐标”下重新描述同一个信号信息量一点没丢。想通这一点后面所有公式的神秘感都会消失。1.2 时域里乱成一团的信号频域里可能干净得惊人为什么工程师总喜欢把信号从时域搬到频域因为很多在时域里极其复杂、完全看不出规律的东西到了频域里都是一目了然的。举个例子一个由50Hz、120Hz、200Hz三个正弦波叠加起来的信号你在时域里看它的波形会觉得那是一团乱麻毫无规律。但你把它做一次傅里叶变换画出来的频谱图就是干干净净三根竖线每根线的高度对应这个频率成分的强弱。滤波、调制解调、振动特征提取、音频均衡全是靠这个“换坐标系”的本事。不过我得提醒一句傅里叶变换不是万能的。它把信号当成一堆“从头到尾频率都不变”的纯正弦波的叠加。如果一个信号的频率随时间变化比如语音、音乐、机器转速变化时的振动常规傅里叶变换只能给出一个“平均”结果看不出频率随时间变化的过程。后来发展的短时傅里叶变换、小波分析都是为了解决这个局限。所以理解傅里叶变换的适用边界和会用公式一样重要。2. 从傅里叶级数到傅里叶变换周期与非周期之间只差一个极限2.1 傅里叶级数把周期信号拆成一组“纯音”学傅里叶变换之前最好先把傅里叶级数吃透因为两者是同一件事的两个阶段。1807年傅里叶提出任何一个周期信号都可以分解成一组正弦波的叠加。一个周期为 (T) 的信号可以写成基频 (f_01/T) 的正弦波加上频率为 (2f_0, 3f_0, 4f_0...) 的谐波每一项乘以各自的系数再全部加起来。关键在于为什么必须是基频的整数倍因为只有整数倍频率的正弦波叠加起来才仍然是以 (T) 为周期的周期信号。非整数倍频率往里一加周期就乱了。方波是理解级数最好的例子。一个方波看起来方方正正但傅里叶说它是无数个正弦波的叠加——基频的幅度最大然后是3倍频、5倍频、7倍频幅度依次减小。你叠加的谐波越多合成波形越接近方波。但有个很有意思的现象叫吉布斯现象无论叠加多少项方波的跳变沿附近始终有过冲不会消失。网上有大量展示这个过程的动画建议找来看看比抄十遍公式管用。2.2 推极限周期变无穷大级数就变成了积分周期信号可以用级数分解那非周期信号怎么办比如你拍一下桌子振动一下就没了它不是周期信号没法用傅里叶级数。傅里叶的思路很直接把非周期信号看成“周期无穷大的周期信号”。周期 (T) 趋近于无穷大时基频 (f_01/T) 趋近于0谐波之间的频率间隔趋近于0离散的一条条谱线最终变成了一条连续的频谱曲线。数学上原来的离散求和就变成了连续积分这就是傅里叶变换[ X(f)\int_{-\infty}^{\infty} x(t)e^{-j2\pi ft}dt ]这个积分不用怕它的本质就是“把信号 (x(t)) 往频率为 (f) 的纯正弦波上投影”积分就是做内积。结果 (X(f)) 是一个复数幅度代表该频率成分的强度相位代表该成分的时间偏移。工程上画频谱图往往只看幅度谱把相位丢了这个习惯很容易埋坑——反变换重建信号时没有相位信息重建出来的东西可能完全不对。所以提醒一句幅度谱和相位谱同样重要。2.3 总被误解的一点频谱到底是离散还是连续这个问题我在工作里问过不少同事能立刻说清楚的不多。记住两条就够周期信号 → 离散谱也就是一根根线谱非周期信号 → 连续谱用乐谱和录音类比最直观。一首歌的乐谱上面只有有限的几个音符这是“周期信号”的离散谱。但真实演奏出来的录音包含了乐器的泛音、空间混响、各种细微的连续频率成分它的频谱是连续一片的。所以如果你处理的是一段真实采集的振动信号它几乎不可能是严格周期的你看到频谱一般是连续谱上面某些位置有明显凸起的峰工程上所谓的“谱峰”就是指这些凸起。3. 计算机只认识离散数据从DFT到FFT的门槛3.1 采样不是把信号抄下来而是按快门计算机没法处理连续信号它只能每隔固定时间取一个数值这个动作叫采样。采样率 (f_s) 就是每秒取多少个点。这里有一个整个信号处理领域最重要的定理——奈奎斯特采样定理要无失真地恢复最高频率为 (f_{max}) 的信号采样率必须大于 (2f_{max})。为什么必须是2倍以上道理特别朴素一个正弦波一个周期内至少要有两个采样点才能确定它到底是高频还是低频。如果每个周期只采一个点你看到的可能是直流采样点再少一点高频信号就会伪装成低频信号这就是“混叠”。很多人在监控视频里看到车轮在倒转也是同一回事——时间上的采样率不够视觉上就产生了假频率。工程上通常不会卡着2倍取而是留出裕量。我做振动测试时一般取关注的最高频率的2.56倍到10倍采样率越高抗混叠越容易做但数据量也越大要平衡。3.2 手算一个DFT5个点就够明白公式了离散傅里叶变换简称DFT公式长这样[ X[k]\sum_{n0}^{N-1}x[n]e^{-j2\pi kn/N} ]为了让你彻底明白它咱们手算一个极简例子。取 (N5)采样序列是 (x[1,2,3,4,5])。当 (k0) 时(e^01)所以 [ X[0]1234515 ] 这个就是直流分量也就是信号的平均值乘以 (N)。当 (k1) 时 [ X[1]12e^{-j2\pi/5}3e^{-j4\pi/5}4e^{-j6\pi/5}5e^{-j8\pi/5} ] 用欧拉公式展开把实部和虚部分别加起来结果是 [ X[1]\approx -2.53.44j ]就这么一个简单计算你就能看出来为了算一个 (X[k])要做 (N) 次复数乘加一共有 (N) 个 (k)所以整个DFT的计算量是 (O(N^2))。(N8192) 时需要6700多万次运算做完一次分析人都要等麻了。而快速傅里叶变换也就是FFT只是把计算复杂度降到 (O(N\log N))它不是一个新变换是同一个变换的快速算法。普通计算机普通的FFT处理8192个点基本上是瞬间完成。这就是为什么实际工程里永远在用FFT而不是直接算DFT。3.3 频率轴怎么画k和实际频率的关系用FFT算完的结果是一串复数每个点对应一个频率但这个频率怎么读公式很简单[ f_k k \cdot \frac{f_s}{N} ]第0个点是直流第 (N/2) 个点对应奈奎斯特频率 (f_s/2)再往后的点对应负频率。画频谱图时通常只看前面一半也就是从0到 (f_s/2)这叫单边谱。还有一个关键参数频率分辨率 (\Delta f f_s/N 1/T)(T) 是信号的实际时长。这意味着你想要1Hz的分辨率至少要采1秒数据想要0.1Hz就得采10秒。时长远一点频率分辨得细一点这是物理规律不是算法能绕过去的。4. 我踩过的坑谱泄漏、加窗与补零4.1 完整排查一个50Hz纯正弦波频谱怎么会“糊”了说一个我自己踩过的大坑。有一回在Python里生成一个标准的50Hz正弦波采样率1000Hz采了1秒然后直接FFT看频谱。理论上这应该是一根完美的谱线结果画出来50Hz附近居然有一大片不太整齐的分量像是一座带“裙边”的山峰。我当时第一反应是FFT代码写错了。换MATLAB算结果一样怀疑采样率设置问题检查没错怀疑点数不够改大FFT点数还是那样。折腾了大半天才终于想到问题出在“截断”上。我们只能对一段有限长度的信号做分析这相当于在时域里给信号乘了一个矩形窗——从0秒到1秒的地方信号突然开始、突然结束。时域相乘对应频域卷积而矩形窗的频谱是sinc函数它有个主瓣旁边还带着一串衰减的旁瓣。于是本来应该是一根线的50Hz就卷积出“糊”成一团的效果。那串“裙边”就是矩形窗的旁瓣。一句话总结谱泄漏不是算法错了是“截断”这个操作本身造成的。只要被分析的信号不是完整周期截断泄漏就一定存在。这段排查经历让我以后再也不敢不看时域波形就直接FFT了。下面是一段能复现这个现象的代码你可以自己跑一下看看import numpy as np import matplotlib.pyplot as plt fs 1000 N 1000 t np.arange(N) / fs x np.sin(2 * np.pi * 50 * t) X np.fft.fft(x) freqs np.fft.fftfreq(N, 1 / fs) plt.figure(figsize(10, 4)) plt.plot(freqs[:N // 2], np.abs(X[:N // 2])) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.title(50Hz sine wave spectrum (rectangular window)) plt.grid(True) plt.show()4.2 加窗不同窗函数就是不同滤镜既然矩形窗的旁瓣太高工程上的标准做法是给信号乘一个“两端趋于0”的窗函数让截断变得柔和一些。时域上把信号边缘慢慢“按下去”频域上旁瓣就会被压低。常用窗函数的特点我用一张表总结一下窗函数主瓣宽度旁瓣衰减适用场景矩形窗窄约-13dB整周期截断、瞬态信号汉宁窗中等约-31dB一般频谱分析首选汉明窗中等约-43dB对频域幅值精度要求较高布莱克曼窗较宽约-58dB分离两个频率相近的信号注意加窗是把双刃剑旁瓣压低了主瓣变宽了原本离得很近的两个频率峰可能就分不开了。不存在“万能窗”。我做振动分析默认先用汉宁窗发现频率分离有问题再换别的。做传声器声学测量时又经常需要根据场景选择其他窗。原则就是先搞清楚你要看的是幅值、频率还是近距频率分离能力再选窗。这里补充一个很多教程不会细讲的坑加窗会让频谱幅值变小因为窗函数本身损失了一部分能量。汉宁窗的相干增益约0.5所以加汉宁窗之后算出来的频谱峰值要除以0.5也就是乘以2再按后续的幅度换算方法处理才能恢复真实幅值。这个细节在做定量分析时特别重要否则你测出来的幅值天生就少一半。4.3 补零不是提高分辨率只是让曲线更圆滑还有一个在论坛上看过无数次的误解——频谱不够精细有人就在信号末尾补几千个零再做FFT发现曲线确实变光滑了于是宣称“分辨率提高了”。这是错的。补零确实会让FFT的点数变大频率轴变密但真实分辨率由信号本身的时长决定公式是 (\Delta f 1/T)。补零不会给信号增加任何新信息它只是在已有的频谱曲线上做插值让曲线看起来更平滑。这就像你把一张低分辨率图片拖大了看颗粒感还在细节不会凭空长出来。真正要提高分辨率只有一条路延长采集时间。比如两个频率成分相差0.5Hz你的数据只有1秒分辨率是1Hz补再多的零、FFT点数再大也永远分不开它们。把这句话刻在脑子里能帮你少走很多弯路。5. 验证你的频谱算对了没有三个自查方法5.1 帕塞瓦尔定理能量守恒核对法做信号处理最怕的就是算了个频谱还不敢确定对不对。我每次换新库、换新代码第一个验证手段是帕塞瓦尔定理。它的意思是信号在时域里的总能量等于在频域里的总能量。具体到离散形式[ \sum_{n0}^{N-1}|x[n]|^2 \frac{1}{N}\sum_{k0}^{N-1}|X[k]|^2 ]把原始信号的平方和算出来再把FFT结果的模平方和除以 (N) 算出来两边应该几乎相等。浮点误差范围内比如1e-10量级正常如果对不上或者差很多倍问题一定出在FFT的归一化方式上。不同软件对FFT的默认缩放处理不一样有的库不除有的库除 (N)这个定理能一秒帮你定位。5.2 幅度换算正弦波峰值到底应该是多少频谱图上的峰高度到底代表多少幅值这个坑尤其隐蔽。以标准正弦波 (x(t)A\cos(2\pi f_0 t)) 为例。做 (N) 点FFT如果直接把 (|X[k]|) 除以 (N)得到的是双边谱在 (f_0) 处的峰值是 (A/2)。如果要画单边谱把正频率部分的幅度乘以2直流除外峰值就是 (A)。这就是“单边谱”和“双边谱”的差别。但要注意这是在没有加窗、信号频率正好落在谱线上时才成立。如果加了汉宁窗还得先把相干增益的系数校正回来也就是前面说的乘以2。不同软件画频谱图时默认处理完全不同有的自动做了单边化有的没有。所以最稳的做法是造一个幅值已知的正弦波走一遍自己的完整流程确认输出幅值和你期望的一致再去处理真实数据。我每到一个新环境、新工具链都会先做这个“标定”习惯了就再也没被幅值问题坑过。5.3 反变换重建头尾相接的终极校验第三种校验方法比较暴力但最有效对变换结果做逆FFT把时域信号重建出来再和原始信号对比。如果最大误差在小数点后10位量级说明你这套变换基本没问题。如果对不上优先查两件事相位信息是不是被丢了。很多人做分析时只保留了幅度谱拿去反变换结果当然是错的。处理完后要用复数谱做反变换不要用幅度谱。共轭对称性是不是被破坏了。实信号做FFT后结果有共轭对称性也就是 (X[N-k]\overline{X[k]})。如果你对频谱做了裁剪、只保留前半部分或者修改了某个频点而没做对称处理反变换出来的就是复数甚至是一团乱码。这个方法对写自定义FFT封装函数的人尤其重要可以拦截掉大量隐性bug。6. 一个真实项目用傅里叶变换定位减速机异响6.1 故障现象与现场情况最后讲一个我印象最深的项目。某条产线的一台减速机最近出现了周期性的“咔嗒”声位置大概在输出轴附近。老师傅听了几次凭经验判断可能是齿轮磨损建议直接拆开检查。但拆减速机不是小事停产、起吊、拆装、回装至少要一整天所有计划都会被打乱。我提出先做一次振动检测用傅里叶变换看看频谱里有什么。现场情况在减速机轴承座上吸了一个加速度传感器采样率设为10kHz。为什么选10kHz因为轴承故障特征频率一般落在几百赫兹到几千赫兹10kHz足够覆盖又不会产生过大的数据量。采集时长2秒因为我们需要至少1Hz的分辨率2秒数据真实分辨率能做到0.5Hz还可以做多次平均来稳定频谱。6.2 频谱分析过程原始时域波形重叠在一起毫无规律肉眼根本看不出异常。对数据做FFT(N) 取8192采样率10k加汉宁窗然后绘制频谱。频谱一出来情况就明朗了。首先看到明显的转频及其谐波——24.6Hz、49.2Hz……对应电机输出转速约1476转/分。在约116Hz附近有一组不太起眼的峰而且旁边带着一系列间隔约24Hz的边频带。这里有个非常重要的谱图判读知识边频带意味着调制。某个高频振动被一个低频特征频率调制表现形式就是频谱上出现一个中心峰两侧对称分布着一串间隔为调制频率的边带。间隔约24Hz和转频吻合说明这个振动和转频之间存在调制关系。按滚动轴承外圈故障特征频率估算公式计算如果用常见参数估算外圈故障特征频率大约是转频的4.73倍也就是[ 4.73 \times 24.6Hz \approx 116.4Hz ]这个116Hz的峰和外圈故障特征频率对得上。于是初步判断不是齿轮问题而是轴承外圈故障。6.3 拆机验证与复盘后来检修窗口期拆开减速机轴承外圈滚道果然有一处明显的疲劳剥落。诊断结论和实际故障完全吻合。这个项目让我想明白一个道理傅里叶变换在工程里的价值不是让你计算“数学题”而是把“凭感觉判断”变成“可量化的诊断依据”。没有频谱分析只能整机拆解费时费力有了频谱分析停在现场半小时就能给出一个高置信度的结论。但也要冷静地说一句傅里叶变换给出的是“特征”不是“结论”。它告诉你这里有116Hz的边频带但116Hz这个频率对应轴承外圈还是滚动体还需要结合机械设计参数去计算。工具用得再好最终还是要落到你对系统的理解上。我后来带新人学傅里叶变换都会让他们先做三件事第一造一个已知信号自己走一遍采样、FFT、幅值换算的完整流程第二故意用矩形窗对非整周期截断感受一次谱泄漏长什么样第三用帕塞瓦尔定理去核对自己的结果搞清楚“对不上”到底差在哪里。这三件事做完公式背后的逻辑就慢慢长在身上了。傅里叶变换不是靠背出来的是靠一次次调试、一个个坑踩出来的。