傅里叶变换实战：用Python手把手教你分析音频频谱（附完整代码）

从波形到音符：用Python实战音频频谱分析，解锁声音的“基因图谱”

你有没有想过，一段动听的音乐、一句清晰的语音，在计算机的“耳朵”里究竟是什么样子？我们日常听到的声音是随时间变化的波形，但工程师和分析师们更常与之打交道的，却是另一种形态——频谱。这就像给声音做了一次“基因测序”，将混杂在一起的复杂波形，分解成不同频率、不同强度的“音符”成分。今天，我们不谈枯燥的数学推导，直接上手Python，用几行代码带你亲眼看见声音的“基因图谱”，并学会如何解读它。无论你是想为音频处理应用添加降噪功能，还是好奇于音乐可视化背后的原理，这篇文章都将提供一套可直接复用的工具箱。

1. 准备你的“数字耳朵”：环境与核心库搭建

在开始“听诊”声音之前，我们需要配置好分析环境。对于音频信号处理，Python生态提供了强大而便捷的工具链。核心是三个库：numpy 负责高效的数值计算和数组操作，scipy 提供了现成的信号处理函数，而 matplotlib 则是我们绘制频谱图的“画板”。如果你使用 conda 管理环境，创建和激活环境的命令会非常清晰。

# 创建并激活一个名为audio_analysis的虚拟环境
conda create -n audio_analysis python=3.9
conda activate audio_analysis

# 安装核心库
pip install numpy scipy matplotlib

对于音频文件的读写，我们有两个主流选择。librosa 是音频分析领域的瑞士军刀，功能全面且API友好，特别适合音乐信息检索。而 soundfile 或 pydub 则在基础的读写和格式转换上非常高效。我通常根据项目复杂度来选择：快速读取和分析用 librosa，需要精细控制或批量处理时用 soundfile。

# 安装音频处理库（选择其一或全部安装）
pip install librosa  # 功能全面，适合分析
# 或
pip install soundfile  # 读写高效，依赖简单

安装完成后，我们可以通过一个简单的测试脚本来验证环境是否就绪，并查看关键库的版本。

import numpy as np
import scipy
import matplotlib.pyplot as plt
import librosa

print(f"NumPy版本: {np.__version__}")
print(f"SciPy版本: {scipy.__version__}")
print(f"Librosa版本: {librosa.__version__}")

# 尝试创建一个简单的正弦波并绘图
t = np.linspace(0, 1, 1000)
signal = np.sin(2 * np.pi * 440 * t)  # 440Hz，标准A音
plt.figure(figsize=(10, 4))
plt.plot(t[:100], signal[:100])  # 只绘制前100个点以便观察
plt.title("测试信号：440Hz正弦波")
plt.xlabel("时间 (秒)")
plt.ylabel("振幅")
plt.grid(True)
plt.show()

如果能看到一个清晰的正弦波形图，恭喜你，你的“数字耳朵”已经准备就绪。接下来，我们需要理解将声音从“波形”转化为“频谱”的核心武器——离散傅里叶变换（DFT）。简单来说，DFT就像一台精密的频率分析仪，它能告诉我们一段数字信号中，各个频率成分的强度和相位。而Python中高效实现DFT的算法，就是快速傅里叶变换（FFT），它通过巧妙的计算将复杂度从O(N²)降到O(N log N)，让我们能在瞬间分析长达数小时的音频。

2. 加载与初探：将声音文件变成数字信号

理论准备就绪，现在让我们接触真实的声音数据。我们以最常见的WAV格式为例，它通常是无压缩的PCM编码，能最真实地保留原始音频信息。使用 librosa.load() 函数，我们可以轻松地将音频文件加载到内存中。

注意：librosa.load() 默认会将多声道音频（如立体声）混合成单声道，并将采样率统一为22050 Hz。这对于频谱分析通常是利大于弊的，因为它简化了数据处理。如果你需要保留原始声道信息，请设置 mono=False。

import librosa
import librosa.display

# 加载音频文件，这里以你本地的一个test.wav文件为例
audio_path = 'your_audio_file.wav'  # 请替换为你的音频文件路径
y, sr = librosa.load(audio_path, sr=None)  # sr=None 表示保持原始采样率

print(f"音频时长: {len(y)/sr:.2f} 秒")
print(f"采样率: {sr} Hz")
print(f"总采样点数: {len(y)}")
print(f"音频数据（前10个点）: {y[:10]}")

加载完成后，y 是一个一维（单声道）或二维（多声道）的NumPy数组，包含了音频的振幅序列；sr 是采样率，表示每秒采集多少个样本。采样率决定了我们能分析的最高频率，根据奈奎斯特采样定理，可分析的最高频率是采样率的一半。例如，CD音质的采样率是44100 Hz，那么它能无损表示的最高频率就是22050 Hz，这已经覆盖了人耳的听觉范围（20Hz-20kHz）。

为了对音频有一个直观的时域认识，我们可以先绘制它的波形图。

plt.figure(figsize=(14, 5))
librosa.display.waveshow(y, sr=sr, alpha=0.6)
plt.title(f"音频波形图: {audio_path}")
plt.xlabel("时间 (秒)")
plt.ylabel("振幅")
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

波形图能告诉我们声音的总体振幅变化、静音段和响亮段，但它无法揭示声音的“成分”。为了听到加载的音频，你可以在Jupyter Notebook中使用 IPython.display.Audio 控件。

from IPython.display import Audio
Audio(data=y, rate=sr)

在深入频谱分析前，还有一个关键步骤：预处理。原始音频信号可能包含直流偏移（即整个波形在零轴上方或下方平移），这会在频谱的0Hz处产生一个巨大的峰值，干扰对其他频率的分析。消除直流分量很简单：

y_processed = y - np.mean(y)  # 减去均值，消除直流偏移

另一个常见操作是归一化，将振幅缩放到[-1, 1]的范围，防止在后续计算中溢出或精度问题。

y_normalized = y_processed / np.max(np.abs(y_processed))

现在，我们手头已经有了干净、规整的数字音频信号，是时候请出主角FFT，来一窥其内在的频率奥秘了。

3. 核心魔法：执行FFT并理解频谱图

FFT是频谱分析的基石。在Python中，我们使用 numpy.fft.fft() 函数来计算它。但直接对整段音频做FFT往往不是最佳选择，因为声音的特性是随时间变化的（想想一句话中元音和辅音的交替）。因此，更实用的方法是短时傅里叶变换（STFT），它将长音频分割成许多短片段（帧），对每一帧分别进行FFT，从而得到一个随时间变化的频谱。

首先，让我们理解一下对单帧信号进行FFT的基本流程和结果含义。

import numpy.fft as fft

# 选取一段音频进行分析，例如开头的1秒钟
duration_to_analyze = 1.0  # 分析1秒
num_samples = int(duration_to_analyze * sr)
audio_segment = y_normalized[:num_samples]

# 执行FFT
N = len(audio_segment)  # 信号长度
Y = fft.fft(audio_segment)  # 得到复数数组
freqs = fft.fftfreq(N, d=1/sr)  # 计算对应的频率轴

# FFT结果是对称的，我们通常只取前半部分（正频率）
magnitude = np.abs(Y[:N//2])  # 幅度谱
phase = np.angle(Y[:N//2])    # 相位谱
freqs_pos = freqs[:N//2]

# 绘制幅度谱
plt.figure(figsize=(12, 5))
plt.plot(freqs_pos, magnitude)
plt.title("单帧音频的幅度频谱")
plt.xlabel("频率 (Hz)")
plt.ylabel("幅度")
plt.xlim(0, 5000)  # 聚焦在0-5kHz，人声和大部分乐器的关键频段
plt.grid(True)
plt.show()

这张图显示了在所选1秒时间段内，音频信号各个频率分量的强度。横轴是频率，纵轴是对应的幅度。一个纯音（如440Hz的音叉声）会在此图上显示为一条单一的竖线。而现实中的声音，如人声或音乐，则是许多不同频率线条的叠加。

然而，单帧FFT丢失了时间信息。STFT通过引入“窗函数”和“帧移”解决了这个问题。窗函数（如汉明窗）用于平滑每一帧信号的边缘，减少频谱泄漏；帧移决定了帧与帧之间的重叠程度。

# 使用Librosa计算STFT
n_fft = 2048  # FFT窗口大小，决定频率分辨率（越大分辨率越高）
hop_length = 512  # 帧移，决定时间分辨率（越小分辨率越高）
window = 'hann'  # 汉宁窗，平衡频率分辨率和频谱泄漏

D = librosa.stft(y_normalized, n_fft=n_fft, hop_length=hop_length, window=window)
# D是一个复数矩阵，行代表频率，列代表时间
magnitude_spectrogram = np.abs(D)  # 得到幅度谱图
db_spectrogram = librosa.amplitude_to_db(magnitude_spectrogram, ref=np.max)  # 转换为分贝标度，更符合人耳感知

# 绘制频谱图（Spectrogram）
plt.figure(figsize=(14, 6))
librosa.display.specshow(db_spectrogram, sr=sr, hop_length=hop_length,
                         x_axis='time', y_axis='log', cmap='viridis')
plt.colorbar(format='%+2.0f dB')
plt.title('音频频谱图 (对数频率轴)')
plt.tight_layout()
plt.show()

这张彩色的频谱图就是声音的“基因图谱”可视化。X轴是时间，Y轴是频率（这里用了对数刻度，更符合人耳对音高的感知），颜色深浅代表该时刻该频率成分的强度（分贝值）。从图中，你可以清晰地看到：

• 谐波结构：对于有音高的声音（如人声、乐器），你会看到一系列等间距的亮线，基频是最低的那条，上方是它的整数倍频（泛音）。
• 瞬态噪声：敲击声、辅音（如“t”、“p”）会表现为垂直的短促亮线，因为它们包含很宽的频率成分。
• 持续音：元音或长音会表现为水平的亮带。

为了更深入地操控分析，了解关键参数的影响至关重要：

参数	含义	影响	典型值建议
n_fft	FFT窗口长度	决定频率分辨率。值越大，频率区分越精细，但时间模糊度增加。	2048, 4096
hop_length	帧移长度	决定时间分辨率。值越小，时间定位越精确，但计算量增大，频率分辨率间接受影响。	512 (通常为n_fft的1/4)
window	窗函数类型	减少频谱泄漏。汉明窗、汉宁窗最常用。	‘hamming‘, ‘hann‘
sr	采样率	决定可分析的最高频率（sr/2）。	22050, 44100

你可以尝试调整这些参数，重新生成频谱图，观察图像细节的变化。例如，将 n_fft 增加到4096，你会看到频率轴上的线条变得更细、更清晰；而将 hop_length 减小到256，时间轴上的变化则会更加连续平滑。

4. 从分析到应用：频谱解读与实战案例

看懂频谱图只是第一步，更重要的是利用这些信息解决实际问题。下面我们通过几个典型场景，将频谱分析技术落地。

案例一：检测音频中的主旋律或基频 在音乐分析或语音处理中，找出主导频率（基频，F0）非常有用。librosa 提供了 pyin 或 yin 等算法来估计基频轨迹。

# 使用PYIN算法估计基频（更适合音乐和有调性声音）
f0, voiced_flag, voiced_probs = librosa.pyin(y_normalized,
                                             fmin=librosa.note_to_hz('C2'), # 最低搜索频率
                                             fmax=librosa.note_to_hz('C7'), # 最高搜索频率
                                             sr=sr)

# 将基频轨迹叠加在频谱图上显示
times = librosa.times_like(f0, sr=sr, hop_length=hop_length)
plt.figure(figsize=(14, 8))

# 绘制频谱图作为背景
librosa.display.specshow(db_spectrogram, sr=sr, hop_length=hop_length,
                         x_axis='time', y_axis='log', cmap='gray_r', alpha=0.7)

# 叠加基频轨迹
plt.plot(times, f0, linewidth=2, color='cyan', label='Estimated F0 (Pitch)')
plt.ylim([librosa.note_to_hz('C2'), librosa.note_to_hz('C7')])
plt.title('频谱图与基频轨迹叠加')
plt.legend()
plt.colorbar(format='%+2.0f dB')
plt.tight_layout()
plt.show()

图中那条蜿蜒的亮青色曲线就是估计出的基频。你可以清晰地看到旋律的起伏，或者语音中声调的变化。

案例二：识别和过滤特定频率（如去除嗡嗡声） 环境中的持续嗡嗡声（如50/60Hz的电源噪声）在频谱图上会表现为一条顽固的水平亮线。我们可以设计一个陷波滤波器来削弱它。

from scipy import signal

# 假设我们要消除一个60Hz的电源噪声
notch_freq = 60.0
quality_factor = 30.0  # Q值，决定滤波器带宽

# 设计一个IIR陷波滤波器
b, a = signal.iirnotch(notch_freq, quality_factor, sr)

# 应用滤波器
y_filtered = signal.filtfilt(b, a, y_normalized)  # 使用filtfilt实现零相位滤波

# 对比滤波前后的频谱
D_before = librosa.stft(y_normalized, n_fft=n_fft, hop_length=hop_length)
D_after = librosa.stft(y_filtered, n_fft=n_fft, hop_length=hop_length)

fig, ax = plt.subplots(2, 1, figsize=(14, 10))
librosa.display.specshow(librosa.amplitude_to_db(np.abs(D_before), ref=np.max),
                         sr=sr, hop_length=hop_length, x_axis='time', y_axis='log',
                         ax=ax[0], cmap='viridis')
ax[0].set_title('原始音频频谱图')
ax[0].axhline(y=notch_freq, color='r', linestyle='--', alpha=0.7, label=f'{notch_freq}Hz')

librosa.display.specshow(librosa.amplitude_to_db(np.abs(D_after), ref=np.max),
                         sr=sr, hop_length=hop_length, x_axis='time', y_axis='log',
                         ax=ax[1], cmap='viridis')
ax[1].set_title('去除60Hz噪声后的频谱图')
ax[1].axhline(y=notch_freq, color='r', linestyle='--', alpha=0.7, label=f'{notch_freq}Hz')
ax[0].legend()
ax[1].legend()
plt.tight_layout()
plt.show()

对比上下两图，你应该能看到在60Hz附近，下方图的红色虚线处的亮线明显减弱或消失了。

案例三：基于频谱特征的简单音频分类 不同的声音在频谱上有不同的“指纹”。例如，掌声是宽频带的短促爆发，鸟鸣是集中在高频的短促谐波。我们可以提取一些统计特征来区分它们。

def extract_spectral_features(audio_signal, sr):
    """提取一组简单的频谱特征"""
    # 计算梅尔频谱图（更接近人耳听觉）
    mel_spec = librosa.feature.melspectrogram(y=audio_signal, sr=sr, n_mels=128)
    mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max)

    # 计算一些统计特征
    spectral_centroid = librosa.feature.spectral_centroid(y=audio_signal, sr=sr)[0]
    spectral_bandwidth = librosa.feature.spectral_bandwidth(y=audio_signal, sr=sr)[0]
    spectral_rolloff = librosa.feature.spectral_rolloff(y=audio_signal, sr=sr)[0]

    features = {
        'mel_mean': np.mean(mel_spec_db, axis=1),
        'spectral_centroid_mean': np.mean(spectral_centroid),
        'spectral_bandwidth_mean': np.mean(spectral_bandwidth),
        'spectral_rolloff_mean': np.mean(spectral_rolloff)
    }
    return features

# 示例：对两段不同的音频提取特征
# 假设你有两段音频数据：applause（掌声）和 bird（鸟鸣）
# features_applause = extract_spectral_features(applause_audio, sr)
# features_bird = extract_spectral_features(bird_audio, sr)

# 可以打印或比较这些特征字典，它们会显著不同。
# 在实际分类任务中，你会用这些特征向量去训练一个机器学习模型（如SVM、随机森林）。

通过计算频谱质心（声音的“亮度”指标）、频谱带宽（能量分布的宽度）和频谱滚降点（频率累计能量达到一定比例的频率点），我们可以量化声音的频谱特性。将这些特征输入分类器，就能实现简单的自动音频标签。

5. 性能优化与高级技巧

当你开始处理更长的音频或需要实时分析时，性能就变得关键。以下是一些提升代码效率的实用技巧：

1. 实时流式处理模拟 对于超长音频或实时音频流，一次性加载到内存计算STFT可能不可行。我们可以模拟一个滑动窗口来处理：

def process_audio_stream(audio_data, sr, window_seconds=0.1, hop_seconds=0.05):
    """模拟流式处理，每次处理一小段音频"""
    window_samples = int(window_seconds * sr)
    hop_samples = int(hop_seconds * sr)

    for start in range(0, len(audio_data) - window_samples, hop_samples):
        chunk = audio_data[start:start + window_samples]
        # 对chunk执行FFT或其他分析
        D_chunk = librosa.stft(chunk, n_fft=1024, hop_length=256)
        # 在这里进行你的实时分析逻辑，例如检测特定频率...
        # 示例：计算当前块的平均频谱能量
        current_energy = np.mean(np.abs(D_chunk)**2)
        # print(f"Time: {start/sr:.2f}s, Energy: {current_energy:.4f}")

2. 选择合适的FFT长度 n_fft 并非越大越好。长度最好是2的整数次幂（如256, 512, 1024, 2048），这样FFT算法效率最高。同时，需要权衡频率分辨率和时间分辨率。一个经验法则是：n_fft 应能覆盖你感兴趣的最短时间特征。

3. 利用对数频率和梅尔尺度 人耳对频率的感知是对数型的，对100Hz到200Hz的变化感知，与对1000Hz到2000Hz的变化感知类似。因此，在分析音乐或语音时，使用对数频率轴或梅尔频谱图往往更有意义。

# 计算梅尔频谱图
mel_spec = librosa.feature.melspectrogram(y=y_normalized, sr=sr, n_mels=128)
mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max)

plt.figure(figsize=(14, 5))
librosa.display.specshow(mel_spec_db, sr=sr, x_axis='time', y_axis='mel', cmap='magma')
plt.colorbar(format='%+2.0f dB')
plt.title('梅尔频谱图 (更符合人耳听觉)')
plt.tight_layout()
plt.show()

梅尔频谱图将线性频率轴映射到梅尔刻度，在低频部分有更高的分辨率，在高频部分分辨率较低，这与人耳的听觉特性相匹配，在语音识别和音乐分类任务中是标准的特征输入。

4. 处理静音和端点检测 在实际应用中，音频常包含大量静音段。直接分析会浪费算力并稀释有效特征。librosa.effects.split 或基于能量的VAD（语音活动检测）可以帮助我们定位有效音频段。

# 基于能量进行静音区间检测和分割
intervals = librosa.effects.split(y_normalized, top_db=20)  # 低于峰值20dB以下视为静音

print(f"检测到 {len(intervals)} 个有效音频区间")
for idx, (start, end) in enumerate(intervals):
    print(f"  区间 {idx+1}: {start/sr:.2f}s 到 {end/sr:.2f}s")
    # 可以针对每个区间单独进行频谱分析

掌握这些优化技巧后，你的音频分析代码将更加高效和健壮。最后，别忘了将处理后的频谱数据或特征保存下来，供后续机器学习模型使用。常用的保存格式有NPZ（NumPy压缩格式）、CSV或HDF5。

# 保存梅尔频谱特征为.npz文件
np.savez_compressed('audio_features.npz',
                    mel_spectrogram=mel_spec_db,
                    spectral_centroid=spectral_centroid,
                    sample_rate=sr)

频谱分析的世界远不止于此，从声纹识别、乐器分类到环境声音监测，其应用无处不在。我最初在处理一段带有严重背景风扇噪声的访谈录音时，正是通过仔细观察频谱图，定位到噪声的精确频带，然后用一个精心设计的带阻滤波器将其滤除，才让对话内容清晰可辨。这个过程让我深刻体会到，将抽象的数学工具与直观的可视化结合，是解决实际工程问题最有力的方式。