利用FFT计算卷积

姓名：高铭遥 班级：16131701 学号：1120171450 成绩：实验二 DFT/FFT 的应用-利用FFT 实现快速卷积[实验目的]1．深刻理解DFT/FFT 的概念和性质，进一步掌握圆周卷积和线性卷积两者之间的关系。

2．掌握DFT/FFT 的应用。

理解FFT 在实现数字滤波（或快速卷积）中的重要作用，更好地利用FFT 进行数字信号处理。

[实验内容及要求]1．给定两个序列()[]2,1,1,2x n =，()[]1,1,1,1h n =--。

首先直接在时域计算两者的线性卷积；然后用FFT 快速计算二者的线性卷积，验证结果。

（1）线性卷积 程序代码：figure(1);N1=4; N2=4; xn=[2,1,1,2]; hn=[1,-1,-1,1];N=N1+N2-1;%卷积后的序列长度 yn=conv(xn,hn);%线性卷积 x=0:N-1;stem(x,yn);title('线性卷积'); 运行结果：（2）FFT 卷积快速卷积 程序代码： figure(1); n=0:1:3; m=0:1:3;N1=length(n);%xn 的序列长度 N2=length(m);%hn 的序列长度 xn=[2,1,1,2]; hn=[1,-1,-1,1];姓名：高铭遥 班级：16131701 学号：1120171450 成绩：N=N1+N2-1;%卷积后的序列长度XK=fft(xn,N);%xn 的离散傅里叶变换 HK=fft(hn,N);%hn 的离散傅里叶变换 YK=XK.*HK;yn=ifft(YK,N);%逆变换if all(imag(xn)==0)&&(all(imag(hn)==0))%实序列的循环卷积仍为实序列 yn=real(yn); endx=0:N-1;stem(x,yn);title('FFT 卷积'); 运行结果：结果分析：对比（1）和（2）直接线性卷积和FFT 快速卷积的结果可以验证，用FFT 线性卷积的结果是与直接卷积的结果相同的，FFT 可以实现快速卷积，提高运算速度。

快速傅里叶变换实现卷积
快速傅里叶变换(FFT)是一种高效计算离散傅里叶变换(DFT)及其逆变换的方法，可以用于实现卷积的计算。

在实现卷积时，可以使用FFT将输入信号和系统函数都转换到频域，然后对应相乘得到输出信号的频谱，最后再通过逆FFT转换回时域得到输出信号。

这种通过FFT实现卷积的方法具有运算速度快、精度高等优点。

具体实现过程如下：
1. 将输入信号x(t)和系统函数h(t)分别转换成离散时间序列x[n]和h[n]。

2. 对x[n]和h[n]分别进行FFT变换，得到X(f)和H(f)。

3. 将X(f)和H(f)对应相乘得到Y(f)。

4. 对Y(f)进行逆FFT变换，得到输出信号y(t)。

需要注意的是，使用FFT实现卷积时，需要对输入信号和系统函数进行采样，并且采样频率应该满足采样定理的要求，否则会产生失真。

此外，还需要注意FFT的点数应该与输入信号和系统函数的长度相同，否则会产生误差。

1.FFT的运算量。

直接计算DFT共需N2次复数乘法和N(N-1)次复数加法。

而FFT仅需计算0.5N㏒2N次复数乘法和N㏒2N次复数加法。

由于在计算机上计算乘法所需的时间比计算加法多得多，所以FFT的运算量比DFT要少的少。

2.FFT及其反变换的MATLAB函数。

y=fft(x);y=fft(x,N).y=ifft(x);y=ifft(x,N)3.利用FFT计算线性卷积。

步骤：将x1(n)和x2(n)都延长到N=N1+N2-1。

计算x1(n)的N点FFT。

X1=FFT（x1）计算x2(n)的N点FFT。

X2=FFT（x2）计算Y= X1 .X2计算Y的反变换y=ifft(Y)在MATLAB中，可用nextpow2来确定呢最接近的值是2的多少次幂。

、4.利用FFT对信号进行频谱分析。

有X（k）群殴相位谱，振幅谱，功率谱。

常用的函数有，计算模值的abs函数和计算相角的angle函数。

因为实信号以fs为采样速率的信号在fs/2 处混叠，所以实信号fft的结果中前半部分对应[0, fs/2],后半部分对应[ -fs/2, 0]1）实信号fft的结果前半部分对应[0, fs/2]是正频率的结果,后半部分对应[ -fs/2, 0]是负频率的结果。

大于fs/2的部分的频谱实际上是实信号的负频率加fs的结果。

故要得到正确的结果，只需将视在频率减去fs即可得到频谱对应的真实负频率2）如果要让实信号fft的结果与[-fs/2, fs/2]对应，则要fft后fftshift一下即可，fftshift的操作是将fft结果以fs/2为中心左右互换3）如果实信号fft的绘图频率f从[-fs/2, fs/2]，并且没有fftshift，则fft正频谱对应f在[0, fs/2]的结果将混叠到(f - fs/2)的位置;fft负频谱对应f在[-fs/2, 0]的结果混叠到f + fs - fs/2 的位置，注意这里f为负值，也就是说此种情况下fft负频谱对应的视在频率减去fs/2即可得到频谱对应的真实负频率二. 复信号情况1）复信号没有负频率，以fs为采样速率的信号，fft的频谱结果是从[0, fs]的。

实验1 利用DFT 分析信号频谱一、实验目的1、加深对DFT 原理的理解。

2、应用DFT 分析信号的频谱。

3、深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。

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

三、实验基础理论1.DFT 与DTFT 的关系：有限长序列的离散时间傅里叶变换(e )j X ω 在频率区间(02)ωπ≤≤ 的N 个等间隔分布的点2(0k N 1)kk N πω=≤≤-上的N 个取样值可以有下式表示：2120(e )|(n)e(k)(0k N 1)N jkn j Nkk NX x X πωπω--====≤≤-∑由上式可知，序列(n)x 的N 点DFT (k)X ,实际上就是(n)x 序列的DTFT 在N 个等间隔频率点2(0k N 1)kk N πω=≤≤-上样本(k)X 。

2.利用DFT 求DTFT方法1：由(k)X 恢复出(e )j X ω的方法如下：由流程知：11(e )(n)e[(k)W]e N j j nkn j nNn n k X x X Nωωω∞∞----=-∞=-∞===∑∑∑继续整理可得到：12()(k)()Ni k kx e X N ωπφω==-∑其中(x)φ为内插函数：sin()2()sin()2N N ωφωω=方法2：实际在MATLAB 计算中，上述插值运算不见得是最好的办法。

由于DFT 是DTFT 的取样值，其相邻两个频率样本点的间距为2N π，所以如果我们增加数据的长度N ，使得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样就可以利用DFT 计算DTFT 。

如果没有更多的数据，可以通过补零来增加数据长度。

3.利用DFT 分析连续信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续信号离散化，这里需要进行两个操作：一是采样，二是截断。

对于连续时间非周期信号(t)a x ,按采样间隔T 进行采样，阶段长度M ，那么：1(j )(t)e(nT)e M j tj nTa a a n X x dt T x -∞-Ω-Ω-∞=Ω==∑⎰对(j )a X Ω 进行N 点频域采样，得到：2120(j )|(nT)e(k)M jkn Na a M kn NTX T x TX ππ--Ω==Ω==∑采用上述方法计算信号(t)a x 的频谱需要注意如下三个问题：（1）频谱混叠；（2）栅栏效应和频谱分辨率； （3）频谱泄露。

实验一 利用傅立叶变换计算线性卷积一、实验目的1． 掌握MATLAB 的使用。

2． 掌握用直接法计算线性卷积的原理和方法3． 掌握利用FFT 及IFFT 计算线性卷积的原理和方法二、实验原理及方法1、线性卷积的定义序列)1N n 0(),n (x -≤≤和序列)1M n 0(),n (h -≤≤的线性卷积y(n)=x(n)*h(n)定义为：10),()()(10-+≤≤-⨯=∑-=M N n m n h m x n y N m 利用直接法计算线性卷积即用线性卷积的定义计算。

2、利用FFT 及IFFT 计算线性卷积的原理和方法如果将序列x(n)和h(n) 补零，使其成为长度为L 的序列(L>=N+M-1), 则x(n)与h(n)的线性卷积y(n)=x(n)*h(n)与L 点圆周卷积相等，而圆周卷积可采用FFT 及IFFT 完成，即求y(n)=x(n)*h(n)可转化为：对上式两端取FFT 得： Y(k)=X(k)H(k)其中：X(k)=FFT[x(n)], H(k)=FFT[h(n)]则：y(n)=IFFT[Y(k)]三、实验仪器及材料⒈ 计算机，并装有MATLAB 程序⒉ 打印机四、实验步骤1、已知两序列： ⎩⎨⎧>≤≤=3n ;03n 0;)5/3()n (h n 用Matlab 随机生成输入信号X （n ），范围为0~2；2、得出用直接法（定义）计算线性卷积y(n)=x(n)*h(n)的结果；3、用Matlab 编制利用FFT 和IFFT （圆周卷积）计算线性卷积y(n)=x(n)*h(n)的程序； 分别令圆周卷积的点数为L=5，7，8，10，打印结果。

4、对比直接法和圆周卷积法所得的结果。

五、实验说明:1、实验前复习线性卷积,圆周卷积及FFT 内容。

2、利用FFT 计算线性卷积是将x(n)、h(n)用补零的方法延长到N+M-1，再用圆周卷积完成，因此要求x(n)、h(n)延长后的长度满足L>=N+M-1，才能保证用圆周卷积计算结果与直接法计算结果相同。

利用FFT和IFFT计算线性卷积利用FFT和IFFT计算线性卷积/**************************************************************** ********说明:用快速傅立叶变换计算两个有限长序列的线性卷积void dftCOMPLEX *x, COMPLEX *y,int n,int flagx:指向复数变量的指针变量,存放要变换的数据,长度n。

y:指向复数变量的指针变量,存放变换的结果,长度n。

n:整型变量,变换数据的长度。

flag: 整型变量,标志符,flag=1时,做DFT, flag-1时,做IDFT。

***************************************************************** *******/说明:用快速傅立叶变换计算两个有限长序列的线性卷积void convoldouble *x,double *y,int lenx,int leny,int lenx:双精度实型一维数组,长度为1en。

开始时存放实序列xi,程序结束时,存放线性卷积的结果。

y:双精度实型一维数组,长度为n。

存放实序列yi。

lenx:整型变量。

序列xi的长度。

leny:整型变量。

序列yi的长度。

len:整型变量,线性卷积的长度。

1enm+n-1,且必须是2的整数次幂/**************************************************************** ********#include “rfft.c”#include “irfft.c”#in clude “math.h”#include “graphics.h”#include “stdlib.h”#include “string.h”#include “draw.h”void signalxhdouble *x,int ch;void convolutiondouble *x,double *h,double *y,int lenx,int lenhvoid fconvoldouble *x,double *y,int lenx,int leny,int len;void mainint i,choice,select,lenx,leny,flag;double *x;// 输入信号double *h; //单位脉冲响应存放变换的结果double *y;//输出信号printf“\nchoice singal : choice??xn select??hn\n”;printf“1.矩形序列,长度10 2.单位冲击序列\n”;printf“3.单位阶跃序列,长度54.矩形序列,长度6\n”;printf“choice :选择输入信号 select:选择输出信号\n”; printf“choice”;sc anf“%d”,&choice;ifchoice1 lenx10;else ifchoice2 lenx1;else ifchoice3 lenx5;else lenx6;xdouble *calloclenx,sizeofdouble;printf“\nselect”;scanf“%d”,&select;ifselect 1 lenh10;else ifselect 2 lenh1;else ifselect 3 lenh5;else lenh6;leny lenx+lenh-1;hdouble *calloclenh,sizeofdouble;y double *callocleny,sizeofdouble;printf“\nchoice :1.continue 2.over\n”; dosingalxhx,choice;singalxhh,select;convolutionx,h,y,lenx,lenh;printf“Do you want to input length of Convolution:Y/N?”;ifgetchar’Y’||’y’printf“input length:len”;scanf“%d”,len;else len intlogleny/log2+0.9999;printf(”\nDirect Caculation of Linear Convolution\n”);fori0;ileny;i++printf“%10.1lf”,y[i];ifi%43 printf“\n”;printf“\n”;fconvolx,h,lenx,lenh,len;printf(”Fast Caculation of Linear Convolution\n”);fori0;ileny;i++printf“%10.1lf”,x[i];ifi%43 printf“\n”;printf“\n”;scanf“&d”,&flag;while!flag;void singalxhdouble *x,int chint i;double al,a,f,t;printf“1.矩形序列,长度10 2.单位冲击序列\n”;printf“3.1,2.4,-2.4,2.4,长度44.矩形序列,长度6\n”;switchchcase 1: fori0;i10;i++ x[i]1;break;case 2: x[0]1;break;case 3: x[0]1;x[1]2.4;x[2]-2.4;x[3]2.4;break;default: fori0;i5;i++ x[i]1;break;void convolutiondouble *x,double *h,double *y,int lenx,int lenh int len;int i,j,k;lenlenx+lenh-1;fork0;klen;k++y[k]0.0;fori0;ik;i++y[k]y[k]+x[i]*h[k-i];void fconvoldouble *x,double *y,int lenx,int leny,int lenint i,len2;double t;forforilenx;ilen;i++ x[i]0.0;forilenh;ilen;i++ y[i]0.0;rfftx,len;rffty,len;len2len/2;x[0]x[0]*y[0];x[len2]x[len2]*y[len2];fori1;ilen2;i++tx[i]*y[i]-x[len-i]*y[len-i];x[len-i] x[i]* y[len-i]+ x[len-i]* y[i]; x[i]t;irfftx,len;void rfftdouble *x,int nint i,j,k,m,i1,i2,i3,i4,nl,n2,n4;double a,e,cc,ss,xt,t1,t2;forj1,i1;i16;i++m=i;j=2*j;ifjn break;n1=n-1;forj0,i0;i<n;i++=ifi<jxtx[j];x[j]x[i];x[i]=x[j];kn/2;whilekj十1j=j?k;kk/2;j=j+k;fori=0;i<n;i+2=xtx[i];x[i]xt+x[i+1];x[i+1]=xt-x[i+1j]; n21;for k=2;km;k++n4=n2;n2=2*n4;n12*n2;e=6.28318530718/n1;fori=0;in;i+n1xtx[i];x[i]xt+x[i+n2];x[i+n2]=xt-x[i+n2];x[i+n2+n4]-x[i+n2+n4]; ae;forj1;jn4?1;j++i1=i+j;i2i-j+n2;i3=i+j+n2;i4i-j+n1;cccosa;ss=sina;aa+e;t1=cc*x[i3]+ss*x[i4];t2=ss*x[i3]-cc*x[i4];x[i4]=x[i2]-t2;x[i3]=-x[i2]-t2;x[i2]=x[i1]-t1;x[i1]=x[i1]+t;void irfftdouble *x,int nint i,j,k,m,i1,i2,i3,i4,i5,i6,i7,i8,n2,n4,n8,id,is; double a,e,a3,t1,t2,t3,t4,t5,cc1,cc3,ss1,ss3;forj1,i1;i16;i++m=i;j=2*j;ifjn break;n2=2*n;fork1;km;k++is0;id=n2;n2=n2/2;R4n2/4;n8n2/8;e=6.28318530718/n2; doforiis;in;i+idi1=i;i2i1+n4;i3i2+n4;i4i3+n4;tl=x[i1]-x[i3];x[i1]x[i1]+x[i3]; x[i2]2*x[i2];x[i3]t1-2*x[i4];x[i4]t1+2*x[i4];ifn41 continue;i1+n8;i3+n8;i4+n8;t1x[i2]-x[i1]/sqrt2.0; t2x[i4]+x[i3]/sqrt2.0; x[i1]x[i1]+x[i2];x[i2]x[i4]-x[i3];x[i3]2*-t2-t1;x[i4]2*-t2+t1;is2*id-n2;id4*id;whileisn-1;ae;forj=1;jn8;j++a3=3*a;cc1cosa;ss1sina;cc3=cosa3;ss3=sina3;a=j+1*e;id2*n2;doforiis;in-1;ii+idi1i-j;i2=i1+n4;i3=i2+n4;i4=i3+n4;15=i+n4-j;i6=i5+n4;i7=i6+n4;i8=i7+n4;t1=x[i1]-x[i6];x[i1]=x[i1]+x[i6]; t2=x[i5]-x[i2];x[i5]=x[i2]+x[i5]; t3=x[i8]+x[i3];x[i6]=x[i8]-x[i3]; t4=x[i4]+x[i7];x[i2]=x[i4]-x[i7];t5t1-t4;t1t1+t4;t4=t2-t3;t2=t2+t3;x[i3]t5*cc1+t4*ss1; x[i7]-t4*cc1+t5*ss1; x[i4]t1*cc3-t2*ss3; x[i8]t2*cc3+t1*ss3;is2*id-n2;id=4*id;whileisn-1;is=0;id=4;doforiis;in;i+idi1=i+1;t1=x[i];x[i]=t1+x[i1];x[i1]=t1-x[i1];is2*id-2;id4*id;whileisn-1forj=0,i0;in-1;i++ ifijt1=x[j];x[j]=X[i];x[i]=t1;k=n/2;whi1ekj+1j=j-k;k=k/2;j=j十k;fori=0;in;i++ x[i]x[i]/n;。