Python实战：用NumPy和Matplotlib实现FFT频谱分析（附完整代码）-优快云博客

Python实战：用NumPy和Matplotlib实现FFT频谱分析（附完整代码）

你是否曾经面对一段音频波形、一组传感器数据，或者任何随时间变化的信号，感到无从下手？那些看似杂乱无章的波动背后，往往隐藏着决定性的规律——频率。在数字信号处理的世界里，快速傅里叶变换（FFT）就是一把打开这扇大门的万能钥匙。它能够将时间域的信号“翻译”成频率域的语言，让我们清晰地看到构成复杂信号的各个频率成分及其强度。对于数据分析师、算法工程师，甚至是硬件开发者而言，掌握FFT的实战应用，就如同掌握了一门洞察数据本质的“内功心法”。

今天，我们不打算深陷于复杂的数学公式推导，而是直接切入Python编程的实战场景。我们将手把手地，从零开始，使用NumPy和Matplotlib这两个强大的库，构建一个完整的FFT频谱分析流程。无论你是刚刚接触信号处理的初学者，还是希望快速将理论转化为代码的实践者，这篇文章都将为你提供可直接运行、可修改、可扩展的代码范例和清晰的思路。我们将从生成一个模拟信号开始，一步步完成FFT计算、结果解读，并最终实现专业级的可视化图表，让你在实战中真正理解频谱分析的魅力。

1. 环境准备与核心库概览

在开始编码之前，确保你的Python环境已经就绪。我强烈建议使用Anaconda来管理环境，它能避免很多依赖库的冲突问题。当然，使用纯净的Python配合pip安装也是完全可行的。

首先，我们需要安装两个核心库：NumPy 和 Matplotlib。如果你使用pip，打开终端或命令提示符，执行以下命令：

pip install numpy matplotlib

对于追求更佳绘图效果和交互体验的用户，也可以考虑安装seaborn库来美化图表，但这并非必需：

pip install seaborn

现在，让我们快速认识一下今天的主角：

NumPy：Python科学计算的基石。它提供了高效的多维数组对象ndarray，以及一系列用于数组操作的函数。在FFT分析中，我们几乎所有的数学运算——从信号生成到复数计算——都将依赖它。其numpy.fft模块封装了高效的FFT算法实现，是我们今天工作的核心引擎。
Matplotlib：Python绘图库的“事实标准”。我们将主要使用其pyplot子模块（通常别名为plt）来创建各种静态、交互式的图表，将计算得到的频谱数据直观地呈现出来。

一个常见的误区是认为FFT分析非常复杂，需要自己实现算法。实际上，得益于NumPy等成熟库，我们只需调用几个函数即可完成核心计算，从而将精力集中在信号的理解和结果的分析上。下面是一个最简单的导入代码块，我们后续的所有代码都将基于此：

import numpy as np
import matplotlib.pyplot as plt

# 为了在Jupyter Notebook中直接显示图表，可以取消下面这行的注释
# %matplotlib inline

print("NumPy版本:", np.__version__)
print("Matplotlib版本:", plt.matplotlib.__version__)

运行这段代码，如果没有报错并成功打印出版本号，那么恭喜你，环境已经准备就绪。

2. 构建你的第一个测试信号

在进行频谱分析之前，我们首先需要一个分析对象——信号。直接从真实世界采集的数据（如音频、振动）往往包含噪声和复杂的成分，不利于初学者理解。因此，我们从一个完全由自己控制的“纯净”合成信号开始。这能让我们清晰地验证FFT分析的正确性。

一个最经典的合成信号是多个正弦波的叠加。正弦波是频率成分最单一的信号，多个不同频率、不同幅度的正弦波叠加，就能模拟出许多复杂信号的特性。

假设我们想模拟一个采样频率为1000 Hz（即每秒采集1000个数据点），持续时间为1秒的信号。其中包含两个主要的频率成分：一个强力的50 Hz信号和一个稍弱的120 Hz信号。此外，为了模拟真实环境，我们还可以加入一点随机噪声。以下是生成该信号的代码：

# 信号参数设置
fs = 1000               # 采样频率 (Hz)
T = 1                   # 信号总时长 (秒)
N = int(fs * T)         # 总采样点数
t = np.linspace(0, T, N, endpoint=False) # 时间轴，不包含终点以避免重复

# 信号的频率成分
freq1 = 50    # 第一个频率成分 (Hz)
amp1 = 0.7    # 第一个频率的振幅
freq2 = 120   # 第二个频率成分 (Hz)
amp2 = 0.3    # 第二个频率的振幅

# 生成纯净的合成信号
signal_pure = amp1 * np.sin(2 * np.pi * freq1 * t) + amp2 * np.sin(2 * np.pi * freq2 * t)

# 添加一些高斯白噪声，模拟真实情况
noise_amplitude = 0.05
noise = noise_amplitude * np.random.randn(N)
signal_noisy = signal_pure + noise

# 可视化原始信号
fig, axes = plt.subplots(2, 1, figsize=(10, 6))

# 绘制纯净信号
axes[0].plot(t, signal_pure, linewidth=1.0, color='blue', alpha=0.7)
axes[0].set_title('纯净合成信号 (50Hz & 120Hz)')
axes[0].set_xlabel('时间 [秒]')
axes[0].set_ylabel('振幅')
axes[0].grid(True, linestyle='--', alpha=0.5)
axes[0].set_xlim(0, 0.1)  # 只看前0.1秒，细节更清晰

# 绘制含噪信号
axes[1].plot(t, signal_noisy, linewidth=1.0, color='red', alpha=0.7)
axes[1].set_title('含噪声的合成信号')
axes[1].set_xlabel('时间 [秒]')
axes[1].set_ylabel('振幅')
axes[1].grid(True, linestyle='--', alpha=0.5)
axes[1].set_xlim(0, 0.1)

plt.tight_layout()
plt.show()

运行这段代码，你会看到两个并排的子图，分别展示了纯净信号和加噪信号在前0.1秒内的波形。即使加入了噪声，你依然能看出信号周期性波动的“轮廓”。这个时域波形图告诉我们信号如何随时间变化，但它无法直接告诉我们信号由哪些频率组成。这就是我们需要FFT的原因。

提示：采样频率fs的选择至关重要。根据奈奎斯特采样定理，要无失真地还原一个信号，采样频率必须大于信号中最高频率成分的两倍。在我们的例子中，最高频率是120Hz，因此fs=1000Hz是绰绰有余的。如果fs设置过低，会导致高频成分“混叠”到低频中，造成分析错误。

3. 执行FFT计算与频谱解读

有了信号数据，接下来就是施展FFT魔法的时刻。使用NumPy进行FFT计算异常简单，核心函数就是np.fft.fft()。但是，理解其输入输出背后的意义，才是正确解读结果的关键。

3.1 执行FFT计算

我们直接对上一节生成的signal_noisy进行分析。

# 对含噪信号执行FFT
fft_values = np.fft.fft(signal_noisy)

# 计算对应的频率轴
n = len(signal_noisy) # 信号长度
freqs = np.fft.fftfreq(n, d=1/fs) # d是采样间隔，等于1/fs

# 取单边频谱（因为频谱是关于0Hz对称的）
half_n = n // 2
fft_magnitude = np.abs(fft_values[:half_n]) * 2 / n # 求幅度，并归一化
single_sided_freqs = freqs[:half_n]

短短几行代码，计算就完成了。我们来拆解一下：

np.fft.fft(signal_noisy)：返回一个复数数组，长度与输入信号n相同。每个复数对应一个频率分量的信息，包含了该分量的幅度和相位。
np.fft.fftfreq(n, d=1/fs)：生成一个与FFT结果对应的频率数组。其范围从-fs/2到+fs/2（实际上略小于）。参数d是采样时间间隔。
由于实数信号的频谱是共轭对称的（即负频率部分是正频率部分的镜像），我们通常只关心正频率部分（从0Hz到fs/2）。half_n = n // 2取得了正频率部分的索引。
np.abs(fft_values)计算复数的模，即该频率分量的幅度。* 2 / n是一个常见的归一化操作，使得频谱图中幅度的刻度与实际信号振幅相匹配（对于除直流分量外的频率）。直流分量（0Hz）通常需要单独处理。

3.2 关键概念与结果解读

FFT的输出直接看是难以理解的。下表总结了几个核心概念，帮助你建立直观认识：

概念	解释	在代码中的体现	注意事项
频率分辨率	频谱图上相邻频率点之间的间隔。它决定了你能区分多近的两个频率。	`df = fs / n`	信号总时长`T`越长（`n`越大），分辨率`df`越小，频谱越精细。
奈奎斯特频率	能够被无混叠分析的最高频率，等于采样频率`fs`的一半。	`f_nyquist = fs / 2`	所有高于此频率的信号成分都会“折叠”回来，造成混叠失真。
幅度谱	各频率分量振幅的大小。是我们最常观察的频谱图。	`np.abs(fft_values)`	反映了该频率成分对原始信号“能量”的贡献大小。
相位谱	各频率分量在时间轴上的起始位置。	`np.angle(fft_values)`	在需要重建信号或分析信号间时序关系时非常重要。
频谱泄漏	由于信号长度有限（非整数个周期），导致能量“泄漏”到旁瓣频率的现象。	表现为主峰周围有拖尾。	可通过加窗函数（如汉宁窗）来抑制。

现在，让我们可视化幅度谱，看看能否从中找到我们预设的50Hz和120Hz。

# 绘制幅度频谱图
plt.figure(figsize=(12, 5))

# 绘制线性坐标下的频谱
plt.subplot(1, 2, 1)
plt.plot(single_sided_freqs, fft_magnitude, color='green', linewidth=1.5)
plt.axvline(x=freq1, color='blue', linestyle='--', alpha=0.7, label=f'{freq1} Hz')
plt.axvline(x=freq2, color='red', linestyle='--', alpha=0.7, label=f'{freq2} Hz')
plt.title('幅度频谱 (线性坐标)')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.xlim(0, 200) # 聚焦在0-200Hz范围
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)

# 绘制对数坐标下的频谱（更能看清小信号）
plt.subplot(1, 2, 2)
plt.plot(single_sided_freqs, 20 * np.log10(fft_magnitude + 1e-10), color='purple', linewidth=1.5) # 加极小值避免log(0)
plt.axvline(x=freq1, color='blue', linestyle='--', alpha=0.7, label=f'{freq1} Hz')
plt.axvline(x=freq2, color='red', linestyle='--', alpha=0.7, label=f'{freq2} Hz')
plt.title('幅度频谱 (对数坐标，单位: dB)')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度 [dB]')
plt.xlim(0, 200)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

在这张图中，你应该能清晰地看到在两个虚线标记的位置（50Hz和120Hz），出现了尖锐的谱峰。它们的幅度比例大约为0.7:0.3，与我们生成信号时设定的振幅比例相符。对数坐标图将幅度转换成了分贝（dB），这使得较小的频率成分（如我们添加的噪声）也变得可见，噪声在整个频带内呈现为一条抬高的“基底”。

4. 高级技巧与实战问题处理

掌握了基础分析后，我们面对真实数据时总会遇到各种挑战。这一节，我们将探讨几个提升分析质量的关键技巧。

4.1 使用窗函数抑制频谱泄漏

前面我们提到过频谱泄漏。当你分析的信号片段不是其周期的整数倍时，就会发生泄漏。这会导致本应集中在一个频率点上的能量，扩散到周围的频率上，使谱峰变宽、旁边出现虚假的旁瓣，影响频率识别和幅度测量的准确性。

加窗是抑制泄漏的标准方法。窗函数在时域上与信号相乘，使信号的起始和结束平滑地过渡到零，从而减少截断带来的突变。常用的窗函数有汉宁窗（Hanning）、汉明窗（Hamming）、布莱克曼窗（Blackman）等。

# 生成一个非整数周期的信号，演示泄漏
freq_leak = 50.5  # 非整数周期频率
t_leak = np.linspace(0, 1, fs, endpoint=False)
signal_leak = np.sin(2 * np.pi * freq_leak * t_leak)

# 应用汉宁窗
window = np.hanning(len(signal_leak))
signal_windowed = signal_leak * window

# 计算加窗前后的FFT
fft_raw = np.fft.fft(signal_leak)
fft_win = np.fft.fft(signal_windowed)
freqs = np.fft.fftfreq(len(signal_leak), d=1/fs)
half = len(signal_leak)//2

# 可视化对比
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 时域信号对比
axes[0, 0].plot(t_leak[:100], signal_leak[:100], label='原始信号')
axes[0, 0].set_title('时域原始信号 (前100点)')
axes[0, 0].set_xlabel('时间 [秒]')
axes[0, 0].set_ylabel('振幅')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(t_leak[:100], signal_windowed[:100], color='orange', label='加窗后信号')
axes[0, 1].plot(t_leak[:100], window[:100], color='red', linestyle='--', alpha=0.5, label='汉宁窗')
axes[0, 1].set_title('时域加窗信号 (前100点)')
axes[0, 1].set_xlabel('时间 [秒]')
axes[0, 1].set_ylabel('振幅')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 频域频谱对比 (对数坐标)
axes[1, 0].plot(freqs[:half], 20*np.log10(np.abs(fft_raw[:half]) + 1e-10))
axes[1, 0].set_title('原始信号频谱 (泄漏严重)')
axes[1, 0].set_xlabel('频率 [Hz]')
axes[1, 0].set_ylabel('幅度 [dB]')
axes[1, 0].set_xlim(40, 60)
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(freqs[:half], 20*np.log10(np.abs(fft_win[:half]) + 1e-10), color='orange')
axes[1, 1].set_title('加窗后信号频谱 (泄漏抑制)')
axes[1, 1].set_xlabel('频率 [Hz]')
axes[1, 1].set_ylabel('幅度 [dB]')
axes[1, 1].set_xlim(40, 60)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

通过对比图可以明显看到，加窗前（左下），50.5Hz的主峰周围有很多高的旁瓣，能量泄漏严重。加汉宁窗后（右下），旁瓣被显著抑制，主峰更“干净”，但代价是主峰略有展宽。这是一个典型的权衡：用频率分辨率换取频谱纯度。

4.2 处理实信号与提高频率估计精度

有时我们不仅要知道频率的大致位置，还需要更精确地估计频率值，尤其是在谱峰不在频率分辨率的整数倍位置上时。一种简单有效的方法是零填充。

零填充是在进行FFT之前，在信号末尾补充若干个0。它并不增加信号的真实信息，因此不能提高真实的频率分辨率，但可以对频谱进行插值，让谱线看起来更平滑，有助于通过观察找到谱峰的顶点，从而更精确地估计频率。

# 对一个接近但不等于频率分辨率整倍数的信号进行分析
f_precise = 52.3  # 精确频率
signal_precise = np.sin(2 * np.pi * f_precise * t)

# 方法1：无零填充
fft_no_pad = np.fft.fft(signal_precise)
freqs_no_pad = np.fft.fftfreq(len(signal_precise), d=1/fs)
half_no_pad = len(signal_precise)//2

# 方法2：进行零填充（将信号长度扩展到原来的4倍）
zero_pad_factor = 4
signal_padded = np.pad(signal_precise, (0, len(signal_precise)*(zero_pad_factor-1)), 'constant')
fft_padded = np.fft.fft(signal_padded)
freqs_padded = np.fft.fftfreq(len(signal_padded), d=1/fs)
half_padded = len(signal_padded)//2

# 绘制对比
plt.figure(figsize=(14, 5))

plt.subplot(1, 2, 1)
plt.stem(freqs_no_pad[:half_no_pad], np.abs(fft_no_pad[:half_no_pad]), basefmt=" ", use_line_collection=True)
plt.title(f'无零填充频谱 (分辨率: {fs/len(signal_precise):.2f} Hz)')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.xlim(50, 55)
plt.grid(True, alpha=0.3)
# 找到峰值位置
peak_idx_no_pad = np.argmax(np.abs(fft_no_pad[:half_no_pad]))
estimated_freq_no_pad = freqs_no_pad[peak_idx_no_pad]
plt.axvline(x=estimated_freq_no_pad, color='red', linestyle='--', label=f'估计: {estimated_freq_no_pad:.2f} Hz')
plt.axvline(x=f_precise, color='green', linestyle=':', label=f'真实: {f_precise} Hz')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(freqs_padded[:half_padded], np.abs(fft_padded[:half_padded]), linewidth=1.5)
plt.title(f'4倍零填充后频谱 (插值后)')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.xlim(50, 55)
plt.grid(True, alpha=0.3)
# 在插值后的频谱上找峰值
peak_idx_padded = np.argmax(np.abs(fft_padded[:half_padded]))
estimated_freq_padded = freqs_padded[peak_idx_padded]
plt.axvline(x=estimated_freq_padded, color='red', linestyle='--', label=f'估计: {estimated_freq_padded:.2f} Hz')
plt.axvline(x=f_precise, color='green', linestyle=':', label=f'真实: {f_precise} Hz')
plt.legend()

plt.tight_layout()
plt.show()

print(f"真实频率: {f_precise} Hz")
print(f"无零填充估计频率: {estimated_freq_no_pad:.2f} Hz (误差: {abs(estimated_freq_no_pad - f_precise):.2f} Hz)")
print(f"零填充后估计频率: {estimated_freq_padded:.2f} Hz (误差: {abs(estimated_freq_padded - f_precise):.2f} Hz)")

左图使用stem绘制离散的谱线，由于频率分辨率是1Hz（1000Hz/1000点），谱峰落在52Hz的格点上，与真实频率52.3Hz有误差。右图经过零填充插值后，频谱曲线变得连续平滑，我们可以更准确地定位到谱峰顶点在52.3Hz附近，估计精度大大提高。

5. 综合实战：分析一段真实音频片段

最后，我们将所有知识融会贯通，处理一个更接近真实世界的任务：分析一段简短音频文件的频谱。这里我们将使用scipy库来读取音频文件，如果你没有安装，请先运行pip install scipy。

注意：由于我们无法在代码中直接包含音频文件，以下示例假设你有一个名为sample_audio.wav的WAV格式文件。你可以录制一段包含固定频率（如用手机软件生成一个440Hz的A4标准音）的音频，或者使用任何简短的音效文件。

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

# 1. 读取音频文件
# 请将 'sample_audio.wav' 替换为你的音频文件路径
try:
    sample_rate, audio_data = wavfile.read('sample_audio.wav')
except FileNotFoundError:
    print("未找到音频文件，将使用模拟音频数据代替。")
    # 生成模拟音频数据：一个A4标准音 (440Hz) 加一点背景噪声
    sample_rate = 44100
    duration = 2.0  # 2秒
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    A4_freq = 440.0
    audio_data = 0.5 * np.sin(2 * np.pi * A4_freq * t) + 0.1 * np.sin(2 * np.pi * 880 * t) # 加一个倍频
    audio_data = audio_data + 0.02 * np.random.randn(len(t)) # 加噪声
    # 将浮点数数据转换为16位PCM格式模拟
    audio_data = (audio_data * 32767).astype(np.int16)

# 如果音频是双声道，取左声道或转换为单声道
if audio_data.ndim > 1:
    print(f"音频为 {audio_data.shape[1]} 声道，取左声道进行分析。")
    audio_data = audio_data[:, 0]

# 2. 截取一段进行分析（例如前0.5秒），避免计算量过大
analysis_duration = 0.5  # 分析0.5秒
num_samples = int(sample_rate * analysis_duration)
audio_segment = audio_data[:num_samples].astype(float) # 转换为浮点型便于计算
time_axis = np.arange(num_samples) / sample_rate

# 3. 加窗并计算FFT
window = np.hanning(num_samples)
windowed_audio = audio_segment * window

fft_result = np.fft.fft(windowed_audio)
freq_axis = np.fft.fftfreq(num_samples, d=1/sample_rate)

# 取单边频谱
half_n = num_samples // 2
magnitude = np.abs(fft_result[:half_n]) * 2 / np.sum(window) # 使用窗函数能量进行归一化
single_sided_freq = freq_axis[:half_n]

# 4. 绘制时域波形和频谱图
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 时域图
axes[0].plot(time_axis, audio_segment, linewidth=0.8, alpha=0.7)
axes[0].set_title('音频信号时域波形 (片段)')
axes[0].set_xlabel('时间 [秒]')
axes[0].set_ylabel('振幅')
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, analysis_duration)

# 频域图 (对数坐标)
axes[1].plot(single_sided_freq, 20 * np.log10(magnitude + 1e-10), linewidth=1.0)
axes[1].set_title('音频信号幅度频谱 (加汉宁窗)')
axes[1].set_xlabel('频率 [Hz]')
axes[1].set_ylabel('幅度 [dB]')
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(20, sample_rate/2)  # 通常关注20Hz以上
axes[1].set_ylim(-80, 0)  # 设置动态范围

# 标记可能的基频和泛音
# 这里简单演示：找到前5个峰值
from scipy.signal import find_peaks
peaks, properties = find_peaks(20*np.log10(magnitude+1e-10), height=-50, distance=50) # 设置最小高度和峰间距离
for i, peak in enumerate(peaks[:5]):
    freq_peak = single_sided_freq[peak]
    db_peak = 20 * np.log10(magnitude[peak] + 1e-10)
    axes[1].plot(freq_peak, db_peak, "rx")
    axes[1].annotate(f'{freq_peak:.1f} Hz', xy=(freq_peak, db_peak), xytext=(5, 5),
                     textcoords='offset points', fontsize=9, color='red')

plt.tight_layout()
plt.show()

# 打印一些基本信息
print(f"采样率: {sample_rate} Hz")
print(f"分析片段长度: {analysis_duration} 秒")
print(f"频率分辨率: {sample_rate/num_samples:.2f} Hz")
print(f"检测到的主要峰值频率: {single_sided_freq[peaks[:5]]}")

这段代码完成了一个简易的音频频谱分析仪的核心功能。它读取（或模拟）音频数据，截取片段，加窗，计算FFT并绘制出专业的频谱图。通过scipy.signal.find_peaks函数，我们还能自动识别并标记出频谱中的主要峰值，这对于分析乐器音高、机器故障特征频率等应用非常有用。

在实际项目中，你可能还需要考虑平均多个时间段的频谱（短时傅里叶变换，STFT）来观察频率随时间的变化，或者使用更精确的峰值检测算法。但无论如何，你现在已经拥有了从信号生成、FFT计算、结果可视化和解读，到应用高级技巧处理实际问题的完整工具箱。