实验二用FFT作谱分析Word文档格式.docx
《实验二用FFT作谱分析Word文档格式.docx》由会员分享,可在线阅读,更多相关《实验二用FFT作谱分析Word文档格式.docx(18页珍藏版)》请在冰豆网上搜索。
intL,B,J,P,k,i;
doublerPartKB,iPartKB;
doublerCf[128],iCf[128]
/*计算旋转因子*/
doublePI2=8.0*atan(1.0);
for(i=0;
i〈N;
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<
J)
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'
);
N=8;
N=16;
elseN=32;
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:
x=cos(n*pi/4)+j*sin(n*pi/8);
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),’。
’);
k’);
ylabel('
|X5(k)|'
title(’恢复后的X5(k)'
ifb==8
k
(1)=conj(f
(1));
form=2:
N
k(m)=conj(f(N-m+2));
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));
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'
b=0;
closeall;
(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的傅里叶变换,这在上面两图可以明显看出。
所以,离散傅里叶变换实际上是对该序列在频域范围内以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
通过结合x4(n)和x5(n)的频谱分析,从以上的图可以看出,将原信号的DFT变换分为共轭对称部分合共个反对称部分,则可以得出原信号的实部对应离散傅里叶变换的共轭对称部分,原信号的虚部对应离散信号的共轭反对称部分。
4。
思考题
(1)在N=8时,x2(n)和x3(n)的幅频特性会相同吗?
为什么?
N=16呢?
N=8时一样,N=16时不一样.因为DFT变换可以看成是将该序列进行周期延拓后的傅里叶级数变换的主值序列。
当N=8时,两序列进行周期延拓后序列相同,所以其傅里叶级数变换的主值序列也相同,进而DFT变换也相同。
而当N=16时,两序列进行周期延拓后序列不相同,所以其傅里叶级数变换的主值序列也相同,进而DFT变换也不相同
(2)如果周期信号的周期预先不知道,如何用FFT进行谱分析?
通过N的不同取值,然后分析频谱,从而不断缩小N的猜测范围,最后得出N值。
5。
实验结论
FFT变换即快速傅里叶变换的性质同DFT即离散傅里叶变换相同。
离散傅里叶变换有两个物理意义,一是,是对该序列的傅里叶变换w的抽样或者说对Z变换单位圆内的抽样。
二是,将该序列进行周期延拓后的傅里叶级数变换的主值序列.