DFT和FFT实验上传.docx
- 文档编号:25524007
- 上传时间:2023-06-09
- 格式:DOCX
- 页数:21
- 大小:170.50KB
DFT和FFT实验上传.docx
《DFT和FFT实验上传.docx》由会员分享,可在线阅读,更多相关《DFT和FFT实验上传.docx(21页珍藏版)》请在冰豆网上搜索。
DFT和FFT实验上传
DFT和FFT实验
一、实验目的和要求
1、掌握DFT变换
2、掌握DFT性质
3、掌握快速傅立叶变换(FFT)
二、实验内容和原理
1、实验内容
1)求有限长离散时间信号的离散时间傅立叶变换
并绘图。
∙已知
∙已知
2)已知序列
,
,绘制
及其离散傅立叶变换
的幅度、相位图。
3)设
,
其中,randn(n)为高斯白噪声。
求出
,m=2,3,4的matlab采用不同算法的执行时间。
4)研究高密度频谱和高分辨率频谱。
设有连续信号
∙以采样频率
对信号x(t)采样,分析下列三种情况的幅频特性。
∙采集数据长度N=16点,做N=16点的FFT,并画出幅频特性。
∙采集数据长度N=16点,补零到256点,做N=256点的FFT,并画出幅频特性。
∙采集数据长度N=256点,做N=256点的FFT,并画出幅频特性。
观察三种不同频率特性图,分析和比较它们的特点以及形成的原因。
2、实验原理
1)DFT
序列x(n)的离散时间傅里叶变换(DTFT)表示为
,
如果x(n)为因果有限长序列,n=0,1,...,N-1,则x(n)的DTFT表示为
x(n)的离散傅里叶变换(DFT)表达式为
序列的N点DFT是序列DTFT在频率区间[0,2π]上的N点灯间隔采样,采样间隔为2π/N。
通过DFT,可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。
X(k)的幅度谱为
,其中下标R和I分别表示取实部、虚部的运算。
X(k)的相位谱为
。
离散傅里叶反变换(IDFT)定义为
。
2)FFT
快速傅里叶变换(FFT)是DFT的快速算法,并不是一个新的映射。
FFT利用了
函数的周期性和对称性以及一些特殊值来减少DFT的运算量,可使DFT的运算量下降几个数量级,从而使数字信号处理的速度大大提高。
若信号是连续信号,用FFT进行谱分析时,首先必须对信号进行采样,使之变成离散信号,然后就可以用FFT来对连续信号进行谱分析。
为了满足采样定理,一般在采样之前要设置一个抗混叠低通滤波器,且抗混叠滤波器的截止频率不得高于与采样频率的一半。
比较DFT和IDFT的定义,两者的区别仅在于指数因子的指数部分的符号差异和幅度尺度变换,因此可用FFT算法来计算IDFT。
三、主要仪器设备
Matlab
四、操作方法和实验步骤
1、认真分析原函数,取点
2、用matlab编写程序,运行程序得出结果
五、实验数据记录、处理和分析
1、求有限长离散时间信号的离散时间傅立叶变换
并绘图。
∙已知
∙已知
【解答】
思路:
这是一道DFT的题,按照题目要求只需要取11个点即可。
第
(1)小题
M文件源代码
N=11;%取点个数为11个
j=sqrt(-1);%定义j为复数单位
f=inline('(0.9*exp(j*pi/3))^n','n');%定义一个函数f(n)
W=0:
2*pi/1000:
2*pi;%定义离散域的基本频率W为数组,间距为2*pi/1000
Xw=zeros(size(W));%定义一个与W位数相等的数组
forn=0:
N-1
Xw=Xw+f(n)*exp(-j*W*n);
end%对f(n)函数做DFT变换
xn=[];
forn=0:
N-1
xn(n+1)=f(n);
end%将f(n)的值放进数组xn里面,便于最后画出xn的图像
magXw=abs(Xw);%定义一个数组magXw,将abw(Xw)的值赋给它
angleXw=angle(Xw);%定义数组angleXw,将angle(Xw)的值赋给它
figure
(1);
plot(xn,'.-');
xlabel('n');ylabel('x(n)');%画出xn的图
figure
(2);
k=0:
1:
N-1;
plot(W,magXw,'-')
xlabel('W');ylabel('|X(W)|');%画出magXw的图像
figure(3);
plot(W,angleXw,'-');
xlabel('W');ylabel('angle(X(W))');%画出angleXw的图像
运行结果
xn图像:
X(W)的幅度图
X(W)的相位图
【分析】
可见
的幅度频谱有11-1=10个极大,11-1=10个极小。
而
的相位则有11-1=10个极大,11-1=10个极小,并且相位在-
和
之间摆动。
第
(2)小题
M文件源代码
N=10;%取点个数为11个
j=sqrt(-1);%定义j为复数单位
f=inline('2^n','n');%定义一个函数f(n)
W=0:
2*pi/1000:
2*pi;%定义离散域的基本频率,将其设置为间距为2*pi/1000的数组
Xw=zeros(size(W));%定义一个数组Xw,位数与W相等
forn=-N:
N
Xw=Xw+f(n)*exp(-j*W*n);
end%对f(n)函数做DTFT变换
xn=[];%定义一个数组xn
forn=-N:
N
xn(n+1+N)=f(n);
end%将f(n)的值放进数组xn里面,便于最后画出xn的图像
magXw=abs(Xw);%将X(W)的模值放进数组magXw
angleXw=angle(Xw);%将X(W)的相位放进数组magXw
figure
(1);
plot(xn,'.-');
xlabel('n');ylabel('x(n)');%画出xn的图
figure
(2);
plot(W,magXw,'-')
xlabel('W');ylabel('|X(w)|');%画出magXw的图像
figure(3);
plot(W,angleXw,'-');
xlabel('W');ylabel('angle(X(w))');%画出angleXw的图像
xn图像
X(W)的幅度图
X(W)的相位图
【分析】
可见
的幅度有一个极大值,一个极小值。
的相位在
之间来回振动,并且中间出现突变的情况。
2)已知序列
,
,绘制
及其离散傅立叶变换
的幅度、相位图。
【解答】
思路:
这是一道DFT的题,按照题目要求只需要取51个点即可。
M文件源代码
N=51;%取点个数为50个
j=sqrt(-1);%定义j为复数单位
f=inline('cos(0.82*pi*n)+2*sin(pi*n)','n');%定义一个函数f(n)
Xk=[];%定义一个数组Xk
W=2*pi/N;%定义离散域的基本频率
fork=1:
N
Xk(k)=0;
forn=0:
N-1
Xk(k)=Xk(k)+f(n)*exp(-j*(k-1)*W*n);
end
end%对f(n)函数做DFT变换
xn=[];%定义一个数组xn
forn=0:
N-1
xn(n+1)=f(n);
end%将f(n)的值放进数组xn里面,便于最后画出xn的图像
magXk=[];%定义一个数组magXk
fork=1:
N
magXk(k)=abs(Xk(k));
end%将X(kW)的模值放进数组magXk
angleXk=[];%定义数组angleXk
fork=1:
N
angleXk(k)=angle(Xk(k));
end%将X(kW)的相位放进数组magXk
figure
(1);
plot(xn,'.-');
xlabel('n');ylabel('x(n)');%画出xn的图
figure
(2);
k=0:
1:
N-1;
plot(k,magXk,'+-')
xlabel('k');ylabel('|X(k)|');%画出magXk的图像
figure(3);
plot(k,angleXk,'x-');
xlabel('k');ylabel('angle(X(k))');%画出angleXk的图像
命令窗口中的运行及其结果:
xn图像
Xk的幅度图
Xk的相位图
【分析】
可见
的幅度频谱拥有两个峰值,
的相位频谱在
之间来回振动,且中间存在3个台阶式的向下跳变,一个台阶式的向上跳变。
3.设
,
其中,randn(n)为高斯白噪声。
求出
,m=2,3,4的matlab采用不同算法的执行时间。
【解答】
思路:
计算DFT算法和FFT算法的运行时间可以使用的etime函数
M文件源代码
m=input('m=:
');%输入m值
N=4^m;%求出所取得x(n)的点数
j=sqrt(-1);%定义j为复数单位
arr=[];%定义一个数组arr
W=2*pi/N;%定义离散域的基本频率
dft_time=0;%定义dft_time为0
t1=clock;%此处为dft计算的时间起点
fork=1:
N
arr(k)=0;
forn=0:
N-1
arr(k)=arr(k)+(sin(0.2*pi*n)+rand
(1))*exp(-j*(k-1)*W*n);
end
end%对x(n)做dft变换
dft_time=etime(clock,t1)%得出dft变换所花的时间
Wn=exp(-j*2*pi/N);%求出旋转因子
xn=zeros(1,N);%定义一个N位数组
fft_time=0;%定义fft_time为0
t2=clock;%此处为fft计算的时间起点
forn=0:
N-1
xn(n+1)=sin(0.2*pi*n)+randn
(1);%将x(n)的值放入数组xn中
end
n1=fliplr(dec2bin([0:
N-1]));%码位倒置步骤1:
将码位转换为二进制,再进行倒序
n2=[bin2dec(n1)];%码位倒置步骤2:
将码位转换为十进制后翻转
x=zeros(2*m+1,N);%定义一个(2m+1)×N的矩阵
fori=1:
N
x(1,i)=xn(n2(i)+1);
end%将码位倒序后的值重新赋值进入矩阵第一列
fori=1:
2*m%进行第v级计算
Number=2^(2*m-i);
Interval_of_Unit=2^(i-1);%每组中每个计算单元的间距
Interval_of_Group=2^i;%每组之间的间距
Wnr=[];%定义一个新数组Wnr
forr=1:
2^(i-1)
Wnr(r)=Wn^((r-1)*N/2^i);
end%将每一级运算的指数因子赋值给Wnr
fork=0:
Number-1
forl=1:
2^(i-1)
x(i+1,l+k*2^i)=x(i,l+k*2^i)+Wnr(l)*x(i,l+k*2^i+2^(i-1));
x(i+1,l+k*2^i+2^(i-1))=x(i,l+k*2^i)-Wnr(l)*x(i,l+k*2^i+2^(i-1));
end
end
end%对x(n)做fft变换
fft_time=etime(clock,t2)%得出dft变换所需要的时间
命令窗口中的运行及其结果
m=:
2
dft_time=
0
fft_time=
0
m=:
3
dft_time=
.022*********
fft_time=
0.00400000000000
m=:
4
dft_time=
0.31900000000000
fft_time=
0.01800000000000
【分析】
从实验结果可以看出,对于同样的x(n),FFT变换的计算时间小于DFT变换的计算时间,并且随着取样点数的增多,这种差距越来越明显。
这是由于FFT算法利用了旋转因子的对称性和周期性,将长序列分解为短序列,分级进行蝶形运算使得复数乘法的次数减少到分解前的一半。
通过对短序列的计算进行适当的组合,从而达到了删除重复运算,提供运算速度的目的。
4)研究高密度频谱和高分辨率频谱。
设有连续信号
∙以采样频率
对信号x(t)采样,分析下列三种情况的幅频特性。
∙采集数据长度N=16点,做N=16点的FFT,并画出幅频特性。
∙采集数据长度N=16点,补零到256点,做N=256点的FFT,并画出幅频特性。
∙采集数据长度N=256点,做N=256点的FFT,并画出幅频特性。
观察三种不同频率特性图,分析和比较它们的特点以及形成的原因。
【解答】
思路:
此处为FFT变换,要注意采样数据长度N和FFT变换点数M不同,因此当M>N时,应该将x(n)中非0-N-1范围内的点补零。
M文件源代码
N=input('N=:
');%输入采样数据长度N
M=input('M=:
');%输入FFT变换点数M
fs=32000;v=log2(M);%给采样频率fs赋值为32000,v为FFT变换的级数
Wn=exp(-j*2*pi/M);%定义旋转因子Wn
j=sqrt(-1);%定义复数单位j
xn=zeros(1,M);%定义一个M位数组xn
forn=0:
N-1
xn(n+1)=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);
end%将x(n)的值赋值进入xn中
n1=fliplr(dec2bin([0:
M-1]));%码位倒置步骤1:
将码位转换为二进制,再进行倒序
n2=[bin2dec(n1)];%码位倒置步骤2:
将码位转换为十进制后翻转
x=zeros(v+1,M);%定义一个(v+1)×M的矩阵x
fori=1:
M
x(1,i)=xn(n2(i)+1);
end%将码位倒序后的x(n)值赋值进入矩阵x的第一行
fori=1:
v%进行第v级计算
Number=2^(v-i);%每一级计算中的“群”数
Interval_of_Unit=2^(i-1);%每组中每个计算单元的间距
Interval_of_Group=2^i;%每组之间的间距
Wnr=[];%定义一个数组Wnr
forr=1:
2^(i-1)
Wnr(r)=Wn^((r-1)*M/2^i);
end%将每一级运算的指数因子赋值给Wnr
fork=0:
Number-1
forl=1:
2^(i-1)
x(i+1,l+k*2^i)=x(i,l+k*2^i)+Wnr(l)*x(i,l+k*2^i+2^(i-1));
x(i+1,l+k*2^i+2^(i-1))=x(i,l+k*2^i)-Wnr(l)*x(i,l+k*2^i+2^(i-1));
end
end
end%对x(n)做fft变换
Xk=[];%定义一个数组Xk
fork=1:
M
Xk(k)=x(v+1,k);
end%将矩阵x第v+1行数赋值给数组
figure
(1);
k=0:
1:
M-1;
stem(k,Xk,'.');%画出X(k)的频率特性图
命令窗口中的运行及其结果
N=16,M=16时
N=16,M=256时
N=16,M=256时
【分析】
N=16,M=16时,频谱密度与分辨率明显小于N=16,M=256和N=256,M=256时的频谱。
这是由于采样的点数过少,而且FFT所取的点数太少。
N=16,M=256时,频谱密度比N=16,M=16时的大,和N=256,M=256时的频谱密度相同。
这是由于采样的点数过少,FFT所取的点数较多。
N=256,M=256时,频谱密度比N=16,M=16时的大,与N=16,M=256时的频谱密度相同,分辨率较高。
这是由于采样点数较多,FFT所取的点数较多。
在采样频率不变的情况下,增加采样点数为原采样点数的整数倍可以较好的还原得到原图像。
而依靠补零的方式进行计算,可能会导致严重的频率泄露,破坏原信号的周期性。
六、讨论、心得
1、通过本次实验,我学会了matlab的基本使用方法,掌握DFT运算和FFT蝶形运算的基本原理。
2、在进行蝶形运算时,涉及到分级运算,以及在每一级运算中,需要对每一个运算单元进行分组,并完成每一组的运算。
这需要十分巧妙地设置循环机构和循环变量。
一开始没有找到合适的循环方式,导致无法用程序完成蝶形运算。
3、在进行蝶形运算时,需要先对原输入序列进行倒序码位运算,这需要十分复杂的程序。
后来通过dec2bin和bin2dec函数完成了十进制和二进制之间的相互转换,又通过fliplr函数完成了二进制码位倒置。
(注:
可编辑下载,若有不当之处,请指正,谢谢!
)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- DFT FFT 实验 上传