快速傅里叶变换实验报告.docx
《快速傅里叶变换实验报告.docx》由会员分享，可在线阅读，更多相关《快速傅里叶变换实验报告.docx（38页珍藏版）》请在冰豆网上搜索。
快速傅里叶变换实验报告
快速傅里叶变换
实验报告
班级：
姓名：
学号:
快速傅里叶变换
一．实验目的
1.在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解；
2.熟悉并掌握按时间抽取FFT算法的程序；
3.了解应用FFT进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT。
二．实验内容
1.仔细分析教材第六章‘时间抽取法FFT’的算法结构，编制出相应的用FFT进行信号分析的C语言（或MATLAB语言）程序；
2.用FFT程序分析正弦信号
y（t）sin（2ft）[u（t）u（tN*T）]t,设u（0）1
分别在以下情况进行分析并讨论所得的结果：
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运算，给出的是正弦信号，需要先对信号进行采样，得到有限
长序列xn，n0,1,2......N
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后补充0endN=N+ZeroNum;
3.码位倒置
本实验做FFT变换的级数为M，Mlog2N
做序列数对应的二进制数的码位倒置，dec2bin（）函数将十进制数转换为二进制数，fliplr（）将二进制数进行码位倒置，bin2dec（）将二进制数转换为十
进制数，并将按码位倒置得到的序列赋值为An，n0,1,2......N
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循环用来控制进行蝶形运算的两点之间的距离。
最终得到的Ak即为FFT变换的结果。
Matlab实现：
forL=1:
1:
M
forJ=0:
1:
（2^（L-1）-1）
fork=（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运算结果序列Ak均取绝对值
得序列Bk，并取出绝对值中最大值m，序列Bk中所有元素均除以m，即得
到归一化处理后的序列。
用stem函数即可实现作图。
Matlab实现：
%%归一化处理
B=abs（A）;%%将矩阵A中元素均取绝对值，得矩阵B
m=max（B）;%%取矩阵B中的最大值
X=B/m;%%A（k）的幅值归一化处理之后的结果
%%作图
fori=1:
1:
N
stem（i-1,X（i））;
%%stem（A,B
）表示以矩阵
A中元素为纵坐标，
B中元素为横坐标
（一一对应）作图
hold
on
%%采样时间点值与元素序数相差
1，故
end
axis（[0N01]）;
%%axis
限定横，纵坐标范围
五．实验结果及分析
本实验时域上加时窗，对应于频域上与sinc函数做卷积，当采样为整数倍周期时，时窗对频谱图无影响，当采样是非整数个周期时，时窗对频谱图影响较大。
采样频率fs对应数字域的2。
a）信号频率f＝50Hz，采样点数N=32,采样间隔T=0.000625s
（1）X（k）值如下表：
X（0）
X
（1）
X
（2）
X（3）
X（4）
X（5）
X（6）
X（7）
0
0-16i
0
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+16i
（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
（2）频谱图如下：
（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.15
10.2519
3.4960
2.1750
1.6339
1.3568
1.2050
1.1273
（2）频谱图如下：
（3）分析：
对于本题，若采样个数改为N64，不补零，则有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.95110.9867-
1.105-
1.3526-1.8670-
3.1952-
11.383-
-7.844+
0.0854i
0.1829i
0.3125i
0.5220i
0.9911i
3.6858i
2.5301i
X（8）
X（9）
X（10）
X（11）
X（12）
X（13）
X（14）
X（15）
-3.077+
-2.000+
-1.537+
-1.288+
-1.140+
-1.048+
-0.991+
-0.961+
0.9511i
0.5718i
0.3925i
0.2826i
0.2045i
0.1432i
0.0912i
0.0445i
X（16）
X（17）
X（18）
X（19）
X（20）
X（21）
X（22）
X（23）
-0.9511-0.9608-
-0.9916-
-1.0482-
-1.1405-
-1.2889-
-1.5376-
-2.0004-
0.0445i
0.0912i
0.1432i
0.2045i
0.2826i
0.3925i
0.5718i
X（24）
X（25）
X（26）
X（27）
X（28）
X（29）
X（30）
X（31）
-3.0777--7.8447-
11.383+
3.1952+
1.8670+
1.3526+
1.1052+
0.9867+
0.9511i
2.5301i
3.6858i
0.9911i
0.5220i
0.3125i
0.1829i
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.1033
0
1.1273
0
1.2050
0
1.3568
0
X（8）
X（9）
X（10）
X（11）
X（12）
X（13）
X（14）
X（15）
1.6339
0
2.1750
0
3.4960
0
10.2519
-16.000i
X（16）
X（17）
X（18）
X（19）
X（20）
X（21）
X（22）
X（23）
-10.153
0
-3.3953
0
-2.0703
0
-1.5226
0
X（24）
X（25）
X（26）
X（27）
X（28）
X（29）
X（30）
X（31）
-1.2361
0
-1.0707
0
-0.9739
0
-0.9225
0
X（32）
X（33）
X（34）
X（35）
X（36）
X（37）
X（38）
X（39）
-0.9063
0
-0.9225
0
-0.9739
0
-1.0707
0
X（40）X（41）X（42）X（43）X（44）X（45）X（46）X（47）
-1.2361
0
-1.5226
0
-2.0703
0
-3.3953
0
X（48）
X（49）
X（50）
X（51）
X（52）
X（53）
X（54）
X（55）
-10.153
16.0000i
10.2519
0
3.4960
0
2.1750
0
X（56）
X（57）
X（58）
X（59）
X（60）
X（61）
X（62）
X（63）
1.6339
0
1.3568
0
1.2050
0
1.1273
0
（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
forJ=0:
1:
（2^（L-1）-1）
fork=（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）;
X=B/m;
%%取矩阵B中的最大值
%%A（k）的幅值归一化处理之后的结果
%%作图
fori=1:
1:
N
stem（i-1,X（i））;%%stem（A,B）表示以矩阵A中元素为纵坐标，B中元素为横坐标
（一一对应）作图
hold
on
%%采样时间点值与元素序数相差
1，故
end
axis（[0N01]）;
%%axis
限定横，纵坐标范围
七．实验总结
通过本次快速傅里叶变换实验，使我对FFT运算有了更深入的了解，让我
认识到课堂上的学习仅限于理论知识的学习，而在工程实践中则会面临各种各样
的问题，通过编程实现FFT运算，更是对我们编程能力的考验。
从本实验中学习到，对正弦信号进行采样，对于采得样本为整数倍周期时，频谱图仍为正弦信号的频谱，有2根谱线，而对于非整数倍周期的，频谱图会发生泄漏，实际上是不能做FFT的。
此外，本实验要求我们不仅能通过Matlab实现FFT运算，做出实验结果，更重要的是能对实验结果进行分析，作出合理的解释，从而对理论知识有更深入的理解。
在此，还要感谢各位老师的热心帮助，老师答疑解惑，认真耐心地讲解，帮助我们完成实验，让我们学习到更多的知识，感谢师恩。
升级会员 