分段卷积算法

分段卷积算法

实时音频处理中，分段卷积算法尤其重要。由于卷积的时域计算量非常大，而 FFT 乘法相当于时域圆周卷积，因此主流方法是使用FFT计算线性卷积。将圆周卷积转换为线性卷积用到的方法包括重叠相加法和重叠保留法。
然而，目前大多数的卷积算法介绍是基于对输入音频流分段处理，而滤波器不分段处理。假如滤波器非常长，这种方法仍然有较大计算量，影响实时性。
本篇文章实现了对输入信号和滤波器均分段的利用重叠保留法实现的频域卷积算法。
本篇文章是基于 F. Wefers, “Partitioned convolution algorithms for real-time auralization,” 2014 一书中的第5章第一节。

重叠保留法

频域乘积相当于时域圆周卷积。为了利用频域乘积获得时域线性卷积，必须对时域信号进行处理。
在对输入信号分段而滤波器不分段的卷积算法中，常用的适用于重叠保留法的处理为：

1. 将信号 $x(n)$ 分割为与滤波器 $h(n)$ 等长的子段 $x^{(i)}(n)$；设滤波器长度为 $L$，为 2 的整数指数倍；
2. 当前帧输入为前后两子段的组合 $x_i=[x^{(i-1)}(n)x^{(i)}(n)]$；对 $x_i$ 做 $2L$ 点 FFT 得到 $X_i$；
3. 将滤波器 $h(n)$ 后补 $L$ 个零到 $2L$ 点，然后做 $2L$ 点 FFT 得到 $H$;
4. $X_i$ 与 $H$ 相乘，获得 $Y_i$；
5. 对 $Y_i$ 做 $2L$ 点 IFFT，取后 $L$ 点作为当前帧输出 $y_i$。

均匀分割卷积算法

本方法为对输入信号和滤波器均进行均匀分割的卷积算法。

• 先处理滤波器：

  1. 长度为 $N$ 的滤波器 $h(n)$ 均等分割为长度为 $B>0$ 的子滤波器 $h^{(i)}(n)$，子滤波器的数量 $P\in \mathbb{N}$ 满足：$P=\lceil \frac{N}{B}\rceil$
  2. 将每个子滤波器后补 $B$ 个零到 $2B$ 点，然后做 $2B$ 点 FFT 转换到频域，对于实数滤波器，每个子滤波器保留前 $B+1$ 点即可；
  3. 将所有子滤波器频域的前 $B+1$ 点保存为维度为 $[P,B+1]$ 的复数数组 $H$。

• 再处理输入信号：

  1. 申请一个与复数数组 $H$ 维度相同的数组 $X$，初始化为 0，作为输入数据的频域缓存；
  2. 将信号 $x(n)$ 分割为长度为 $B$ 的子段 $x^{(i)}(n)$；
  3. 当前帧输入为前后两子段的组合 $x_i=[x^{(i-1)}(n)x^{(i)}(n)]$；对 $x_i$ 做 $2B$ 点 FFT 得到 $X_i$；
  4. 在频域上，将 $X$ 中的每行向后移动一行，将 $X_i$ 的前 $B+1$ 点存入 $X$ 的第一行中。

• 滤波过程：

  1. 数组 $H$ 与 数组 $X$ 的对应行相乘，其结果按列求和，获得 $Y_i$;
  2. 对 $Y_i$ 做 $2L$ 点 IFFT，取后 $L$ 点作为当前帧输出 $y_i$。

Python 实现

def ir_convolve_freq(x, h, block_size):
	B = block_size  # 分段长度
	N = int(len(x))  # 输入信号总长度，本处未做实时处理
	L = int(len(h))  # 滤波器总长度

	P = int(np.ceil(L / B))  # 滤波器分段数
	h = np.pad(h, (0, P * B - L))  # 对滤波器补零至 P*B
	h = h.reshape((P, B))  # 将滤波器重构为维度为 [P,B] 的数组
	h = np.pad(h, ((0, 0), (0, B)))  # 在滤波器数组的每行后补 B 个 0
	H = np.fft.rfft(h)  # 转换到频域
	
	X = np.zeros_like(H)  # X 为与 H 相同维度的数组，并初始化为 0
	buffer = np.zeros(2 * B)  # 当前输入帧的缓冲区
	y = np.array([])  # 输出
	
	flag = 0  
	while flag < N + B:  
	# flag < N + B 的原因是，当 flag>=N 时，输入x中的数据已用尽，但 buffer 中仍有数据
		buffer[:B] = buffer[B:]  # buffer 的后半移至前半，即为上一帧数据
		buffer[B:] = np.zeros(B)  
		
		if flag + B < N:
			buffer[B:] = x[flag:flag + B] # 新一帧存入 buffer 的后半
		else:
			for i in range(N - flag):  # 对最后一帧单独处理，因为 x 长度可能不是 B 的整数倍
				buffer[B + i] = x[flag + i]
				
		X[1:] = X[:-1]  # X 向后移动一行
		X[0] = np.fft.rfft(buffer)  # 当前帧的 FFT 存入 X 第一行
		
		Y = np.zeros(B + 1, dtype=complex)
		for i in range(P):  # X 与 H 按行相乘，按列求和，获得 Y
			Y += X[i] * H[i]
			
		y = np.concatenate((y, np.fft.irfft(Y)[B:]))  # 对 Y 做 IFFT，取后 B 点作为当前帧输出
		
		flag += B
		
	for j in range(1, P):  # x 和 buffer 中的数据均已用尽，无新输入，但 X 中仍有缓存，继续计算
		Y = np.zeros(B + 1, dtype=complex)
		for i in range(j, P):
			Y += X[i - j] * H[i]
			
		y = np.concatenate((y, np.fft.irfft(Y)[B:]))
		
	return y[:N + L - 1]  # 返回卷积结果