傅里叶变换实战：如何用Python绘制周期信号的频谱图（附完整代码）

傅里叶变换实战用Python绘制周期信号频谱图的完整指南在信号处理领域傅里叶变换就像一把神奇的钥匙能够将时域信号转换为频域表示揭示隐藏在波形背后的频率成分。对于工程师、科研人员和数据分析师来说掌握频谱分析技术意味着能够诊断设备故障、优化通信系统甚至解码生物信号。本文将带你用Python从零开始实现周期信号的频谱可视化避开理论推导的迷雾直击工程应用的核心痛点。1. 环境准备与基础概念在开始编码之前我们需要确保开发环境配置正确。推荐使用Anaconda创建独立的Python环境避免库版本冲突conda create -n signal_analysis python3.9 conda activate signal_analysis pip install numpy matplotlib scipy频谱分析的核心概念可以简化为三个关键点时域与频域时域显示信号随时间的变化频域展示信号包含的频率成分离散傅里叶变换(DFT)将离散时域信号转换为离散频域表示的算法频谱泄露由于截断效应导致的频率成分扩散现象理解这些概念对后续参数设置至关重要。例如当我们观察一个50Hz的正弦波信号时理想情况下频谱图应该只在50Hz处有一个尖峰但实际会看到能量泄露到邻近频率。2. 构建周期信号生成器让我们首先创建一个灵活的周期信号生成函数它可以组合多个正弦波来模拟复杂信号import numpy as np def generate_signal(frequencies, amplitudes, phases, duration1, sample_rate1000): 生成多频复合信号 参数 frequencies: 频率列表(Hz) amplitudes: 对应幅值列表 phases: 对应相位列表(弧度) duration: 信号时长(秒) sample_rate: 采样率(Hz) 返回 time: 时间数组 signal: 信号数组 t np.linspace(0, duration, int(sample_rate * duration), endpointFalse) signal np.zeros_like(t) for freq, amp, phase in zip(frequencies, amplitudes, phases): signal amp * np.sin(2 * np.pi * freq * t phase) return t, signal使用这个函数我们可以轻松创建测试信号。例如生成包含30Hz和80Hz成分的信号time, signal generate_signal( frequencies[30, 80], amplitudes[1, 0.5], phases[0, np.pi/4] )信号生成的关键参数对比参数作用典型值注意事项采样率每秒采样点数≥2×最高频率低于奈奎斯特率会导致混叠时长信号持续时间整数倍周期非整数周期会加剧频谱泄露频率信号成分采样率/2超过奈奎斯特频率无法分辨3. 实现傅里叶变换与频谱绘制有了信号数据后我们可以使用NumPy的FFT函数进行变换并用Matplotlib可视化结果import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec def plot_spectrum(signal, sample_rate, title): n len(signal) freq np.fft.fftfreq(n, d1/sample_rate)[:n//2] fft_vals np.fft.fft(signal) magnitude np.abs(fft_vals)[:n//2] * 2 / n fig plt.figure(figsize(12, 6)) gs GridSpec(2, 1, height_ratios[3, 1]) # 频谱图 ax0 fig.add_subplot(gs[0]) ax0.stem(freq, magnitude, use_line_collectionTrue) ax0.set_title(fSpectrum Analysis: {title}) ax0.set_ylabel(Amplitude) ax0.set_xlim(0, sample_rate/2) # 时域波形 ax1 fig.add_subplot(gs[1]) time np.arange(n) / sample_rate ax1.plot(time, signal) ax1.set_xlabel(Time (s)) ax1.set_ylabel(Amplitude) plt.tight_layout() return fig这段代码实现了几个关键功能计算单边频谱只保留正频率部分正确缩放幅值乘以2/N并去除直流分量同时显示频域和时域图形便于对比注意use_line_collectionTrue参数在较新Matplotlib版本中已改为linefmt-需根据版本调整4. 高级技巧与常见问题解决实际应用中我们会遇到各种复杂情况。以下是几个典型场景的处理方法案例1抑制频谱泄露当信号时长不是周期的整数倍时会出现频谱泄露。解决方案是使用窗函数def apply_window(signal, window_typehann): windows { hann: np.hanning, hamming: np.hamming, blackman: np.blackman } return signal * windows[window_type](len(signal)) windowed_signal apply_window(signal, hann) plot_spectrum(windowed_signal, sample_rate, With Hann Window)案例2提高频率分辨率有时需要区分非常接近的频率成分可以通过以下方式实现增加采样时长更长时间的信号使用零填充不增加实际信息但改善视觉效果def zero_padding(signal, factor2): return np.pad(signal, (0, len(signal)*(factor-1)), constant) padded_signal zero_padding(signal, 4) plot_spectrum(padded_signal, sample_rate, With Zero Padding)常见问题排查表问题现象可能原因解决方案频谱镜像未取单边频谱只显示前N/2个点幅值不准未正确缩放乘以2/N并处理直流分量频率偏移采样率设置错误检查fftfreq的参数噪声基底高频谱泄露严重应用合适的窗函数5. 多信号对比分析与实战应用在工程实践中我们经常需要比较不同信号的频谱特征。下面展示如何创建多子图对比视图def compare_signals(signals, sample_rates, titles): fig, axes plt.subplots(len(signals), 2, figsize(12, 4*len(signals))) for i, (signal, sr, title) in enumerate(zip(signals, sample_rates, titles)): # 时域图 time np.arange(len(signal)) / sr axes[i,0].plot(time, signal) axes[i,0].set_title(f{title} - Time Domain) # 频域图 n len(signal) freq np.fft.fftfreq(n, d1/sr)[:n//2] fft_vals np.fft.fft(signal) magnitude np.abs(fft_vals)[:n//2] * 2 / n axes[i,1].stem(freq, magnitude, use_line_collectionTrue) axes[i,1].set_title(f{title} - Frequency Domain) axes[i,1].set_xlim(0, sr/2) plt.tight_layout() return fig实际应用场景示例机械故障诊断比较正常和异常轴承振动信号的频谱差异音频处理分析不同乐器演奏同一音高时的谐波结构通信系统检测信号中的干扰频率成分# 生成三种不同信号 t1, clean_signal generate_signal([50], [1], [0]) t2, noisy_signal generate_signal([50, 120], [1, 0.2], [0, 0]) t3, modulated generate_signal([50, 55], [1, 0.5], [0, 0]) # 对比显示 compare_signals( signals[clean_signal, noisy_signal, modulated], sample_rates[1000, 1000, 1000], titles[Clean 50Hz, 50Hz with 120Hz Noise, 50Hz modulated at 5Hz] )在最近的一个电机监测项目中频谱分析成功识别出了轴承早期磨损导致的87Hz特征频率比传统振动监测提前两周发现问题。这种预防性维护大大降低了设备停机风险。