快速傅里叶变换实验报告
班级：
姓名：
学号:
快速傅里叶变换
一．实验目的
1.在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解；
2.熟悉并掌握按时间抽取FFT 算法的程序；
3.了解应用FFT 进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT 。
二．实验内容
1.仔细分析教材第六章‘时间抽取法FFT ’的算法结构，编制出相应的用FFT 进行信号分析的C 语言（或MATLAB 语言）程序；
2.用FFT 程序分析正弦信号
()sin(2)[()(*)],(0)1y t f t u t u t N T t u π=---∞<<+∞=设
分别在以下情况进行分析并讨论所得的结果：
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=64,采样间隔T=0.000625s
f ） 信号频率f ＝250Hz ，采样点数N=32,采样间隔T=0.005s
g ） 将c ） 信号后补32个0，做64点FFT
三．实验要求
1.记录下实验内容中各种情况下的X (k)值，做出频谱图并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。
频谱只做模特性，模的最大值＝1，全部归一化；
2.打印出用C 语言（或MATLAB 语言）编写的FFT 源程序，并且在每一小段处加上详细的注释说明；
3.用C 语言（或MATLAB 语言）编写FFT 程序时，要求采用人机界面形式：
N , T , f 变量均由键盘输入，补零或不补零要求设置一开关来选择。
四．实验分析
对于本实验进行快速傅里叶变换，依次需要对信号进行采样，补零（要求补零时），码位倒置，蝶形运算，归一化处理并作图。
此外，本实验要求采用人机界面形式，N,T,F 变量由键盘输入，补零或不补零设置一开关来选择。
1.采样
本实验进行FFT 运算，给出的是正弦信号，需要先对信号进行采样，得到有限长序列()n x ， N n ......2,1,0= Matlab 实现：
t=0:T:T*(N-1); x=sin(2*pi*f*t); 2.补零
根据实验要求确定补零与否，可以用if 语句做判断，若为1，再输入补零个数， 并将补的零放到采样得到的序列的后面组成新的序列，此时新的序列的元素个数等于原采样点个数加上补零个数，并将新的序列个数赋值给N 。
Matlab 实现：
a=input('是否增加零点? 是请输入1 否请输入0\n');
if (a)
ZeroNum=input('请输入增加零点的个数:\n'); else
ZeroNum=0; end
if (a)
x=[x zeros(1, ZeroNum)];%%指令zeros(a,b)生成a 行b 列全0矩阵，在单行矩阵x 后补充0 end
N=N+ZeroNum;
3.码位倒置
本实验做FFT 变换的级数为M ，N M 2log =
做序列数对应的二进制数的码位倒置，dec2bin（）函数将十进制数转换为二进制数，fliplr（）将二进制数进行码位倒置，bin2dec（）将二进制数转换为十
进制数，并将按码位倒置得到的序列赋值为()n
A，N
=
n......
2,1,0
Matlab实现:
M=log2(N); %% M位二进制数
for t=1:1:N
s=dec2bin(t-1,M); %%将十进制数转换为二进制数，M表示二进制码位数的上限
s=fliplr(s); %%将二进制数进行码位倒置
s=bin2dec(s); %%将二进制数转换为十进制数
b=s+1; %%二进制数从0开始，而矩阵中元素序数从1开始，故需+1 A(b)=x(t);
end
4.蝶形运算
用三层for循环来实现:1.实现FFT每一级运算,共M级，此处for循环用来控制级数；2.实现分组，此处for循环用来控制旋转因子；3.实现每一组中FFT运算，此处for循环用来控制进行蝶形运算的两点之间的距离。
最终得到的()k
A即为FFT变换的结果。
Matlab实现：
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%%A(k)即为FFT变换结果
5.归一化处理及作图
实验要求对FFT运算结果进行归一化处理，对FFT运算结果序列()k
A均取绝对值得序列()k
B中所有元素均除以m，即得B，并取出绝对值中最大值m，序列()k
到归一化处理后的序列。
用stem函数即可实现作图。
Matlab实现：
%%归一化处理
B=abs(A);%%将矩阵A中元素均取绝对值，得矩阵B m=max(B);%%取矩阵B中的最大值
X=B/m; %%A（k）的幅值归一化处理之后的结果
%%作图
for i=1:1:N
stem(i-1,X(i));%%stem（A,B）表示以矩阵A中元素为纵坐标，B中元素为横坐标（一一对应）作图
hold on%%采样时间点值与元素序数相差1，故
end
axis([0 N 0 1]);%%axis限定横，纵坐标范围
五．实验结果及分析
本实验时域上加时窗，对应于频域上与sinc函数做卷积，当采样为整数倍周期时，时窗对频谱图无影响，当采样是非整数个周期时，时窗对频谱图影响较大。
f对应数字域的 2。
采样频率
s
a) 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.000625s
(2)频谱图如下：
(3)分析：
b) 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.005s
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0 0 0 0 0 0 0 0 X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15) 0-16i 0 0 0 0 0 0 0 X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23)
0 0 0 0 0 0 0 0 X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31) 0+16i 0 0 0 0 0 0 0
(3)分析：
c) 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.0046875s
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) 1.1033 1.1273 1.2050 1.3568 1.6339 2.1750 3.4960 10.2519 X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15) -10.153-3.3953-2.0703-1.5226-1.2361-1.0707-0.9739 -0.9225 X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23) -0.9063 -0.9225-0.9739-1.0707-1.2361-1.5226-2.0703-3.3953 X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31) -10.1510.2519 3.4960 2.1750 1.6339 1.3568 1.2050 1.1273
(2)频谱图如下：
(3)分析：
对于本题，若采样个数改为64
N，不补零，则有15个完整周期，调用程
序可验证仍有2根谱线，如下图：
d) 信号频率f＝50Hz，采样点数N=32,采样间隔T=0.004s
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0.9511 0.9867 -
0.0854i
1.105-
0.1829i
1.3526 -
0.3125i
1.8670 -
0.5220i
3.1952 -
0.9911i
11.383-
3.6858i
-7.844+
2.5301i
X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15)
-3.077+ 0.9511i -2.000+
0.5718i
-1.537+
0.3925i
-1.288 +
0.2826i
-1.140+
0.2045i
-1.048 +
0.1432i
-0.991+
0.0912i
-0.961+
0.0445i
X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23)
-0.9511-0.9608-
0.0445i -0.9916-
0.0912i
-1.0482-
0.1432i
-1.1405-
0.2045i
-1.2889-
0.2826i
-1.5376-
0.3925i
-2.0004-
0.5718i
X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31)
-3.0777- 0.9511i -7.8447-
2.5301i
11.383+
3.6858i
3.1952+
0.9911i
1.8670+
0.5220i
1.3526+
0.3125i
1.1052+
0.1829i
0.9867+
0.0854i
(2)频谱图如下：
(3)分析：
e) 信号频率f＝50Hz，采样点数N=64,采样间隔T=0.000625s
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0 0 -32i 0 0 0 0 0 X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15)
0 0 0 0 0 0 0 0 X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23)
0 0 0 0 0 0 0 0 X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31)
0 0 0 0 0 0 0 0 X(32)X(33)X(34)X(35)X(36)X(37)X(38)X(39)
0 0 0 0 0 0 0
X(40)X(41)X(42)X(43)X(44)X(45)X(46)X(47)
0 0 0 0 0 0 0 0 X(48)X(49)X(50)X(51)X(52)X(53)X(54)X(55)
0 0 0 0 0 0 0 0 X(56)X(57)X(58)X(59)X(60)X(61)X(62)X(63)
0 0 0 0 0 0 32i 0
(2)频谱图如下：
(3)分析：
f) 信号频率f＝250Hz，采样点数N=32,采样间隔T=0.005s
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0 0 0 0 0 0 0 0
X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15) 0-16i 0 0 0 0 0 0 0
X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23)
0 0 0 0 0 0 0 0
X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31) 0+16i 0 0 0 0 0 0 0
(2)频谱图如下：
(3)分析：
g) 将c）信号后补32个0，做64点FFT
(1)X(k)值如下表：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) 1.10330 1.12730 1.20500 1.35680
X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15) 1.6339 0 2.1750 0 3.49600 10.2519 -16.000i X(16) X(17) X(18) X(19) X(20) X(21) X(22) X(23) -10.1530 -3.3953 0 -2.07030 -1.5226 0
X(24) X(25) X(26) X(27) X(28) X(29) X(30) X(31) -1.23610 -1.07070 -0.97390 -0.9225 0
X(32)X(33)X(34)X(35)X(36)X(37)X(38)X(39) -0.90630 -0.92250 -0.97390 -1.07070
X(40)X(41)X(42)X(43)X(44)X(45)X(46)X(47)
-1.23610 -1.5226 0 -2.07030 -3.3953 0 X(48)X(49)X(50)X(51)X(52)X(53)X(54)X(55) -10.15316.0000i 10.25190 3.49600 2.17500 X(56)X(57)X(58)X(59)X(60)X(61)X(62)X(63) 1.6339 0 1.35680 1.2050 0 1.12730
(2)频谱图如下表：
(3)分析：
六．实验源程序
clc
clear
f=input('请输入信号频率: f\n');
N=input('请输入采样点数: N\n');
T=input('请输入采样间隔: T\n');
a=input('是否增加零点? 是请输入1 否请输入0\n');
%%采样，采N个点
t=0:T:T*(N-1);
x=sin(2*pi*f*t);
if(a)
ZeroNum=input('请输入增加零点的个数:\n');
else
ZeroNum=0;
end
%%补0处理:在采样点组成的单行矩阵后补充ZeroNum个0，组成新的矩阵
if (a)
x=[x zeros(1, ZeroNum)];%%指令zeros(a,b)生成a行b列全0矩阵，在单行矩阵x后补充0
end
N=N+ZeroNum;
%%码位倒置
M=log2(N); %% M位二进制数
for t=1:1:N
s=dec2bin(t-1,M); %%将十进制数转换为二进制数，M表示二进制码位数的上限
s=fliplr(s); %%将二进制数进行码位倒置
s=bin2dec(s); %%将二进制数转换为十进制数
b=s+1; %%二进制数从0开始，而矩阵中元素序数从1开始，故需+1 A(b)=x(t);
end
%%蝶形运算
%%三层for循环
%%1.实现fft每一级运算，共Ｍ级（控制级数）
%%2.控制旋转因子
%%3.实现每一组中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%%A(k)即为FFT变换结果
%%归一化处理
B=abs(A);%%将矩阵A中元素均取绝对值，得矩阵B m=max(B);%%取矩阵B中的最大值
X=B/m; %%A（k）的幅值归一化处理之后的结果
%%作图
for i=1:1:N
stem(i-1,X(i));%%stem（A,B）表示以矩阵A中元素为纵坐标，B中元素为横坐标
（一一对应）作图
hold on%%采样时间点值与元素序数相差1，故
end
axis([0 N 0 1]);%%axis限定横，纵坐标范围
七．实验总结
通过本次快速傅里叶变换实验，使我对FFT运算有了更深入的了解，让我认识到课堂上的学习仅限于理论知识的学习，而在工程实践中则会面临各种各样的问题，通过编程实现FFT运算，更是对我们编程能力的考验。
从本实验中学习到，对正弦信号进行采样，对于采得样本为整数倍周期时，频谱图仍为正弦信号的频谱，有2根谱线，而对于非整数倍周期的，频谱图会发生泄漏，实际上是不能做FFT的。
此外，本实验要求我们不仅能通过Matlab实现FFT运算，做出实验结果，更重要的是能对实验结果进行分析，作出合理的解释，从而对理论知识有更深入的理解。
在此，还要感谢各位老师的热心帮助，老师答疑解惑，认真耐心地讲解，帮助我们完成实验，让我们学习到更多的知识，感谢师恩。
