数字信号处理实验实验五汇总

数字信号处理实验报告
实验名称：应用FFT实现信号频谱分析
学生姓名：z
学生学号：
学生班级：
上课时间：周二上午
指导老师：
一、 实验目的
（1） 能够熟练掌握快速离散傅里叶变换的原理及应用
FFT 进行频谱分析的基本方法。

（2） 对离散傅里叶变换的主要性质及FFT 在数字信号处
理中的重要作用有进一步的了解。

二、 实验原理
1、离散傅里叶变换（DFT ）及其主要性质
DFT 表示离散信号的离散频谱，DFT 的主要性质中有奇偶对称性、虚实特性等。

通过实验可以加深理解。

实序列的DFT 具有偶对称的实部和奇对称的虚部，这可以证明如下： 由定义，可得
X(k)=∑-=1
)(N n kn N W n x
=)2sin()()2cos()(1
1
0kn N n x j kn N n x N n N n ∑∑-=-=-ππ
X(N-k)=∑-=-1
0)()(N n n k N N
W
n x =∑-=-1
)(N n kn n
Nn
W
W
n x =∑-=-1
)(N n kn N W n x
=)2sin()()2cos()(1
1
0kn N n x j kn N n x N n N n ∑∑-=-=+ππ
所以，X(k)=X *(N-k)
实序列DFT 的这个特性，在本实验中可以通过实指数序列及三角序列看出来。

对于单一频率的三角序列来说，它的DFT 谱线也是单一的，这个物理意义可以从实验中得到验证，在理论上可以推导如下： 设
)()2s i n ()(n R n N
n x N π
= 其DFT 为 X(k)=∑-=-1
02)(N n kn N
e
n x π
=kn N
j N n e n N π
π21
)2sin(--=∑=kn N j n N j N n N j e e e j π
π
π
22102)(21---=-∑
=)(21)
1(210
)1(2+--=---∑k n N j N n k n N j e e j π
π
从而，X(0)=0)(21210
2=---=∑n N j N n N
j e e j π
π
X(1)=2
2)1(21104N
j j N e j N n n N j -==-∑-=-π
X(2)=0 …… X(N-2)=0
X(N-1)=2
2)(21210)2(2N
j j N e e j n j N n n N N j =-=-∑-=--ππ
以上这串式中X(0)反映了x(n)的直流分量，X(1)是x(n)的一次谐波，又根据虚实特性X *(N-1)=X(1),而其他分量均为零。

当周期减少时显然)3
2sin(
n N π
R N (n)的谱只应该在k=3及k=N-3才有
分量，实验者可以通过和上述相同的步骤加以理论证明。

由于)(2cos
n R N N ∙π与)(2sin n R n N
N ∙π
相位差2π，所以它的DFT 只包括
实部而没有虚部，以上这些性质可在本实验中得到验证。

2、利用DFT 对信号进行频谱分析
DFT 的重要应用之一是对时域连续信号的频谱进行分析，称为傅里叶分析，时域连续信号离散傅里叶分析的基本步骤如图5-1所示。

LPF t s c →)()()()
()()(k V DFT D A j H n v n x t x a c
→−−→−⊗−−→−−−
→−Ω
图5-1 时域连续信号离散傅里叶分析的处理步骤
其中消混叠低通滤波器LPF （预滤波器）的引入，是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。

实际工作中，时域离散信号x(n)的时宽是很长的甚至是无限长的（例如语音或音乐信号）。

由于DFT 需要，必须把x(n)限制在一定时间间隔之内，即进行数据截断。

数据的截断相当于加窗处理。

因此，在计算机x(n)de DFT 之前，用一个时域有限的窗函数ω(n)加到x(n)上是非常必要的。

X c (t)通过D A 变换器转换成取样序列x(n)—时域离散序列。

其频谱用)(ωj e X 表示，它的频率ω的周期函数，即
)2(1)(T
m
j T j X T e X m c j πωω
+=∑∞-∞=
其中X c (j Ω)或)(T
j X c ω
为x c (t)的频谱。

在实际应用中，消混叠低通滤波器的阻带不可能是无限衰减的，
↑)
(n w
故由X c (j Ω)周期延拓得到的)(ωj e X 有非零重叠，即出现混叠现象。

由于进行DFT 的需要，必须对序列x(n)进行加窗处理，即v(n)=x(n)·ω(n)。

加窗对频域的影响，用周期卷积表示。

θπ
θωπ
πθ
ω
d e W e
X e V j j j )()(21
)()(--
⎰=
最后是进行DFT 运算，加窗后的DFT 为 10,)()(1
02-≤≤=∑-=-N k e
n v k V N n nk N
j π
其中假设窗函数长L 小于或等于DFT 长度N 。

有限长序列v(n)=x(n)·ω(n)的DFT 相当于v(n)傅里叶变换的等间隔取样。

k N
j e V k V πωω2|
)()(=
=
V(k)便是s c (t)的离散频率函数。

因为DFT 频率间隔为
N
π
2，且模拟频率Ω和数字频率ω间的关系ω=ΩT ，所以离散频率点对应的模拟频率为 NT
k
π2=
Ω 显然频率分辨率Δf 为 Δf =
NT
1 利用DFT 计算频谱，只会给频谱k N k πω2=
或k NT
k π2=Ω的频率分量，即频率的取样值，而不可能得到连续的频谱函数。

如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。

为了保持原来频谱形状不变的情况下，使谱线加密，即使频域取样点数增加，从而使原来看不到的频谱分量变得可以看到，可以通过
在信号数据的末端补加一些零值点，使DFT 计算周期内点数增加，但又不改变原有的记录数据的方法实现。

3、快速离散傅里叶变换（FFT ）
快速离散傅里叶变换是计算离散傅里叶变换的一种快速算法，为了提高运算速度，FFT 将DFT 的计算逐次分解成较小点数的DFT 。

按时间抽取FFT 算法把输入序列x(n)按其n 值为偶数或是奇数分解成越来越短的序列。

按频域抽取FFT 算法是把输出序列X(k)按其k 值是偶数或是奇数来分解成越来越短的序列。

三、实验内容
（1）实验前学生应认真学习《数字信号处理》中有关章节的内
容，···········掌握DFT 的基本理论和应用FFT 计算信号频谱的原理与方法；
（2）编写一个调用FFT 函数的通用程序，可计算下列三种序列的离
散频谱；
指数序列：v 1(n)=(0.9)n u(n); 周期为N 的余弦序列：v 2(n)=cos(N
π
2n)，且0≤ n ≤ N-1； 符合函数序列：v 3(n)=0.9sin(
N π2n)+0.6sin(3
2N πn) （3）计算实指数序列v 1(n)的N 点离散频谱V 1（k ），记录N 为不同 的2的幂次方时的V 1（k ）值，并与理论V 1(e j ωk )进行比较。

（4）计算周期为N 的余弦序列v 2(n)的N 点FFT 、2N 点FFT 及（N+2）点FFT ，记录结果并分析说明。

（5）已知信号x(t)=0.15sin(2πf1t)+sin(2πf2t)-0.1sin(2πf3t)，其中f1=1Hz，f2=2Hz，f3=3Hz，取样频率为32Hz。

编程实现：
① 32点FFT，画出其幅度谱。

② 64点FFT，画出其幅度谱，比较两者间的差异，思考实际频率与
离散频谱图中横坐标k的对应关系。

1,书中例5.1 MATLAB程序如下：
clear all
N=100;
n=0:N-1;
xn=2*sin(0.48*pi*n)+cos(0.52*pi*n);
XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('x(n)');
title('x(n) N=100');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|X(n)|');
title('X(k) N=100');
执行结果如下图形所示：
2，实验内容第（2）小题
ⅰ.指数序列：v1(n)=(0.9)n u(n) MATLAB程序如下：
clear all
N=100;
n=0:N-1;
xn=0.9.^n;
XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v1(n)');
title('v1(n) N=100');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V1(n)|');
title('V1(k) N=100');
执行结果如下图形所示：
第（2）小题
2n)，且0≤ n ≤ N-1ⅱ. 周期为N的余弦序列：v2(n)=cos(
N
MATLAB程序如下：
clear all
N=100;
n=0:N-1;
xn=cos(2*pi/N*n);
XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v2(n)');
title('v2(n) N=100');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V2(n)|');
title('V2(k) N=100');
执行结果如下图形所示：
ⅲ. 符合函数序列：v 3(n)=0.9sin(
N π2n)+0.6sin(3
2N πn) MATLAB 程序如下：
clear all
N=100; n=0:N-1;
xn=0.9*sin(2*pi/N*n)+0.6*sin(2*pi/N*3*n); XK=fft(xn,N); magXK=abs(XK); phaXK=angle(XK); subplot(1,2,1); plot(n,xn);
xlabel('n');ylabel('v3(n)'); title('v3(n) N=100'); subplot(1,2,2);
k=0:length(magXK)-1; stem(k,magXK,'.');
xlabel('k');ylabel('|V3(n)|'); title('V3k) N=100');
第（3）小题N=128
MATLAB程序如下：
clear all
N=128;
n=0:N-1;
xn=0.9.^n;
XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v1(n)');
title('v1(n) N=128');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V1(n)|'); title('V1(k) N=128');
N=64 MATLAB程序如下：
clear all
N=64;
n=0:N-1;
xn=0.9.^n;
XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v1(n)');
title('v1(n) N=64');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V1(n)|');
title('V1(k) N=64');
执行结果如下图形所示：
N=32 MATLAB程序如下：
clear all
N=32;
n=0:N-1;
xn=0.9.^n; t(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v1(n)');
title('v1(n) N=32');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V1(n)|');
title('V1(k) N=32');
执行结果如下图形所示：
第（4）小题ⅰ. N点FFT MATLAB程序如下：
clear all
N=100;
M=N;
n=0:N-1;
xn=cos(2*pi/N*n);
XK=fft(xn,M);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v2(n)');
title('v2(n) M=N');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V2(n)|');
title('V2(k) M=N');
执行结果如下图形所示：
ⅱ.2N点FFT MATLAB程序如下：
clear all
N=100;
M=2*N;
n=0:N-1;
xn=cos(2*pi/N*n);
XK=fft(xn,M);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v2(n)');
title('v2(n) M=2*N');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V2(n)|');
title('V2(k) M=2*N');
执行结果如下图形所示：
ⅲ.N+2点FFT MATLAB程序如下：
clear all
N=100;
M=N+2;
n=0:N-1;
xn=cos(2*pi/N*n);
XK=fft(xn,M);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v2(n)');
title('v2(n) M=N+2');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V2(n)|');
title('V2(k) M=N+2');
执行结果如下图形所示：
第（5）小题ⅰ.32点FFT
MATLAB程序如下：
clear all
N=32;
n=0:N-1;
xn=0.15*sin(2*pi/32*n)+sin(2*pi*2/32*n)-0.1*sin(2*pi*3/32*n); XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v3(n)');
title('v3(n) N=32');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V3(n)|');
title('V3k) N=32');
执行结果如下图形所示：
ⅱ.64点FFT MATLAB程序如下：
clear all
N=64;
n=0:N-1;
xn=0.15*sin(2*pi/32*n)+sin(2*pi*2/32*n)-0.1*sin(2*pi*3/32*n); XK=fft(xn,N);
magXK=abs(XK);
phaXK=angle(XK);
subplot(1,2,1);
plot(n,xn);
xlabel('n');ylabel('v3(n)');
title('v3(n) N=32');
subplot(1,2,2);
k=0:length(magXK)-1;
stem(k,magXK,'.');
xlabel('k');ylabel('|V3(n)|');
title('V3k) N=32');
执行结果如下图形所示：
四、实验总结
通过本次实验，我能够熟练掌握快速离散傅里叶变换的原理，以及理解了频谱分析的意义，学会了使用求序列的DFT函数fft和求IDFT函数ifft。

通过对书中例子的分析，掌握每个Matlab语句的意义和作用，实验内容中的题目都是与之相关的，只要稍作修改即可。