数字信号实验5报告

广州大学学生实验报告

开课学院及实验室：电子楼317EDA

年月日 学院 机械与电气

工程学院

年级、

专业、班

姓名

学号

实验课

程名称

数字信号处理实验 成绩 实验项

目名称

用DFT （FFT ）对信号进行频谱分析

指导教师

张承云

一、 实验目的

学习DFT 的基本性质及对时域离散信号进行频谱分析的方法，进一步加深对频域概念和数字频率的理解，掌握MATLAB 函数中FFT 函数的应用。

二、 实验原理

离散傅里叶变换（DFT ）对有限长时域离散信号的频谱进行等间隔采样，频域函数被离 散化了，便于信号的计算机处理。

设x (n )是一个长度为M 的有限长序列，x (n )的N 点傅立叶变换：

其中N

j

N e

W π

2-=，它的反变换定义为：

∑-=-=

1

)(1

)(N k nk N

W

k X N

n x

令

k N

W

z -=，则有：

∑-=-==1

0)()

(N n nk

N k N

W z W n x z X

可以得到，

k N

W z z X k X -==)

()(，

k

N

W

z -=是Z 平面单位圆上幅角为

k

N πω2=

的点，就是将单位圆进

行N 等分以后第K 个点。所以，X(K)是Z 变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足Nyquist 定理时，就不会发生频谱混叠。DFT 是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。

如果用FFT 对模拟信号进行谱分析，首先要把模拟信号转换成数字信号，转换时要求知道模拟信号的最高截至频率，以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的3~4倍。另外要选择对模拟信号的观测时间，如果采样频率和观测时间确定，则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关，最小的观测时间和分辨率成倒数关系。最小的采样点数用教材相关公式确定。要求选择的采样点数和观测时间大于它的最小值。

如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的周期，要尽量选择观测时间长一些，以减少阶段效应的影响。

用FFT 对模拟信号作谱分析是一种近似的谱分析，首先一般模拟信号（除周期信号以外）的频谱是连续谱，而用FFT 作谱分析得到的是数字谱，因此应该取FFT 的点数多一些，用它的包络作为模拟信号的近似谱。另外，如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行滤波，或者尽量采样频率取高一些。

一般频率混叠发生在折叠频率附近，分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的，分析时要截断，截断的长度和分辨率有关，但也要尽量取长一些，取得太短会截断引起的误差会很大。举一个极端的例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不是周期的倍数，这种截断效应也会小一些。

如同理论课教材所讨论的，在运用DFT 进行频谱分析的时候可能有三种误差，即： （1）混叠现象

当采样率不满足Nyquist 定理，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实的反映原信号的频谱。所以，在利用DFT 分析连续信号频谱的时候，必须注意这一问题。避免混叠现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。 （2）泄漏现象

实际中的信号序列往往很长，甚至是无限长。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT 来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。值得一提的是，泄漏是不能和混叠完全分离开的，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。 （3）栅栏效应

因为DFT 是对单位圆上Z 变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度看， 用DFT 来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真是的频谱。这样的话就会有一些频谱的峰点或谷点被“栅栏”挡住，不能被我们观察到。减小栅栏效应的一个方法是在源序列的末端补一些零值，从而变动DFT 的点数。这种方法的实质是改变了真是频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些频谱的峰点或谷点显露出来。注意，这时候每根谱线所对应的频和原来的已经不相同了。

从上面的分析过程可以看出，DFT 可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。

DFT 运算量较大，快速离散傅里叶变换算法FFT 是解决方案。FFT 并不是DFT 不相同

的另一种变换，而是为了减少DFT 运算次数的一种快速算法。它是对变换式进行一次次的分解，使其成

为若干小点数DFT 的组合，从而减小运算量。常用的FFT 是以2为基数的，其长度为M

N 2=。它的运

算效率高，程序比较简单，使用也十分的方便。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT ，可以用末尾补零的方法，使其长度延长至2的整数次方。IFFT 一般也可以通过FFT 程序来完成。

三、 实验内容和步骤

(一) 离散信号

给定参考实验信号如下：

21j 0

()[()]()e 01N kn

N N n X k DFT x n x n k N π--===≤≤-∑

3()

x n ：用

14()()

x n R n =以8为周期进行周期性延拓形成的周期序列。

(1) 分别以变换区间N ＝8，16，32，对

14()()

x n R n =进行DFT(FFT)，画出相应的幅频特性曲线；

(2) 分别以变换区间N ＝4，8，16，对x 2(n )分别进行DFT(FFT)，画出相应的幅频特性曲线； (3) 对x 3(n )进行频谱分析，并选择变换区间，画出幅频特性曲线。 程序：

主函数（脚本文件）：

draw=@subfunction1;%产生函数句柄

n0=input('please input n0');%定义n0的值 n1=input('please input n1');%定义n1的值 n=n0:n1;%设置函数自变量的范围 x1=n>=0&n<4;%产生x1函数 x2=cos(pi*n/4);%产生x2函数 x3=0.*n;%产生容器数组

for i=0:2:n1%计算周期化次数 x3_mid=n>=0+i*4&n<=0+i*4+3; x3=x3+x3_mid;%正半周期化 end

for i=0:-2:n0%计算周期化次数

x3_mid=n>=-5+i*4-3&n<=-5+i*4; x3=x3+x3_mid;%负半周期化 end

%------------x1的三种幅频特性------------ draw(x1,8,n0,1,n);%x1的8点

title('x_{1}的8点幅频特性曲线');%设置标题 draw(x1,16,n0,2,n);%x1的16点

title('x_{1}的16点幅频特性曲线');%设置标题 draw(x1,32,n0,3,n);%x1的32点

title('x_{1}的32点幅频特性曲线');%设置标题 %------------x2的三种幅频特性------------ draw(x2,8,n0,4,n);%x2的8点

title('x_{2}的8点幅频特性曲线');%设置标题 draw(x2,16,n0,5,n);%x2的16点

title('x_{2}的16点幅频特性曲线');%设置标题 draw(x2,32,n0,6,n);%x2的32点

title('x_{2}的32点幅频特性曲线');%设置标题 %------------x3的三种幅频特性------------ draw(x3,8,n0,7,n);%x3的8点

title('x_{3}的8点幅频特性曲线');%设置标题 draw(x3,16,n0,8,n);%x3的16点

title('x_{3}的16点幅频特性曲线');%设置标题 draw(x3,32,n0,9,n);%x3的32点 子函数：

function subfunction1(x,n,n0,p,N)

%---------------计算不同fft 的子函数--------------- %x-----------------------需要进行fft 的序列 %n-----------------------fft 的点数

%n0----------------------x 序列对应的n 轴的起始位置 %p-----------------------子图选取位置

%N-----------------------fft 对应的横轴数组 Xfft=fft(x,n);%计算n 点的fft X_fft=N.*0;%产生容器数组

for i=0-n0+1:n+abs(n0) %产生与N 对应长度作为n 点FFT X_fft(i)=Xfft(i-abs(n0)); end

subplot(3,3,p);stem(2*N/n,abs(X_fft),'.');%画图 axis([0,2,0,inf]);%调整图形大小 xlabel('{\omega}/{\pi}');%设置x 轴 ylabel('|Xe^{j{\omega}}|');%确定Y 轴

（二）连续信号 实验信号：

1()()

x t R t τ=

选择 1.5ms τ=，式中()

R t τ的波形以及幅度特性如图7.1所示。 2()sin(2/8)

x t ft ππ=+

式中频率f 自己选择。

3()cos8cos16cos 20x t t t t

πππ=++

分别对三种模拟信号选择采样频率和采样点数。 对1()x t ()

R t τ=，选择采样频率

4s f kHz

=，8kHz ，16kHz ，采样点数用τ.

s

f 计算。

对

2()sin(2/8)

x t ft ππ=+，周期1/T f =，频率f 自己选择，采样频率4s f f

=，观测时间0.5p T T =，

T ，2T ，采样点数用p s T f 计算。

图5.1 R(t)的波形及其幅度特性 对

3()cos8cos16cos 20x t t t t

πππ=++，选择采用频率

64s f Hz

=，采样点数为16，32，64。

3. 分别对它们转换成序列，按顺序用

123(),(),()

x n x n x n 表示。

4. 分别对它们进行FFT 。如果采样点数不满足2的整数幂，可以通过序列尾部加0满足。

5. 计算幅度特性并进行打印。

以上各序列的DFT(FFT)均可调用MATLAB 信号处理工具箱函数fft 计算。

采样部分

主函数（脚本文件）：

ADtransfer=@subfunction2;%产生采样函数句柄

fftADtransfer=@subfunction3;%产生采样函数句柄 n0=input('please input n0');%定义n0的值 n1=input('please input n1');%定义n1的值 n=n0:n1;%确定离散采样范围 syms t;%定义时间变量

x1=heaviside(t)-heaviside(t-1.5*10^-3);%产生x1函数 x2=sin(2*pi*t+pi/8);%产生x2函数

x3=cos(8*pi*t+16*pi*t+20*pi*t);%产生x3函数

%---------------------x1的采样--------------------- figure(1);ADtransfer(x1,4000,t,n,1)%x1的4kHz 采样

figure(2);fftADtransfer(x1,4000,t,n,1,n0,n1)%x1的4kHz 采样幅频特性 figure(1);title('x_{1}的4kHz 采样');%确定标题 figure(1);ylabel('x_{1}(n)');%确定y 轴

figure(2);title('x_{1}的4kHz 采样幅频特性曲线');%确定标题 figure(2);ylabel('|X_1(e^{j{\omega}n})|');%确定y 轴 figure(1);ADtransfer(x1,8000,t,n,2)%x1的8kHz 采样

figure(2);fftADtransfer(x1,8000,t,n,2,n0,n1)%x1的8kHz 采样幅频特性 figure(1);title('x_{1}的8kHz 采样');%确定标题 figure(1);ylabel('x_{1}(n)');%确定y 轴

figure(2);title('x_{1}的8kHz 采样幅频特性曲线');%确定标题 figure(2);ylabel('|X_1(e^{j{\omega}n})|');%确定y 轴 figure(1);ADtransfer(x1,16000,t,n,3)%x1的16kHz 采样

figure(2);fftADtransfer(x1,16000,t,n,3,n0,n1)%x1的16kHz 采样幅频特性 figure(1);title('x_{1}的16kHz 采样');%确定标题 figure(1);ylabel('x_{1}(n)');%确定y 轴

figure(2);title('x_{1}的16kHz 采样幅频特性曲线');%确定标题 figure(2);ylabel('|X_1(e^{j{\omega}n})|');%确定y 轴 %---------------------x2的采样--------------------- figure(1);ADtransfer(x2,10,t,n,4)%x2的10Hz 采样

figure(2);fftADtransfer(x2,10,t,n,4,n0,n1)%x2的10Hz 采样幅频特性 figure(1);title('x_{2}的10Hz 采样');%确定标题 figure(1);ylabel('x_{2}(n)');%确定y 轴

figure(2);title('x_{2}的10Hz 采样幅频特性曲线');%确定标题 figure(2);ylabel('|X_2(e^{j{\omega}n})|');%确定y 轴 figure(1);ADtransfer(x2,20,t,n,5)%x2的20Hz 采样

figure(2);fftADtransfer(x2,20,t,n,5,n0,n1)%x2的20Hz 采样幅频特性 figure(1);title('x_{2}的20Hz 采样');%确定标题 figure(1);ylabel('x_{2}(n)');%确定y 轴

figure(2);title('x_{2}的20Hz 采样幅频特性曲线');%确定标题 figure(2);ylabel('|X_2(e^{j{\omega}n})|');%确定y 轴 figure(1);ADtransfer(x2,40,t,n,6)%x2的40Hz 采样

figure(2);fftADtransfer(x2,40,t,n,6,n0,n1)%x2的40Hz采样幅频特性figure(1);title('x_{2}的40Hz采样');%确定标题

figure(1);ylabel('x_{2}(n)');%确定y轴

figure(2);title('x_{2}的40Hz采样幅频特性曲线');%确定标题

figure(2);ylabel('|X_2(e^{j{\omega}n})|');%确定y轴

%---------------------x3的采样---------------------

figure(1);ADtransfer(x3,10,t,n,7)%x3的10Hz采样

figure(2);fftADtransfer(x3,10,t,n,7,n0,n1)%x3的10Hz采样幅频特性figure(1);title('x_{3}的10Hz采样');%确定标题

figure(1);ylabel('x_{3}(n)');%确定y轴

figure(2);title('x_{3}的10Hz 采样幅频特性曲线');%确定标题

figure(2);ylabel('|X_3(e^{j{\omega}n})|');%确定y轴

figure(1);ADtransfer(x3,20,t,n,8)%x3的20Hz采样

figure(2);fftADtransfer(x3,20,t,n,8,n0,n1)%x3的20Hz采样幅频特性figure(1);title('x_{3}的20Hz采样');%确定标题

figure(1);ylabel('x_{3}(n)');%确定y轴

figure(2);title('x_{3}的20Hz采样幅频特性曲线');%确定标题

figure(2);ylabel('|X_3(e^{j{\omega}n})|');%确定y轴

figure(1);ADtransfer(x3,40,t,n,9)%x3的40Hz采样

figure(2);fftADtransfer(x3,40,t,n,9,n0,n1)%x3的40Hz采样幅频特性figure(1);title('x_{3}的40Hz采样');%确定标题

figure(1);ylabel('x_{3}(n)');%确定y轴

figure(2);title('x_{3}的40Hz采样幅频特性曲线');%确定标题

figure(2);ylabel('|X_3(e^{j{\omega}n})|');%确定y轴

子函数1：

function subfunction2(x,fs,t,n,p)

%------------采样函数------------

%x-----------需要采样的信号

%fs----------采样频率

%t-----------原来模拟域连续自变量（符号变量）

%n-----------采样数字域离散自变量（数组）

%p-----------子图序号

digitalizex=subs(x,t,1/fs*n);%fsHz采样

subplot(3,3,p);stem(n,digitalizex,'.');%画图

axis tight;%控制显示范围

xlabel('n');%确定x轴

end

子函数2：

function subfunction3(x,fs,t,n,p,n0,n1)

%------------显示幅频特性函数------------

%x-----------需要采样的信号

%fs----------采样频率 %t-----------原来模拟域连续自变量（符号变量）

%n-----------采样数字域离散自变量（数组）

%p-----------子图序号

%n0----------n区间左值

%n1----------n区间右值

digitalizex=subs(x,t,1/fs*n);%fsHz采样

subplot(3,3,p);stem(2*n/(n1-n0+1),abs(fft(double(digitalizex),length(n))),'.');%画图axis tight;%控制显示范围

xlabel('\omega/\pi');%确定x轴

end