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

数字信号处理实验——快速傅里叶变换学院：班级：姓名：学号：-1-快速傅里叶变换一、实验目的1.在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解；2.熟悉并掌握按时间抽取FFT算法的程序；3.了解应用FFT进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT。二、实验要求1.记录下实验内容中各种情况下的X(k)值，做出频谱图并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。频谱只做模特性，模的最大值＝1，全部归一化；2.打印出用C语言（或MATLAB语言）编写的FFT源程序，并且在每一小段处加上详细的注释说明；3.用C语言（或MATLAB语言）编写FFT程序时，要求采用人机界面形式：N,T,f变量均由键盘输入，补零或不补零要求设置一开关来选择。三、实验内容1.仔细分析教材第六章‘时间抽取法FFT’的算法结构，编制出相应的用FFT进行信号分析的C语言（或MATLAB语言）程序；2.用FFT程序计算有限长度正弦信号()sin(2**)[()()],0ytftututNTtNT分别在以下情况下所得的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循环来将采样点逐个存入数组x(j)中。forj=0:1:N-1x(j+1)=sin(2*pi*f*j*T);end2.补零采用insert语句对是否补零来进行判断，如果补零，就将补零数加到采样点-2-数，且补零的那些点在数组中都存为0。insert=input('是否补0？[0-否；1-是]:');ifinsert==1ZERO=input('补零数：');forj=N:1:N+ZERO-1x(j+1)=0;endN=N+ZERO;end3.码位倒置M为fft运算的级数，在这里被用为表示数组中数据下表的二进制代码位数上限。用for循环来实现对数组中的每个数据进行处理。其中dec2bin（）函数将十进制代码转换为二进制代码，fliplr（）将二进制的矩阵进行左右对称的翻转，即实现了码位倒置，bin2dec（）在讲二进制代码转换为十进制代码。M=log2(N);fort=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与分组和级数有关。forL=1:1:M%级数forJ=0:1:(2^(L-1)-1)%分组fork=(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（）规定图形的-3-纵轴范围（0,1）和横轴范围（0，N）。forj=1:1:Nstem(j-1,X(j));holdonendaxis([0N01]);7.完整程序f=input('输入信号频率f/Hz：');%输入信号频率f%N=input('输入采样点数N：');%输入采样点数N%T=input('输入采样周期T/s：');%输入采样间隔T%%采样forj=0:1:N-1x(j+1)=sin(2*pi*f*j*T);end%补0insert=input('是否补0？[0-否；1-是]:');ifinsert==1ZERO=input('补零数：');forj=N:1:N+ZERO-1x(j+1)=0;endN=N+ZERO;end%码位倒置%M=log2(N);fort=1:1:Ns=dec2bin(t-1,M);s=fliplr(s);s=bin2dec(s);A(s+1)=x(t);end%按时间抽取的FFT蝶形运算%forL=1:1:MforJ=0:1:(2^(L-1)-1)fork=(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-4-%模值归一化%x=abs(A);y=max(x);X=x/y;%绘图forj=1:1:Nstem(j-1,X(j));holdonendaxis([0N01]);五、实验结果分析在本次试验中，进行DFT运算的序列是我们对一个在时域为连续周期的信号进行采样得到的，这个有限长序列（在本题中为32个点）可以看作是对无限长的采样得到的序列进行加窗得到的，这个窗是一个长度为32的矩形窗。由FT的性质知，在时域连续周期的信号的频域是离散非周期的。对时域信号的采样会导致频域发生周期延拖，而延拖的周期为采样频率，如果对采样后的序列进行加窗，那么会在频谱上卷积一个sinc函数。最终得到的DFT的结果是对这个已经卷积过的频谱进行频域采样得到的。a.信号频率f＝50Hz，采样点数N=32,采样间隔T=0.000625s信号的频率为50Hz，采样周期为0.000625s，则采样频率为1600Hz，在数字域，采样频率对应的是2*pi。-5-由于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-6-采样周期为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-7-采样周期为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-8-信号频率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，其它均没有变化，就是-9-说，在对频谱采样前，本题的频谱和c）的是相同的。在c中采样的间隔是pi/16,而在本题中采样的间隔是pi/32。所以在2*n点上的值与c）对应的n点的值相同。在2*n+1点上除了在第15和第49点上有峰值外，其余点为0。六、思考题1）在实验a).b).c)和d)中，正弦信号的初始相位对频谱图中的幅度特性是否有影响？为什么？由各个图象对比可以得出：正弦信号的初始相位对频谱图中的幅度特性没有影响，原因主要在于采样之后进行了周期延拓。2）信号补零后做FFT是否可以提高信号频谱的分辨率？为什么？不能，若采样频率为fs,FFT长度为N，则频谱分辨率为fs/N,设时间长度为DT，频谱分辨率＝1/DT，只要DT不变，频谱分辨率就不会变。数据后面补零可以克服栅栏效应。