Python scipy.signal 频谱分析实战：3种频谱图（线性/对数/功率）生成与对比

Python scipy.signal 频谱分析实战3种频谱图生成与深度对比在数字信号处理领域频谱分析是揭示信号频率成分的核心技术。当我们面对一段音频、振动信号或任何时域波形时频谱图就像是一把解构信号频率组成的显微镜。本文将带您使用Python的scipy.signal库通过完整代码示例生成三种最具工程价值的频谱图线性振幅谱、对数振幅谱和自功率谱并深入分析它们的应用场景与差异。1. 频谱分析基础与环境准备信号分析的本质是在时域和频域之间架起桥梁。时域信号告诉我们振幅如何随时间变化而频域则揭示这些变化背后的频率成分。傅里叶变换是连接这两个世界的数学工具它将时间函数分解为不同频率的正弦波组合。1.1 必备工具安装确保您的Python环境已安装以下核心科学计算库pip install numpy scipy matplotlib1.2 生成测试信号我们将创建一个包含多个频率成分的复合信号作为分析对象import numpy as np from scipy import signal import matplotlib.pyplot as plt # 采样参数 fs 1000 # 采样率1000Hz T 1.0 # 信号时长1秒 n int(fs * T) # 采样点数 t np.linspace(0, T, n, endpointFalse) # 生成测试信号50Hz正弦波 120Hz正弦波 随机噪声 sig 0.6 * np.sin(2*np.pi*50*t) 1.0 * np.sin(2*np.pi*120*t) noise 0.4 * np.random.normal(sizelen(t)) signal_with_noise sig noise # 绘制时域波形 plt.figure(figsize(12, 4)) plt.plot(t, signal_with_noise) plt.title(时域信号 (含噪声)) plt.xlabel(时间 [秒]) plt.ylabel(振幅) plt.grid() plt.show()这段代码生成了一个包含50Hz和120Hz正弦波分量并叠加高斯白噪声的测试信号。时域图显示信号随时间变化的振幅但无法直观看出包含哪些频率成分——这正是频谱分析要解决的问题。2. 线性振幅谱分析线性振幅谱是最直接的频谱表示方式它显示各频率成分的原始振幅大小。2.1 快速傅里叶变换实现# 计算FFT frequencies np.fft.fftfreq(n, d1/fs)[:n//2] fft_vals np.fft.fft(signal_with_noise)[:n//2] amplitude_spectrum np.abs(fft_vals) / (n//2) # 归一化 # 绘制线性振幅谱 plt.figure(figsize(12, 4)) plt.plot(frequencies, amplitude_spectrum) plt.title(线性振幅谱) plt.xlabel(频率 [Hz]) plt.ylabel(振幅) plt.xlim(0, 200) # 聚焦0-200Hz范围 plt.grid() plt.show()2.2 线性谱的特点与应用从输出结果可以清晰看到在50Hz和120Hz处有两个明显的峰值对应我们信号中的两个正弦波分量。线性振幅谱的优势在于保持原始振幅比例关系计算量小实时性好适合分析强信号成分但它的局限性也很明显当信号中存在大幅值和小幅值成分时小信号可能被噪声淹没而难以观察。这正是对数振幅谱要解决的问题。3. 对数振幅谱分析对数振幅谱通过对振幅取对数压缩动态范围使小信号得以显现。3.1 对数变换实现# 计算对数振幅谱 (dB) log_amplitude_spectrum 20 * np.log10(amplitude_spectrum 1e-12) # 避免log(0) # 绘制对数振幅谱 plt.figure(figsize(12, 4)) plt.plot(frequencies, log_amplitude_spectrum) plt.title(对数振幅谱 (dB)) plt.xlabel(频率 [Hz]) plt.ylabel(振幅 [dB]) plt.xlim(0, 200) plt.grid() plt.show()3.2 对数谱的优势与适用场景对数变换后我们可以看到50Hz分量约 -4.5dB120Hz分量约 0dB最大分量噪声基底约 -40dB这种表示方式特别适合以下场景分析微弱信号与强信号的共存系统音频处理人耳对声音的感知近似对数关系需要观察宽动态范围信号的场合注意dB计算中的1e-12是为了避免对零取对数相当于给信号加了一个极小的偏置。4. 自功率谱密度分析功率谱描述信号功率在频域的分布情况在随机信号分析和噪声研究中尤为重要。4.1 使用Welch方法计算功率谱# 使用scipy.signal.welch计算功率谱密度 freqs, psd signal.welch(signal_with_noise, fs, nperseg1024) # 绘制功率谱 plt.figure(figsize(12, 4)) plt.plot(freqs, psd) plt.title(自功率谱密度) plt.xlabel(频率 [Hz]) plt.ylabel(功率谱密度 [V²/Hz]) plt.xlim(0, 200) plt.grid() plt.show()4.2 功率谱的工程价值功率谱与振幅谱的主要区别特性振幅谱功率谱物理量振幅功率单位V或dBV²/Hz相位信息保留丢失适用信号确定性信号随机信号估计方法直接FFT平均周期图法功率谱特别适用于噪声分析与测量随机振动研究系统频率响应识别特征频率提取即使信号非周期5. 三种频谱图的参数化对比与选择指南在实际工程中频谱类型的选择取决于具体应用需求。下面通过一组对比实验展示不同参数对频谱图的影响。5.1 窗函数影响对比windows [boxcar, hann, hamming, blackman] plt.figure(figsize(12, 8)) for i, window in enumerate(windows, 1): freqs, psd signal.welch(signal_with_noise, fs, windowwindow, nperseg256) plt.subplot(2, 2, i) plt.plot(freqs, 10*np.log10(psd)) plt.title(f窗函数: {window}) plt.xlabel(频率 [Hz]) plt.ylabel(功率 [dB]) plt.xlim(0, 200) plt.grid() plt.tight_layout() plt.show()5.2 频谱分析参数选择建议根据实际测试结果我们总结出以下实用建议窗函数选择矩形窗boxcar频率分辨率最高但频谱泄漏严重汉宁窗hann平衡频率分辨率和泄漏通用推荐海明窗hamming类似汉宁窗但第一旁瓣更低布莱克曼窗blackman旁瓣抑制最好但主瓣最宽分段长度nperseg较长分段提高频率分辨率适合稳态信号较短分段提高时间分辨率适合瞬变信号重叠率noverlap通常设为分段长度的50-75%减少方差6. 实战案例电机振动信号分析让我们将这些技术应用于一个真实场景——电机振动监测。假设我们采集到以下振动信号# 模拟电机振动信号 motor_rpm 1800 # 转速1800转/分钟 shaft_freq motor_rpm / 60 # 轴转频30Hz bearing_freq shaft_freq * 3.7 # 假设轴承故障特征频率111Hz motor_vibration ( 2.0 * np.sin(2*np.pi*shaft_freq*t) 0.3 * np.sin(2*np.pi*bearing_freq*t) 0.5 * np.random.normal(sizelen(t)) ) # 计算功率谱 freqs, psd signal.welch(motor_vibration, fs, nperseg1024, windowhann) # 绘制结果 plt.figure(figsize(12, 4)) plt.plot(freqs, 10*np.log10(psd)) plt.axvline(shaft_freq, colorr, linestyle--, alpha0.5) plt.axvline(bearing_freq, colorg, linestyle--, alpha0.5) plt.title(电机振动功率谱 (红色:轴转频, 绿色:轴承故障频率)) plt.xlabel(频率 [Hz]) plt.ylabel(功率 [dB]) plt.xlim(0, 150) plt.grid() plt.show()在这个案例中对数功率谱清晰地显示了30Hz处的轴转频分量强111Hz处的轴承故障特征频率较弱宽频带的随机振动这种分析可以帮助工程师早期发现设备异常避免 catastrophic failure。7. 高级技巧与性能优化当处理长时间信号或需要实时分析时频谱计算的效率变得至关重要。7.1 快速卷积技巧# 使用FFT卷积加速滤波 def bandpass_filter(sig, lowcut, highcut, fs, order5): nyq 0.5 * fs low lowcut / nyq high highcut / nyq b, a signal.butter(order, [low, high], btypeband) return signal.filtfilt(b, a, sig) # 比较直接卷积和FFT卷积的速度 %timeit signal.lfilter(b, a, motor_vibration) # 传统时域滤波 %timeit np.convolve(motor_vibration, kernel, modesame) # 直接卷积 %timeit signal.fftconvolve(motor_vibration, kernel, modesame) # FFT卷积7.2 频谱分析参数优化表应用场景推荐窗函数分段长度重叠率频谱类型音频分析汉宁窗1024-409650-75%对数振幅谱振动监测海明窗512-204850%功率谱瞬态信号矩形窗256-10240%线性振幅谱随机噪声布莱克曼窗2048-819275%功率谱7.3 内存优化策略处理超长信号时可采用分段处理策略def chunked_psd(signal, fs, chunk_size1024, overlap0.5): step int(chunk_size * (1 - overlap)) n_chunks (len(signal) - chunk_size) // step 1 psd_sum np.zeros(chunk_size//2 1) for i in range(n_chunks): chunk signal[i*step : i*step chunk_size] freqs, psd signal.welch(chunk, fs, npersegchunk_size) psd_sum psd return freqs, psd_sum / n_chunks这种方法可以处理任意长度的信号同时保持统计可靠性。