数字信号处理 用FFT作谱分析 第二次实验
数字信号实验报告
实验项目名称：用FFT作谱分析
所属课程名称：数字信号处理
实验类型：综合型
指导教师：
实验日期：2013.11.26 班级：
学号：
姓名：
目录
一、实验目的 (1)
二、实验步骤 (1)
三、上机实验内容 (2)
1.程序
2.运行截图
3.注释
四、思考题 (11)
一、 实验目的：
（1） 进一步加深DFT 算法原理和基本性质的理解（因为FFT 只是DFT 的一种快
速算法，所以FFT 的运算结果必然满足DFT 的基本性质）。 （2） 熟悉FFT 算法原理和FFT 子程序的应用。
（3） 学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现
的分析误差及其原因，以便在实际中正确应用FFT 。
二、 实验步骤：
（1） 复习DFT 的定义、性质和用DFT 作谱分析的有关内容。
（2） 复习FFT 算法原理与编程思想，并对照DIT —FFT 运算流图和程序框图，
读懂本实验提供的FFT 子程序。
（3） 编制信号产生子程序，产生以下典型信号供谱分析用：
()1x n =R4(n)
()21,038,470,n n x n n n +<=<=??
=-<=<=???其它n
()34,033,470,n n x n n n -<=<=??
=-<=<=???其它n
()4x n =cos(pi/4*n) ()5x n =sin(pi/8*n)
()6x t =cos(pi*8*t)+ sin(pi*16*t)+cos(20*pi*t)
应当注意，如果给出的是连续信号xa(t),则首先要根据其最高频率确定采样速率fs 以及由频率选择采样点数N ，然后对其进行软件采样（即计算x(n)=xa(nT),0<=n<=N-1）,产生对应序列x(n)。对信号x6(t)，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行分析，否则有可能产生较大的分析误差。请实验者根据DFT 的隐含周期性思考这个问题。 （4） 编写主程序。
下图给出了主程序框图，供参考。本实验提供FFT 子程序和通用绘图子程序。
三、上机实验内容
（1） 对2中所给出的信号逐个进行谱分析。下面给出针对各信号的FFT 变换区
间N 以及对连续信号()6x t 的采样频率fs 。
()1x n ,()2x n ,()3x n ,()4x n ,()5x n ,：N=8，16 ()6x t ：fs=64（hz ）, N=16,32,64
实验结果： 1．()1x n =R4(n) 原程序： n=[0:7];
x=[1 1 1 1 0 0 0 0] f1=fft(x,8) f2=fft(x,16) subplot(2,2,1) stem(n,x);
axis([0 8 0 2]) xlabel('n')
ylabel('x1(n)') title('x1的波形') subplot(2,2,4) k=[0:15] stem(k,abs(f2)); axis([0 16 0 5]) xlabel('k')
ylabel('|x1(k)|')
title('x1(n)的8点fft') subplot(2,2,3) k=[0:7]
stem(k,abs(f1)); axis([0 10 0 5]) xlabel('k')
ylabel('|x1(k)|')
title('x1(n)的8点fft')
得到的波形图如下：
分析：对8 点和16 点FFT 变换分析如下： 离散傅里叶变换的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、7π /8 的傅里叶变换，这在上面两图可以明显看出。所以，离散傅里叶变换实际上是对该序列在频域范围内以2π /N 的间隔进行抽样。
2．()21,038,470,n n x n n n +<=<=??
=-<=<=???
其它n
原程序：n=[0:7];
x=[1 2 3 4 4 3 2 1] f1=fft(x,8) f2=fft(x,16) subplot(2,2,1) stem(n,x);
axis([0 8 0 4]) xlabel('n')
ylabel('x2(n)') title('x2的波形') subplot(2,2,4) k=[0:15] stem(k,abs(f2)); axis([0 16 0 20]) xlabel('k')
ylabel('|x2(k)|')
title('x2(n)的8点fft') subplot(2,2,3) k=[0:7]
stem(k,abs(f1)); axis([0 10 0 20]) xlabel('k')
ylabel('|x2(k)|')
title('x2(n)的8点fft')
分析：对8 点和16 点FFT 变换分析如下： 离散傅里叶变换的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、7π /8 的傅里叶变换，这在上面两图可以明显看出。所以，离散傅里叶变换实际上是对该序列在频域范围内以2π /N 的间隔进行抽样。
3．()34,033,470,n n x n n n -<=<=??
=-<=<=???
其它n
原程序：n=[0:7];
x=[4 3 2 1 1 2 3 4] f1=fft(x,8) f2=fft(x,16) subplot(2,2,1) stem(n,x);
axis([0 8 0 4]) xlabel('n')
ylabel('x3(n)') title('x3的波形') subplot(2,2,4) k=[0:15]
stem(k,abs(f2)); axis([0 16 0 20]) xlabel('k')
ylabel('|x3(k)|')
title('x3(n)的8点fft') subplot(2,2,3) k=[0:7]
stem(k,abs(f1)); axis([0 8 0 20]) xlabel('k')
ylabel('|x3(k)|')
title('x3(n)的8点fft')
分析：对8 点和16 点FFT 变换分析如下：离散傅里叶变换的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、7π/8 的傅里叶变换，这在上面两图可以明显看出。所以，离散傅里叶变换实际上是对该序列在频域范围内以2π/N 的间隔进行抽样。
4．()
4
x n=cos(pi/4*n)
原程序：
n=[0:7];
x=cos(0.25*pi*n) f1=fft(x,8) subplot(2,2,1) stem(n,x);
axis([0 8 -4 4]) xlabel('n')
ylabel('x4(n)') title('x4的波形') n=[0:15]
x=cos(0.25*pi*n) f2=fft(x,16) subplot(2,2,2) stem(n,x);
axis([0 16 -4 4]) xlabel('n') ylabel('x4(n)')
title('x4的波形') subplot(2,2,4)
k=[0:15]
stem(k,abs(f2));
axis([0 16 0 20]) xlabel('k')
ylabel('|x4(k)|')
title('x4(n)的16点fft') subplot(2,2,3)
k=[0:7]
stem(k,abs(f1));
axis([0 8 0 20]) xlabel('k')
ylabel('|x4(k)|')
title('x4(n)的8点fft')
波形图：
分析：对这两组图的分析除了同x4(n)类似的内容外，还可以得出以下分析：如上图所示，N 取值越大，越接近余弦曲线，即越接近原输入连续信号，且时域周期频域离散。该信号的周期为8，对其取8 点进行离散傅里叶变换的结果实际上是对原输入的周期信号进行傅里叶级数变换所截取的主值序列，从上图可以看出周期函数的傅里叶变换在以2π/8 为间隔的变换值即K=0,1,2,3,4,5,6,7 的所对应的离散傅里叶变换值和其傅里叶级数是一样的。 5．()5x n =sin(pi/8*n) 原程序：n=[0:7];
x=sin((pi*n)/8) f1=fft(x,8) subplot(2,2,1) stem(n,x);
axis([0 8 -4 4]) xlabel('n')
ylabel('x5(n)') title('x5的波形') n=[0:15]
x=sin(0.125*pi*n) f2=fft(x,16) subplot(2,2,2) stem(n,x);
axis([0 16 -4 4]) xlabel('n') ylabel('x5(n)') title('x5的波形') subplot(2,2,4) k=[0:15]
stem(k,abs(f2)); axis([0 16 0 20]) xlabel('k')
ylabel('|x5(k)|')
title('x5(n)的16点fft') subplot(2,2,3) k=[0:7]
stem(k,abs(f1)); axis([0 8 0 20]) xlabel('k')
ylabel('|x5(k)|')
title('x5(n)的8点fft')
波形图：
分析：从以上两组图可以看出，若按X1(n)进行分析，则明显不对。现分析如下：原信号周期为16，所以当N=8 时，未能取完一个周期的值，N=16 则取完了一个周期的值，所以这是两个不同的序列，所以按照X1(n)的分析方式是不对的，因为本身它们的傅里叶变换就是不一样的。由于离散傅里叶变换是该序列周期延拓后所对应的傅里叶级数变换的主值序列，所以，当N=16 时，所得的DFT 值与X5(n)的傅里叶级数变换的主值序列是一致的，而N=8 时是X5(n)的部分序列的周期延拓后的傅里叶级数变换的主值序列，因此两者的值是不同的。
6．()
6
x t=cos(pi*8*t)+ sin(pi*16*t)+cos(20*pi*t)
原程序：
Ts=1/64;
n=0:15;
Xa=cos(8*n*Ts*pi)+cos(16*n*Ts*pi )+cos(20*n*Ts*pi);
f1=fft(Xa,16);
subplot(3,2,1);
stem(n,Xa);
axis([0 15 -2 3]);
xlabel('n');
ylabel('X6(n)');
title('X6(n) N=16');%显示x6（n）N=16
k=0:15
subplot(3,2,2);
stem(k,abs(f1));
axis([0 16 0 15]);
xlabel('k');
ylabel('|X6(k)|');
title('X6(n) N=16 的16点FFT');%显示X6(n)的16点FFT
n=0:31;
Xb=cos(8*n*Ts*pi)+cos(16*n*Ts*pi )+cos(20*n*Ts*pi);
f2=fft(Xb,32);
subplot(3,2,3);
stem(n,Xb);
axis([0 32 -2 3]);
xlabel('n');
ylabel('X6(n)');
title('X6(n) N=32');%显示x6（n）N=32
subplot(3,2,4);
stem(abs(f2));
axis([0 32 0 20]);
xlabel('k');
ylabel('|X6(k)|');
title('X6(n) N=32 的32点FFT');%显示X6(n)的32点FFT
n=0:63;
Xc=cos(8*n*Ts*pi)+cos(16*n*Ts*pi )+cos(20*n*Ts*pi); f3=fft(Xc,64); subplot(3,2,5); stem(n,Xc);
axis([0 64 -2 3]); xlabel('n');
ylabel('X6(n)'); title('X6(n) N=64');%显示x6（n ）N=64
subplot(3,2,6); stem(abs(f3)); axis([0 64 0 40]); xlabel('k');
ylabel('|X6(k)|');
title('X6(n) N=64 的64点FFT');%显示X6(n)的64点FFT
波形图：
分析:原连续信号的周期为0.5，当采样频率Fs=64Hz 时，所形成的序列周期为0.5*64=32。所以只有N>32,才能取完一个周期的序列。这一点，从上面三个图可以清晰看出。其中N=32 和N=16 的图形分析，可以参考x4(n)的分析。N=32 和 N=64 的图形分析，可以参考x5(n)的分析。 （2）令
45()()()
x n x n x n =+，用FFT 计算8点和16点离散傅立叶变换，
[]
()()x k DFT x n =，并根据DFT 的对称性，由()x k 求出
[]
44()()x k DFT x n =和
[]
55()()x k DFT x n =并与（1）中所得结果比较。提示（取N=16时，44()()
x n x N n =-，
55()()
x n x N n =--
实验结果： n=[0:7];
x=cos(0.25*pi*n)+sin(0.125*pi*n) f1=fft(x,8) subplot(2,2,1) stem(n,x);
axis([0 8 -4 4]) xlabel('n')
ylabel('x7(n)')
title('x7的波形') n=[0:15]
x=cos(0.25*pi*n)+sin(0.125*pi*n) f2=fft(x,16) subplot(2,2,2) stem(n,x);
axis([0 16 -4 4])
xlabel('n')
ylabel('x7(n)') title('x7的波形') subplot(2,2,4) k=[0:15]
stem(k,abs(f2)); axis([0 16 0 20]) xlabel('k')
ylabel('|x7(k)|') title('x7(n)的16点fft') subplot(2,2,3) k=[0:7]
stem(k,abs(f1)); axis([0 8 0 20]) xlabel('k')
ylabel('|x7(k)|')
title('x7(n)的8点fft')
波形图：
n=[0:15];
x=cos(0.25*pi*n)+sin(0.125*pi*n) f1=fft(x,16)
Re=(f1+conj(f1))/2 Im=(f1-conj(f1))/2 subplot(2,2,1) stem(n,abs(Re)); axis([0 16 0 20]) xlabel('k') 波形图：
ylabel('|Re(x7(k))|') title('恢复后的x4(k)') subplot(2,2,2) stem(abs(Im));
axis([0 16 0 20]) xlabel('k')
ylabel('|Im(x7(k))|') title('恢复后的x5(k)')
（3）令45()()()x n x n jx n =+，重复（2） n=[0:15];
x=cos(0.25*pi*n)+j*sin(0.125*pi*n)
f1=fft(x,16) subplot(2,2,2) stem(n,abs(f1)); axis([0 16 0 10]) xlabel('k')
ylabel('|x8(k)|') title('x8(n)的16点fft') subplot(2,2,1) k=0:7
f2=fft(x,8)
stem(k,abs(f2)); axis([0 8 0 10]) xlabel('k')
ylabel('|x8(k)|')
title('x8(n)的8点fft')
波形图：
n=[0:15];
x=cos(0.25*pi*n)+j*sin(0.125*pi*n)
f=fft(x,16)
k(1)=conj(f(1)) m=2:16
k(m)=conj(f(16-m+2)) fe=(f+k)/2 fo=(f-k)/2
xr=ifft(fe,16) xo=ifft(fo,16)/j subplot(2,2,1) stem(n,xr);
axis([0 16 -1 1]) xlabel('n')
ylabel('|x4(n)|') title('x4(n)的波形') subplot(2,2,2) stem(n,abs(fe)); axis([0 16 0 10]) xlabel('k')
ylabel('|x8e(k)|')
title('x4(n)的16点fft') subplot(2,2,3) stem(n,xo);
axis([0 16 -1 1]) xlabel('n')
ylabel('|x5(n)|') title('x5(n)的波形') subplot(2,2,4) stem(n,abs(fo)); axis([0 16 0 10]) xlabel('k')
ylabel('|x8o(k)|')
title('x5(n)的16点fft')
波形图：
4. 思考题
(1) 在N=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?
答：N=8是x2(n)和x3(n)的幅频特性相同；N=16是x2(n)和x3(n)的幅频特性不相同。N=8时造成了栅栏效应漏掉大的频谱分量。
(2) 如果周期信号的周期预先不知道，如何用FFT进行谱分析?
答：如果周期预先不知道，可先截取M点惊醒FFT，再将截取长度扩大一倍，比较两次变换，如果而至的主谱差别满足分析误差要求，则用其近似表示周期信号的频谱，否则继续将截取长度加倍，直至前后两次分析所得主谱频率差别满足误差要求。