语音信号处理实战：用Matlab和Python彻底搞懂频谱、相位谱与幅度谱
第一次接触语音信号处理时，我被各种"谱"搞得晕头转向——频谱、相位谱、幅度谱、功率谱，它们到底有什么区别？为什么FFT后得到的是复数？实部和虚部又代表什么？如果你也有类似的困惑，不妨跟着我用Matlab和Python两种工具，通过实际代码和可视化对比来彻底理清这些概念。
1. 准备工作：理解基础概念与工具选择
在开始编码之前，我们需要明确几个核心概念。频谱(Spectrum)是信号在频域中的完整表示，包含幅度和相位信息；幅度谱(Magnitude Spectrum)只关注信号在不同频率上的强度；相位谱(Phase Spectrum)则记录各频率成分的相位角。而功率谱(Power Spectrum)可以理解为幅度谱的平方，反映信号能量分布。
为什么选择Matlab和Python进行对比？Matlab在信号处理领域有着悠久历史，其内置函数和工具箱非常成熟；Python则凭借丰富的科学计算库（如NumPy、SciPy）和更友好的开源生态，正成为越来越多研究者的首选。通过两种工具的实现对比，我们能更深入地理解这些概念的本质。
首先准备一段测试音频。我建议使用16kHz采样率、16位量化的单声道WAV文件，时长约3-5秒。可以使用自己录制的语音，或者从开放音频数据集中获取。将文件命名为test.wav并放在项目目录下。
2. 信号加载与FFT基础实现
2.1 Matlab实现
在Matlab中加载音频文件非常简单：
[x, fs] = audioread('test.wav'); duration = length(x)/fs; % 计算音频时长(秒) t = (0:length(x)-1)/fs; % 时间轴进行快速傅里叶变换(FFT)：
N = length(x); % 信号长度 X = fft(x); % 执行FFT f = (0:N-1)*(fs/N); % 频率轴2.2 Python实现
Python中使用librosa或soundfile加载音频：
import numpy as np import librosa x, fs = librosa.load('test.wav', sr=None) # sr=None保持原始采样率 duration = len(x)/fs t = np.arange(len(x))/fsFFT计算使用NumPy：
N = len(x) X = np.fft.fft(x) # 执行FFT f = np.fft.fftfreq(N, 1/fs) # 频率轴2.3 FFT结果解析
无论Matlab还是Python，FFT后得到的X都是一个复数数组，这是理解各种"谱"的关键：
- 实部(Real part): 表示各频率成分的余弦幅度
- 虚部(Imaginary part): 表示各频率成分的正弦幅度
- 模(Magnitude):
abs(X)，即实部和虚部的平方和开方 - 相位角(Phase angle):
angle(X)，即虚部与实部的反正切
下表对比了两种工具的关键操作：
| 操作 | Matlab | Python (NumPy) |
|---|---|---|
| FFT | fft(x) | np.fft.fft(x) |
| 频率轴 | (0:N-1)*(fs/N) | np.fft.fftfreq(N, 1/fs) |
| 取模 | abs(X) | np.abs(X) |
| 相位角 | angle(X) | np.angle(X) |
3. 幅度谱与相位谱的可视化对比
3.1 幅度谱绘制
幅度谱展示了信号在不同频率上的能量分布。由于FFT结果是对称的（对于实数信号），我们通常只显示前半部分。
Matlab实现：
magnitude = abs(X(1:N/2+1)); % 取前半部分模值 f_plot = f(1:N/2+1); % 对应的频率轴 figure; plot(f_plot, 20*log10(magnitude)); % 用dB表示 xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)'); title('Magnitude Spectrum'); grid on;Python实现：
import matplotlib.pyplot as plt magnitude = np.abs(X[:N//2+1]) f_plot = f[:N//2+1] plt.figure() plt.plot(f_plot, 20*np.log10(magnitude)) plt.xlabel('Frequency (Hz)') plt.ylabel('Magnitude (dB)') plt.title('Magnitude Spectrum') plt.grid() plt.show()3.2 相位谱绘制
相位谱反映了各频率成分的相位信息，虽然人耳对相位不敏感，但在某些应用中（如语音合成）非常重要。
Matlab实现：
phase = angle(X(1:N/2+1)); % 取前半部分相位 figure; plot(f_plot, phase*180/pi); % 转换为角度 xlabel('Frequency (Hz)'); ylabel('Phase (degrees)'); title('Phase Spectrum'); grid on;Python实现：
phase = np.angle(X[:N//2+1]) plt.figure() plt.plot(f_plot, phase*180/np.pi) plt.xlabel('Frequency (Hz)') plt.ylabel('Phase (degrees)') plt.title('Phase Spectrum') plt.grid() plt.show()3.3 结果对比分析
观察幅度谱和相位谱，我们可以发现：
- 幅度谱通常有明显的峰值，对应语音的共振峰
- 相位谱看起来往往"杂乱无章"，这是因为：
- 语音信号的相位本身就很复杂
- 相位对时间偏移非常敏感
- 人耳对相位变化不敏感，进化过程中更关注幅度信息
提示：在实际应用中，幅度谱通常比相位谱更有用，这也是为什么很多语音处理算法（如MFCC）主要基于幅度信息。
4. 功率谱与语谱图的进阶分析
4.1 功率谱计算
功率谱是幅度谱的平方，表示信号功率在频域的分布。
Matlab实现：
power_spectrum = abs(X(1:N/2+1)).^2; figure; plot(f_plot, 10*log10(power_spectrum)); xlabel('Frequency (Hz)'); ylabel('Power (dB)'); title('Power Spectrum'); grid on;Python实现：
power_spectrum = np.abs(X[:N//2+1])**2 plt.figure() plt.plot(f_plot, 10*np.log10(power_spectrum)) plt.xlabel('Frequency (Hz)') plt.ylabel('Power (dB)') plt.title('Power Spectrum') plt.grid() plt.show()4.2 语谱图实现
语谱图是时频分析的重要工具，展示信号能量随时间变化的频率分布。
Matlab实现：
window = hann(256); % 汉宁窗 noverlap = 128; % 重叠样本数 nfft = 1024; % FFT点数 spectrogram(x, window, noverlap, nfft, fs, 'yaxis'); title('Spectrogram');Python实现：
plt.figure() plt.specgram(x, NFFT=1024, Fs=fs, window=np.hanning(256), noverlap=128, cmap='jet') plt.xlabel('Time (s)') plt.ylabel('Frequency (Hz)') plt.title('Spectrogram') plt.colorbar() plt.show()4.3 应用场景对比
不同的谱图适用于不同场景：
| 谱图类型 | 主要信息 | 典型应用场景 |
|---|---|---|
| 幅度谱 | 频率强度分布 | 音高检测、共振峰分析 |
| 相位谱 | 频率相位关系 | 语音合成、声源定位 |
| 功率谱 | 能量分布 | 语音活动检测、噪声估计 |
| 语谱图 | 时频能量变化 | 语音识别、发音分析 |
5. 常见问题与实用技巧
在实际应用中，有几个关键点需要注意：
频谱泄露问题：
- 原因：信号截断导致的边界不连续
- 解决方案：使用窗函数（如汉宁窗）
window = np.hanning(len(x)) X_windowed = np.fft.fft(x * window)频率分辨率与采样率的关系：
- 频率分辨率 = 采样率/FFT点数
- 要提高分辨率，可以增加FFT点数（补零）
对数尺度的重要性：
- 人耳对声音的感知是对数式的
- 用dB表示能更好反映听觉特性
实数信号的对称性：
- 对于实数信号，FFT结果是对称的
- 只需分析前半部分即可
两种工具的选择建议：
- 研究原型开发：Python更灵活
- 工程实现：Matlab可能更稳定
- 性能关键应用：考虑Python + Numba加速
注意：进行频域分析时，FFT点数最好选择2的幂次方，这样计算效率最高。如果信号长度不足，可以通过补零(zero-padding)来实现。
6. 从理论到实践：一个完整的语音分析案例
让我们通过一个完整的例子，分析一段包含两个不同频率正弦波的合成语音：
Matlab实现：
fs = 16000; % 采样率16kHz t = 0:1/fs:1-1/fs; % 1秒时间轴 f1 = 500; f2 = 1500; % 两个频率成分 x = 0.5*sin(2*pi*f1*t) + 0.2*sin(2*pi*f2*t + pi/4); % pi/4相位偏移 % 添加汉宁窗 window = hann(length(x))'; x_windowed = x .* window; % 计算FFT N = 4096; % 补零到4096点 X = fft(x_windowed, N); f = (0:N-1)*(fs/N); % 绘制幅度谱 figure; plot(f(1:N/2), 20*log10(abs(X(1:N/2)))); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)'); title('Two-tone Signal Spectrum'); grid on;Python实现：
import numpy as np import matplotlib.pyplot as plt fs = 16000 t = np.arange(0, 1, 1/fs) f1, f2 = 500, 1500 x = 0.5*np.sin(2*np.pi*f1*t) + 0.2*np.sin(2*np.pi*f2*t + np.pi/4) # 加窗 window = np.hanning(len(x)) x_windowed = x * window # FFT N = 4096 X = np.fft.fft(x_windowed, N) f = np.fft.fftfreq(N, 1/fs) # 绘图 plt.figure() plt.plot(f[:N//2], 20*np.log10(np.abs(X[:N//2]))) plt.xlabel('Frequency (Hz)') plt.ylabel('Magnitude (dB)') plt.title('Two-tone Signal Spectrum') plt.grid() plt.show()这个例子清晰地展示了两个频率成分（500Hz和1500Hz）在幅度谱中的峰值，以及它们之间的幅度差异。通过这样的实际案例，我们能更直观地理解频谱分析的实际意义。
