实验二用fft作谱分析.docx

实验二用fft作谱分析.docx

《实验二用fft作谱分析.docx》由会员分享，可在线阅读，更多相关《实验二用fft作谱分析.docx（20页珍藏版）》请在冰点文库上搜索。

实验二用fft作谱分析

实验二用FFT作谱分析

1.实验目的

（1）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。

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

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

2.实验步骤

（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。

（2）复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。

（3）编制信号产生子程序，产生以下典型信号供谱分析用：

（4）编写主程序。

/*DIT-FFT函数（C语言）*/

fft—基2DIT—FFT函数

要求：

指向复数数组指针X，FFT长度为2m，m为正整数

FFT输出结果放在输入复数数组中。

/*计算N点FFT子程序*/

/*xr：

=信号序列实部，xi：

=信号序列虚部，N：

=FFT变换区间长度N=2^M*/

/*如果信号长度小于N，应该给xr，xi后面补0*/

/*计算如果X（K）的实部和虚部分别储存在数组xr和xi中*/

VoidFft（doublexr[],doublexi[],intN,intM）

{

intL,B,J,P,k,i;

doublerPartKB,iPartKB;

doublerCf[128],iCf[128]

/*计算旋转因子*/

doublePI2=8.0*atan（1.0）;

for（i=0;i

{

rCf[i]=cos（i*PI2/N）;

iCf[i]=sin（i*PI2/N）;

}

ChangeOrder（xr,xi,N）;

/*计算各级蝶形*/

for（L=1;L<=M;L++）

{

B=（int）（pow（2,（L-1））+0.5）;

for（J=0;J<=B-1;J++）

{

P=J*（（int）（pow（2,（M-L））+0.5））;

for（k=J;k<=N-1;k+=（int）（pow（2,L）+0.5））

{

rPartKB=xr[k+B]*rCf[P]-xi[k+B]*iCf[P];

iPartKB=xi[k+B]*rCf[P]+xr[k+B]*iCf[P]

xr[k+B]=xr[k]-rPartKB;

xi[k+B]=xi[k]-iPartKB;

xr[k]=xr[k]+rPartKB;

xi[k]=xi[k]+iPartKB;

}

}

}

}

/*倒序子程序*/

voidChangeOrdor（doublexr[],doublexi[],intN）

{

intLH,N1,I,J,K;

doubleT;

LH=N/2;J=LH;N1=N–2;

for（I=1;I<=N1;I++）

{

if（I

{

T=xr[I];xr[I]=xr[J];xr[J]=T;

T=xi[I];xi[I]=xi[J];xi[J]=T;

}

K=LH;

while（J>=K）

{

J=J-K;

K=（int）（K/2+0.5）;

}

J=J+K;

}

}

（1）按实验内容要求，上机实验并写出实验报告。

本实验采用的是MATLAB语言，因此FFT子程序直接调用MATLAB语言中的FFT函数就可以实现。

下面给出完整的MATLAB程序

%实验二，用FFT做谱分析

b=menu（'请选择信号x1（n）--x8（n）','x1（n）','x2（n）','x3（n）','x4（n）','x5（n）','x6（n）','x7=x4+x5','x8=x4+jx5','Exit'）;

ifb==9

b=0;

end

i=0;

closeall;

while（b）

ifb==6

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

iftemp==1

N=16;

elseiftemp==2

N=32;

elseN=64;

end

fs=64;

n=0:

N-1;

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

else

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

iftemp==1

N=8;

elseiftemp==2

N=16;

elseN=32;

end

ifb==1

x=[11110000];

elseifb==2

x=[12344321];

elseifb==3

x=[43211234];

elseifb==4

n=0:

N-1;

x=cos（0.25*pi*n）;

elseifb==5

n=0:

N-1;

x=sin（（pi*n）/8）;

elseifb==7

n=0:

N-1;

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

elseifb==8

n=0:

N-1;

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

end

end

end

end

end

end

end

end

%%TOCalculateFFT

f=fft（x,N）;

i=i+1;

figure（i）;

printf（x,abs（f）,abs（N）,abs（b））;

ifN==16

ifb==7

k=conj（f）;

x4=（f+k）/2;%Re[X7（k）=x4（k）

figure（i+2）;

subplot（2,2,1）;

stem（abs（x4）,'.'）;

xlabel（'k'）;

ylabel（'|X4（k）|'）;

title（'恢复后的X4（k）'）;

x5=（f-k）/2;%jIm[X7（k）=X5（k）

subplot（2,2,3）;

Stem（abs（x5）,'.'）;

xlabel（'k'）;

ylabel（'|X5（k）|'）;

title（'恢复后的X5（k）'）;

end

ifb==8

k

（1）=conj（f

（1））;

form=2:

N

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

end

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

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

xr=ifft（fe,N）;%xr=x4（n）

b=4;

figure（i+1）

printf（xr,abs（fe）,abs（N）,abs（b））;

xi=ifft（fo,N）/j;%xi=x5（n）

b=5;

figure（i+2）

printf（xi,abs（f）,abs（N）,abs（b））;

end

end

b=menu（'请选择信号x1（n）--x8（n）','x1（n）','x2（n）','x3（n）','x4（n）','x5（n）','x6（n）','x7=x4+x5','x8=x4+jx5','Exit'）;

ifb==9

b=0;

end

closeall;

end

（5）按实验内容要求，上机实验，并写出实验报告。

3.上机实验内容

（1）对2中所给出的信号逐个进行谱分析。

下面给出针对各信号的FFT变换区间N以及对连续信号x6（t）的采样频率Fs，供实验时参考。

x1（n），x2（n），x3（n），x4（n），x5（n）：

N=8，16

x6（t）：

Fs=64Hz，N=16，32，64

x1（n）

N=8N=16

由以上两个图分析如下：

离散傅里叶变换的N点变换在频域范围内表现为对傅里叶变换即Z变换在单位圆上的抽样。

所以N取8点时，k=0,1,2,3,4,5,6,7与N取16点时，k=0,2,4,6,8,10,12,14的离散傅里叶变换值对应相等，即它们都等于原信号在w=0、

/8、

/4、3

/8、4

/8、5

/8、6

/8、6

/8的傅里叶变换，这在上面两图可以明显看出。

所以，离散傅里叶变换实际上是对该序列在频域范围内以2

/N的间隔进行抽样。

x2（n）N=8N=16

x3（n）N=8N=16

X4（n）N=8N=16

从以上两组图可以看出，若按X1（n）进行分析，则明显不对。

现分析如下：

原信号周期为16，所以当N=8时，未能取完一个周期的值，N=16则取完了一个周期的值，所以这是两个不同的序列，所以按照X1（n）的分析方式是不对的，因为本身它们的傅里叶变换就是不一样的。

由于离散傅里叶变换是该序列周期延拓后所对应的傅里叶级数变换的主值序列，所以，当N=16时，所得的DFT值与X5（n）的傅里叶级数变换的主值序列是一致的，而N=8时是X5（n）的部分序列的周期延拓后的傅里叶级数变换的主值序列，因此两者的值是不同的。

X6（n）N=16N=32

N=64

原连续信号的周期为0.5，当采样频率Fs=64Hz时，所形成的序列周期为0.5*64=32。

所以只有N>32,才能取完一个周期的序列。

这一点，从上面三个图可以清晰看出。

其中N=32和N=16的图形分析，可以参考x4（n）的分析。

N=32和N=64的图形分析，可以参考x5（n）的分析。

（2）令x（n）=x4（n）+x5（n），用FFT计算8点和16点离散傅里叶变换。

X（k）=DFT[x（n）]

X7N=8

N=16

两者都不能显示一个周期内的所有序列，尽管N=8时的序列是N=16时序列的一部分，但是它们确属两个不同的序列。

所以它们的傅里叶变换不同，即不能按照x1（n）的进行分析而且周期延拓后所取得傅里叶级数的主值序列不同，即DFT变换值不同。

（3）令x（n）=x4（n）+jx5（n），重复

（2）。

X8N=8

N=16

通过结合x4（n）和x5（n）的频谱分析，从以上的图可以看出，将原信号的DFT变换分为共轭对称部分合共个反对称部分，则可以得出原信号的实部对应离散傅里叶变换的共轭对称部分，原信号的虚部对应离散信号的共轭反对称部分。

4.思考题

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

为什么？

N=16呢？

　　Ｎ＝８时一样，Ｎ＝１６时不一样。

因为DFT变换可以看成是将该序列进行周期延拓后的傅里叶级数变换的主值序列。

当N=8时，两序列进行周期延拓后序列相同，所以其傅里叶级数变换的主值序列也相同，进而DFT变换也相同。

而当N=16时，两序列进行周期延拓后序列不相同，所以其傅里叶级数变换的主值序列也相同，进而DFT变换也不相同

（2）如果周期信号的周期预先不知道，如何用FFT进行谱分析？

通过Ｎ的不同取值，然后分析频谱，从而不断缩小Ｎ的猜测范围，最后得出Ｎ值。

5.实验结论

FFT变换即快速傅里叶变换的性质同DFT即离散傅里叶变换相同。

离散傅里叶变换有两个物理意义，一是，是对该序列的傅里叶变换w的抽样或者说对Z变换单位圆内的抽样。

二是，将该序列进行周期延拓后的傅里叶级数变换的主值序列。