Python scipy.signal 频谱图实战3种振幅谱线性/对数/功率生成与对比信号处理工程师常常需要在时域和频域之间切换视角而频谱图正是连接这两个世界的桥梁。当你第一次看到杂乱无章的时域信号在频域中展现出清晰的频率成分时那种豁然开朗的感觉令人难忘。本文将带你用Python的scipy.signal库亲手实现三种最具代表性的频谱图——线性振幅谱、对数振幅谱和自功率谱并通过实际代码演示它们的差异与应用场景。1. 信号生成与基础准备在开始频谱分析之前我们需要先创建一个包含多个频率成分的测试信号。这个信号将作为我们后续分析的基础。import numpy as np from scipy import signal import matplotlib.pyplot as plt # 设置采样参数 sample_rate 44100 # 采样率44.1kHz接近CD音质 duration 2.0 # 信号持续时间2秒 n_samples int(sample_rate * duration) t np.linspace(0, duration, n_samples, endpointFalse) # 生成包含三个频率成分的测试信号 freq1 100 # 100Hz低频 freq2 1000 # 1kHz中频 freq3 5000 # 5kHz高频 # 各频率成分的振幅 amp1 0.8 amp2 0.3 amp3 0.1 # 生成正弦波信号 signal_clean amp1 * np.sin(2 * np.pi * freq1 * t) \ amp2 * np.sin(2 * np.pi * freq2 * t) \ amp3 * np.sin(2 * np.pi * freq3 * t) # 添加高斯白噪声 noise 0.05 * np.random.normal(0, 1, n_samples) signal_noisy signal_clean noise这个测试信号由三个正弦波组成频率分别为100Hz、1kHz和5kHz振幅依次递减。我们还特意添加了少量高斯白噪声来模拟真实环境中的信号。这样的设计让我们能够清晰观察频谱图对不同频率和振幅成分的表现。信号可视化是理解的第一步plt.figure(figsize(12, 6)) plt.plot(t[:1000], signal_noisy[:1000]) # 只显示前1000个采样点 plt.title(时域信号前1000个采样点) plt.xlabel(时间 (秒)) plt.ylabel(振幅) plt.grid() plt.show()这段代码会显示信号起始部分的时域波形。在时域中我们只能看到信号振幅随时间的变化而无法直接识别出其中包含的频率成分。这正是我们需要频域分析的原因。2. 快速傅里叶变换(FFT)基础傅里叶变换是时域和频域之间的数学桥梁。在计算机上我们使用快速傅里叶变换(FFT)算法高效地实现这一转换。from scipy.fft import fft, fftfreq # 计算FFT n_fft 4096 # FFT点数 yf fft(signal_noisy[:n_fft]) xf fftfreq(n_fft, 1/sample_rate)[:n_fft//2] # 计算幅度谱 magnitude_spectrum 2/n_fft * np.abs(yf[:n_fft//2])这里有几个关键点需要注意n_fft决定了频率分辨率值越大分辨率越高但计算量也越大FFT结果是对称的我们只需要前半部分幅度需要乘以2/N进行归一化直流分量除外FFT结果可视化plt.figure(figsize(12, 6)) plt.plot(xf, magnitude_spectrum) plt.title(基本FFT幅度谱) plt.xlabel(频率 (Hz)) plt.ylabel(幅度) plt.xlim(0, 6000) # 限制显示范围到6kHz plt.grid() plt.show()这个基本的FFT幅度谱已经能够显示出信号的三个主要频率成分但它的实用价值有限。接下来我们将使用scipy.signal中更专业的频谱分析工具。3. 线性振幅谱分析线性振幅谱是最直观的频谱表示方式它直接显示信号中各频率成分的振幅大小。# 使用scipy.signal.spectrogram计算线性振幅谱 f, t, Sxx signal.spectrogram(signal_noisy, fssample_rate, windowhann, nperseg1024, noverlap512, modemagnitude) # 绘制线性振幅谱 plt.figure(figsize(12, 6)) plt.pcolormesh(t, f, Sxx, shadinggouraud) plt.title(线性振幅谱图) plt.ylabel(频率 [Hz]) plt.xlabel(时间 [秒]) plt.colorbar(label幅度) plt.ylim(0, 6000) # 限制频率显示范围 plt.show()关键参数说明nperseg每个段的长度决定频率分辨率noverlap段之间的重叠点数影响时间分辨率window窗函数类型这里使用汉宁窗减少频谱泄漏线性振幅谱特别适合观察信号中主要频率成分的相对强度。在我们的测试信号中可以清晰看到100Hz、1kHz和5kHz三个成分它们的振幅差异与我们在信号生成时设置的完全一致。提示当需要精确测量信号中各频率成分的绝对振幅时线性振幅谱是最佳选择。但在分析包含大动态范围即同时存在很强和很弱的信号成分的信号时它可能无法清晰显示弱信号。4. 对数振幅谱分析对数振幅谱通过对幅度取对数通常以10为底单位为分贝dB来压缩动态范围使得弱信号更容易被观察到。# 计算对数振幅谱dB f_log, t_log, Sxx_log signal.spectrogram(signal_noisy, fssample_rate, windowhann, nperseg1024, noverlap512, modemagnitude) Sxx_db 20 * np.log10(Sxx_log) # 转换为分贝 # 绘制对数振幅谱 plt.figure(figsize(12, 6)) plt.pcolormesh(t_log, f_log, Sxx_db, shadinggouraud) plt.title(对数振幅谱图 (dB)) plt.ylabel(频率 [Hz]) plt.xlabel(时间 [秒]) plt.colorbar(label幅度 (dB)) plt.ylim(0, 6000) plt.show()对数变换的公式为dB 20 * log10(amplitude)。这种表示方式有几个显著优势能够同时显示大信号和小信号动态范围更广更符合人类听觉的感知特性人耳对声音强度的感知近似对数关系在音频处理等领域是标准表示方法在我们的测试信号中虽然5kHz成分的绝对振幅只有100Hz成分的1/8但在对数谱中仍然清晰可见。这种特性使得对数振幅谱成为分析包含噪声或弱信号成分的理想选择。5. 自功率谱分析自功率谱简称功率谱显示的是信号中各频率成分的功率能量分布它是幅度谱的平方。# 计算功率谱 f_power, t_power, Sxx_power signal.spectrogram(signal_noisy, fssample_rate, windowhann, nperseg1024, noverlap512, modepsd) # 功率谱密度模式 # 绘制功率谱 plt.figure(figsize(12, 6)) plt.pcolormesh(t_power, f_power, 10 * np.log10(Sxx_power), shadinggouraud) plt.title(功率谱图 (dB/Hz)) plt.ylabel(频率 [Hz]) plt.xlabel(时间 [秒]) plt.colorbar(label功率谱密度 (dB/Hz)) plt.ylim(0, 6000) plt.show()功率谱的特点和应用表示信号能量在频域的分布对随机信号分析特别有用如噪声在通信系统中用于分析信号带宽单位通常是dB/Hz功率谱的一个关键特性是它丢弃了相位信息只保留幅度信息。这意味着不同信号可能具有相同的功率谱。在我们的测试信号中功率谱清晰地显示了三个频率成分的能量分布。6. 三种频谱图的对比分析现在我们已经实现了三种频谱图让我们通过表格对比它们的主要特点频谱类型数学表示单位优点缺点典型应用场景线性振幅谱V或任意单位直观保持绝对振幅关系动态范围有限精确振幅测量谐波分析对数振幅谱20log₁₀dB大动态范围显示弱信号非线性刻度音频分析噪声检测自功率谱dB/Hz显示能量分布适合随机信号丢失相位信息噪声分析带宽测量为了更直观地比较让我们将三种频谱图并排显示# 计算三种频谱图 f, t, Sxx_lin signal.spectrogram(signal_noisy, fssample_rate, nperseg1024, modemagnitude) _, _, Sxx_log signal.spectrogram(signal_noisy, fssample_rate, nperseg1024, modemagnitude) Sxx_db 20 * np.log10(Sxx_log) _, _, Sxx_power signal.spectrogram(signal_noisy, fssample_rate, nperseg1024, modepsd) # 创建对比图 fig, (ax1, ax2, ax3) plt.subplots(1, 3, figsize(18, 6)) # 线性振幅谱 im1 ax1.pcolormesh(t, f, Sxx_lin, shadinggouraud) ax1.set_title(线性振幅谱) ax1.set_ylabel(频率 [Hz]) ax1.set_xlabel(时间 [秒]) fig.colorbar(im1, axax1, label幅度) # 对数振幅谱 im2 ax2.pcolormesh(t, f, Sxx_db, shadinggouraud) ax2.set_title(对数振幅谱 (dB)) ax2.set_xlabel(时间 [秒]) fig.colorbar(im2, axax2, label幅度 (dB)) # 功率谱 im3 ax3.pcolormesh(t, f, 10 * np.log10(Sxx_power), shadinggouraud) ax3.set_title(功率谱 (dB/Hz)) ax3.set_xlabel(时间 [秒]) fig.colorbar(im3, axax3, label功率谱密度 (dB/Hz)) for ax in [ax1, ax2, ax3]: ax.set_ylim(0, 6000) plt.tight_layout() plt.show()从对比图中可以明显看出线性振幅谱直接反映了各频率成分的振幅比例关系对数振幅谱增强了小信号的可见性功率谱提供了能量分布的视角特别适合分析噪声成分7. 高级应用与参数调优在实际应用中频谱分析的参数选择会显著影响结果。让我们探讨几个关键参数的优化策略。窗函数选择不同的窗函数在频率分辨率和频谱泄漏之间有不同的权衡。windows [boxcar, hann, hamming, blackman] plt.figure(figsize(12, 8)) for i, window in enumerate(windows, 1): f, t, Sxx signal.spectrogram(signal_noisy, fssample_rate, windowwindow, nperseg1024) Sxx_db 20 * np.log10(Sxx) plt.subplot(2, 2, i) plt.pcolormesh(t, f, Sxx_db, shadinggouraud) plt.title(f窗函数: {window}) plt.ylabel(频率 [Hz]) plt.xlabel(时间 [秒]) plt.ylim(0, 6000) plt.tight_layout() plt.show()段长度与重叠nperseg和noverlap参数影响时间和频率分辨率。segment_lengths [256, 512, 1024, 2048] plt.figure(figsize(12, 8)) for i, nseg in enumerate(segment_lengths, 1): f, t, Sxx signal.spectrogram(signal_noisy, fssample_rate, npersegnseg, noverlapnseg//2) Sxx_db 20 * np.log10(Sxx) plt.subplot(2, 2, i) plt.pcolormesh(t, f, Sxx_db, shadinggouraud) plt.title(f段长度: {nseg}点) plt.ylabel(频率 [Hz]) plt.xlabel(时间 [秒]) plt.ylim(0, 6000) plt.tight_layout() plt.show()通过调整这些参数我们可以在时间分辨率和频率分辨率之间找到最佳平衡点具体取决于应用需求。例如语音分析通常需要较高的时间分辨率而机械振动分析则更关注频率分辨率。