哈工大数字信号处理实验报告

H a r b i n I n s t i t u t e o f T e c h n o l o g y

实验报告

课程名称：数字信号处理

实验题目：用FFT作谱分析、用窗函数法设计数字滤波器

院系：电子与信息工程学院

班级：

哈尔滨工业大学

一、实验目的

(1) 进一步加深DFT算法原理和基本性质的理解。

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

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

二、实验步骤

(1) 复习DFT的定义、性质和用DFT作谱分析的有关内容。

(2) 复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图。

(3) 编制信号产生程序，产生以下典型信号供谱分析用：

1423()()

1,03()847

403()347

0x n R n n n x n n n n n x n n n =?+≤≤?

=-≤≤???-≤≤??

=-≤≤???

456()cos 4

()sin

8

()cos8cos16cos 20x n n x n n

x t t t t π

π

πππ=≤≤=≤≤=++，0n 19，0n 19

(4) 按实验内容要求，上机实验，并写出实验报告。 三、实验原理

1、FFT 产生背景。

离散傅立叶变换从理论上解决了傅立叶变换应用于实际的可能性，但若直接按DFT 公式计算，运算量太大（与N 2成比例）。快速傅立叶变换（FFT ）是离散傅立叶变换的快速算法，它大大减少了离散傅立叶变换的运算次量，一般可缩短一、二故数量级，且N 越大改善越明显，从而使DFT 的运算在实际工作中才真正得到广泛应用。

2、FFT 原理。

减少离散傅立叶变换运算次数的方法基于W nk N

的对称性和周期性。

（1）对称性

W W kn

N

kn N -==*n)-k (N N

)(W

（2）周期性

W W W nk

N

k

N n N

kn

N N -+===)()(mod k n N

W

）（

由此可得

W W

W W W k

N

N k N

N N nk

N

k

n N N -=-===+--)

2/(2/)(k)-n(N N

1W

这样可以将DFT 运算中有些项进行合并，利用W nk

N 的对称性和周期性，将长序列的DFT

分解为短序列的DFT 。

3、FFT 运算量。

在采用按时间抽取法的情况下，当N=2M 时共有M 级蝶形，每级都有N/2个蝶形运算。计算每个蝶形需要一次复数乘法和两次复数加法，因而每级运算需N/2次复数乘法和N 次复数加法，这样M 级运算总共需要：

复数乘法次数：

m F =N/2.M=N/2.log 2N 复数加法次数：

a F =NM=Nlog 2N 在一般情况下，复数乘法所需时间比复数加法多，因此以复数乘法为例将DFT 运算量与FFT 运算量进行对比。直接进行DFT 运算复数乘法次数为N 2次，利用FFT 运算复数乘法为N/2.log 2N 次。

计算量之比为：

N N N log log N 222

22N FFT DFT ==

复数乘法次数复数乘法次数

四、实验结果、源程序及误差分析

1、

实验程序：

N1=8;

N2=16;

x1=[1,1,1,1];

y11=fft(x1,8);

y12=fft(x1,16)

subplot(2,2,1);

stem(0:3,x1);

title('函数x1时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y11));

title('N=8，x1的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y12));

title('N=16，x1的DFT图像');box off

实验误差分析：理论的DFT函数应为抽样函数，由实验结果可以看出，随着FFT点数的

增加，FFT图像包络函数越趋近于抽样函数。

2、

实验结果：

实验程序：

N1=8;

N2=16;

xa=1:1:4;

xb=4:-1:1;

x2=[xa,xb];

y21=fft(x2,8);

y22=fft(x2,16);

subplot(2,2,1);

stem(0:7,x2);

title('函数x2时域图像');

subplot(2,2,2);box off

stem(0:(N1-1),abs(y21));

title('N=8，x2的DFT图像');box off

subplot(2,2,4);

stem(0:(N2-1),abs(y22));

title('N=16，x2的DFT图像');box off

实验误差分析：理论的DFT函数应为形如抽样函数的平方的函数，由实验结果可以看出，随着FFT点数的增加，FFT图像包络函数接近抽样函数的平方。

3、

实验结果：

实验程序：N1=8;

N2=16;

xa=4:-1:1;

xb=1:1:4;

x3=[xa,xb];

y31=fft(x3,8);

y32=fft(x3,16);

subplot(2,2,1);

stem(0:7,x3);

title('函数x3时域图像');

subplot(2,2,2);box off

stem(0:(N1-1),abs(y31));

title('N=8时，x3的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y32));

title('N=16时，x3的DFT图像');box off

实验误差分析：随着N的增加，曲线接近理论值，具体的分析和2中的信号相同

4、

实验结果：

实验程序：

N1=8;

N2=16;

n=0:1:19

x4=cos(0.25*pi*n);

n=0:1:7;

x411=cos(0.25*pi*n); n=8:1:15;

x412=cos(0.25*pi*n); n=16:1:19;

x413=cos(0.25*pi*n); y411=fft(x411,8);

y412=fft(x412,8);

y413=fft(x413,8);

y41=y411+y412+y413; n=0:1:15;

x421=cos(0.25*pi*n);

n=16:1:19;

x422=cos(0.25*pi*n);

y421=fft(x421,16);

y422=fft(x422,16);

y42=y421+y422;

subplot(2,2,1);

stem(0:19,x4);

title('函数x4时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y41));

title('N=8,x4的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y42));

title('N=16,x4的DFT图像');box off

实验误差分析：理论的DFT函数为冲激函数，从实验结果来看，有两个明显的峰值可以看到，与理论值接近。

5、

实验结果：

实验程序：

n=0:1:19;

N1=8;

N2=16;

x5=sin(0.125*pi*n); n=0:1:7;

x511=sin(0.125*pi*n); n=8:1:15;

x512=sin(0.125*pi*n); n=16:1:19;

x513=sin(0.125*pi*n); y511=fft(x511,8);

y512=fft(x512,8);

y513=fft(x513,8);

y51=y511+y512+y513; n=0:1:15;

x521=sin(0.125*pi*n);

n=16:1:19;

x522=sin(0.125*pi*n);

y521=fft(x521,16);

y522=fft(x522,16);

y52=y521+y522;

subplot(2,2,1);

stem(0:19,x5);

title('函数x5时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y51));

title('N=8,x5的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y52));

title('N=16,x5的DFT图像');box off

实验误差分析:理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真，当N=8时，满足抽样定理，因此可得到与理论曲线一致的图形。

6、x4(n)+x5(n)

实验结果：

实验程序：

n=0:1:19;

N1=8;

N2=16;

x4=cos(0.25*pi*n);

x5=sin(0.125*pi*n);

x45=x4+x5;

n=0:1:7;

x4511=cos(0.25*pi*n)+sin(0.125*pi*n); n=8:1:15;

x4512=cos(0.25*pi*n)+sin(0.125*pi*n); n=16:1:19;

x4513=cos(0.25*pi*n)+sin(0.125*pi*n); y4511=fft(x4511,8);

y4512=fft(x4512,8);

y4513=fft(x4513,8);

y451=y4511+y4512+y4513;

n=0:1:15;

x4521=cos(0.25*pi*n)+sin(0.125*pi*n);

n=16:1:19;

x4522=cos(0.25*pi*n)+sin(0.125*pi*n);

y4521=fft(x4521,16);

y4522=fft(x4522,16);

y452=y4521+y4522;

subplot(2,2,1);

stem(0:19,x45);

title('函数x45时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y451));

title('N=8,x45的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y452));

title('N=16,x45的DFT图像');box off

实验误差分析：这个函数的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真。

7、x4(n)+j*x5(n)

实验结果：

实验程序：

n=0:1:19;

N1=8;

N2=16;

x4=cos(0.25*pi*n);

x5=sin(0.125*pi*n);

x45=x4+x5;

n=0:1:7;

x4511=cos(0.25*pi*n)+j*sin(0.125*pi*n); n=8:1:15;

x4512=cos(0.25*pi*n)+j*sin(0.125*pi*n); n=16:1:19;

x4513=cos(0.25*pi*n)+j*sin(0.125*pi*n); y4511=fft(x4511,8);

y4512=fft(x4512,8);

y4513=fft(x4513,8);

y451=y4511+y4512+y4513;

n=0:1:15;

x4521=cos(0.25*pi*n)+j*sin(0.125*pi*n);

n=16:1:19;

x4522=cos(0.25*pi*n)+j*sin(0.125*pi*n);

y4521=fft(x4521,16);

y4522=fft(x4522,16);

y452=y4521+y4522;

subplot(2,2,1);

stem(0:19,x45);

title('函数x45时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y451));

title('N=8,x45的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y452));

title('N=16,x45的DFT图像');box off

实验误差分析：理论FFT抽样频谱为两个单位冲激函数，因此只要满足抽样定理，FFT 都可以与理论图形一样的曲线。当N=8时，不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。

8、

实验结果：

实验程序：n=0:1:69;

N1=16;

N2=32;

N3=64;

fs=64;

x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=0:1:15;

x611=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=16:1:31;

x612=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=32:1:47;

x613=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=48:1:63;

x614=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69;

x615=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y611=fft(x611,N1);

y612=fft(x612,N1);

y613=fft(x613,N1);

y614=fft(x614,N1);

y615=fft(x615,N1);

y61=y611+y612+y613+y614+y615;

n=0:1:31;

x621=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=32:1:63;

x622=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69;

x623=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y621=fft(x621,N2);

y622=fft(x622,N2);

y623=fft(x623,N2);

y62=y621+y622;

n=0:1:63;

x631=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69;

x632=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y631=fft(x631,N3);

y632=fft(x632,N3);

y63=y631+y632;

subplot(2,2,1);

stem(0:69,x6);

title('函数x6时域图像');box off

subplot(2,2,2);

stem(0:(N1-1),abs(y61));

title('N=16,x6的DFT图像');box off

subplot(2,2,3);

stem(0:(N2-1),abs(y62));

title('N=32,x6的DFT图像');box off

subplot(2,2,4);

stem(0:(N3-1),abs(y63));

title('N=64,x6的DFT图像');box off

实验误差分析：理论FFT抽样频谱为两个单位冲激函数，因此只要满足抽样定理，FFT 都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真，当N=32,64时，满足抽样定理，因此可得到与理论曲线一致的图形。

五、思考题

(1) 在N=8时，x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?

答：从实验结果看出，N=8时，x2(n)和x3(n)的幅频特性相同。而N=16时，两者不同。原因是N=8时，x3(n)只是以长度8为周期，将其延拓称周期序列然后加以移位，最后截取主值区间的序列值。因此x3(n)是x2(n)进行循环位移后的结果。移位后只影响DFT相频特性而不影响幅频特性。而当N=16时，两个序列不满足循环位移的性质，因此幅频特性发生变化。

(2) 如果周期信号的周期预先不知道，如何用FFT进行谱分析?

答：若该信号是模拟信号，用FFT进行频谱分析时，首先对信号进行采样，使之变成离散信号，然后用FFT对离散信号进行频谱分析。由采样定理，采样频率fs应当大于信号最高频率的2倍。如果信号的周期预先不知道，可以选择一个较大的采样频率，用试验法逐步增大采样频率，观察频谱特性。

一、实验目的

(1)熟悉矩形窗、汉宁窗、海明窗和布莱克曼窗。

(2) 掌握用上述窗函数法设计FIR 数字滤波器的原理和方法。 (3) 熟悉线性相位FIR 数字滤波器特性。 (4) 了解各种窗函数对滤波特性的影响。 二、实验原理与方法

如果所希望的滤波器的理想频率响应函数为Hd(e j ω)，则其对应的单位脉冲响应为

1

()()2j j n d d h n H e e d π

ωωπ

ω

π

-

=

?

用窗函数w(n)将hd(n)截断， 并进行加权处理， 得到:

()()()

d h n h n n ω=

h(n)作为实际设计的FIR 数字滤波器的单位脉冲响应序列， 其频率响应函数H(ej ω)

为

1

()()N j j n

n H e h n e ωω--==∑

如果要求线性相位特性， 则h(n)还必须满足：

()(1)h n h N n =±--

根据上式中的正、 负号和长度N 的奇偶性又将线性相位FIR 滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性， 可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。 三、实验内容与步骤

(1)复习用窗函数法设计FIR 数字滤波器一节内容，阅读本实验原理，掌握设计步骤。 (2)编写程序，其中幅度特性要求用dB 表示。 备注：不许用 freqz() 这个函数

设