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

          快速傅里叶变换实验

          姓名        

            班级       

            学号      

                   日期 2011.10.25

实验名称               快速傅里叶变换

一．实验目的

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程序计算有限长度正弦信号

分别在以下情况下所得的DFT结果并进行分析和讨论：

a）信号频率f＝50Hz，采样点数N=32,采样间隔T=0.000625s

b）信号频率f＝50Hz，采样点数N=32,采样间隔T=0.005s

c）信号频率f＝50Hz，采样点数N=32,采样间隔T=0.0046875s

d）信号频率f＝50Hz，采样点数N=32,采样间隔T=0.004s

e）信号频率f＝50Hz，采样点数N=,采样间隔T=0.000625s

f）信号频率f＝250Hz，采样点数N=32,采样间隔T=0.005s

g）  将c） 信号后补32个0，做点FFT

四．内容分析

FFT运算时，依次要进行采样，码位倒置，蝶形运算，并且要用三层循环来控制，最内层循环控制做蝶形算法，中层循环完成每一级中相同系数的蝶形算法，最外层循环控制级数递推。

1数组采样：

  for j=0:1:N-1

    x(j+1)=sin(2*pi*f*j*T);               

end

2.判断是否补零点

补零的点在数组中都存为0

insert=input('是否补0？ 1(yes)/0(no):');       %补零

if insert==1                              

    ZERO=input('补零数：');    

    for j=N:1:N+ZERO-1

        x(j+1)=0;                      

    end

    N=N+ZERO;                           

end

3.码位倒置

  M为fft运算的级数，在这里被用为表示数组中数据下表的二进制代码位数上限。用for循环来实现对数组中的每个数据进行处理。其中dec2bin（）函数将十进制代码转换为二进制代码，fliplr（）将二进制的矩阵进行左右对称的翻转，即实现了码位倒置，bin2dec（）在讲二进制代码转换为十进制代码。

  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

4.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;             

         end

     end

end

5.归一化

对采样点的幅值去绝对值，再用每一个绝对值除以其中的最大值，这样就可以得到归一化后的图形

x=abs(A);                           

y=max(x);                            

X=x/y;

6.画图

 用for循环来实现采样点的幅值与采样点的对应绘图。Axis（）规定图形的纵轴范围（0,1）和横轴范围（0，N）。

 for j=1:1:N

    stem(j-1,X(j));                 

    hold on

end

axis([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=,采样间隔T=0.000625s

信号的频率为50Hz，采样周期为0.000625s，所以采样频率为1600Hz，在数字域，采样频率对应的是2*pi，由于sin信号的频谱为在+1/16pi

和-1/16pi上的冲激，所以时域采样后，在频域以2*pi为周期进行周期延拖。加窗后，卷积一个sinc函数，由于窗的长度为，所以所以其主瓣的宽度为1/16pi，因为进行DFT有个点，所以第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，做点FFT

本题与c）题的区别在于最终做了一个点的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-1

    x(j+1)=sin(2*pi*f*j*T);               

end

%补0

insert=input('是否补0？ 1(yes)/0(no):');    

if insert==1                              

    ZERO=input('补零数：');    

    for j=N:1:N+ZERO-1

        x(j+1)=0;                      

    end

    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):2^L: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;             

        end

    end

end

%模值归一化%

x=abs(A);                           

y=max(x);                            

X=x/y;

%绘图

for j=1:1:N

    stem(j-1,X(j));                 

    hold on

end

axis([0 N 0 1]);