ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

傅里叶变换实战指南:从频谱分析到Python信号与图像处理

傅里叶变换实战指南:从频谱分析到Python信号与图像处理 先别急着翻书我给你讲个更直白的场景你在录音棚里录了一段吉他结果混进了空调的低频轰鸣和隔壁房间的狗叫。时域里这段波形乱成一团你想把狗叫抠掉鼠标拖半天也选不准位置。但如果你把这段声音扔进傅里叶变换里事情立刻变清楚——低频喘气声待在左边狗叫的高频刺耳声待在右边你只需要在频域里把对应频段“静音”一下再变回时域干净的吉他声就回来了。这就是傅里叶变换最典型的应用逻辑不在原来的坐标系里死磕换一个视角把问题变简单。这几年我不管是做音频处理、振动分析还是机器视觉里的图像滤波绕来绕去最后都在跟傅里叶变换打交道。今天的系列第一篇就把信号处理、图像处理这两个最常碰到的应用场景摊开讲配合Python代码直接演示看完你就能在自己项目里上手。1. 先从直觉说起傅里叶变换到底在干什么1.1 一个生活化的类比把波形拆成“配料表”很多教材一上来就丢公式把人劝退。我换个说法任何一个复杂的波形都可以看成是一堆频率不同、幅度不同、相位不同的正弦波叠加出来的。傅里叶变换干的事情就是把这堆混合在一起的配料逐个挑出来列一张配料清单哪个频率分量有多少克幅度、什么时候下锅相位。你每天早上喝的那杯拿铁咖啡豆、牛奶、糖浆混在一起你喝到的是一杯整体。但如果你想搞清楚这杯咖啡到底放了多少糖就得想办法把成分分离出来。傅里叶变换对信号干的正是这件事。时域信号是“咖啡”频域谱就是“配料表”。这个视角一换很多问题就豁然开朗了。时域里你看到的是纠缠在一起的波形频域里你能清清楚楚看到每一个周期成分的贡献大小。比如一段50Hz工频干扰叠加在有用信号上时域里波形只是多了些毛刺频域里50Hz处会竖起一根明显的尖峰用频域处理就是精准拔掉这根刺。1.2 为什么工程师离不开频域视角我最早接触傅里叶变换是在大学信号与系统课上当时只觉得是数学折磨。真正让我改观的是工作后做设备振动诊断电机轴承磨损时时域波形几乎看不出异常但频谱里会出现特定频率的边带族。那一刻我才意识到频域不是数学家的玩具而是工程诊断的通用语言。在音频领域均衡器就是调整不同频段的增益在通信领域OFDM把高速数据流拆到多个正交子载波上传输在图像处理领域JPEG压缩的核心就是丢弃高频细节。这些看似完全不同的技术底层全是傅里叶变换的变体。你不需要一次性把所有分支都吃透但抓住“时域和频域是同一事物的两个视角”“卷积在时域复杂、在频域是乘法”这两条主线后面学任何应用都会顺畅很多。2. 公式与常用变换对看懂它们才敢用2.1 连续公式与离散公式先分清你处理的是什么傅里叶变换的连续形式长这样$$X(f) \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft} dt$$看着吓人但拆开就一句话把信号x(t)和不同频率的复指数e^{-j2πft}做内积看每个频率分量有多大。逆变换则是把这些分量加回去重建原始信号。但在计算机里你处理的永远是离散采样序列所以真正用的是离散傅里叶变换DFT$$X[k] \sum_{n0}^{N-1} x[n] e^{-j2\pi kn/N}$$DFT把N个时域采样点变成N个频域采样点第k个点对应频率k×fs/N其中fs是采样率。这里必须注意DFT计算量是O(N^2)几千个点还凑合几百万个点就废了。所以实际工程里99%用的是快速傅里叶变换FFT它是DFT的高效算法将复杂度降到O(N log N)。Python里直接用numpy.fft.fft就能调用内部已经做了优化。2.2 常用变换对速查我实际项目里翻来覆去就那几个表和解题不一样工程上常用的变换对其实非常固定。我把这些年用得最多的列成一张表贴在工位上那种时域信号傅里叶变换频域典型场景直流常数冲激函数位于0Hz去直流分量时定位基线单频正弦对称的两根谱线识别工频干扰、振动基频矩形脉冲sinc函数sa函数理解窗函数频谱泄漏高斯函数高斯函数图像高斯滤波的理论基础冲激串冲激串采样定理、频谱周期性单位阶跃1/jπf加冲激分析阶跃响应的频域特征比如你处理音频时发现频谱在0Hz处有个巨大的尖峰那是有直流偏置直接在频域把这个点置零就能去掉。又比如图像滤波时用高斯核正是因为高斯函数的傅里叶变换还是高斯空域平滑对应频域低通不会引入振铃。这些变换对不用死记用多了自然熟。但有一个必须印在脑子里矩形窗对应sinc函数而sinc函数有正负旁瓣——这就是频谱泄漏和图像振铃的理论根源。后面排查问题时会反复提到。3. 实例一Python音频信号去噪实战3.1 构造一个含噪信号先学会看频谱演示之前先解释为什么选音频场景。音频信号是典型的一维信号采样率固定、数据量大、降噪需求明确是理解FFT滤波最好的练兵场。我建议你自己动手生成测试数据而不是直接录一段真实音频因为合成信号你知道“正确答案”方便验证。下面这段代码生成一个5Hz工频干扰叠加两个低频成分的信号再加入随机噪声import numpy as np fs 1000 # 采样率 1000Hz t np.arange(0, 1, 1/fs) # 1秒时间轴 # 有用信号5Hz 50Hz 两个正弦叠加 x 1.5*np.sin(2*np.pi*5*t) 0.8*np.sin(2*np.pi*50*t) # 加入随机噪声 noise 0.6*np.random.randn(len(t)) signal x noise你画出时域波形会发现5Hz和50Hz的正弦成分被噪声彻底淹没肉眼根本分辨不出来。这时候做FFT一切就会明朗X np.fft.fft(signal) # 快速傅里叶变换 freqs np.fft.fftfreq(len(X), 1/fs) # 频率轴 # 只取正频率部分画幅度谱 half len(X)//2 plt.plot(freqs[:half], np.abs(X[:half]))频谱图上你会在5Hz和50Hz处看到两根清晰的尖峰其余全是噪声底。这就是时域和频域视角的差异。顺便强调一个惯用套路fft结果零频在第一个元素用fftfreq生成频率轴时前半是正频率、后半是负频率工程上常常用fftshift把零频移到中心看起来更直观。用不用shift不影响数据处理只影响可视化。3.2 完整的FFT滤波流程频率置零不是唯一办法最朴素的做法是找到噪声所在频点直接置零然后逆变换X_filter X.copy() X_filter[np.abs(freqs) 1] 0 # 去掉直流偏置 X_filter[(np.abs(freqs) 8) (np.abs(freqs) 40)] 0 # 去掉低频段噪声 X_filter[np.abs(freqs) 200] 0 # 去掉超高频噪声 y np.real(np.fft.ifft(X_filter))这里有几个关键点。第一正频率和负频率是对称的处理时务必把对称位置一起置零否则你逆变换回来的信号会有虚部相位也全乱套。第二硬性置零等于在频域乘了一个矩形窗时域会引入振铃效应Gibbs现象表现为信号在跳变处出现振荡尾巴。第三更好的方案是用平滑过渡的滤波器比如用scipy.signal里的butterworth滤波器设计巴特沃斯带通在频域乘一个平滑的传递函数能有效减轻振铃。我实际团队里的做法是能用现成滤波器尽量别手动抠频点只有做这种教学演示或频谱分析时才直接用置零法展示原理。抠频点一步到位但代价是信号质量下降。3.3 参数选择的门道采样率、窗函数、频谱泄漏这一步最容易翻车。很多人以为只要FFT完事就结束了实际上FFT的效果依赖几个前置参数。首先是采样率。根据奈奎斯特采样定理采样率至少是信号最高频率的2倍否则高频成分会混叠到低频段。工程上我习惯至少留出2.5到3倍余量后面加抗混叠滤波器也有空间。其次是分析时长。FFT的频率分辨率是fs/NN是参与FFT的点数。想分辨间隔1Hz的两个频率峰你的采样时长至少需要1秒。靠补零不能提高分辨率只是把频谱插值得更平滑这一点很多教程没讲透。再次是窗函数。直接对一段截断信号做FFT相当于给信号乘了矩形窗会产生频谱泄漏让单频成分的能量散到旁边频点去。解决办法是加汉宁窗或汉明窗window np.hanning(len(signal)) X_windowed np.fft.fft(signal * window)用窗函数测频率泄漏更小测幅度更准。代价是主瓣变宽两个很近的频率峰更难区分。所以做纯幅度测量我用平顶窗做窄带信号频率识别我用汉宁窗这属于工程经验书本上不会教这么细。4. 实例二机器视觉里的图像傅里叶变换4.1 二维FFT与频谱图怎么看图像的傅里叶变换是把二维DFT分别沿行和列各做一次变换公式记不住没关系关键是要理解图像频谱图的读法。用numpy实现非常直接import cv2 import numpy as np img cv2.imread(texture.jpg, cv2.IMREAD_GRAYSCALE) f np.fft.fft2(img.astype(np.float32)) f_shift np.fft.fftshift(f) # 把零频挪到中心 magnitude np.log(np.abs(f_shift) 1e-6) # 取对数放大动态范围读频谱图的几个核心看板中心亮点是直流分量对应图像的平均灰度越亮说明整体越亮。中心向外辐射的亮点代表不同频率的成分离中心越远频率越高对应图像中越细微的结构边缘、纹理、噪点。频谱的方向对应图像中纹理的方向如果图片里有大量横向条纹频谱会在横轴上出现一条明亮的线。我做过织物瑕疵检测布纹的周期性纹理在频谱里看非常明显有瑕疵时对应频点就会出现异常的峰值。4.2 高通、低通、带通滤波器实操频域滤波的流程是FFT - 设计频域掩膜 - 相乘 - 逆FFT。掩膜只能作用于幅度谱还是可以作用于复数频谱严格说掩膜可以直接乘在复频谱上但为了不破坏相位往往只处理幅度。图像中相位承载了大量结构信息乱动相位会得到面目全非的图像。低通滤波保留中心低频滤掉外围高频效果是图像变模糊去噪点rows, cols img.shape crow, ccol rows//2, cols//2 # 构造低通掩膜半径r内的中心区域为1其余为0 mask_low np.zeros((rows, cols), np.float32) cv2.circle(mask_low, (ccol, crow), 30, 1, -1) f_shift_filtered f_shift * mask_low f_filtered np.fft.ifftshift(f_shift_filtered) img_filtered np.real(np.fft.ifft2(f_filtered))高通滤波反过来滤掉低频保留高频图像会变成只剩边缘和细节的锐化效果常用来做边缘提取预处理。带通滤波则是保留某个频带比如提取特定尺度的纹理在检测布匹纹理缺陷时非常有用。需要特别提醒二维FFT后频谱角落的四个角是最高频区域不是噪声区。很多人第一次看到频谱图把四个角的亮点误当成噪声去滤结果图像整体模糊。高频成分对应锐利边缘滤过头细节全丢。4.3 频域做模板匹配与图像配准的另一种思路频域在机器视觉里不只是做滤波。我举一个实际经历给两张有重叠区域的航拍图像做拼接传统做法是找特征点匹配但如果图像纹理稀疏特征点经常找不到。这时候用傅里叶相位相关法就有效果。相位相关的核心是傅里叶变换的平移性质空域的平移对应频域的线性相位移。两张图像的相对位移可以通过互功率谱的逆变换峰值来确定。用Python做大概三十行def phase_correlation(img1, img2): f1 np.fft.fft2(img1) f2 np.fft.fft2(img2) cross f1 * np.conj(f2) normalized cross / (np.abs(cross) 1e-6) result np.fft.ifft2(normalized) return np.unravel_index(np.argmax(np.abs(result)), img1.shape)返回的坐标就是两幅图像间的位移量。这里面有个陷阱计算结果受图像边缘影响大如果两幅图光照不一致或内容不完全重叠峰值会不明显得先用高通或边缘增强预处理突出结构信息。图像配准之后还能通过FFT做旋转角度的估计原理涉及傅里叶变换的旋转性质这块本篇先不展开后续系列里我会单独讲。5. 常见问题与排查技巧实录5.1 频谱图为什么左右对称实数信号的DFT结果具有共轭对称性第k个频点和第N-k个频点的幅度相等、相位相反。这导致你看到的频谱图总是左右对称。处理时只关心正频率部分即可做谱分析画图时取前一半。做滤波时手动置零记得正负频率一起处理否则逆变换信号是复数相位也出错。还有一个衍生问题如果你把频谱左侧当“低频”右侧当“高频”那方向就踏错了。正确理解是频率轴从负到正排列对称点是零频。5.2 频谱泄漏永远消不掉吗彻底消除频谱泄漏不可能因为处理信号必然截断截断就等效于乘窗而窗函数的频谱不可能无限窄。只能缓解。缓解手段有选择旁瓣更低的窗函数比如布莱克曼窗或凯泽窗保证采样时长尽量是信号周期的整数倍对周期性信号有用但真实信号很难满足用更长的时间窗口提高频率分辨率让泄漏能量集中到更窄的频带里。我实际处理振动数据时常常先用汉宁窗加多次平均叠加平均法把随机噪声打下去周期性振动信号就能看得更清楚。这里没有银弹只能适配场景。5.3 图像频谱图中心的“亮十字”是啥很多人在图像频谱图里看到穿过中心的横竖亮线以为是信号有问题。其实这是图像边缘不连续造成的。FFT默认图像在边界处首尾相接如果图像左右边缘灰度不连续频谱横轴方向就会出现一条亮线竖线同理。消除办法是给图像乘以一个二维窗函数让边缘慢慢衰减到零。或者做边缘零填充后再FFT处理完再裁回来。工业检测项目里我经常用边缘填充来压制这种伪影效果立竿见影。还有一种情况是亮十字特别集中说明原图中有固定频率的条纹噪声这往往来自传感器行扫描噪声或电源纹波可以用陷波滤波器针对性地扣掉对应频率的尖峰。5.4 时域卷积对应频域乘法这个特性怎么用卷积定理大概是傅里叶变换里最有工程价值的性质两个信号在时域做卷积等于各自傅里叶变换后在频域相乘再逆变换回去。这意味着卷积运算的复杂度可以从O(N^2)降到O(N log N)。我第一次实际用这个性质是在做大尺寸图像的模板匹配。直接做卷积模板在图像上滑窗计算2560x1440的图加一个100x100的模板就要跑十几秒。而利用卷积定理用FFT计算频域乘积同样尺寸只需几十毫秒。代码模式固定f_img np.fft.fft2(img) f_tmpl np.fft.fft2(template, img.shape) # 模板补零到图像尺寸 result np.real(np.fft.ifft2(f_img * np.conj(f_tmpl)))这其实就是相位相关法的“亲戚”本质上都是频域匹配。但要注意大模板用FFT才划算模板尺寸小的时候直接滑窗反而更快因为FFT有固定开销。判断阈值我一般按模板像素数占图像像素数的5%来划边界小于这个比例直接滑动窗口大于就上FFT。6. 一些实操体会做了这么多年信号和图像处理我最大的感受是傅里叶变换的价值不在公式推导在于建立“换个坐标系看问题”的思维习惯。很多在时域里剪不断理还乱的现象拖到频域里就是一根刺或一坨噪声思路立刻简单了。如果你刚入门给你三条建议。第一先在自己熟悉的领域找一个具体问题比如音频降噪或者图像去模糊拿着代码跑通一遍比读十篇理论文章都有用。第二一定要把fftshift、fftfreq这种API的手感和边界条件弄清楚这决定了你能不能把数学作用到真实数据上。第三永远先画频谱图再说话。数据到手第一步不是写滤波逻辑而是看频谱长什么样很多时候问题自己就暴露出来了。这一篇里的一维FFT和二维FFT只是第一梯队工具后面还有短时傅里叶变换、小波变换、相位相关匹配这些延伸后续继续更新。你要是有具体应用场景带着问题来留言交流效果最好。
返回列表