实验二 用FFT做谱分析.docx

实验二 用FFT做谱分析.docx

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

实验二用FFT做谱分析

实验_二_题目_用FFT做谱分析_________第_14__周星期_三_第_6,7__节

一、实验目的

（1）进一步加深DFT算法原理和基本性质的理解

（2）熟悉FFT算法原理和FFT子程序的应用

（3）学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。

二、实验内容

分别对下列信号进行谱分析（采样频率为fs=64HZ）

x1（n）=R4（n）

x2（n）=[12344321]

FFT变换区间长度NN=8N=16N=32

x3（n）=[43211234]

x4（n）=cos（npi/4）

x5（n）=sin（npi/8）

FFT变换区间长NN=16N=32N=64

x6（t）=cos（8pit）+cos（16pit）+cos（20pit）

FFT变换区间长度NN=8N=16N=32

x7=x4+x5

x8=x4+jx5

三、实验程序

clear;%清除变量

closeall;%关闭全部绘图窗口

b=menu（'请选择信号x1（n）--x8（n）','x1（n）=R4（n）','x2（n）=[12344321]',...

'x3（n）=[43211234]','x4（n）=cos（npi/4）','x5（n）=sin（npi/8）',...

'x6（n）=cos（8pit）+cos（16pit）+cos（20pit）','x7=x4+x5','x8=x4+jx5','Exit'）;

i=0;

A=[8,16,32,64];

while（b~=9）%当选择EXIT时，返回值为9，则退出循环

ifb==6

m=menu（'请选择FFT变换区间长度N','N=16','N=32','N=64'）;

N=A（m+1）;

fs=64;

n=0:

（N-1）;

x=cos（8*pi*n/fs）+cos（16*pi*n/fs）+cos（20*pi*n/fs）;

else

m=menu（'请选择FFT变换区间长度N','N=8','N=16','N=32'）;

N=A（m）;

n=0:

（N-1）;

ifb==1

x=[1,1,1,1,0,0,0,0,zeros（1,N-8）];

elseifb==2

x=[1,2,3,4,4,3,2,1,zeros（1,N-8）];

elseifb==3

x=[4,3,2,1,1,2,3,4,zeros（1,N-8）];

elseifb==4

x=cos（n*pi/4）;

elseifb==5

x=sin（n*pi/8）;

elseifb==7

x=cos（n*pi/4）+sin（n*pi/8）;

elseifb==8

x=cos（n*pi/4）+j*sin（n*pi/8）;

end

end

%先画出信号源图

i=i+1;

figure（i）;%创建绘图窗口

subplot（2,2,1）;%指定1号子图

xlabel（'n'）;%标记X坐标

ifb==8

stem（n,abs（x）,'.m'）;

ylabel（'|x（n）|'）;

title（['|x8（n）|的波形']）;

else

stem（n,x,'.r'）;

ylabel（'x（n）'）;

title（['x',num2str（b）,'（n）的波形']）;

end

%进行FFT

f=fft（x,N）;

%再画出FFT波形

subplot（2,2,3）;

stem（n,abs（f）,'.b'）;

xlabel（'k'）;

ylabel（'|X（k）|'）;

title（['x',num2str（b）,'（n）的N=',num2str（N）,'点FFT']）;

ifN==16

ifb==7

subplot（2,2,2）;

stem（n,abs（real（f））,'.b'）;%Re[X7（k）]

xlabel（'k'）;

ylabel（'|Re[X7（k）]|'）;

title（'|[X7（k）]的实部|,可以与|X4（k）|对比'）;

subplot（2,2,4）;

stem（n,abs（imag（f））,'.b'）;%Im[X7（k）]

xlabel（'k'）;

ylabel（'|jIm[X7（k）]|'）;

title（'|[X7（k）]的虚部|,可以与|X5（k）|对比'）;

end

end

ifb==8

k

（1）=conj（f

（1））;

form=2:

N

k（m）=conj（f（N-m+2））;

end

fe=（f+k）/2;%求X8（k）的共轭对称分量

fo=（f-k）/2;%求X8（k）的共轭反对称分量

x8real=ifft（fe,N）;

x8imag=ifft（fo,N）/j;

i=i+1;

figure（i）;

subplot（2,2,1）;

stem（n,x8real,'.r'）;

xlabel（'n'）;

ylabel（'x8real（n）'）;

title（['x8（n）实部的波形']）;

subplot（2,2,3）;

stem（n,abs（fe）,'.b'）;

xlabel（'k'）;

ylabel（'|X8e（k）|'）;

title（['|X8（k）的共轭对称分量|波形']）;

subplot（2,2,2）;

stem（n,x8imag,'.r'）;

xlabel（'n'）;

ylabel（'x8imag（n）'）;

title（['x8（n）虚部的波形']）;

subplot（2,2,4）;

stem（n,abs（fo）,'.b'）;

xlabel（'k'）;

ylabel（'|X8o（k）|'）;

title（['|X8（k）的共轭反对称分量|波形']）;

end

b=menu（'请选择信号x1（n）--x8（n）','x1（n）=R4（n）','x2（n）=[12344321]',...

'x3（n）=[43211234]','x4（n）=cos（npi/4）','x5（n）=sin（npi/8）',...

'x6（n）=cos（8pit）+cos（16pit）+cos（20pit）','x7=x4+x5','x8=x4+jx5','Exit'）;

end

closeall;

四、实验结果

x1（n）的波形:

N=8

N=16

N=32

由3张图可知道,N值越大，频率分辨率越高。

x2（n）的波形:

N=8

N=16

N=32

由3张图可知道,N值越大，频率分辨率越高。

x3（n）的波形

N=8

N=16

N=32

由3张图可知道,N值越大，频率分辨率越高。

x4（n）的波形

N=8

根据参数可得出X4（t）的频率f=8Hz,

由于N=8,频率分辨率F0=fs/N=8Hz

因此在FFT图里的点N=1有高幅值

因为截取的为周期序列的整数倍，

所以所得出的谱正确。

N=16

根据参数可得出X4（t）的频率f=8Hz,

由于N=16,频率分辨率F0=fs/N=4Hz

因此在FFT图里的点N=2有高幅值

因为截取的为周期序列的整数倍，

所以所得出的谱正确。

N=32

根据参数可得出X4（t）的频率f=8Hz,

由于N=32,频率分辨率F0=fs/N=2Hz

因此在FFT图里的点N=4有高幅值

因为截取的为周期序列的整数倍，

所以所得出的谱正确。

x5（n）的波形

N=8

根据参数可得出X5（t）的频率f=4Hz,

由于N=8,频率分辨率F0=fs/N=8Hz

因为截取的不是为周期序列的整数倍，而且频率分辨率不够，所得出的谱有较大的误差，所以FFT图包含一些频率分量，不能清楚看清原信号的频率f.

N=16

根据参数可得出X5（t）的频率f=4Hz,

由于N=16,频率分辨率F0=fs/N=4Hz

因此在FFT图里的点N=1有高幅值

因为截取的是周期序列的整数倍，所以所得的频谱正确。

N=32

根据参数可得出X5（t）的频率f=4Hz,

由于N=32,频率分辨率F0=fs/N=2Hz

因此在FFT图里的点N=2有高幅值

因为截取的为周期序列的整数倍，所以所得出的谱正确

x6（n）的波形

N=16

根据参数可得出X6（t）里包含3个频率分别为f1=4,f2=8，f3=10

X6（n）=cos（npi/8）+cos（npi/4）+cos（5npi/32）

由于N=16,频率分辨率F0=fs/N=4Hz

因为截取的不是x6里各周期序列的整数倍，所得出的谱有频谱泄漏，FFT图里可以看出信号cos（8pit）和cos（16pit）的频率f1=4,f2=8（在点N=1，2处有较大的幅值），而且频率分辨率不够高，不能分辨开第三个信号cos（20pit）的频率f3

N=32

根据参数可得出X6（t）里包含3个频率分别为f1=4,f2=8，f3=10

由于N=32,频率分辨率F0=fs/N=2Hz

因此在FFT图里的点N=2有高幅值

在N=4有高幅值N=5也有高幅值

因为截取的为周期序列的整数倍，所以所得出的谱正确。

N=64

根据参数可得出X6（t）里包含3个频率分别为f1=4,f2=8，f3=10

由于N=64,频率分辨率F0=fs/N=1Hz

因此在FFT图里的点N=4有高幅值

在N=8有高幅值N=10也有高幅值

因为截取的为周期序列的整数倍，所以所得出的谱正确。

x7（n）的波形

N=8

N=16

N=32

x8（n）的波形

N=8

N=16

N=32

五、实验结果分析

误差产生的原因：

（1）对周期序列的截取不当，造成频谱泄漏

（2）抽样点数N太少，频率分辨率不够

用FFT做谱分析时参数的选择:

（1）抽样频率要满足奈奎斯特准则，不小于信号最高频率的2倍

（2）在抽样频率一定的情况下，抽样点数N要适当。

太小会造成频率分辨

力不够，太大会造成数据冗余。

对周期序列，最好截取周期的整数倍进行谱分析

实验结果表明:

（1）如上图Re[X7（k）]和jIm[X7（k）]的模,它们正是16点|X4（k）|和|X5（k）|

（2）|X8e（k）|=（1/2）|X8（k）+X*（N-k）|及x4（n）=x8r（n）=IDFT[X8e（k）]，如图所示，正好x4（n）16点的X4（k）及x4（n）相同。

|X8o（k）|=（1/2）|X8（k）+X*（N-k）|及x4（n）=x8i（n）=IDFT[X8o（k）]/j,如图所示,正好与x5（n）16点的X5（k）及x5（n）相同。

六、思考题

（1）在N=8时，x2（n）和x3（n）的幅频特性相同。

在N=16时，则不相同。

因为对x2（n）和x3（n）进行DFT时，原信号进行周期延拓，周期延拓后的x2（n）和x3（n）的波形相同，只是相位不一样而已。

所以他们的幅频特性相同。

而在N=16时，由于在时域数据末端添加了一些零点值。

所以在周期延拓后x2（n）和x3（n）的波形不相同，因而幅频特性也不一样。

（2）周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求。