用Python+NumPy可视化理解FS/FT/DFT/FFT的物理意义

信号处理工程师最常遇到的困惑之一，就是那些看起来相似的缩写：FS、FT、DFS、DTFT、DFT、FFT。它们都带着"傅里叶"的基因，却各自有着不同的数学表达和适用场景。本文将通过Python代码和可视化手段，带你从波形和频谱的视角重新认识这些概念。

1. 从周期信号到傅里叶级数(FS)

让我们从一个简单的周期方波开始。方波是理解傅里叶分析的理想起点，因为它包含丰富的谐波成分。

import numpy as np
import matplotlib.pyplot as plt

def square_wave(t, T):
    """生成周期为T的方波信号"""
    return np.where(np.sin(2*np.pi*t/T) > 0, 1, -1)

T = 2  # 周期(秒)
t = np.linspace(0, 3*T, 1000)  # 3个周期的时间序列
signal = square_wave(t, T)

plt.figure(figsize=(10,4))
plt.plot(t, signal)
plt.title("周期方波信号(时域)")
plt.xlabel("时间(s)")
plt.ylabel("幅值")
plt.grid()
plt.show()

傅里叶级数告诉我们，任何周期信号都可以表示为正弦波的叠加。对于方波，其傅里叶级数展开为：

$$ f(t) = \frac{4}{\pi}\sum_{k=1,3,5...}^{\infty}\frac{1}{k}\sin(2\pi k f_0 t) $$

其中$f_0=1/T$是基频。让我们用Python实现这个级数的前N项：

def fourier_series(t, T, N):
    """计算方波的傅里叶级数前N项和"""
    result = np.zeros_like(t)
    f0 = 1/T
    for k in range(1, N+1, 2):  # 只考虑奇数次谐波
        result += (4/(np.pi*k)) * np.sin(2*np.pi*k*f0*t)
    return result

# 绘制不同项数下的近似效果
plt.figure(figsize=(12,8))
for i, N in enumerate([1, 3, 5, 10, 20, 50], 1):
    approx = fourier_series(t, T, N)
    plt.subplot(3,2,i)
    plt.plot(t, signal, label='原始信号')
    plt.plot(t, approx, label=f'{N}项近似')
    plt.legend()
    plt.grid()
plt.tight_layout()
plt.show()

随着项数的增加，我们看到方波的形状逐渐清晰。这就是FS的物理意义： 用不同频率的正弦波合成任意周期信号 。

2. 非周期信号的傅里叶变换(FT)

当信号周期T趋近于无穷大时，周期信号就变成了非周期信号，傅里叶级数也演变为傅里叶变换。让我们通过高斯脉冲信号来理解这一转变。

def gaussian_pulse(t, sigma):
    """生成高斯脉冲信号"""
    return np.exp(-t**2/(2*sigma**2))

sigma = 0.5  # 脉冲宽度参数
t_ft = np.linspace(-5, 5, 1000)
signal_ft = gaussian_pulse(t_ft, sigma)

plt.figure(figsize=(10,4))
plt.plot(t_ft, signal_ft)
plt.title("高斯脉冲信号(时域)")
plt.xlabel("时间(s)")
plt.ylabel("幅值")
plt.grid()
plt.show()

高斯脉冲的傅里叶变换解析解为：

$$ F(\omega) = \sigma\sqrt{2\pi}e^{-\sigma^2\omega^2/2} $$

我们同时用NumPy的FFT计算数值解：

# 计算傅里叶变换
dt = t_ft[1] - t_ft[0]  # 采样间隔
freq = np.fft.fftfreq(len(t_ft), dt)
fft_result = np.fft.fft(signal_ft) * dt  # 乘以dt近似连续FT

# 只取正频率部分
positive_freq = freq > 0
freq_pos = freq[positive_freq]
fft_pos = fft_result[positive_freq]

# 解析解
analytic = sigma * np.sqrt(2*np.pi) * np.exp(-sigma**2*(2*np.pi*freq_pos)**2/2)

plt.figure(figsize=(10,4))
plt.plot(2*np.pi*freq_pos, np.abs(fft_pos), label='FFT数值解')
plt.plot(2*np.pi*freq_pos, analytic, '--', label='解析解')
plt.title("高斯脉冲的频谱")
plt.xlabel("角频率(rad/s)")
plt.ylabel("幅值")
plt.legend()
plt.grid()
plt.show()

FT的物理意义在于： 非周期信号可以看作无限多个连续频率的正弦波的叠加 ，其频谱是连续的。

3. 离散时间傅里叶变换(DTFT)与离散傅里叶变换(DFT)

在实际应用中，我们处理的是离散采样信号。DTFT和DFT就是为离散信号设计的频域分析工具。

考虑一个简单的离散指数衰减信号：

n = np.arange(0, 50)  # 离散时间点
a = 0.9  # 衰减系数
x = a**n  # 离散指数衰减信号

plt.figure(figsize=(10,4))
plt.stem(n, x)
plt.title("离散指数衰减信号")
plt.xlabel("采样点n")
plt.ylabel("幅值")
plt.grid()
plt.show()

DTFT的数学定义为：

$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n]e^{-j\omega n} $$

对于我们的有限长信号，可以计算其DTFT：

def dtft(x, omega):
    """计算离散时间傅里叶变换"""
    n = np.arange(len(x))
    return np.sum(x * np.exp(-1j * omega * n))

omega = np.linspace(-np.pi, np.pi, 500)
X = np.array([dtft(x, w) for w in omega])

plt.figure(figsize=(12,4))
plt.subplot(121)
plt.plot(omega, np.abs(X))
plt.title("DTFT幅度谱")
plt.xlabel("数字频率ω")
plt.ylabel("|X(e^jω)|")
plt.grid()

plt.subplot(122)
plt.plot(omega, np.angle(X))
plt.title("DTFT相位谱")
plt.xlabel("数字频率ω")
plt.ylabel("∠X(e^jω)")
plt.grid()
plt.tight_layout()
plt.show()

DTFT的频谱是连续的，但在计算机中我们只能计算有限个频率点，这就是DFT：

N = len(x)  # DFT点数
X_dft = np.fft.fft(x, N)  # 计算N点DFT
omega_dft = 2*np.pi/N * np.arange(N)  # 对应的数字频率

plt.figure(figsize=(12,4))
plt.subplot(121)
plt.stem(omega_dft, np.abs(X_dft))
plt.title(f"{N}点DFT幅度谱")
plt.xlabel("数字频率ω")
plt.ylabel("|X[k]|")
plt.grid()

plt.subplot(122)
plt.stem(omega_dft, np.angle(X_dft))
plt.title(f"{N}点DFT相位谱")
plt.xlabel("数字频率ω")
plt.ylabel("∠X[k]")
plt.grid()
plt.tight_layout()
plt.show()

DFT可以理解为 对DTFT频谱的均匀采样 。当DFT点数增加时，DFT结果会更接近DTFT：

N_large = 500  # 更大的DFT点数
X_dft_large = np.fft.fft(x, N_large)
omega_dft_large = 2*np.pi/N_large * np.arange(N_large)

plt.figure(figsize=(10,4))
plt.plot(omega, np.abs(X), label='DTFT')
plt.stem(omega_dft_large, np.abs(X_dft_large), 'r', markerfmt='ro', basefmt=" ", label=f'{N_large}点DFT')
plt.title("DFT对DTFT的采样效果")
plt.xlabel("数字频率ω")
plt.ylabel("幅值")
plt.legend()
plt.grid()
plt.show()

4. 快速傅里叶变换(FFT)的工程实践

FFT是DFT的高效算法实现，在实际工程中应用广泛。让我们通过一个多频信号的处理案例来理解FFT的优势。

# 生成包含多个频率成分的信号
fs = 1000  # 采样率(Hz)
T = 1.0/fs
L = 1000  # 信号长度
t_fft = np.arange(L)*T

# 信号包含50Hz、120Hz和随机噪声
f1, f2 = 50, 120
signal_fft = 0.7*np.sin(2*np.pi*f1*t_fft) + np.sin(2*np.pi*f2*t_fft)
signal_fft += 0.1*np.random.randn(len(t_fft))  # 添加噪声

plt.figure(figsize=(10,4))
plt.plot(t_fft[:100], signal_fft[:100])  # 只显示前100个点
plt.title("多频信号(时域)")
plt.xlabel("时间(s)")
plt.ylabel("幅值")
plt.grid()
plt.show()

使用FFT分析其频谱：

# 计算FFT
fft_result = np.fft.fft(signal_fft)
freq = np.fft.fftfreq(L, T)[:L//2]  # 只取正频率部分

# 计算幅度谱
magnitude = 2/L * np.abs(fft_result[:L//2])

plt.figure(figsize=(10,4))
plt.plot(freq, magnitude)
plt.title("信号的单边幅度谱")
plt.xlabel("频率(Hz)")
plt.ylabel("幅值")
plt.grid()
plt.show()

从频谱中可以清晰地看到50Hz和120Hz两个频率成分。FFT的高效性使得实时频谱分析成为可能：

import time

# 比较DFT和FFT的计算时间
N = len(signal_fft)

# DFT直接实现
def dft(x):
    N = len(x)
    n = np.arange(N)
    k = n.reshape((N,1))
    return np.sum(x * np.exp(-2j*np.pi*k*n/N), axis=1)

start = time.time()
dft_result = dft(signal_fft)
dft_time = time.time() - start

# FFT计算
start = time.time()
fft_result = np.fft.fft(signal_fft)
fft_time = time.time() - start

print(f"DFT计算时间: {dft_time:.4f}秒")
print(f"FFT计算时间: {fft_time:.6f}秒")
print(f"速度提升倍数: {dft_time/fft_time:.1f}倍")

对于N=1000点的变换，FFT通常比直接计算DFT快100倍以上。这正是FFT在工程实践中不可或缺的原因。

5. 频谱泄漏与窗函数应用

在实际信号处理中，我们常常需要处理有限长度的信号段，这会导致频谱泄漏现象。让我们通过一个例子来理解这个问题及其解决方案。

考虑一个单一频率的正弦信号：

fs = 1000  # 采样率
T = 1/fs
f = 100.4  # 不是整数倍的频率
t_leak = np.arange(0, 0.1, T)  # 0.1秒时长
signal_leak = np.sin(2*np.pi*f*t_leak)

plt.figure(figsize=(10,4))
plt.plot(t_leak, signal_leak)
plt.title("100.4Hz正弦信号")
plt.xlabel("时间(s)")
plt.ylabel("幅值")
plt.grid()
plt.show()

计算其FFT频谱：

N = len(signal_leak)
fft_leak = np.fft.fft(signal_leak)
freq_leak = np.fft.fftfreq(N, T)[:N//2]
magnitude_leak = 2/N * np.abs(fft_leak[:N//2])

plt.figure(figsize=(10,4))
plt.plot(freq_leak, magnitude_leak)
plt.title("频谱泄漏现象")
plt.xlabel("频率(Hz)")
plt.ylabel("幅值")
plt.grid()
plt.show()

可以看到，虽然信号只有一个频率成分，但频谱却散布在很宽的频率范围内。这是因为信号截断相当于乘以矩形窗，导致频谱卷积。

为了减轻泄漏，我们可以使用窗函数：

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

# 计算窗化后的FFT
fft_windowed = np.fft.fft(signal_windowed)
magnitude_windowed = 2/np.sum(window) * np.abs(fft_windowed[:N//2])

plt.figure(figsize=(10,4))
plt.plot(freq_leak, magnitude_leak, label='无窗')
plt.plot(freq_leak, magnitude_windowed, label='汉宁窗')
plt.title("窗函数对频谱泄漏的影响")
plt.xlabel("频率(Hz)")
plt.ylabel("幅值")
plt.legend()
plt.grid()
plt.show()

窗函数虽然不能完全消除泄漏，但可以显著减少频谱扩散，提高频率分辨率。常用的窗函数还有汉明窗、布莱克曼窗等，各有不同的主瓣宽度和旁瓣衰减特性。

6. 从理论到实践：音频频谱分析案例

让我们用一个实际的音频处理案例来综合应用前面学到的知识。我们将分析一段包含多个频率成分的合成音频信号。

from scipy.io import wavfile

# 生成多频音频信号
fs_audio = 44100  # 音频采样率
duration = 2.0  # 持续时间(秒)
t_audio = np.arange(0, duration, 1/fs_audio)

# 三个频率成分
freqs = [440, 880, 1320]  # A4, A5, E6
audio_signal = sum(np.sin(2*np.pi*f*t_audio) for f in freqs)
audio_signal /= np.max(np.abs(audio_signal))  # 归一化

# 保存为WAV文件(用于后续分析)
wavfile.write('multi_tone.wav', fs_audio, (audio_signal*32767).astype(np.int16))

# 绘制前50ms的波形
plt.figure(figsize=(10,4))
plt.plot(t_audio[:2200], audio_signal[:2200])  # 2200 samples ≈ 50ms
plt.title("多频音频信号(前50ms)")
plt.xlabel("时间(s)")
plt.ylabel("幅值")
plt.grid()
plt.show()

现在我们来分析这个音频信号的频谱特性：

# 计算短时傅里叶变换(STFT)
n_fft = 2048  # FFT点数
frequencies, times, stft = plt.mlab.specgram(
    audio_signal, 
    NFFT=n_fft, 
    Fs=fs_audio,
    noverlap=n_fft//2,
    window=np.hanning(n_fft)
)

# 绘制频谱图
plt.figure(figsize=(12,6))
plt.pcolormesh(times, frequencies, 10*np.log10(stft), shading='auto')
plt.colorbar(label='强度(dB)')
plt.title("音频信号的频谱图")
plt.xlabel("时间(s)")
plt.ylabel("频率(Hz)")
plt.ylim(0, 2000)
plt.show()

从频谱图中可以清晰地看到三个稳定的频率成分。这种时频分析在音乐处理、语音识别等领域应用广泛。

7. 傅里叶变换在图像处理中的应用

傅里叶分析不仅适用于一维信号，也可以扩展到二维图像处理。让我们看看图像的频域表示。

from scipy import misc

# 加载示例图像
image = misc.ascent()
plt.figure(figsize=(6,6))
plt.imshow(image, cmap='gray')
plt.title("原始图像")
plt.axis('off')
plt.show()

# 计算二维傅里叶变换
fft_image = np.fft.fft2(image)
fft_shifted = np.fft.fftshift(fft_image)  # 将低频移到中心
magnitude = 20*np.log(np.abs(fft_shifted))

plt.figure(figsize=(12,6))
plt.subplot(121)
plt.imshow(np.abs(fft_shifted), cmap='gray')
plt.title("频域幅度谱")
plt.axis('off')

plt.subplot(122)
plt.imshow(np.angle(fft_shifted), cmap='gray')
plt.title("频域相位谱")
plt.axis('off')
plt.show()

图像的低频成分对应于平滑区域，高频成分对应于边缘和细节。我们可以通过频域滤波来增强或抑制特定特征：

# 创建低通滤波器
rows, cols = image.shape
crow, ccol = rows//2, cols//2
mask = np.zeros((rows, cols), np.uint8)
r = 30  # 滤波半径
mask[crow-r:crow+r, ccol-r:ccol+r] = 1

# 应用滤波
filtered = fft_shifted * mask
filtered = np.fft.ifftshift(filtered)
reconstructed = np.fft.ifft2(filtered).real

plt.figure(figsize=(12,6))
plt.subplot(121)
plt.imshow(mask, cmap='gray')
plt.title("低通滤波器")
plt.axis('off')

plt.subplot(122)
plt.imshow(reconstructed, cmap='gray')
plt.title("滤波后图像")
plt.axis('off')
plt.show()

这种频域处理技术在图像压缩、去噪、增强等方面有广泛应用。