Python实战：用NumPy和Matplotlib实现时域信号到频域信号的转换

你是否曾经盯着一个音频波形图，看着它上下起伏，却想知道这个声音里到底藏着哪些“音符”？或者，在处理传感器数据时，面对一串随时间变化的数字，感觉无从下手，无法一眼看出其中主导的振动频率？这就像只看到了一个复杂机器的外壳，却不知道内部有哪些齿轮在转动。时域分析给了我们信号在时间轴上的“外貌”，而频域分析则像一台X光机，能让我们透视信号内在的频率“骨架”。

对于开发者、数据分析师和硬件工程师来说，掌握时域到频域的转换，是一项从“看现象”到“究本质”的关键技能。它能帮你从嘈杂的脑电图中分离出特定节律，从工厂设备的振动数据中提前发现故障特征，甚至优化无线通信的信号质量。今天，我们不谈深奥的数学证明，而是直接上手Python，用NumPy和Matplotlib这两把利器，亲手将时间序列“翻译”成频率图谱。我会带你一步步搭建代码，可视化每一个中间结果，并分享几个我实际调试信号时踩过的坑和总结的技巧。无论你是想处理音频、分析振动，还是单纯对信号处理感到好奇，这篇实战指南都将为你提供一个坚实、可操作的起点。

1. 理解核心：傅里叶变换的工程视角

在敲代码之前，我们得先统一一下思想。傅里叶变换听起来很高深，但从工程应用的角度看，你可以把它理解为一个极其高效的“成分分析器”。

想象一下，你面前有一杯混合果汁。时域观点就是描述这杯果汁每一口的味道（随时间变化的整体口感）。而频域观点，则是通过一台精密仪器，分析出这杯果汁里含有30%的橙子、50%的苹果和20%的芒果（各频率成分的强度）。傅里叶变换就是这台“果汁成分分析仪”。它不关心你什么时候喝，只关心里面到底有什么。

在数字世界，我们处理的是离散的信号样本，所以实际使用的是离散傅里叶变换（DFT）。而NumPy提供的np.fft.fft函数，实现的是更为高效的**快速傅里叶变换（FFT）**算法。你只需要记住一点：FFT是DFT的一种快速计算方法，对于大多数应用，我们直接调用FFT就足够了。

这里有一个关键概念叫频谱泄漏。这好比你的分析仪分辨率不够高，把橙子的味道错误地“泄漏”到了旁边苹果的味道区间里。在实际计算中，如果信号的长度不是信号周期的整数倍，就会发生这种情况。后文我们会用加窗的方法来缓解它。

注意：FFT计算得到的结果是复数，包含了每个频率成分的幅度和相位信息。我们通常更关心幅度，所以会对结果取绝对值（np.abs）。

2. 环境搭建与基础信号生成

工欲善其事，必先利其器。我们首先确保环境就绪，并创建一些用于实验的“标准信号”。

2.1 安装与导入

如果你使用的是Anaconda，那么NumPy和Matplotlib通常已经预装好了。如果没有，可以通过pip快速安装：

pip install numpy matplotlib

接下来，在Python脚本或Jupyter Notebook的开头，我们导入必要的库：

import numpy as np
import matplotlib.pyplot as plt
# 设置Matplotlib在Notebook内嵌显示（如果是Jupyter环境）
%matplotlib inline
# 设置全局字体和图像清晰度
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
plt.rcParams['figure.dpi'] = 150  # 提高图形分辨率

2.2 构造一个合成信号

为了清晰地看到转换效果，我们手动混合两个正弦波，制造一个简单的时域信号。这就像我们先知道果汁的配方，然后再用仪器去验证分析结果是否正确。

# 定义信号参数
sample_rate = 1000  # 采样率，每秒1000个点
duration = 1.0      # 信号持续时间，1秒
t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)  # 时间轴

# 定义两个频率成分
freq1 = 50   # 50 Hz
freq2 = 120  # 120 Hz

# 生成两个正弦波并叠加
signal_1 = 0.7 * np.sin(2 * np.pi * freq1 * t)  # 幅度0.7，频率50Hz
signal_2 = 1.0 * np.sin(2 * np.pi * freq2 * t)  # 幅度1.0，频率120Hz
composite_signal = signal_1 + signal_2  # 合成信号

# 添加一些高斯白噪声，让信号更接近真实情况
noise = 0.2 * np.random.randn(len(t))
noisy_signal = composite_signal + noise

现在，我们有了三个信号：两个纯净的单频信号signal_1和signal_2，一个理想的合成信号composite_signal，以及一个更真实的带噪信号noisy_signal。让我们先看看它们在时域里的样子：

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)

# 绘制纯净合成信号
axes[0].plot(t, composite_signal, linewidth=0.8)
axes[0].set_ylabel('幅度')
axes[0].set_title('纯净合成信号 (50Hz + 120Hz)')
axes[0].grid(True, alpha=0.3)

# 绘制噪声信号
axes[1].plot(t, noisy_signal, linewidth=0.8, color='orange')
axes[1].set_ylabel('幅度')
axes[1].set_title('带噪声的合成信号')
axes[1].grid(True, alpha=0.3)

# 绘制前50ms的细节，以便观察波形
zoom_t = t[:50]  # 前50个采样点，对应50ms
axes[2].plot(zoom_t, noisy_signal[:50], '.-', linewidth=1, markersize=3, color='red')
axes[2].set_xlabel('时间 [秒]')
axes[2].set_ylabel('幅度')
axes[2].set_title('信号前50毫秒细节（带噪声）')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

通过时域图，我们能直观感受到信号的波动，但很难准确说出其中包含哪几个主要频率，更别说各自的强度了。这正是频域分析要解决的问题。

3. 执行FFT与频谱可视化

现在进入核心环节：对时域信号应用FFT，并将结果以频谱图的形式呈现出来。

3.1 执行FFT计算

我们对noisy_signal进行操作，因为它更接近真实场景。

# 执行快速傅里叶变换
fft_result = np.fft.fft(noisy_signal)

# 计算频率轴
n = len(noisy_signal)  # 信号长度
freqs = np.fft.fftfreq(n, d=1/sample_rate)  # 生成对应的频率坐标

# 计算幅度谱 (取绝对值)
magnitude_spectrum = np.abs(fft_result)

# 通常我们只关心正频率部分（对于实数信号，频谱是共轭对称的）
half_n = n // 2
positive_freqs = freqs[:half_n]
positive_magnitude = magnitude_spectrum[:half_n]

这里有几个技术细节需要解释：

• np.fft.fftfreq 根据数据点数量和采样间隔，自动生成正确的频率坐标轴。
• 对于由实数组成的时域信号（绝大多数情况），其频谱关于零点对称。因此，我们通常只绘制正频率部分，它已经包含了全部信息。
• 幅度谱 magnitude_spectrum 是复数 fft_result 的模，代表了每个频率成分的能量大小。

3.2 绘制双边与单边频谱

为了理解全貌，我们先看看完整的双边频谱，再聚焦于更有用的单边频谱。

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 绘制双边频谱
axes[0].stem(freqs, magnitude_spectrum, linefmt='b-', markerfmt=' ', basefmt='k-', use_line_collection=True)
axes[0].set_xlim(-sample_rate/2, sample_rate/2)  # 显示范围限制在采样率一半以内
axes[0].set_xlabel('频率 [Hz]')
axes[0].set_ylabel('幅度')
axes[0].set_title('双边幅度频谱 (包含正负频率)')
axes[0].grid(True, alpha=0.3)
axes[0].axvline(0, color='red', linestyle='--', linewidth=0.8)  # 标记零频

# 绘制单边频谱 (仅正频率)
axes[1].stem(positive_freqs, positive_magnitude, linefmt='g-', markerfmt='go', basefmt='k-', use_line_collection=True)
axes[1].set_xlim(0, sample_rate/2)  # 奈奎斯特频率
axes[1].set_xlabel('频率 [Hz]')
axes[1].set_ylabel('幅度')
axes[1].set_title('单边幅度频谱 (仅正频率)')
axes[1].grid(True, alpha=0.3)
# 标记我们预设的频率成分
axes[1].axvline(freq1, color='red', linestyle=':', linewidth=1, alpha=0.7, label=f'{freq1} Hz')
axes[1].axvline(freq2, color='blue', linestyle=':', linewidth=1, alpha=0.7, label=f'{freq2} Hz')
axes[1].legend()

plt.tight_layout()
plt.show()

在单边频谱图上，你应该能清晰地看到在50Hz和120Hz附近出现了明显的尖峰。50Hz处的峰值幅度大约在0.7左右，120Hz处的峰值幅度大约在1.0左右，这与我们合成信号时设定的幅度比例基本吻合。周围的“毛刺”则是我们添加的高斯白噪声的贡献，它广泛分布于所有频率上。

3.3 功率谱密度：更专业的视角

在工程中，功率谱密度（PSD） 比简单的幅度谱更常用，因为它能更好地反映信号功率在频率上的分布，并且其纵坐标单位更有物理意义（如V²/Hz）。

# 计算功率谱密度 (PSD)
psd = (np.abs(fft_result) ** 2) / (sample_rate * n)  # 周期图法估计
# 对于单边谱，除了0和奈奎斯特频率点，其他点功率需要乘以2
psd_single_sided = psd[:half_n].copy()
psd_single_sided[1:-1] *= 2

# 绘制功率谱密度
plt.figure(figsize=(10, 5))
plt.plot(positive_freqs, 10 * np.log10(psd_single_sided), linewidth=1.5)  # 转换为分贝(dB)刻度
plt.xlim(0, sample_rate/2)
plt.xlabel('频率 [Hz]')
plt.ylabel('功率谱密度 [dB/Hz]')
plt.title('信号的单边功率谱密度 (PSD)')
plt.grid(True, alpha=0.3)
plt.axvline(freq1, color='red', linestyle=':', alpha=0.6, label='50 Hz 主峰')
plt.axvline(freq2, color='blue', linestyle=':', alpha=0.6, label='120 Hz 主峰')
plt.legend()
plt.tight_layout()
plt.show()

使用分贝刻度后，不同频率成分的功率对比更加直观，背景噪声的“基底”也清晰可见。这种表示方式在通信、声学等领域是标准做法。

4. 处理实际问题：频谱泄漏与加窗函数

之前我们用的是理想情况：信号恰好是1秒整，且50Hz和120Hz成分在1秒内都是完整的周期。但现实往往没那么完美。

4.1 制造频谱泄漏

让我们创建一个频率为53Hz的信号，但只采样0.8秒。53Hz在0.8秒内不是整数个周期。

# 制造非整周期信号
freq_leak = 53  # Hz
duration_leak = 0.8  # 秒
t_leak = np.linspace(0, duration_leak, int(sample_rate * duration_leak), endpoint=False)
signal_leak = np.sin(2 * np.pi * freq_leak * t_leak)

# 计算其频谱
fft_leak = np.fft.fft(signal_leak)
n_leak = len(signal_leak)
freqs_leak = np.fft.fftfreq(n_leak, d=1/sample_rate)
magnitude_leak = np.abs(fft_leak)[:n_leak//2]
freqs_leak_pos = freqs_leak[:n_leak//2]

# 绘制
plt.figure(figsize=(10, 5))
plt.stem(freqs_leak_pos, magnitude_leak, linefmt='b-', markerfmt='bo', basefmt='k-', use_line_collection=True)
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.title(f'频谱泄漏示例：{freq_leak}Hz信号，持续{duration_leak}秒（非整周期）')
plt.grid(True, alpha=0.3)
plt.axvline(freq_leak, color='red', linestyle='--', linewidth=1, label='真实频率')
plt.legend()
plt.tight_layout()
plt.show()

你会发现，频谱的尖峰变宽了，并且在主峰周围出现了许多不应该存在的“旁瓣”。能量从主频“泄漏”到了其他频率上，这就是频谱泄漏。它会降低频率分辨率，让临近的小信号被淹没。

4.2 应用加窗函数

为了抑制泄漏，我们可以在做FFT之前，给时域信号乘以一个窗函数。窗函数的两端平滑地衰减到零，强制让信号的首尾相接变得连续，从而减少截断带来的突变。最常用的窗函数是汉宁窗（Hanning Window）。

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

# 分别计算加窗前后的频谱
fft_raw = np.fft.fft(signal_leak)
fft_win = np.fft.fft(signal_windowed)
mag_raw = np.abs(fft_raw)[:n_leak//2]
mag_win = np.abs(fft_win)[:n_leak//2]

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

# 时域信号对比
axes[0, 0].plot(t_leak, signal_leak, label='原始信号')
axes[0, 0].plot(t_leak, window, label='汉宁窗', alpha=0.7)
axes[0, 0].set_xlabel('时间 [秒]')
axes[0, 0].set_ylabel('幅度')
axes[0, 0].set_title('时域：原始信号与窗函数')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(t_leak, signal_windowed, color='orange')
axes[0, 1].set_xlabel('时间 [秒]')
axes[0, 1].set_ylabel('幅度')
axes[0, 1].set_title('时域：加窗后的信号')
axes[0, 1].grid(True, alpha=0.3)

# 频域频谱对比 (线性坐标)
axes[1, 0].plot(freqs_leak_pos, mag_raw, label='不加窗', linewidth=1)
axes[1, 0].set_xlabel('频率 [Hz]')
axes[1, 0].set_ylabel('幅度')
axes[1, 0].set_title('频域：不加窗频谱 (线性)')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].axvline(freq_leak, color='red', linestyle=':', linewidth=0.8)

axes[1, 1].plot(freqs_leak_pos, mag_win, color='green', label='加汉宁窗', linewidth=1)
axes[1, 1].set_xlabel('频率 [Hz]')
axes[1, 1].set_ylabel('幅度')
axes[1, 1].set_title('频域：加汉宁窗频谱 (线性)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].axvline(freq_leak, color='red', linestyle=':', linewidth=0.8)

plt.tight_layout()
plt.show()

通过对比可以明显看到，加窗后频谱的旁瓣（两侧的“小鼓包”）被显著抑制了，主峰看起来更“干净”。但代价是主峰有轻微的展宽。这是一个典型的权衡：加窗降低了频谱泄漏，但牺牲了一些频率分辨率。不同的窗函数（如Hamming, Blackman, Kaiser）有不同的旁瓣抑制和主瓣宽度特性，可以根据具体应用选择。

提示：在测量未知信号的精确频率和幅度时，加窗是必不可少的步骤。但对于像我们最初例子那样已知是纯正弦波组合的情况，如果采样时长恰好是信号周期的整数倍，则可以不加窗，以获得最尖锐的谱线。

5. 实战案例：从音频文件中提取频谱

让我们把学到的知识用在一个更贴近实际的场景：分析一段真实的音频信号。我们将使用scipy库来读取一个WAV文件。

# 首先安装scipy（如果尚未安装）
# pip install scipy

from scipy.io import wavfile

# 假设我们有一个名为‘example_audio.wav’的音频文件
# 这里我们生成一个模拟的音频信号作为示例，因为无法提供真实文件
# 生成一段包含440Hz（标准音A）和880Hz（高八度A）的模拟音频
audio_sample_rate = 44100  # 标准音频采样率
audio_duration = 3.0       # 3秒
t_audio = np.linspace(0, audio_duration, int(audio_sample_rate * audio_duration), endpoint=False)

# 生成两个频率的音频，并加入衰减包络使其更自然
freq_a = 440.0
freq_a_high = 880.0
envelope = np.exp(-0.5 * t_audio)  # 指数衰减包络
audio_signal = envelope * (0.5 * np.sin(2 * np.pi * freq_a * t_audio) +
                           0.3 * np.sin(2 * np.pi * freq_a_high * t_audio) +
                           0.05 * np.random.randn(len(t_audio)))  # 少量噪声

# 归一化到16位PCM范围（模拟WAV文件）
audio_signal_int16 = np.int16(audio_signal / np.max(np.abs(audio_signal)) * 32767)

# 绘制音频波形（时域）
plt.figure(figsize=(12, 4))
plt.plot(t_audio[:10000], audio_signal[:10000])  # 只绘制前10000个点以便观察细节
plt.xlabel('时间 [秒]')
plt.ylabel('振幅')
plt.title('模拟音频信号波形 (时域，前0.23秒)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

现在，我们对这段“音频”进行频谱分析。由于音频信号很长，我们通常使用**短时傅里叶变换（STFT）**来观察频谱随时间的变化，这也就是声谱图（Spectrogram）的原理。但作为入门，我们先分析其中一小段稳定部分。

# 截取中间一段相对稳定的信号进行分析（避开开头和结尾的瞬态）
start_idx = int(audio_sample_rate * 1.0)  # 从第1秒开始
end_idx = start_idx + 4096  # 分析4096个采样点，约93毫秒
audio_segment = audio_signal[start_idx:end_idx]

# 应用汉宁窗后计算FFT
window_audio = np.hanning(len(audio_segment))
audio_segment_windowed = audio_segment * window_audio
fft_audio = np.fft.fft(audio_segment_windowed)
freqs_audio = np.fft.fftfreq(len(audio_segment), d=1/audio_sample_rate)
magnitude_audio = np.abs(fft_audio)[:len(audio_segment)//2]
freqs_audio_pos = freqs_audio[:len(audio_segment)//2]

# 绘制音频段频谱
plt.figure(figsize=(10, 5))
plt.plot(freqs_audio_pos, 20 * np.log10(magnitude_audio), linewidth=1)  # 分贝刻度
plt.xlim(0, 5000)  # 聚焦于0-5kHz范围
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度 [dB]')
plt.title('模拟音频信号片段频谱 (加汉宁窗)')
plt.grid(True, alpha=0.3)
plt.axvline(freq_a, color='red', linestyle='--', alpha=0.7, label=f'{int(freq_a)} Hz (A4)')
plt.axvline(freq_a_high, color='blue', linestyle='--', alpha=0.7, label=f'{int(freq_a_high)} Hz (A5)')
plt.legend()
plt.tight_layout()
plt.show()

在频谱图中，你应该能看到在440Hz和880Hz处有明显的峰值。由于我们加了衰减包络和噪声，谱线底部会有一些噪声基底，并且主峰可能不是无限尖锐，这正是真实音频分析中会遇到的情况。

为了更全面地观察，我们还可以计算并绘制整个音频文件的声谱图，它展示了频率成分随时间的变化。

from scipy import signal

# 计算声谱图
frequencies, times, Sxx = signal.spectrogram(audio_signal,
                                             fs=audio_sample_rate,
                                             window='hann',
                                             nperseg=1024,   # 每个段的长度
                                             noverlap=512)   # 段之间重叠的样本数

# 绘制声谱图
plt.figure(figsize=(12, 6))
plt.pcolormesh(times, frequencies, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis')
plt.colorbar(label='强度 [dB]')
plt.ylabel('频率 [Hz]')
plt.xlabel('时间 [秒]')
plt.title('模拟音频信号声谱图 (Spectrogram)')
plt.ylim(0, 2000)  # 限制频率显示范围
plt.tight_layout()
plt.show()

在声谱图中，纵轴是频率，横轴是时间，颜色深浅代表该频率在对应时间的能量强度。你可以看到两条明亮的水平线，分别对应440Hz和880Hz，并且它们的亮度随时间逐渐减弱（因为我们的衰减包络）。声谱图是分析非平稳信号（如音乐、语音）的利器。

6. 性能优化与高级技巧

当你处理超长信号（例如数小时的高采样率传感器数据）时，直接对整个数组做FFT可能会耗尽内存。这时，你需要一些策略。

策略一：分段处理 将长信号分割成重叠的片段，对每段做FFT，然后对结果进行平均（称为Welch方法）。这不仅能降低内存需求，还能通过平均来抑制随机噪声，获得更平滑的频谱估计。scipy.signal中的welch函数就是干这个的。

from scipy.signal import welch

# 使用Welch方法估计PSD
freqs_welch, psd_welch = welch(audio_signal,
                               fs=audio_sample_rate,
                               window='hann',
                               nperseg=1024,
                               noverlap=512,
                               scaling='density')

plt.figure(figsize=(10, 5))
plt.semilogy(freqs_welch, psd_welch)  # 使用对数纵坐标
plt.xlim(0, 2000)
plt.xlabel('频率 [Hz]')
plt.ylabel('功率谱密度 [V²/Hz]')
plt.title('使用Welch方法估计的功率谱密度 (更平滑)')
plt.grid(True, alpha=0.3)
plt.axvline(freq_a, color='red', linestyle=':', alpha=0.7)
plt.axvline(freq_a_high, color='blue', linestyle=':', alpha=0.7)
plt.tight_layout()
plt.show()

策略二：使用实数FFT 如果你的输入信号肯定是实数（没有虚部），可以使用np.fft.rfft和np.fft.rfftfreq。它们只计算正频率部分，速度更快，输出数组大小减半。

# 使用实数FFT (更快，内存更省)
fft_real = np.fft.rfft(audio_segment_windowed)
freqs_real = np.fft.rfftfreq(len(audio_segment), d=1/audio_sample_rate)
magnitude_real = np.abs(fft_real)

# 绘制对比（应与之前结果一致）
plt.figure(figsize=(10, 5))
plt.plot(freqs_real, 20 * np.log10(magnitude_real))
plt.xlim(0, 5000)
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度 [dB]')
plt.title('使用 np.fft.rfft 计算的频谱 (仅正频率)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

一个常见陷阱：幅度校正 当你对信号加窗后，时域信号的幅度被改变了，这会导致FFT计算出的频谱幅度失真。为了得到正确的幅度，需要对结果进行校正。对于汉宁窗，其相干增益约为0.5，校正因子约为2.0（1/0.5）。更严谨的做法是计算窗函数的能量，然后进行归一化。

# 幅度校正示例 (以汉宁窗为例)
window = np.hanning(N)  # N为窗口长度
window_power = np.sum(window**2)  # 窗函数的能量
scale_factor = np.sqrt(window_power)  # 幅度校正因子的一种计算方式

# 更常用的简单校正（适用于幅度谱）：
corrected_magnitude = magnitude_spectrum / (np.sum(window) / 2.0)
# 对于功率谱，校正因子不同

处理真实数据时，我习惯在计算完频谱后，先用一个已知幅度和频率的正弦波信号测试一下整个流程，确保幅度标定是正确的。这能避免后续分析中出现系统性的定量错误。

从时域到频域的转换，就像为你的数据打开了一扇新的观察窗口。最初面对密密麻麻的时域波形可能感到无从下手，但一旦掌握了FFT这个工具，并理解了窗函数、PSD、分段平均这些技巧，你就能从看似混乱的数据中提取出清晰的特征频率。无论是调试异常的机械振动，还是分析一段音频的构成，这个流程都是相通的。我自己的经验是，在编写分析脚本时，一定要把绘图部分做好，将时域波形、频谱、声谱图并列展示，这样能最快地帮你建立直觉，发现数据中的规律和问题。