快速卷积与快速相关——MATLAB

一、实验目的

1、学会FFT算法程序(或函数）的使用方法;
2、了解序列的线性卷积和圆周卷积之间的关系;
3、验证有限长FFT算法实现线性相关运算的快速计算方法;
4、解FFT的点数对快速卷积与快速相关运算结果所产生的影响;
5、了学会利用FFT算法进行有限长序列的线性卷积的快速计算;
6、掌握基-2快速傅立叶变换（Fast FourierTransform,FFT)的算法原理及其程序实现方法.

二、实验原理

1、有限长序列的线性卷积和圆周卷积线性卷积和圆周卷对于有限长序列，存在两种形式的卷积,即积。设x(n)是长度光JM的有限长序列，y(n)是长度为N的有限长序列,则二者的线性卷积可表示为:
M-1

线性卷积的结果序列f(n)是一个长度为L=N +M -1的有限有限长序列，且长序列。将x(n)及y(n)均补零增长为L点的二者的L点的圆周卷积可表示为:

圆周卷积的结果序列f(n)是一个长度为L的有限长序列，由圆周卷积的点数所决定。有限长序列的线性卷积和圆周卷积之间的关系可用公式表示如下:

即:圆周卷积是线性卷积以圆周卷积的点数几为周期进行周期延拓后所取的主值序列。因而,在圆周卷积的点数大于或等于线性卷积的长度时、圆周卷积结果和线性卷积结果相等，这也是快速卷积算法的理论基础之一。
2、离散傅里叶变换的卷积性质
离散傅里叶变换的卷积性质也是快速卷积算法的另一理论基础。若f(n)是有限长序列x(n)和有限长序列y(n)的L点圆周卷积,即公式（5-2），则 f(n)的L点离散傅里叶变换为: F_c (k)=X(k)Y(k)
3、Matlab中FFT与IFFT的实现
离散傅立叶变换（Discrete Fourier Transform, DFT)实现了频域的离散化,方便了计算机处理,在数字信号处理中有着非常重要的作用。但直接计算DFT的运算量与变换长度N的平方成正比,计算量太大。而快速傅立叶变换FFT则是快速计算DFT的有效算法，大大提高了DFT的运算效率，在信号频谱的分析、滤波器频率响应的计算，以及线性卷积的快速计算等方面起着非常重要的作用。FFT 采用分组计算的方式进行DFT的快速计算,具体算法原理参看教材，在附录B中也给出了常用的基-2时间抽取FFT算法和分裂基FFT 算法的C语言程序。相应的,IFFT 则为离散傅里叶反变换，即 IDFT 的快速计算方法。在Matlab中，提供了f(t)和 ifi(t)两个函数来分别实现快速傅立叶变换的正变换和反变换。Ft(t)和if(t)两个函数是用机器语言而不是Matlab 指令写成的，执行速度很快。除了输入、输出参数的含义不同之外，这两个函数的调用方法完全相同，因此以ff()函数为例说明二者的使用方法。ff()函数常用调用格式有两种:
(1）Xk=fft(xn)
其中，xn为输入时域序列x(n)，返回结果xk为x(n)的离散傅里叶变换X(k)。当xn是矩阵时（对应于多通道信号)，计算xn中每一列信号的离散傅里叶变换。当xn的长度是2的整数幂,采
用基2快速算法计算,否则采用较慢的混合基算法进行计算。
（2）Xk=fft(xn，NFFT)
这种调用格式相比较于上一种调用格式,多了一个输入参数NFFT,用于指定FFT的点数。当NFFT的值是2的整数幂,采用基-2快速算法计算，否则采用较慢的混合基算法进行计算。当xn的长度大于NFFT时,对 xn进行自动截断;当xn的长度小于NFFT时,在xn后自动进行补零。
4、快速卷积基本算法的原理
利用FFT进行有限长序列的线性卷积的快速计算即为快速卷积算法。按照上述相关原理,利用FFT计算线性卷积,即计算公式中的f(n)的算法步骤可用下图来表示:

其中，FFT运算的点数L应满足L≥N+ M-1。为了采用基-2的算法，常常需要L还应满足L=2M
5、快速相关算法原理
利用 FFT 讲行有限长序列相关运算的快速实现则称为快速相关算法。若(n)是长度头M的有限长序列，y(n)是长度为N的有限长序列，则一者的线性百相关序列R(m)与二者的线性卷积运算之间的关系可表示

由于线性卷积可以用FFT讲行快速计算，则按公式，相关运算也可以利用fft进行快速计算，在此不再赘述其原理。需要时,根据傅立叶变换的反折和共轭特性可以减少一次FFT的使用提高计算效率。

三、实验步骤、数据记录及处理

本实验利用Matlab中提供的fft()和 ifft()函数进行快速卷积算法和快速相关算法性卷积和圆周卷积之间的关系进行验证。具体的实验内容和实验的实现，并对有限长序列线
步骤如下所示:
其中:
1、用 Matlab生成两个有限长序列x(n) y(n)

(1）基于fft()和 ifft()函数，编程利用4点快速卷积算法计算有限长序列x(n)与y(n)的卷积,结果令为c1(n)。
(2）基于fft()和 ifft()函数,编程利用速卷积算法计算有限长序列x(n)与y(n)的卷积，结果令为c2(n)
(3）调用conv()函数计算有限长序列x(n)与y(n)的卷积,结果令为c3(n)。分别绘制序列x(n)、y(n)、c1(n)、c2(n)和c3(n)的图形。对结果进行分析，并通过实验结果验证有限长序列线生卷积和圆周卷积之间的关系。
2、设两个有限长序列x(n)和h(n)分别为:
(1）x(n)=(sin 0.4n)·R,(n+1);
(2）h(n)=(0.9)"R,o(n+1)
按图所示的快速卷积算法原理编写完整的快速卷积算法程序计算y(n)=x(n)*h(n)，绘制序列图形，并通过 conv()函数对结果进行验证。
3、将实验1中实验内容4所给定的信号利用快速相关算法进行自相关序列的计算,并将结果与实验1中的结果进行对比验证。
实验程序：

clc;clear all;close all;
nx=0:1:2;x=[1,1,1];%确定x(n)与y(n)的自变量取值范围
ny=0:1:2;y=[2,2,2];%生成x(n)与h(n)
L=pow2(nextpow2(length(x)+length(y)-1));%确定FFT快速卷积的点数
xk=fft(x,L);%计算x的L点FFT,结果为x(k)
yk=fft(y,L);%计算y的L点FFT,结果为y(k)
YK=xk.*yk;%计算y(k)
c1n=ifft(YK,4);%四点
c2n=ifft(YK,L);%八点
c3n=conv(x,y);%线性卷积
n1=0:1:3;%确定n1的取值范围
n2=0:1:7;%确定n2的取值范围
n3=0:1:4;%确定n3的取值范围
%绘制x(n)与y(n)波形
figure('name','x(n)与y(n)序列');
subplot(2,1,1);stem(nx,x);grid on;xlabel('n');ylabel('x(n)');title('序列x(n)');
subplot(2,1,2);stem(ny,y);grid on;xlabel('n');ylabel('y(n)');title('序列y(n)');
%绘制x(n)与y(n)序列卷积波形
figure('name','x(n)与y(n)序列卷积');
subplot(3,1,1);stem(n1,c1n);grid on;xlabel('n');ylabel('c1(n)');title('x(n)与y(n)的4点fft');
subplot(3,1,2);stem(n2,c2n);grid on;xlabel('n');ylabel('c2(n)');title('x(n)与y(n)的8点fft');
subplot(3,1,3);stem(n3,c3n);grid on;xlabel('n');ylabel('c3(n)');title('x(n)与y(n)的线性卷积');

clc;clear all;close all;
nxn=1:15;nhn=1:20;%确定x(n)与y(n)的自变量取值范围
xn=sin(0.4*nxn);hn=0.9.^nhn;%生成x(n)与h(n)
L=pow2(nextpow2(length(xn)+length(hn)-1));%确定FFT快速卷积的点数
Xk=fft(xn,L);%计算x的L点FFT,结果为x(k)
Hk=fft(hn,L);%计算y的L点FFT,结果为y(k)
Yk=Xk.*Hk;%计算YK
y1n=ifft(Yk,L);%对YK调用IFFT，求得y1(n)
y2n=conv(xn,hn);%计算y2(n)的卷积
n1=(nxn(1)+nhn(1)):1:(L+nxn(1)+nhn(1)-1);%确定n1的自变量取值范围
n2=(nxn(1)+nhn(1)):1:(nxn(15)+nhn(20));%确定n2的自变量取值范围
%绘制有限长序列想x(n)与h(n)
figure('name','x(n)与h(n)序列');
subplot(2,2,1);stem(nxn,xn);grid on;xlabel('n');ylabel('x(n)');title('x(n)=(sin0.4n)*R15（n+1）');
subplot(2,2,2);stem(nhn,hn);grid on;xlabel('n');ylabel('h(n)');title('h(n)=(（0.9)^n)R20(n+1)');
subplot(2,2,3);stem(n1,y1n);grid on;xlabel('n');ylabel('y(n)');title('快速卷积法计算y(n)');
subplot(2,2,4);stem(n2,y2n);grid on;xlabel('n');ylabel('y(n)');title('conv()函数计算y(n)');

clc;clear all;close all;
n1=1:100;n2=1:100;%确定n1,n2取值范围
x=[101,82,66,35,31,7,20,92,154,125,85,68,38,23,10,24,83,132,131,118,90,67,60,47,41,21,16,6,4,7,14,34,45,43,48,42,28,10,8,2,0,1,5,12,14,35,46,41,30,24,16,7,4,2,8,17,36,50,62,67,71,48,28,8,13,57,122,138,103,86,63,37,24,11,15,40,62,98,124,96,66,64,54,39,21,7,4,23,55,94,96,77,59,44,47,30,16,7,37,74];
mean(x(:))
x=x-mean(x);
Rx=xcorr(x);
y=x(end:-1:1);
subplot(2,1,1);stem(Rx);grid on;xlabel('t');ylabel('Rx');title('xcorr()函数计算的自相关序列');
L=pow2(nextpow2(length(n1)+length(n2)-1));%确定FFT快速卷积的点数
Xk=fft(x,L);%计算x的L点FFT,结果为x(k)
Yk=fft(y,L);%计算y的L点FFT,结果为y(k)
Rk=Xk.*Yk;%计算YK
Rn=ifft(Rk,L);%对RK调用IFFT，求得y1(n)
n=(n1(1)+n2(1)):(L+n1(1)+n2(1)-1);%确定n的自变量取值范围
subplot(2,1,2);stem(1:200,Rn(1:200));grid on;xlabel('t');ylabel('R(n)');title('fft计算的自相关序列');

四、思考题

（1）对实验内容2中快速卷积基本算法的运算量进行分析,即乘法和加法运算次数。

（2）利用实例说明快速卷积基本算法的适用条件，即在什么情况下效率最高。

（3）实验内容3中,如何通过DFT的性质，减少一次ftt（）函数的调用,提高计算效率?

五、总结

通过此次实验的练习，加深理解FFT 在实现数字滤波（或快速卷积）中的重要作用，更好的利用FFT进行数字信号处理，并且掌握了循环卷积和线性卷积两者之间的关系。