使用FFT进行快速FIR滤波

摘自<Understanding Digital Signal Processing>第三版，13.10 Fast FIR Filtering Using the FFT一节

基于的理论：频域上的乘积等效于时域上的卷积。

基本的计算流程，如下图所示。

将输入信号x(n)x(n)x(n)和滤波器参数h(k)h(k)h(k)分别进行FFT，得到X(m)X(m)X(m)和H(m)H(m)H(m)，在频域上进行乘积，然后，进行IFFT。

对于Q−tapQ-tapQ−tap的FIR滤波器，其标准的卷积方程为
y(n)=∑k=0Q−1h(k)x(n−k)=h(k)∗x(n)y(n)=\sum ^{Q-1} _{k=0} {h(k)x(n-k)}=h(k)*x(n) y(n)=k=0∑Q−1h(k)x(n−k)=h(k)∗x(n)
假设h(k)h(k)h(k)的长度为QQQ，x(n)x(n)x(n)的长度为PPP，则最终输出的y(n)y(n)y(n)的长度为L=Q+P−1L=Q+P-1L=Q+P−1.

为了使快速卷积技术能得到有效的结果，前向和逆FFT的尺寸必须等于或大于LLL，因此，N−pointN-pointN−point的FFT尺寸为N≥LN\ge LN≥L，并对h(k)h(k)h(k)和x(n)x(n)x(n)进行zero−padzero-padzero−pad，使其长度等于NNN。所需要的输出y(n)y(n)y(n)为逆FFT的前LLL个样本的实数部分。关于滤波器参数的FFT只需计算一次并存储下来即可。

当输出序列长度太大，必须进行分段处理时，有两种处理方法：overlap−and−saveoverlap-and-saveoverlap−and−save和overlap−and−addoverlap-and-addoverlap−and−add。

Overlap-and-save

具体流程，如下图所示。

对于给定的FIR滤波器的脉冲响应h(k)h(k)h(k)的长度为QQQ，输入序列x(n)x(n)x(n)的长度为PPP，具体步骤：

选择FFT的尺寸为NNN，其中NNN为2的指数，约等于4倍QQQ；
在h(k)h(k)h(k)的末端增加(N−Q)(N-Q)(N−Q)个零值样本，进行N−pointN-pointN−point的FFT，生成复序列H(m)H(m)H(m);
计算整数M，M=N−(Q−1)M=N-(Q-1)M=N−(Q−1)；
在x(n)x(n)x(n)的前MMM样本之前插入(Q−1)(Q-1)(Q−1)个零值样本，创建N−pointN-pointN−point的FFT的输入序列x1(n)x_1(n)x1(n)；
对x1(n)x_1(n)x1(n)进行N−pointN-pointN−point的FFT计算，并乘以H(m)H(m)H(m)，然后，对乘积进行N−pointN-pointN−point的逆FFT。去掉逆FFT结果的前(Q−1)(Q-1)(Q−1)个样本，生成M−pointM-pointM−point的输出y1(n)y_1(n)y1(n)；
将x1(n)x_1(n)x1(n)的前(Q−1)(Q-1)(Q−1)个零值样本增加到第二个原始长度为M的x(n)x(n)x(n)的开头部分，创建第二个N−pointN-pointN−point的FFT输入序列x2(n)x_2(n)x2(n)；
对x2(n)x_2(n)x2(n)进行N−pointN-pointN−point的FFT计算，并乘以H(m)H(m)H(m)，然后，进行逆FFT，去除前(Q−1)(Q-1)(Q−1)个样本，生成第二个M点的输出数据y2(n)y_2(n)y2(n)；
重复步骤6-7，直至处理完整个输入序列。
将生成的y1(n),y2(n),y3(n),....y_1(n),y_2(n),y_3(n),....y1(n),y2(n),y3(n),....进行串联合并，生成最终的线性卷积滤波器的输出序列y(n)y(n)y(n)。
可以试验不同的NNN值，看是否存在一个最优的NNN使得你的硬件和软件实现的计算载荷最小。不过，在任何情况下，NNN必须不小于(M+Q−1)(M+Q-1)(M+Q−1)。

Overlap-and-add

具体流程，如下图所示。

对于这种方法，输入序列x(n)x(n)x(n)被分段成长度为MMM的数据片段，数据的overlap发生在逆FFT计算。对于给定的FIR滤波器的脉冲响应h(k)h(k)h(k)的长度为QQQ，输入序列x(n)x(n)x(n)的长度为PPP，具体步骤：

选择FFT尺寸为NNN，其中NNN为2的指数，约等于2倍的Q；
在h(k)h(k)h(k)的末端增加(N−Q)(N-Q)(N−Q)个零值样本，进行N−pointN-pointN−point的FFT，生成复序列H(m)H(m)H(m);
计算整数M，M=N−(Q−1)M=N-(Q-1)M=N−(Q−1)；
在原始序列x1(n)x_1(n)x1(n)的前MMM样本末端增加(Q−1)(Q-1)(Q−1)个零值样本，创建N−pointN-pointN−point的FFT的输入序列x1(n)x_1(n)x1(n)；
对x1(n)x_1(n)x1(n)进行N−pointN-pointN−point的FFT计算，并乘以H(m)H(m)H(m)，然后，对乘积进行N−pointN-pointN−point的逆FFT。并保留前MMM个样本，生成MMM样本的输出数据y1(n)y_1(n)y1(n)；
在第二个原始序列x2(n)x_2(n)x2(n)的前MMM样本末端增加(Q−1)(Q-1)(Q−1)个零值样本，创建N−pointN-pointN−point的FFT的输入序列x2(n)x_2(n)x2(n)；
对第二个N−pointN-pointN−point的输入序列进行FFT计算，并乘以H(m)H(m)H(m)，然后，对乘积进行N−pointN-pointN−point的逆FFT。将上一步的逆FFT的最后Q−1Q-1Q−1个样本增加到当前逆FFT序列的开头的Q−1Q-1Q−1个样本上。保留Q−1Q-1Q−1个样本相加过程的结果序列中的开头的MMM个样本作为输出数据y2(n)y_2(n)y2(n);
重复步骤6-7，直至处理完整个数据序列。
将生成的y1(n),y2(n),y3(n),....y_1(n),y_2(n),y_3(n),....y1(n),y2(n),y3(n),....进行串联合并，生成最终的线性卷积滤波器的输出序列y(n)y(n)y(n)。
可以试验不同的NNN值，看是否存在一个最优的NNN使得你的硬件和软件实现的计算载荷最小。不过，在任何情况下，NNN必须不小于$