实验3 利用FFT计算线性卷积

实验3 利用FFT计算线性卷积
一、实验目的
1.掌握利用FFT计算线性卷积的原理及具体实现方法。

2.加深理解重叠相加法和重叠保留法。

3.考察利用FFT计算线性卷积各种方法的适用围。

二、实验设备与环境
计算机、matlab软件环境。

三、实验基础理论
1.线型卷积和圆周卷积
设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为
的长度为L+M-1。

x(n)和h(n)的N点圆周卷积为
圆周卷积与线性卷积相等而不产生交叠的必要条件为
圆周卷积定理：根据DFT的性质，x(n)和h(n)的N点圆周卷积的DFT 等于它们DFT的乘积
2.快速卷积
快速卷积算法用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。

可以将快速卷积的步骤归纳如下：
(1)为了使线性卷积可以用圆周卷积来计算，必须选择
;同时为了能使用基-2FFT完成卷积运算，要求。

采用补零的办法使x(n)和h(n)的长度均为N。

(2)计算x(n)和h(n)的N点FFT
(3)组成乘积
(4)利用IFFT计算Y(K)的IDFT，得到线性卷积y(n)
3.分段卷积
我们考察单位取样响应为h(n)的线性系统，输入为x(n),输出为
y(n)，则
当输入序列x(n)极长时，如果要等x(n)全部集齐时再开始进行卷积，会使输出相对输入有较大的延时，再者如果序列太长，需要大量存贮单元。

为此我们把x(n)分段，分别求出每段的卷积，合在一起得到最后总的输出。

这种方法称为分段卷积分段卷积可细分为重叠保留法和重叠相加法。

重叠保留法：设x(n)的长度为，h(n)的长度为M。

我们把序列x(n)分成多段N点序列(n),每段与前一段重叠M-1个样本。

由于第一段没有前一段保留信号，为了修正，我们在第一个输入端前面填充M-1个零。

计算每一段与h(n)的圆周卷积，则其每段卷积结果的前M-1个样本不等于线性卷积值，不是正确的样本值。

所以我们将每段卷积结果的前M-1个样本社区，只保留后面的N-M+1个正确输出样本，把这些输出样本合起来，得到总输出。

利用FFT实现重叠保留发的步骤如下：
(1)在x(n)前面填充M-1个零，扩大以后的序列为
(2)将x(n)分为若干N点子段，设L=N-M+1为每一段的有效数据
长度，则第i段(n)()的数据为
(3)计算每一段与h(n)的N点圆周卷积，利用FFT计算圆周卷积：
(4)舍去每一段卷积结果的前M-1个样本，连接剩下样本，得到卷
积结果y(n)。

重叠相加法：设h(n)长度为M，将信号x(n)分解成长为L的子段，建议L选择与的M数量级相同，以(n)表示没每段信号，则
每一段卷积的长度为L+M-1,所以在做求和时，相邻两段序列有M-1个样本重叠，即前一段的最后M-1个样本喝下一段的前M-1个序列重叠，这个重叠部分相加，再与不重叠部分共同组成输出y(n)。

利用FFT实现重叠保留法的步骤如下：
(1)将x(n)分为若干L点子段(n)。

(2)计算每一段与h(n)的卷积，根据快速卷积算法利用FFT 计算卷积。

(3)将各段相加，得到输出y(n)
四、实验容
假设要计算序列x（n）=u（n）-u（n-L），0<=n<=L和h(n)=cos(0.2n),0<=n<=M的线性卷积，完成以下实验容
1.设L=M，根据线性卷积的表达式和快速卷积的原理，分别编程实现计算两个序列线性卷积的方法，比较当序列长度分别为8,16,32,64,256,512,1024时，两种方法计算线性卷积所需时间。

实验代码及结果：
程序Exp3_1代码
L = input('L=');
M = input('M=');
xn = linspace(1, 1, L);
hn = cos(0.2*pi*[0:M-1]);
yn1 = zeros(1, L+M-1);
tic
for k = 0:L+M-1
for n = max([0 k-M+1]):min([k L-1])
yn1(k+1) = yn1(k+1)+xn(n+1)*hn(k-n+1);
end
end
toc
tic
Xk = fft(xn, L+M-1);
Hk = fft(hn, L+M-1);
Yk = Xk.*Hk;
yn2 = ifft(Yk);
toc
>> Exp3_1
L=8
M=8
Elapsed time is 0.000136 seconds.
Elapsed time is 0.000033 seconds.
>> Exp3_1
M=16
Elapsed time is 0.000205 seconds. Elapsed time is 0.000030 seconds.
>> Exp3_1
L=32
M=32
Elapsed time is 0.000419 seconds. Elapsed time is 0.000056 seconds.
>> Exp3_1
L=64
M=64
Elapsed time is 0.000911 seconds. Elapsed time is 0.000117 seconds.
>> Exp3_1
L=256
M=256
Elapsed time is 0.005773 seconds. Elapsed time is 0.000261 seconds.
>> Exp3_1
L=512
Elapsed time is 0.017060 seconds.
Elapsed time is 0.000190 seconds.
>> Exp3_1
L=1024
M=1024
Elapsed time is 0.055688 seconds.
Elapsed time is 0.000524 seconds.
由此可见，线性卷积慢得多
2.当L=2048且M=256时，比较直接计算线性卷积和快速卷积所需时间，进一步考察当L=4096且M=256时两种算法所需的时间。

实验代码及结果：
>> Exp3_1
L=2048
M=256
Elapsed time is 0.035421 seconds.
Elapsed time is 0.000448 seconds.
由此可见，快速卷积快
>> Exp3_1
M=256
Elapsed time is 0.069497 seconds.
Elapsed time is 0.001463 seconds.
由此可见，快速卷积快
3.编程实现利用重叠相加法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与第二题的结果进行比较。

实验代码及结果：
程序Exp3_3代码
L = input('L=');
M = input('M=');
xn = linspace(1, 1, L);
hn = cos(0.2*pi*[0:M-1]);
tic
Xk = fft(xn, L+M-1);
Hk = fft(hn, L+M-1);
Yk = Xk.*Hk;
yn2 = ifft(Yk);
toc
yn = [];
tic
x1n = [linspace(0,0,M-1) xn];
Hk = fft(hn, 2*M);
iend = floor(length(x1n)/(M+1))+1;
for i = 0: iend
xin = x1n(i*(M+1)+1 : min([length(x1n) i*(M+1)+2*M]));
Xik = fft(xin, 2*M);
Yik = Xik.*Hk;
yin = ifft(Yik);
yn = [yn yin(M:2*M)];
end
yn = yn(1: L+M-1);
toc
>> Exp3_3
L=2048
M=256
Elapsed time is 0.000340 seconds.
Elapsed time is 0.000513 seconds.
>> Exp3_3
L=8192
M=256
Elapsed time is 0.003971 seconds.
Elapsed time is 0.002105 seconds.
分析：当L = 2048 M=256时，重叠相加比较慢，原因是程序控制流耗时较多。

但是，随着序列L的增长，重叠相加变得较快
4.编程实现利用重叠保留法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与第二题的结果进行比较。

实验代码及结果：
程序Exp3_4代码
L = input('L=');
M = input('M=');
xn = linspace(1, 1, L);
hn = cos(0.2*pi*[0:M-1]);
tic
Xk = fft(xn, L+M-1);
Hk = fft(hn, L+M-1);
Yk = Xk.*Hk;
yn2 = ifft(Yk);
toc %Ö±½Ó¼ÆËã¾í»ý
tic
Hk = fft(hn, 2*M-1);
iend = floor(length(xn)/M);
i = 0;
xin = xn(i*M+1 : min([(i+1)*M length(xn)]));
Xik = fft(xin, 2*M-1);
Yik = Xik.*Hk;
yin = ifft(Yik);
yn = yin; %µÚÒ»´Î¼ÆËã
for i = 1: iend
xin = xn(i*M+1 : min([(i+1)*M length(xn)]));
Xik = fft(xin, 2*M-1);
Yik = Xik.*Hk;
yin = ifft(Yik);
yn = [yn(1:(length(yn)-M+1))
yn((length(yn)-(M-2)):length(yn)) + yin(1:M-1) yin(M:2*M-1)]; %Æ´½ÓÐòÁÐ
end
yn = yn(1:L+M-1);
toc
>> Exp3_4
L=2048
M=256
Elapsed time is 0.000441 seconds.
Elapsed time is 0.001055 seconds.
>> Exp3_4
L=32768
M=1024
Elapsed time is 0.018384 seconds.
Elapsed time is 0.016179 seconds.
分析：当L = 2048 M=256时，重叠保留比较慢，原因是程序控制流耗时较多
但是，随着序列L及M的增长，重叠相加变得较快。

五、收获与体会
正式由于DFT存在着快速算法——FFT,所以DFT才成为了实用的工具，这也使得实时信号处理成为可能。

本实验通过快速算法和朴素算法的对比，进一步体会到了快速算法的优越性。

由于实验者本人MATLAB编程水平很低，程序控制语句的应用并不出色，所以分段卷积在某些情况下出现了相反的结果。

然而，其趋势是明显的，如果用专用硬件实现这些算法，其潜力是巨大的。