信号处理实战：5个用Python实现傅里叶变换的经典案例（附完整代码）-CSDN博客

信号处理实战：5个用Python实现傅里叶变换的经典案例（附完整代码）

傅里叶变换，这个听起来有些高深的数学工具，其实早已渗透到我们数字生活的方方面面。从你手机播放的音乐，到社交媒体上分享的照片，再到智能手表监测的心率，背后都有它的身影。对于工程师和开发者而言，理解其理论固然重要，但更重要的是能将它转化为解决实际问题的代码。如果你曾对着一堆公式感到无从下手，或者觉得理论推导与工程实现之间隔着一道鸿沟，那么这篇文章正是为你准备的。我们将彻底抛开繁琐的数学证明，直接切入Python的numpy.fft和scipy.fft模块，通过五个从易到难、覆盖不同领域的实战案例，手把手带你将傅里叶变换从理论公式变成可运行、可调试、可优化的工程代码。你会发现，掌握了这些核心技巧，处理音频降噪、图像压缩、信号滤波等问题，将变得有章可循。

1. 基础热身：从时域到频域，理解FFT的核心输出

在开始复杂的应用之前，我们必须先和numpy.fft.fft这个函数成为“老朋友”。很多初学者直接调用fft后，对着输出的一堆复数数组发懵，不知道如何解读。这一步没走稳，后面的所有应用都是空中楼阁。

首先，我们要建立一个清晰的认知：快速傅里叶变换（FFT）计算的是离散傅里叶变换（DFT），它假设我们的信号是周期性的。当我们对一个长度为N的时域信号序列进行FFT时，我们会得到一个同样长度为N的复数数组。这个数组包含了信号的频率成分信息。

关键点在于如何正确解读这个复数数组的索引（下标）与真实频率的对应关系。 这里有一个经典的“顺序”问题。numpy.fft.fft输出的默认顺序是[0, 1, 2, ..., N/2, -N/2+1, ..., -1]，其中索引0对应直流分量（0 Hz），索引1到N/2对应正频率，索引N/2+1到N-1对应负频率。为了绘制出符合人类阅读习惯的频谱图（从0Hz到最高频率），我们需要使用numpy.fft.fftfreq来获取正确的频率轴，或者用numpy.fft.fftshift将零频率分量移动到频谱中心。

让我们通过一个合成信号来直观感受一下。假设我们采样率为1000 Hz，采集了1秒钟的数据，那么我们就有了1000个数据点。我们合成一个包含50Hz和120Hz正弦波的信号，并加上一些随机噪声。

import numpy as np
import matplotlib.pyplot as plt

# 设置参数
fs = 1000  # 采样频率 (Hz)
T = 1.0    # 信号时长 (秒)
N = int(fs * T)  # 采样点数
t = np.linspace(0.0, T, N, endpoint=False) # 时间轴

# 合成信号：50Hz + 120Hz 正弦波，加噪声
signal = 0.7 * np.sin(2.0 * np.pi * 50.0 * t) + 1.0 * np.sin(2.0 * np.pi * 120.0 * t)
noise = 0.5 * np.random.randn(N)  # 高斯白噪声
y = signal + noise

# 执行FFT
Y = np.fft.fft(y)
# 计算频率轴 (单边频谱，只取正频率部分)
freqs = np.fft.fftfreq(N, 1/fs)[:N//2]
# 计算幅度谱 (取绝对值，并乘以2/N进行缩放，恢复真实幅度，直流分量除外)
magnitude = 2.0/N * np.abs(Y[0:N//2])
# 注意：索引0的直流分量不应该乘以2
magnitude[0] = np.abs(Y[0]) / N

# 绘图
fig, axs = plt.subplots(2, 1, figsize=(10, 6))
axs[0].plot(t, y)
axs[0].set_xlabel('时间 [秒]')
axs[0].set_ylabel('幅度')
axs[0].set_title('含噪声的时域信号')
axs[0].grid(True)

axs[1].stem(freqs, magnitude)
axs[1].set_xlabel('频率 [Hz]')
axs[1].set_ylabel('幅度')
axs[1].set_title('单边幅度频谱')
axs[1].set_xlim([0, 200])  # 聚焦在0-200Hz范围
axs[1].grid(True)

plt.tight_layout()
plt.show()

运行这段代码，你会清晰地看到时域波形一片混乱，但频域频谱中，在50Hz和120Hz处出现了尖锐的峰值，这正是我们合成的两个正弦波的频率。噪声则表现为整个频带上的低幅度“基底”。这个例子揭示了傅里叶变换最根本的能力：将混杂在时间中的不同频率成分分离并展示出来。

注意：fftfreq返回的频率值可能是负数（对应负频率），在绘制单边谱时我们通常只取前N//2个点（正频率部分）。幅度需要乘以2/N进行缩放，以反映原始信号中该频率成分的真实振幅（对于实数信号，其频谱是共轭对称的，能量平均分布在正负频率上）。

2. 音频降噪实战：从嘈杂录音中提取清晰人声

音频处理是傅里叶变换的经典战场。设想一个场景：你在会议室用手机录制了一段重要的访谈，但背景有持续的空调嗡鸣声（比如100Hz的低频嗡嗡声）和一些高频嘶嘶声。我们的目标就是滤除这些特定噪声。

思路很直接：将音频信号转换到频域，找到噪声对应的频率区间，将这些频率成分的幅度置零或衰减，然后再转换回时域。这个过程称为频域滤波。这里我们使用scipy库来读写音频文件，因为它比numpy提供了更完整的音频处理生态。

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.fft import fft, ifft, fftfreq
import sounddevice as sd # 用于播放音频，可选

def denoise_audio(input_path, output_path, low_cutoff=80, high_cutoff=8000, noise_freqs=[]):
    """
    对音频进行降噪处理。
    参数：
        input_path: 输入音频文件路径
        output_path: 输出音频文件路径
        low_cutoff: 高通滤波截止频率，低于此频率的成分将被滤除（去除低频噪声）
        high_cutoff: 低通滤波截止频率，高于此频率的成分将被滤除（去除高频嘶嘶声）
        noise_freqs: 需要陷波滤除的特定频率列表，例如[100, 200]表示滤除100Hz和200Hz
    """
    # 读取音频文件
    sample_rate, data = wavfile.read(input_path)
    # 如果音频是双声道，取第一个声道进行处理（或分别处理）
    if len(data.shape) > 1:
        data = data[:, 0].astype(np.float32)
    else:
        data = data.astype(np.float32)
    
    # 归一化到[-1, 1]范围
    data = data / np.max(np.abs(data))
    
    N = len(data)
    # 执行FFT
    Y = fft(data)
    # 获取频率轴
    freqs = fftfreq(N, 1/sample_rate)
    
    # 创建滤波器（在频域操作）
    filter_mask = np.ones(N, dtype=bool)
    
    # 1. 高通滤波：滤除低频噪声（如风声、震动）
    filter_mask &= (np.abs(freqs) > low_cutoff)
    # 2. 低通滤波：滤除高频噪声（如嘶嘶声）
    filter_mask &= (np.abs(freqs) < high_cutoff)
    # 3. 陷波滤波：滤除特定频率的噪声（如50/60Hz电源工频，空调声）
    for nf in noise_freqs:
        bandwidth = 5  # 陷波带宽，单位Hz
        filter_mask &= ((np.abs(freqs) < (nf - bandwidth)) | (np.abs(freqs) > (nf + bandwidth)))
    
    # 应用滤波器：将滤除的频率成分幅度设为零
    Y_filtered = Y.copy()
    Y_filtered[~filter_mask] = 0
    
    # 执行逆FFT，得到处理后的时域信号
    data_filtered = np.real(ifft(Y_filtered))
    # 重新归一化，防止削波
    data_filtered = data_filtered / np.max(np.abs(data_filtered))
    
    # 保存为16位PCM格式的wav文件