西安交大数字信号处理实验

1、数字信号处理实验快速傅里叶变换学院： 班级： 姓名： 学号：快速傅里叶变换一、实验目的1. 在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解；2. 熟悉并掌握按时间抽取FFT算法的程序；3. 了解应用FFT进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT。二、实验要求1. 记录下实验内容中各种情况下的X (k)值，做出频谱图并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。频谱只做模特性，模的最大值1，全部归一化；2. 打印出用C语言（或MATLAB 语言）编写的FFT源程序，并且在每一小段处加上详细的注释说明；3. 用C语言（或MATL

2、AB 语言）编写FFT程序时，要求采用人机界面形式：N , T , f 变量均由键盘输入，补零或不补零要求设置一开关来选择。三、实验内容1. 仔细分析教材第六章时间抽取法FFT 的算法结构，编制出相应的用FFT进行信号分析的C语言（或MATLAB 语言）程序；2. 用FFT程序计算有限长度正弦信号分别在以下情况下所得的DFT结果并进行分析和讨论：a) 信号频率f50Hz，采样点数N=32,采样间隔T=0.000625sb) 信号频率f50Hz，采样点数N=32,采样间隔T=0.005sc) 信号频率f50Hz，采样点数N=32,采样间隔T=0.0046875sd) 信号频率f50Hz，采样点数

3、N=32,采样间隔T=0.004se) 信号频率f50Hz，采样点数N=64,采样间隔T=0.000625sf) 信号频率f250Hz，采样点数N=32,采样间隔T=0.005sg) 将c）信号后补32个0，做64点FFT四、程序编写FFT运算：先进行采样，后对采样序列进行码位倒置，再进行蝶形算法。1.采样用for循环来将采样点逐个存入数组x(j)中。for j=0:1:N-1x(j+1)=sin(2*pi*f*j*T); end2.补零采用insert语句对是否补零来进行判断，如果补零，就将补零数加到采样点数，且补零的那些点在数组中都存为0。insert=input('是否补0？0-

4、否；1-是:'); if insert=1 ZERO=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:N s=dec2bin(t-1,M); s=fliplr

5、(s); s=bin2dec(s); A(s+1)=x(t); end4.fft蝶形运算 对fft运算采用了三个for循环来实现。第一个循环实现fft每一级的运算，所以它的范围应该为1M,循环一次做一级的fft运算。第二个循环实现每一级运算中的分组，所以其分组与所处的运算级数有关，例如第一级有16的组，第二级有8个组，所以其上限为2(L-1)-1。第三个循环实现每一个组中的fft运算，但每个组中进行的fft运算次数与分组有关，所以运算时对要进行蝶形运算的采样点的下标要用k和级数来控制，而对于系数W与分组和级数有关。for L=1:1:M %级数for J=0:1:(2(L-1)-1) %分组f

6、or k=(J+1):2L: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

7、onendaxis(0 N 0 1);7.完整程序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？0-否；1-是:'); if insert=1 ZERO=input('补零数：'); for j=N:1:N+ZERO-1x(j+1)=0; end

8、N=N+ZERO; end%码位倒置%M=log2(N); for t=1:1:N s=dec2bin(t-1,M); s=fliplr(s); s=bin2dec(s); A(s+1)=x(t); end%按时间抽取的FFT蝶形运算%for L=1:1:M for J=0:1:(2(L-1)-1) for k=(J+1):2L:N 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; endendend%模值归一化%x=abs(A); y=max

9、(x); X=x/y;%绘图for j=1:1:Nstem(j-1,X(j); hold onendaxis(0 N 0 1);五、实验结果分析在本次试验中，进行DFT运算的序列是我们对一个在时域为连续周期的信号进行采样得到的，这个有限长序列（在本题中为32个点）可以看作是对无限长的采样得到的序列进行加窗得到的，这个窗是一个长度为32的矩形窗。由FT的性质知，在时域连续周期的信号的频域是离散非周期的。对时域信号的采样会导致频域发生周期延拖，而延拖的周期为采样频率，如果对采样后的序列进行加窗，那么会在频谱上卷积一个sinc函数。最终得到的DFT的结果是对这个已经卷积过的频谱进行频域采样得到的。a

10、. 信号频率f50Hz，采样点数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的结果只有在

11、第1个点和第31个点有大小相等的值，其它点为0。b. 信号频率f50Hz，采样点数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。所以

12、DFT的结果只有在第8个点和第24个点有大小相等的值，其它点为0.c. 信号频率f50Hz，采样点数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

13、/16,25*pi/16这4个距离sinc函数主瓣最近的频率处会有比较大的幅值，其他点依次衰减。d. 信号频率f50Hz，采样点数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点也为一个峰值，其余点的

14、值一次在幅度上衰减。e. 信号频率f50Hz，采样点数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函数

15、都为0。所以DFT的结果只有在第2个点和第62个点有大小相等的值，其它点为0.f. 信号频率f250Hz，采样点数N=32,采样间隔T=0.005s信号频率f250Hz，采样点数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点上除了在