用FFT做谱分析

实验二用FFT做谱分析

一、实验目的

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

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

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

二、实验原理

用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2π/N≤D。可以根据此时选择FFT的变换区间N。

误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N 要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

对模拟信号的频谱时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

三、

实验内容和步骤

1、对以下典型信号进行谱分析：

⎪⎩

⎪

⎨⎧≤≤-≤≤-=⎪⎩

⎪

⎨⎧≤≤-≤≤+==其它n

n n n n n x 其它n

n n n n n x n R n x ,07

4,

330,4)(,07

4,

830,1)()

()(3241

4()cos

4

x n n π

=

5()cos(/4)cos(/8)x n n n ππ=+

6()cos8cos16cos20x t t t t πππ=++ 2、对于以上信号，x1(n)~x5(n) 选择FFT 的变换区间N 为8和16 两种情况进行频谱分析。分别打印其幅频特性曲线。 并进行对比、分析和讨论;；x6(t)为模拟周期信号，选择 采样频率Hz F s 64=，变换区间N=16,32,64 三种情况进行谱分析。分别打印其幅频特性，并进行分析和讨论。

3、令x7(n)=x4(n)+x5(n),用FFT 计算8点和16点离散傅里叶变换X （k ）=DFT[x(n)],并根据DFT 的对称性，由X(k)求出X4（k ）=DFT[x4(n)]和X5(k)=DFT[x5(n)]。

4、令x8(n)=x4+jx5(n),重复（3）。 四、实验结果及数据分析

1、实验程序：

%实验二，用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');

if b==9

b=0;

end

i=0;

close all;

while(b)

if b==6

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

if temp==1

N=16;

elseif temp==2

N=32;

else N=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'); if temp==1

N=8;

elseif temp==2

N=16;

else N=32;

end

if b==1

x=[1 1 1 1 0 0 0 0];

else if b==2

x=[1 2 3 4 4 3 2 1];

else if b==3

x=[4 3 2 1 1 2 3 4];

else if b==4

n=0:N-1;

x=cos(0.25*pi*n);

else if b==5

n=0:N-1;

x=sin((pi*n)/8);

else if b==7

n=0:N-1;

x=cos(n*pi/4)+sin(n*pi/8);

else if b==8

n=0:N-1;

x=cos(n*pi/4)+j*sin(n*pi/8);

end

end

end

end

end

end

end

end

%%TO Calculate FFT

f=fft(x,N);

i=i+1;

figure(i);

printf(x,abs(f),abs(N),abs(b));

if N==16

if b==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

if b==8

k(1)=conj(f(1));

for m=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)