数字信号处理实验报告 (2)

数字信号处理实验快速傅里叶变换实验姓名班级学号日期 2011.10.25实验名称 快速傅里叶变换一．实验目的1 在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解；2 熟悉并掌握按时间抽取FFT 算法的程序；3 了解应用FFT 进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT 。

二.实验要求1 记录下实验内容中各种情况下的X (k)值，做出频谱图并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。

频谱只做模特性，模的最大值＝1，全部归一化；2 打印出用C 语言（或MA TLAB 语言）编写的FFT 源程序，并且在每一小段处加上详细的注释说明；3 用C 语言（或MATLAB 语言）编写FFT 程序时，要求采用人机界面形式：N , T , f 变量均由键盘输入，补零或不补零要求设置一开关来选择。

三.实验内容1.仔细分析教材第六章‘时间抽取法FFT ’的算法结构，编制出相应的用FFT 进行信号分析的C 语言（或MATLAB 语言）程序；2.用FFT 程序计算有限长度正弦信号()s i n (2**)[()()],0y t f t u t u t N T t N T π=--<< 分别在以下情况下所得的DFT 结果并进行分析和讨论：a ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.000625sb ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.005sc ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.0046875sd ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.004se ） 信号频率f ＝50Hz ，采样点数N=64,采样间隔T=0.000625sf ） 信号频率f ＝250Hz ，采样点数N=32,采样间隔T=0.005sg ） 将c ） 信号后补32个0，做64点FFT四．内容分析FFT 运算时，依次要进行采样，码位倒置，蝶形运算，并且要用三层循环来控制，最内层循环控制做蝶形算法，中层循环完成每一级中相同系数的蝶形算法，最外层循环控制级数递推。

1数组采样：for j=0:1:N-1x(j+1)=sin(2*pi*f*j*T);end2.判断是否补零点补零的点在数组中都存为0insert=input('是否补0？ 1(yes)/0(no):'); %补零if insert==1ZERO=input('补零数：');for j=N:1:N+ZERO-1x(j+1)=0;endN=N+ZERO;end3.码位倒置M为fft运算的级数，在这里被用为表示数组中数据下表的二进制代码位数上限。

用for循环来实现对数组中的每个数据进行处理。

其中dec2bin（）函数将十进制代码转换为二进制代码，fliplr（）将二进制的矩阵进行左右对称的翻转，即实现了码位倒置，bin2dec（）在讲二进制代码转换为十进制代码。

M=log2(N);for t=1:1:Ns=dec2bin(t-1,M);s=fliplr(s);s=bin2dec(s);A(s+1)=x(t);end4.fft蝶形运算对fft运算采用了三个for循环来实现。

第一个循环实现fft每一级的运算，所以它的范围应该为1—M,循环一次做一级的fft运算。

第二个循环实现每一级运算中的分组，所以其分组与所处的运算级数有关，例如第一级有16的组，第二级有8个组，所以其上限为2^(L-1)-1。

第三个循环实现每一个组中的fft运算，但每个组中进行的fft运算次数与分组有关，所以运算时对要进行蝶形运算的采样点的下标要用k和级数来控制，而对于系数W与分组和级数有关。

for L=1:1:M %级数for J=0:1:(2^(L-1)-1) %分组for k=(J+1):2^L:N %组内fft运算T=A(k)+A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N);A(k+2^(L-1))=A(k)-A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N); A(k)=T;endendend5.归一化对采样点的幅值去绝对值，再用每一个绝对值除以其中的最大值，这样就可以得到归一化后的图形x=abs(A);y=max(x);X=x/y;6.画图用for循环来实现采样点的幅值与采样点的对应绘图。

Axis（）规定图形的纵轴范围（0,1）和横轴范围（0，N）。

for j=1:1:Nstem(j-1,X(j));hold onendaxis([0 N 0 1]);五．实验结果分析在本次试验中，进行DFT运算的序列是我们对一个在时域为连续周期的信号进行采样得到的，这个有限长序列（在本题中为32个点）可以看作是对无限长的采样得到的序列进行加窗得到的，这个窗是一个长度为32的矩形窗。

由FT的性质知，在时域连续周期的信号的频域是离散非周期的。

对时域信号的采样会导致频域发生周期延拖，而延拖的周期为采样频率，如果对采样后的序列进行加窗，那么会在频谱上卷积一个sinc函数。

最终得到的DFT的结果是对这个已经卷积过的频谱进行频域采样得到的。

a. 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.000625s信号的频率为50Hz，采样周期为0.000625s，则采样频率为1600Hz，在数字域，采样频率对应的是2*pi，由于sin信号的频谱为在pi/16和-pi/16上的冲激，所以时域采样后，在频域以2*pi为周期进行周期延拖。

加窗后，卷积一个sinc 函数，由于窗的长度为32，所以所以其主瓣的宽度为pi/16，因为进行DFT要有32个点，所以第k个点对应的频谱为（k-1）*pi/16,由上面的sinc函数的主瓣宽度知，除了pi/16和31*pi/16（由于频域的周期延拖），其他的点处，sinc函数都为0，DFT的结果只有在第1个点和第31个点有大小相等的值，其它点为0.b. 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.005s变化的是采样间隔变为0.005s，这样采样频率就变为200Hz，同样采样频率对应数字域的2*pi，所以sin信号的谱在数字域对应的频率为+0.5*pi和-0.5*pi,周期延拖后在1.5*pi处也会有，故在0到2*pi中有两个冲激，窗的宽度还是32，所以sinc函数的主瓣宽度和上题一样，因为进行DFT要有32个点，所以第k 个点对应的频谱为（k-1）*pi/16,由上面的sinc函数的主瓣宽度知，除了pi/2和3*pi/2（由于频域的周期延拖），其他的点处，sinc函数都为0，所以DFT的结果只有在第8个点和第24个点有大小相等的值，其它点为0.c. 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.0046875s采样周期为0.0046875s，所以采样频率为213.33Hz，不是整数，因为在数字域采样频率对应的是2*pi，所以可以知道sin函数在数字域里对应的是0.234*pi，和-0.234*pi，可见不是pi/16的整数倍，所以其卷积sinc函数后，由于sinc函数的主瓣仍然保持为pi/16，所以频率在k*pi/16上不为0而是两个sinc函数的叠加，由于sinc函数从主瓣向两边延伸时，幅度渐渐衰减，所以最终频谱为在7*pi/16,8*pi/16和24*pi/16,25*pi/16这4个距离sinc函数主瓣最近的频率处会有比较大的幅值，其他点依次衰减。

d. 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.004s采样周期为0.004s，所以采样频率为250Hz，采样频率对应数字域的2*pi，所以sin函数的频谱在数字域对应的是2*pi/5和-2*pi/5，sinc函数与上述题目相同，主瓣为pi/16，所以在k*pi/16上不为0而是两个sinc函数的叠加，我们计算可知，sinc函数的主瓣最高点在6.4倍的pi/16处，距离6倍的pi/16近，距离7倍的pi/16远，所以DFT的频率图像在n=6点最高，由于周期性，n=26点也为一个峰值，其余点的值一次在幅度上衰减。

e. 信号频率f＝50Hz，采样点数N=64,采样间隔T=0.000625s信号的频率为50Hz，采样周期为0.000625s，所以采样频率为1600Hz，在数字域，采样频率对应的是2*pi，由于sin信号的频谱为在+1/16pi和-1/16pi上的冲激，所以时域采样后，在频域以2*pi为周期进行周期延拖。

加窗后，卷积一个sinc函数，由于窗的长度为64，所以所以其主瓣的宽度为1/16pi，因为进行DFT有64个点，所以第k个点对应的频谱为（k-1）*pi/32,由上面的sinc函数的主瓣宽度知，除了pi/16和31*pi/16（由于频域的周期延拖），其他的点处，sinc函数都为0，所以DFT的结果只有在第2个点和第62个点有大小相等的值，其它点为0.f. 信号频率f＝250Hz，采样点数N=32,采样间隔T=0.005s信号频率f＝250Hz，采样点数N=32,采样间隔T=0.005s，采样频率为200Hz，采样频率对应数字频域的2*pi，所以sin函数的频谱在2.5*pi和-2.5*pi上有冲激。

由于采样后产生周期性，在0.5*pi和1.5*pi处也有冲激。

sinc函数的主瓣宽度还是pi/16，因为0.5*pi和1.5*pi均是1/16pi的整数倍，所以sinc函数在其余k*pi/16处均为0，所以最终DFT的频谱只在n=8和24有等值的冲激，别的点为0.g. 将c）信号后补32个0，做64点FFT本题与c）题的区别在于最终做了一个64点的DFT，其它均没有变化，就是说，在对频谱采样前，本题的频谱和c）的是相同的。

在c中采样的间隔是pi/16,而在本题中采样的间隔是pi/32。

所以在2*n点上的值与c）对应的n点的值相同；在2*n+1点上除了在第15和第49点上有峰值外，其余点为0。

附件：clc;clear;f=input('信号频率f(Hz)：'); %输入信号频率f%N=input('采样点数N：'); %输入采样点数N%T=input('采样周期T(s)：'); %输入采样间隔T%%采样for j=0:1:N-1x(j+1)=sin(2*pi*f*j*T);end%补0insert=input('是否补0？ 1(yes)/0(no):');if insert==1ZERO=input('补零数：');for j=N:1:N+ZERO-1x(j+1)=0;endN=N+ZERO;end%码位倒置%M=log2(N);for t=1:1:Ns=dec2bin(t-1,M);s=fliplr(s);s=bin2dec(s);A(s+1)=x(t);end%按时间抽取的FFT蝶形运算%for L=1:1:Mfor J=0:1:(2^(L-1)-1)for k=(J+1):2^L:NT=A(k)+A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N);A(k+2^(L-1))=A(k)-A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N); A(k)=T;endendend%模值归一化%x=abs(A);y=max(x);X=x/y;%绘图for j=1:1:Nstem(j-1,X(j));hold onendaxis([0 N 0 1]);。