实验三DFT和FFT频谱分析.docx
《实验三DFT和FFT频谱分析.docx》由会员分享，可在线阅读，更多相关《实验三DFT和FFT频谱分析.docx（13页珍藏版）》请在冰豆网上搜索。
实验三DFT和FFT频谱分析
实验三DFT和FFT频谱分析
一、实验目的
1.掌握DFT频谱分析的原理与编程方法。
2.理解FFT算法的编程思想。
2.熟练掌握利用FFT对信号作频谱分析，包括正确地进行参数选择、画频谱
及读频谱图。
3.利用FFT频谱分析进行快速卷积和太阳黑子周期性检测。
二、实验环境
1.Windowsxp以上操作系统
2.安装MATLAB2007a软件
三、实验原理
1.离散傅里叶变换（DFT）
设序列为x（n）,长度为N,则
X（ejωk）=DFT[x（n）]=
x（n）e-jωkn,
其中ωk=
（k=0,1,2,…,M-1）,通常M>N,以便观察频谱的细节。
|X（ejωk）|----x（n）的幅频谱。
矚慫润厲钐瘗睞枥庑赖賃軔。
矚慫润厲钐瘗睞枥庑赖賃。
2.谱分析参数选择
1）设信号x（t）最高频率为fc,对其进行取样得x（n）,根据取样定理，取样频率fs必须满足：
fs>=2fc。
聞創沟燴鐺險爱氇谴净祸測。
聞創沟燴鐺險爱氇谴净祸。
2）设谱分辨率为F,则最小记录时间tpmin=1/F；取样点数
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=0~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）
2）幅频谱:
|X（k）|=
，由于x（n）为实信号，因此|X（k）|对称，
Matlab语句：
abs（Y）
iii）功率谱:
PSD（k）=|X（k）|2/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和N2。
为使循环卷积等于线性卷积，用补0的方
法使x（n）,h（n）长度均为N，则N须满足N≥N1+N2-1；为用FFT计算DFT,则N
还须满足N=2E（E为整数）。
2）用FFT计算X（k）,H（k）（N点）。
3）Y（k）=ifft;y（n）=ifft[Y（K）]。
四、实验内容
1.根据公式设计DFT原理程序，并计算：
x（n）=[1,1,1,1]的4,16,64点DFT并绘图。
%DFT/IDFT程序DFT.m
clc
clear
xn=input（'x（n）='）;%输入序列x（n）=[1111]
M=length（xn）;%x（n）的长度M
N=input（'变换区间N='）;%变换区间N
xn=[xnzeros（1,N-M）];%补0,使xn长度为N
n=0:
N-1;k=0:
N-1;
nk=n'*k;
wn=exp（-j*2*pi/N）;%旋转因子wn
wnK=wn.^nk;
xk=xn*wnK%作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）|。
%程序DIT.m
clear
clc
x=input（'x='）;%输入序列
N=input（'N='）;%做fft的点数
x（length（x）+1:
N）=zeros（1,N-length（x））;%补0x（1:
N）酽锕极額閉镇桧猪訣锥顧荭。
酽锕极額閉镇桧猪訣锥顧。
l=log2（N）;
x1=zeros（1,N）;
forj1=1:
N%倒序
x1（j1）=x（bin2dec（fliplr（dec2bin（j1-1,l）））+1）;
end
%%FFT（DIT）%%
M=2;
while（M<=N）
W=exp（-2*j*pi/M）;%旋转因子W
V=1;
fork=0:
1:
M/2-1%k为每级蝶形运算旋转因子的个数
fori=0:
M:
N-1%i为各群的首序号
p=k+i;
q=p+M/2;
A=x1（p+1）;
B=x1（q+1）*V;
x1（p+1）=A+B;%本级蝶形运算，x1最终存放X（k）
x1（q+1）=A-B;
end
V=V*W;%旋转因子W的变化
end
M=2*M;%第M级
end
%%%%%%
subplot（211）;stem（x,'.'）;gridon;%画图
title（'x（n）'）;%标题
subplot（212）;stem（abs（x1）,'.'）;gridon;%画图
title（'|X（k）|'）;%标题
3.FFT谱分析
设信号为x（t）=sin（2πf1t）+sin（2πf2t）+随机噪声，f1=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.m
clear
clc
fs=1000;
t=0:
1/fs:
0.25;%时间范围
N=256;%做fft的点数
f1=50;f2=120;%信号频率
s=sin（2*pi*f1*t）+sin（2*pi*f2*t）;%产生x（n）
x=s+randn（size（t））;%信号+噪声x（n）厦礴恳蹒骈時盡继價骚卺癩。
厦礴恳蹒骈時盡继價骚卺。
Y=fft（x,N）;%对x做N点fft
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）））;%画幅度谱（N/2点）鹅娅尽損鹌惨歷茏鴛賴縈诘。
鹅娅尽損鹌惨歷茏鴛賴縈。
subplot（313）;plot（f,PSD（1:
N/2））;%画功率谱（N/2点）籟丛妈羥为贍偾蛏练淨槠挞。
籟丛妈羥为贍偾蛏练淨槠。
画出图形窗口显示的图形，并注名每个图形的含义。
2）回答下列问题：
i）观察幅频谱图，可以发现，信号x（n）含有的两个频率分量分别是
50.8Hz和121.1Hz。
ii）在该程序中的“f=fs/N*（0:
N/2-1）;”下添加“k=0:
N/2-1;”，
“plot（f,abs（Y（1:
N/2）））;”改为“plot（k,abs（Y（1:
N/2）））;”
重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是
13和31。
它们的含义为主频点。
再将该程序中的N改为512,重新运行该程序并观察幅频谱图，这时图中
两峰值对应的下标分别是26和61。
结果是否和上面的相同？
不同
为什么？
N不同。
iii）本例的频谱分辨率F是3.9Hz,改变f2=60Hz,问：
在幅频谱中，能否分辨f1和f2
对应的频率分量？
不能。
为什么？
间隔小于频谱分辨率。
再改变f2=52Hz,问：
在幅频谱中，能否分辨f1和f2对应的频率分量？
不能。
为什么？
间隔小于频谱分辨率。
再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在398.45Hz；
问：
在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？
不能。
为什么？
不符合采样定理。
为了正确检测f2对应的频率分量，则fs至少取多少Hz?
1200Hz。
在该程序中改变fs，验证你的结论。
預頌圣鉉儐歲龈讶骅籴買闥。
預頌圣鉉儐歲龈讶骅籴買。
iv）比较幅频谱和功率谱，可以发现功率谱具有突出主频点的特性。
4.FFT实现任意两个序列的快速卷积。
%程序fftjuanji.m
clear
clc
x1=input（'x1='）;x2=input（'x2='）;%输入序列
N1=length（x1）;N2=length（x2）;%序列x1（n）,x2（n）的长度
E=ceil（log2（N1+N2-1））;%ceil---向+∞方向取整
N=2^E;%做FFT的点数
x1=[x1,zeros（1,N-N1）];%补零使x1,x2长度为N
x2=[x2,zeros（1,N-N2）];
X1=fft（x1,N）;%对x1做N点的fft
X2=fft（x2,N）;
Y=X1.*X2;%数列X1和X2的乘积
y=ifft（Y,N）%对y做N点的ifft
结果分析：
1）回到MATLAB窗口，键入：
x1=[111],x2=[12],回车。
结果:
y=1332
2）问：
可用Matlab中的什么函数验算上述卷积结果？
Y=conv（x1,x2）
5.利用谱分析观察太阳黑子周期性。
以100年中记录到的太阳黑子出现次数为信号x（n）,对x（n）作功率谱，从中观察太阳黑子周期性。
%程序taiyangheizi.m
clear
clc
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点）
填写空格中的画图语句并绘出结果图形。
2）从s的功率谱观察到，其幅度最高处对应的横坐标f=0.086Hz,
则太阳黑子每隔11.62年出现一次最高峰。
3）在对s做FFT时，为何可取fs=1Hz,N=128点？
太阳黑子每一年记录一次，记录100年，但N要取2的指数倍所以取128
升级会员 