如果输入信号特别特别长(比如一段 1 小时的音频),或者信号是实时源源不断进来的(比如直播语音),你就不能等信号全部录完再做一次超大的 FFT。
这就好比你要吃一根无限长的面条,你不能一口气吞下去,必须切成一小段一小段(Block)来吃。
书中介绍的两种“切面条”的方法,就是:
重叠相加法 (Overlap-Add Method)
重叠存储法 (Overlap-Save Method)
核心难点预警:为什么要“重叠”?
在讲具体方法前,必须理解一个物理现象:卷积会让信号变长。
多出来的这部分叫做“拖尾”(Tail)。就像你敲一下钟(输入是短促的),钟声会延续一段时间(输出变长了)。
当我们把长信号切断处理时,每一段都会产生“拖尾”。如果不处理好这个拖尾,拼接的时候声音就会“卡顿”或“失真”。
方法一:重叠相加法 (Overlap-Add)
直观理解:切的时候切断,拼的时候把“拖尾”叠在一起。
根据你提供的图(图 8.5, 8.6)的逻辑:
总结口诀:输入不重叠,输出要叠加。
方法二:重叠存储法 (Overlap-Save)
直观理解:切的时候这就带点重叠,算完后把坏掉的头扔掉,剩下的直接拼。
根据你提供的图(图 8.8)的逻辑:
总结口诀:输入要重叠,输出切两头(扔掉头,保留尾)。
两种方法的对比(帮你理清思路)
| 特性 | 重叠相加法 (Overlap-Add) | 重叠存储法 (Overlap-Save) |
| 输入切分 | 不重叠(干脆利落切断) | 重叠(每段都要包含上一段的尾巴) |
| 补零操作 | 必须补零(防止时域混叠) | 通常不补零(利用循环卷积特性) |
| 输出处理 | 变长了,尾部和下一段头部相加 | 头部是坏的,扔掉头部,剩下直接拼 |
| 计算量 | 两者差不多,都是为了利用 FFT 的高速 | 两者差不多 |
| 书中图示 | 图 8.6:明显的尾部叠加 | 图 8.8:明显的头部打叉(丢弃) |
完整代码
import numpy as np import matplotlib.pyplot as plt import os # 设置字体以支持中文显示 (如果你的系统没有SimHei,可以尝试换成 Microsoft YaHei 或 Arial Unicode MS) plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False def next_power_of_2(x): return 1 if x == 0 else 2**(x - 1).bit_length() def ensure_dir(directory): if not os.path.exists(directory): os.makedirs(directory) # ========================================== # 1. 重叠相加法 (带过程记录) # ========================================== def overlap_add_convolution_debug(x, h, N_fft=None): Lx = len(x) M = len(h) if N_fft is None: N_fft = next_power_of_2(M * 4) L = N_fft - M + 1 H = np.fft.fft(h, N_fft) # 输出总长度 total_len = Lx + M - 1 y = np.zeros(total_len) # 用于记录中间过程以便画图 debug_blocks = [] for start in range(0, Lx, L): # 切分 x_block = x[start : min(start + L, Lx)] # 补零并FFT X = np.fft.fft(x_block, N_fft) # 频域乘法 + IFFT y_block_time = np.fft.ifft(X * H).real # 记录这一段的结果(为了画图) # 我们记录它在总时间轴上的位置 block_full_padded = np.zeros(total_len) end_pos = min(start + N_fft, total_len) valid_len = end_pos - start block_full_padded[start : end_pos] = y_block_time[:valid_len] debug_blocks.append(block_full_padded) # 叠加到总输出 if start + N_fft > total_len: y[start : total_len] += y_block_time[:total_len - start] else: y[start : start + N_fft] += y_block_time return y, debug_blocks, L, M # ========================================== # 2. 重叠存储法 (带过程记录 - 修复版) # ========================================== def overlap_save_convolution_debug(x, h, N_fft=None): Lx = len(x) M = len(h) if N_fft is None: N_fft = next_power_of_2(M * 4) L = N_fft - M + 1 H = np.fft.fft(h, N_fft) overlap = M - 1 # 修复点:为了让输出能够覆盖全长,输入需要补更多的零 # 只要 padded 足够长,循环就能跑完 x_padded = np.concatenate([np.zeros(overlap), x, np.zeros(N_fft * 2)]) y = [] debug_blocks = [] # 记录 (数据, 是否有效掩码) # 修复点:循环的范围必须覆盖到 Lx + M - 1 (输出长度) # 我们以输出索引为基准进行循环 output_len = Lx + M - 1 # range 这里稍微放宽一点,保证最后一段能处理完 for start in range(0, output_len, L): # 取 N_fft 长度,这就自然包含了前一段的 M-1 个重叠点 # 这里的 start 对应的是 x_padded 里的偏移量 x_block = x_padded[start : start + N_fft] # 计算循环卷积 Y_block = np.fft.ifft(np.fft.fft(x_block, N_fft) * H).real # 记录画图数据 debug_blocks.append(Y_block) # 丢弃前 M-1 个点 (Save valid part) y_valid = Y_block[overlap:] y.extend(y_valid) # 裁剪到正确长度 return np.array(y)[:output_len], debug_blocks, overlap # ========================================== # 3. 画图函数 (模仿书中风格) # ========================================== def plot_book_style(x, h, save_dir="result_images"): ensure_dir(save_dir) # 计算一次标准结果作为参考 y_std = np.convolve(x, h, mode='full') # --- 1. 绘制重叠相加法 (Overlap-Add) 原理图 --- y_ola, blocks_ola, L_ola, M = overlap_add_convolution_debug(x, h, N_fft=256) plt.figure(figsize=(12, 8)) # 只画前3个block演示原理,画太多看不清 num_show = min(3, len(blocks_ola)) for i in range(num_show): plt.subplot(num_show + 1, 1, i+1) # 画出当前这一段 plt.plot(blocks_ola[i], 'tab:blue', label=f'段 {i+1} 卷积结果') # 画出重叠区域的提示 if i > 0: # 标出上一段尾巴和这一段头部的重叠区 plt.axvspan((i)*L_ola, (i)*L_ola + (M-1), color='yellow', alpha=0.3, label='重叠相加区域') plt.xlim(0, L_ola * (num_show + 1)) plt.ylabel(f'段 {i+1}') plt.grid(True, linestyle=':') if i == 0: plt.legend(loc='upper right') # 最后画总和 plt.subplot(num_show + 1, 1, num_show + 1) plt.plot(y_ola, 'k', linewidth=1.5, label='最终叠加结果 (y[n])') plt.xlim(0, L_ola * (num_show + 1)) plt.xlabel('时间序列 n') plt.ylabel('总输出') plt.title('图 A: 重叠相加法 (Overlap-Add) - 尾部叠加示意') plt.legend() plt.tight_layout() plt.savefig(os.path.join(save_dir, "overlap_add_demo.png")) print(f"图 A 已保存至: {os.path.join(save_dir, 'overlap_add_demo.png')}") # --- 2. 绘制重叠存储法 (Overlap-Save) 原理图 --- y_ols, blocks_ols, overlap_len = overlap_save_convolution_debug(x, h, N_fft=256) plt.figure(figsize=(12, 8)) N_fft = 256 valid_len = N_fft - overlap_len for i in range(num_show): plt.subplot(num_show + 1, 1, i+1) # 获取当前块的循环卷积结果 current_block = blocks_ols[i] # 1. 画前面的“坏点” (Bad data) - 用红色叉叉或虚线表示 bad_part = current_block[:overlap_len] bad_x = np.arange(overlap_len) plt.plot(bad_x, bad_part, 'r--', alpha=0.5) plt.scatter(bad_x, bad_part, marker='x', color='red', s=20, label='混叠/丢弃部分' if i==0 else "") # 2. 画后面的“好点” (Good data) good_part = current_block[overlap_len:] good_x = np.arange(overlap_len, len(current_block)) plt.plot(good_x, good_part, 'g-', label='有效保留部分' if i==0 else "") plt.title(f'段 {i+1} (循环卷积结果)') plt.ylabel('幅值') plt.grid(True, linestyle=':') if i == 0: plt.legend() # 绘制拼接后的结果 plt.subplot(num_show + 1, 1, num_show + 1) # 为了展示拼接,我们把 valid 部分拼起来画 stitched = [] for i in range(num_show): stitched.extend(blocks_ols[i][overlap_len:]) plt.plot(stitched, 'k', linewidth=1.5, label='拼接结果') plt.title('图 B: 重叠存储法 (Overlap-Save) - 头部丢弃示意') plt.xlabel('相对时间 n') plt.legend() plt.tight_layout() plt.savefig(os.path.join(save_dir, "overlap_save_demo.png")) print(f"图 B 已保存至: {os.path.join(save_dir, 'overlap_save_demo.png')}") # --- 3. 验证正确性并打印 --- print("-" * 30) print(f"标准卷积长度: {len(y_std)}") print(f"重叠相加长度: {len(y_ola)}") print(f"重叠存储长度: {len(y_ols)}") is_ola_ok = np.allclose(y_std, y_ola, atol=1e-8) is_ols_ok = np.allclose(y_std, y_ols, atol=1e-8) print(f"重叠相加法 验证通过? -> {is_ola_ok}") print(f"重叠存储法 验证通过? -> {is_ols_ok}") # ========================================== # Main Execution # ========================================== if __name__ == "__main__": # 生成测试数据 np.random.seed(42) x = np.random.randn(2000) # 输入信号 h = np.random.randn(100) # 滤波器 # 运行并绘图 plot_book_style(x, h)