前言
我们知道,如果使用DFT来计算线性卷积,需要分别将M点数据长度的
x
(
n
)
x(n)
x(n)和N点滤波器长度的
h
(
n
)
h(n)
h(n)分别补零至L点以上,其中
L
=
M
+
N
−
1
L=M+N-1
L=M+N−1然而,当M与N相差非常大时,需要补很多零,另外,由于FFT属于批处理算法,需要等到长序列全部输入完毕后才能利用FFT这一快速算法计算DFT,实时性难以保证。因此,可以采用分段卷积的方法来解决这一问题。分段卷积的具体方法,主要包括重叠相加法和重叠保留法。
重叠相加法
理论分析
重叠相加法的核心思想在于,将长序列均匀分段,每段分别与短序列做线性卷积,而后将卷积结果重叠部分叠加起来。结果可以用分段函数来表示,对给定的n,可以相应求出结果值y(n)。
具体分析如下:
设长序列
x
(
n
)
x(n)
x(n)长到无限长,滤波器
h
(
n
)
h(n)
h(n)的长度为N,将
x
(
n
)
x(n)
x(n)按照每M点一部分进行划分,其中设M>N,每段用
x
k
(
n
)
x_k(n)
xk(n)表示,则有
x
k
(
n
)
=
{
x
(
n
+
k
M
)
,0
≤
n
≤
M-1, k=0,1,2
…
0
,others
x_k(n)=\begin{cases} x(n+kM) & \text{,0 $\le$ n $\le$ M-1, k=0,1,2…} \\ 0 & \text{,others} \end{cases}
xk(n)={x(n+kM)0,0 ≤ n ≤ M-1, k=0,1,2…,others
反过来,
x
(
n
)
x(n)
x(n)可以用
x
k
(
n
)
x_k(n)
xk(n)表示为
x
(
n
)
=
∑
k
=
0
∞
x
k
(
n
−
k
M
)
x(n)=\sum_{k=0}^{\infty}x_k(n-kM)
x(n)=k=0∑∞xk(n−kM)
设卷积结果为
y
(
n
)
y(n)
y(n),有
y
(
n
)
=
x
(
n
)
∗
h
(
n
)
=
∑
k
=
0
∞
x
k
(
n
−
k
M
)
∗
h
(
n
)
=
∑
k
=
0
∞
y
k
(
n
−
k
M
)
y(n)=x(n)*h(n)\\ =\sum_{k=0}^{\infty}x_k(n-kM)*h(n)\\ =\sum_{k=0}^{\infty}y_k(n-kM)
y(n)=x(n)∗h(n)=k=0∑∞xk(n−kM)∗h(n)=k=0∑∞yk(n−kM)
其中,
y
k
(
n
)
=
x
k
(
n
)
∗
h
(
n
)
y_k(n)=x_k(n)*h(n)
yk(n)=xk(n)∗h(n),
0
≤
n
≤
N
+
(
M
−
1
)
−
1
0\le n \le N+(M-1)-1
0≤n≤N+(M−1)−1
则
y
k
(
n
−
k
M
)
y_k(n-kM)
yk(n−kM)的非零部分的自变量取值范围为
0
≤
n
−
k
M
≤
N
+
(
M
−
1
)
−
1
0\le n-kM \le N+(M-1)-1
0≤n−kM≤N+(M−1)−1
即
k
M
≤
n
≤
N
+
(
k
+
1
)
M
−
2
kM\le n \le N+(k+1)M-2
kM≤n≤N+(k+1)M−2
下面讨论重叠部分:
由于
y
k
(
n
−
k
M
)
y_k(n-kM)
yk(n−kM)非零部分的自变量取值范围为
k
M
≤
n
≤
N
+
(
k
+
1
)
M
−
2
kM\le n \le N+(k+1)M-2
kM≤n≤N+(k+1)M−2
则
y
k
+
1
(
n
−
(
k
+
1
)
M
)
y_{k+1}(n-(k+1)M)
yk+1(n−(k+1)M)非零部分的自变量取值范围为
(
k
+
1
)
M
≤
n
≤
N
+
(
k
+
2
)
M
−
2
(k+1)M\le n \le N+(k+2)M-2
(k+1)M≤n≤N+(k+2)M−2
则可以得出,重叠区间为
(
k
+
1
)
M
≤
n
≤
N
+
(
k
+
1
)
M
−
2
(k+1)M\le n \le N+(k+1)M-2
(k+1)M≤n≤N+(k+1)M−2
可见,重叠点数为
N
+
(
k
+
1
)
M
−
2
−
(
k
+
1
)
M
+
1
=
N
−
1
N+(k+1)M-2 - (k+1)M +1 = N-1
N+(k+1)M−2−(k+1)M+1=N−1
所以,只要将结果分段,未重叠部分直接得到最终结果,重叠部分相加后得到结果,就可逐段得到最终卷积结果。
代码验证
下面举个例子:
设
h
(
n
)
=
2
,
4
,
1
,
6
h(n)={2,4,1,6}
h(n)=2,4,1,6,
x
(
n
)
=
n
2
−
4
n
+
1
,
0
≤
n
≤
20
x(n)=n^2-4n+1, 0\le n \le 20
x(n)=n2−4n+1,0≤n≤20,取分段数M为6.
这里使用MATLAB,将采用重叠相加法得到的最终卷积结果与直接使滤波器与长序列线性卷积得到的结果相比较,从而验证算法的正确性。
clc;close all;clear;
hn = [2,4,1,6];
n = 0 : 20;
xn = n.^2-4*n+1;
M = 6;
N = length(hn);
La = length(xn);
L = M + N - 1;
km = floor(La/M);
break_mtx = zeros(km, M);
end_set = zeros(1, La-km*M);
y_mtx = zeros(km, L);
end_y_set = zeros(1, La-km*M+N-1);
for k = 0 : km
if (k~=km)
break_mtx(k+1, :) = xn(1+k*M:(k+1)*M);
y_mtx(k+1, :) = conv(break_mtx(k+1, :), hn);
else
end_set = xn(1+k*M:La);
end_y_set = conv(end_set, hn);
end
end
y_final = zeros(1, La+N-1);
y_mtx = [y_mtx;[end_y_set,zeros(1,L-(La-km*M+N-1))]];
for k = 0 : km
if(k==0)
y_final(1+k*M:(k+1)*M) = y_mtx(k+1,1:M);
y_final((k+1)*M+1:L+k*M) = y_mtx(k+1, M+1:L)+y_mtx(k+2,1:N-1);
elseif(k~=km)
y_final(N+k*M:(k+1)*M) = y_mtx(k+1,N:M);
y_final((k+1)*M+1:L+k*M) = y_mtx(k+1, M+1:L)+y_mtx(k+2,1:N-1);
else
y_final(N+k*M:end) = y_mtx(k+1,N:La-km*M+N-1);
end
end
y_test = conv(xn, hn);
CorrectNum = sum(y_final-y_test==0);
if(CorrectNum==length(y_test)) disp('SUCCESSFUL');end
结果:
SUCCESSFUL
重叠保留法
理论分析
重叠保留法的要点在于,利用循环卷积,通过初始移位使得重叠点为不需要的点,舍弃重叠点拼接得到最终结果。
设滤波器 h ( n ) h(n) h(n)的长度为N,将长序列 x ( n ) x(n) x(n)按照每M个点为一段进行划分。设M>N,我们知道,循环卷积的前N-1个点是重叠点。我们希望重叠点无用被弃置,所以在开始分段之前,先将长序列 x ( n ) x(n) x(n)整体右移N-1个点,空位补零,然后再分段,每段与 h ( n ) h(n) h(n)进行循环卷积。这里要注意,此处的段,不是直接划分,而是在每段的前N-1个点,取上一个段的后N-1个点。将每一个卷积后的段的前N-1个点舍弃,保留并拼接剩下的M-N+1个点,即可得到最终的卷积结果。
代码验证
仍以上面的例子,进行MATLAB验证:
clc;clear;close all;
hn = [2,4,1,6];
n = 0:20;
xn = n.^2-4*n+1;
M = 6;
N = length(hn);
La = length(xn);
km = ceil(La/(M-N+1));
remains = km*(M-N+1)-La;
km = km+1;
x_mtx = zeros(km, M);
y_mtxrow = zeros(1, km*(M-N+1));
for k = 1 : km
if(k==1)
x_mtx(k,:) = [zeros(1,N-1), xn(1+(k-1)*(M-N+1):k*(M-N+1))];
elseif(k==km)
x_mtx(k,:) = [xn((k-1)*(M-N+1)-(N-2):(k-1)*(M-N+1)), zeros(1, M-N+1)];
elseif(k==km-1)
x_mtx(k,:) = [xn((k-1)*(M-N+1)-(N-2):(k-1)*(M-N+1)), xn(1+(k-1)*(M-N+1):La), zeros(1, remains)];
else
x_mtx(k,:) = [xn((k-1)*(M-N+1)-(N-2):(k-1)*(M-N+1)), xn(1+(k-1)*(M-N+1):k*(M-N+1))];
end
temp = ifft(fft(x_mtx(k,:), M).*fft(hn, M), M);
y_mtxrow(1+(k-1)*(M-N+1):k*(M-N+1)) = temp(N:M);
end
y_final = y_mtxrow(1:La+N-1);
y_test = conv(xn, hn);
CorrectNum = sum(y_final-y_test<1e-4);%fft引入微小误差
if(CorrectNum==length(y_test)) disp('SUCCESSFUL');end
结果:
SUCCESSFUL