信号分析与处理实验报告程序+报告内容(非常全)

实验四离散傅里叶变换及其快速算法一、实验目的1、掌握计算序列的离散傅里叶变换（DFT）的方法。2、掌握实现时间抽取快速傅里叶变换（FFT)的编程方法。3、通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，比较DFT与FFT的运算效率。4、加深对DFT与序列的傅里叶变换和Z变换之间关系的理解。二、实验原理2.1离散傅里叶变换（DFT）算法有限列长为N的序列)(nx的DFT变换为nknNnWnxkX10)()(12,1,0Nk其逆变换为10)(1)(NknkNWkXNnx1,1,0Nn2.2快速傅里叶变换基本思路当N很大的时候，求一个N点的DFT要完成NN次复数乘法和(1)NN次复数加法，其计算量相当大。快速傅里叶变换巧妙地利用NW因子的周期性和对称性，构造了一个DFT快速算法。其基本思路如下：1、利用nkNW的对称性使DFT运算中有些项合并)()(knNknNnNkN、利用nkNW的周期性和对称性使长序列的DFT分解为更小点数的DFT)()()(nNkNNnkNnkN为了讨论方便，设2N，其中为整数。如果不满足这个条件，可以认为得加上若干零点来达到。由DFT的定义知：nknNnWnxkX10)()(12,1,0Nk其中()xn是列长为(0,11)NnN的输入序列，把它按n的奇偶分成两个序列：)()12()()2(21rxrxrxrx12,1,0Nr又由于22222jNjNNNWeeW，则111200()()()()()NNnknkknNNnnnnXkxnWxnWXkWXk为偶数为奇数上式表明了一个N点的DFT可以被分解为两个N/2点的DFT。同时，这两个N/2点的DFT按照上式又可以合成为一个N点的DFT。为了要用点数为N/2点的X1(k)、X2(k)来表达N点的X(k)值还必须要用W系数的周期性，即()222NrkrkNNWW这样可得22()21110022()()()2NNNrkrkNNrrNXkxrWxrW即11()()2NXkXk同理可得22()()2NXkXk另外再加上WkN的对称性kNkNNNkNN2)2(就可以将X(k)的表达式分为前后两个部分：前半部分12()()()kNXkXkWXk12,1,0Nr后半部分()212()()()222NkNNNNXkXkWXk)()(21kXWkXkN12,1,0Nr由以上分析可见，只要求出区间0,12N内各个整数k值所对应的X1(k)、X2(k)的值，即可求出0,1N区间内的全部X（k）值，这一点恰恰是FFT能大量节省计算的关键所在。2.3时间抽取的快速傅里叶变换计算机算法1、程序输入序列的元素数目必须为2的整数次幂，即2MN,整个运算需要M级蝶形运算；2、输入序列应该按二进制的码位倒置排列，输出序列按自然序列排列；每个蝶形运算的输出数据占用其他输入数据的存储单元，实现“即位运算”；每级包括N/2个基本蝶形运算，共有M*N/2个基本蝶形运算；3、第L级中有/(2)NL个群，群与群的间隔为2L；4、同一级的各个群的系数Wn分布相同，第L级的群中有12L个系数；5、处于第L级的群的系数是(1)/2LNpNW(p=1,2,3，…….,12L)；6、对于第L级的蝶形运算，两个输入数据的间隔为12L。三、实验程序及注释3.1用C语言实现FFT与IFFT#includestdio.h#includemath.h#includestdlib.h/*清屏命令clrscr()需加载的预处理函数*/#defineswap(a,b)temp=(a);(a)=(b);(b)=temp/*宏定义的交换函数*//*快速傅利叶变换程序，数组A、B分别是带变换序列的实部和虚部ap是变换的类别，ap=1是做FFT;ap=-1是做IFFT*/voidfft(floatA[],floatB[],unsignedM,intap){unsignedlongN,I,J,K,L,LE,LE1,P,Q,R;inti;floatWr,Wi,W1r,W1i,WTr,WTi,theta,Tr,Ti,temp;N=pow(2,M);/*N=2M是序列的总长度*/J=0;for(I=0;IN-1;I++)/*码位倒置*/{if(JI){swap(A[I],A[J]);swap(B[I],B[J]);}K=N1;/*即k=k/2*/while(K=2&&J=K){J-=K;K=1;}J+=K;}for(L=1;L=M;L++)/*外层循环由级数L控制，执行M次*/{LE=1L;/*LE=2L是群间隔*/LE1=LE/2;/*LE1=2L-1是每个群的系数W数目*/Wr=1.0;Wi=0.0;theta=(-1)*ap*3.1415926536/LE1;/*ap用于在IFFT时改变极性*/W1r=cos(theta);W1i=sin(theta);for(R=0;RLE1;R++)/*中层循环由群系数控制，执行次*/{for(P=R;PN-1;P+=LE)/*R是群系数的编号，P、Q是基本蝶形运算两个输入数据在数组中的编号，循环每次完成同一个系数的蝶形运算*/{Q=P+LE1;Tr=Wr*A[Q]-Wi*B[Q];Ti=Wr*B[Q]+Wi*A[Q];/*Tr、Ti是的实部和虚部*/A[Q]=A[P]-Tr;/*即*/B[Q]=B[P]-Ti;A[P]+=Tr;/*即*/B[P]+=Ti;}WTr=Wr;/*Wr、Wi是的实部和虚部*/WTi=Wi;/*下面用和差化积公式求的实虚部Wr、Wi*/Wr=WTr*W1r-WTi*W1i;/**/Wi=WTr*W1i+WTi*W1i;/**/}}if(ap==-1){/*ap=-1是做IFFT，应使变换后的序列除以N*/for(i=0;iN;i++){A[i]/=N;B[i]/=N;}}return;}main(){inti,t,ap;floatA[100],B[100];unsignedM;system(cls);/*清屏*/printf(\npleaseinputbianhuanleibie(FFT:1orIFFT:-1):);scanf(%d,&ap);/*输入序列变换的类别*/printf(\npleaseinputjishuM:);/*M是级数*/scanf(%d,&M);t=pow(2,M);for(i=0;it;i++){printf(\npleaseinputx(%d):(ReIm),i);scanf(%f%f,&A[i],&B[i]);/*输入离散序列的实、虚部*/}fft(A,B,M,ap);/*进行FFT(或IFFT)*/printf(\n*************Result**************);for(i=0;it;i++)printf(\nX(%d):%.5f+j%.5f,i,A[i],B[i]);/*输出变换的结果*/}实例验证：用课本序列计算FFT得：pleaseinputbianhuanleibie(FFT:1orIFFT:-1):1pleaseinputjishuM:2pleaseinputx(0):(ReIm)10pleaseinputx(1):(ReIm)20pleaseinputx(2):(ReIm)-10pleaseinputx(3):(ReIm)40*************Result**************X(0):6.00000+j0.00000X(1):2.00000+j2.00000X(2):-6.00000+j0.00000X(3):2.00000+j-2.00000将上述做IFFT得：pleaseinputbianhuanleibie(FFT:1orIFFT:-1):-1pleaseinputjishuM:2pleaseinputx(0):(ReIm)60pleaseinputx(1):(ReIm)22pleaseinputx(2):(ReIm)-60pleaseinputx(3):(ReIm)2-23.2matlab实现DFT与IDFT（包含频谱混叠）现以一连续时间信号为例，通过采样，得到离散时间信号；再对离散时间信号进行DFT与IDFT。3.2.1实现DFT的子程序function[xy]=dft(x,N)%%该程序完成傅立叶变换%x输入要进行傅立叶变换的数列%N输入数列所包含的元素个数forn=1:Nfork=1:NX(n,k)=x(:,k)*exp(-j*2*pi*(k-1)*(n-1)/N);endendxy=conj(sum(X,2)');3.2.2实现IDFT的子程序function[xy]=idft(x,N)%%该程序完成傅立叶逆变换%x输入要进行傅立叶逆变换的数列%N输入数列所包含的元素个数forn=1:Nfork=1:NX(n,k)=x(:,k)*exp(j*2*pi*(k-1)*(n-1)/N);endendxy=sum(X')/N;3.2.3主程序%%利用傅里叶变换绘制信号的幅频特性clc;clearallTs=1/125;%采样时间间隔N=64;%采样点数T0=N*Ts;%采样时间长度M=N*Ts/T0;%M*T0为采样时间，M=1表示采样为整周期采样,M不等于1会发生频谱泄露%可以观察整周期和非整周期的区别，体会频谱泄露概念forn=1:Nxx(n)=10*sin(2*pi*50*n*Ts+pi/3);%要进行傅立叶变换的连续时间函数endsubplot(2,1,1)t=0:(M*T0)/N:(M*T0*(N-1))/N;stem(t,xx);%采样xlabel('t/s');ylabel('x(t)/(A)');title('时域波形');ticyy=dft(xx,N);%对时域信号做傅立叶变换toc%deltaf=1/(T0*M);%书上（6-78）%forn=1:N/2+1%yy1(n)=yy(n)/N;%书上（6-76）%end%subplot(1,2,2)%f=0:deltaf:N/2*deltaf;%%stem(f,abs(yy1))subplot(2,1,2)yy1=idft(yy,N);%对时域信号做傅立叶逆变换stem(t,yy1)xlabel('t/s');ylabel('x(t)/(A)');title('傅立叶反变换');3.2.4实验结果图（workspace里的各数据结果在此不一一列出）由图可知，离散时间信号经过傅里叶变换和傅里叶反变换后，没有发生改变，说明了程序的正确性。3.2.5频谱混叠%%验证主程序clc;clearallTs=0.005;%采样时间间隔N=64;%采样点数T0=N*Ts;%采样时间长度M=N*Ts/T0;%M*T0为采样时间，M=1表示采样为整周期采样,M不等于1会发生频谱泄露%可以观察整周期和非整周期的区别，体会频谱泄露概念forn=1:Nxx(n)=10*sin(2*pi*50*n*Ts+pi/3);%要进行傅立叶变换的连续时间函数endsubplot(2,1,1)t=0:(M*T0)/N:(M*T0*(N-1))/N;yy=fft(xx,N);%对时域信号做傅立叶变换deltaf=1/(T0*M);%书上（6-78）forn=1:N/2+1yy1(n)=yy(n);%书上（6-76）endsubplot(2,1,1)f=0:deltaf:N/2*deltaf;%stem(f,abs(yy1))xlabel('f/Hz');ylabel('|X(k)|/(A)');title('傅立叶变换');forn=1:Nxx(n)=10*sin(2*pi*50*n*Ts+pi/3);%要进行傅立叶变换的连续时间函数end%%验证频谱混叠Ts1=