试验三DFT和FFT频谱分析.docx
《试验三DFT和FFT频谱分析.docx》由会员分享，可在线阅读，更多相关《试验三DFT和FFT频谱分析.docx（11页珍藏版）》请在冰豆网上搜索。
试验三DFT和FFT频谱分析
实验三DFT和FFT频谱分析
一、实验目的
1.掌握DFT频谱分析的原理与编程方法。
2.理解FFT算法的编程思想。
2.熟练掌握利用FFT对信号作频谱分析，包括正确地进行参数选择、画频谱及读频谱图。
3.利用FFT频谱分析进行快速卷积和太阳黑子周期性检测。
二、实验环境
“.Windowsxp以上操作系统
2.安装MATLAB2007a软件
三、实验原理
1.离散傅里叶变换（DFT）
设序列为x（n）,长度为N,则
NA
X（ej3k）=DFT[x（n）]='x（n）e-jn,
n=0
2n-k
其中业=k（k=0,1,2,…-1）,通常M>N,以便观察频谱的细节。
|X（ejw）|----x（n）的幅频谱。
M
2.谱分析参数选择
1）设信号x（t）最高频率为fc,对其进行取样得x（n）,根据取样定理，取样频率fs必须满足：
fs>=2fc。
2）设谱分辨率为F,则最小记录时间tpmin=MF；取样点数
N>2fc/F;为使用快速傅里叶变换（FFT）进行谱分析，N还须满足：
N=
2E（E为整数）。
3.用FFT计算信号x（n）的频谱。
[设x（n）为实信号]
快速傅里叶变换（FFT）是DFT的一种快速算法，其使得DFT的运算速度大为加快。
1）对信号x（n）作N点FFT,得频谱X（k）（k=O~N-1）
X（k）=XR（k）+jXI（k）（k=0~N/2-1）,XR（k）—X（k）的实部；Xi（k）—X（k）的虚部。
Matlab语句：
Y=fft（x.N）
其中：
x----x（n）;Y----X（k）
/~22
2）幅频谱:
|X（k）|=xR（K）XI（K），由于x（n）为实信号，因此|X（k）|
对称，
Matlab语句：
abs（Y）
2*
iii）功率谱:
PSD（k）=|X（k）|/N=X（k）X（k）/N
Matlab语句：
PSD=Y.*conj（Y）/N
其中：
conj（Y）--X*（k）[X（k）的共轭]
4.读频谱图
频谱图中任意频率点k对应实际频率为:
fk=fs/N*k。
5.用FFT实现线性卷积运算
用FFT实现y（n）=x（n）*h（n）的步骤为：
1）设x（n）及h（n）的长度分别为N1和N?
。
为使循环卷积等于线性卷积，用补0的方
法使x（n）,h（n）长度均为N,则N须满足N》;为用FFT计算DFT,则N
还须满足N=2E（e为整数）。
2）用FFT计算X（k）,H（k）（N点）。
3）Y（k）=if;y（n）=ifft[Y（K）]
四、实验内容
1.根据公式设计DFT原理程序，并计算：
x（n）=[1,1,1,1]的4,16,64点DFT并绘图。
%DFT/IDFT程序DFT.m
clc
clearxn=input（'x（n）='）;M=length（xn）;
N=input（'变换区间N='）;xn=[xnzeros（1,N-M）];n=0:
N-1;k=0:
N-1;nk=n'*k;wn=exp（-j*2*pi/N）;wnK=wn.Ank;
xk=xn*wnK
%输入序列x（n）=[1111]
%x（n）的长度M
%变换区间N
%补0,使xn长度为N
%旋转因子wn
%作x（n）的DFT=xk
subplot（211）;stem（k,abs（xk）,'.'）;gridon;
%显示xk的幅频谱（离散曲线）
subplot（212）;plot（k,abs（xk））;gridon;
%显示xk的幅频谱（连续曲线）运行结果：
问：
由此得出怎样的结论？
答：
n越大越接近原来的dft
2.理解DIT-FFT算法原理程序，并用它计算X（k）=FFT[R4（n）],分别取N=4,8,16
和64,绘出幅频谱|X（k）|。
%程序DlT.m
clear
clcx=input（'x='）;N=input（'N='）;x（length（x）+1:
N）=zeros（1,N-length（x））;l=log2（N）;
x1=zeros（1,N）;
forj1=1:
Nend%%FFT（DIT）%%M=2;
while（Mv=N）
W=exp（-2*j*pi/M）;
V=1;
fork=0:
1:
M/2-1
fori=0:
M:
N-1p=k+i;q=p+M/2;A=x1（p+1）;B=x1（q+1）*V;x1（p+1）=A+B;x1（q+1）=A-B;
end
V=V*W;
end
M=2*M;
end
%%%%%%
subplot（211）;stem（x,'.'）;gridon;title（'x（n）'）;
subplot（212）;stem（abs（x1）,'.'）;gridon;title（'|X（k）|'）;
1.5
2.5
3.5
F
:
x（n）
|X（k）|
x（n）
1
f1
F1
1
4
4
i1
n
8
|X（k）|
1
0.5
0
4
3
2
1
0
1
0.5
0
4
3
2
1
0
x（n）
|X（k）|
x（n）
0.5
10
20
30
40
50
70
f
60
|X（k）|
3.FFT谱分析
设信号为x（t）=sin（2nf+sin（22t）+随机噪声，fi=50Hz,f2=120Hz,以取样频率
fs=1kHz对x（t）进行取样，样本长度tp=0.25s,得x（n）,对x（n）作256点FFT,得频谱X（k）,画原信号x（n）,幅频谱|X（k）|以及功率谱PSD（k）,对信号进行谱分析。
%程序pufenxi.mclear
clcfs=1000;
t=0:
1/fs:
0.25;
N=256;f1=50;f2=120;s=sin（2*pi*f1*t）+sin（2*pi*f2*t）;x=s+randn（size（t））;
Y=fft（x,N）;
PSD=Y.*conj（Y）/N;f=fs/N*（0:
N/2-1）;subplot（311）;plot（x）;subplot（312）;plot（f,abs（Y（1:
N/2）））;subplot（313）;plot（f,PSD（1:
N/2））;
150rr
100-
iii）本例的频谱分辨率F是3.9Hz,改变f2=60Hz,问：
在幅频谱中，能否分辨fi和f?
对应的频率分量？
不能。
为什么？
间隔小于频谱分辨
率。
再改变f2=52Hz，问：
在幅频谱中，能否分辨fi和f2对应的频率分量？
不能。
为什么？
间隔小于频谱分辨率。
再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在398.45Hz；
问：
在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？
不能。
为什么？
不符合采样定理。
为了正确检测f2对应的频率分量，则fs至少取多少Hz?
1200Hz。
在该程序中改变fs，验证你的结论。
突出主频点的特性。
iv）比较幅频谱和功率谱，可以发现功率谱具有
4.FFT实现任意两个序列的快速卷积。
%程序fftjuanji.m
clear
clc
x1=[111],x2=[12],回车。
结果:
y=1332
Y=conv（x1,x2）
2）问：
可用Matlab中的什么函数验算上述卷积结果?
5.利用谱分析观察太阳黑子周期性。
以100年中记录到的太阳黑子出现次数为信号x（n）,对x（n）作功率谱，从中观察太阳黑子周期
性。
%程序taiyangheizi.mclearclc
x=[101826635317209215412585683823102483...
132131118906760474121166471434454348...
4228108201512143546413024167428...
17365062677148288135712213810386633724...11154062981249666645439217423559496...
77594447301673774];%100年中太阳黑子出现的次数subplot（211）;plot（x）%画x（n）
N=128;fs=1;%fs=1Hz,N=128点
s=x-mean（x）;%对x作零均值化处理（去除直流分量）
Y=fft（s,n）;%对S故N点fft
PSD=,Y.*conj（Y）/N;%做功率谱PSD
f=fs/N*（0:
N/2-1）;%将频率定标为实际频率f
subplot（212）;plot（f,PSD（1:
N/2））;%画功率谱（N/2点）
1）填写空格中的画图语句并绘出结果图形。
4
X10
2
1.5
1
0.5
0
00.050.10.150.20.250.30.350.40.450.5
2）从s的功率谱观察到，其幅度最高处对应的横坐标f=0.086Hz,则太阳黑子每隔11.62年出现一次最高峰。
3）在对s做FFT时，为何可取fs=1Hz,N=128点？
太阳黑子每一年记录一次，记录100年，但N要取2的指数倍所以取128
升级会员 