实验三DFT和FFT频谱分析.docx

实验三DFT和FFT频谱分析.docx

《实验三DFT和FFT频谱分析.docx》由会员分享，可在线阅读，更多相关《实验三DFT和FFT频谱分析.docx（12页珍藏版）》请在冰豆网上搜索。

实验三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点）

1）画出图形窗口显示的图形，并注名每个图形的含义。

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点）

1）填写空格中的画图语句并绘出结果图形。

2）从s的功率谱观察到，其幅度最高处对应的横坐标f=0.086Hz,

则太阳黑子每隔11.62年出现一次最高峰。

3）在对s做FFT时，为何可取fs=1Hz,N=128点？

太阳黑子每一年记录一次，记录100年，但N要取2的指数倍所以取128