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

Harbin Institute of Technology

实验报告

课程名称： 数字信号处理 实验题目：用FFT作谱分析、用窗函数法设计数字滤波器 院 系： 电子与信息工程学院 班 级：

哈尔滨工业大学

一、实验目的

(1) 进一步加深DFT算法原理和基本性质的理解。 (2) 熟悉FFT算法原理和FFT子程序的应用。

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

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

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

0

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

x1(n)?R4(n)?n?1,0?n?3?x2(n)??8?n4?n?7??0?4?n0?n?3?x3(n)??n?34?n?7?0?

x4(n)?cosx5(n)?sin?4nn，0?n?19，0?n?19?8x6(t)?cos8?t?cos16?t?cos20?t

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

1、FFT产生背景。

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

2、FFT原理。

减少离散傅立叶变换运算次数的方法基于W（1）对称性

nkN的对称性和周期性。

W（2）周期性

k(N-n)N?(WN)*?WNknkn?kn

W由此可得

（kn）(modN)N?WN?WN(N?n)k(n?N)k?WN?nk

WWW

n(N-k)N?WN?WN?nkN/2NN??1??WNk(k?N/2)

nk这样可以将DFT运算中有些项进行合并，利用WN的对称性和周期性，将长序列的DFT

1

分解为短序列的DFT。

3、FFT运算量。

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

复数乘法次数：

mF=N/2.M=N/2.log2N 复数加法次数：

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

计算量之比为：

DFT复数乘法次数N?2N?NFFT复数乘法次数Nlog2Nlog22

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

1、

2

2

实验程序： 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点数的

3

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

2、

实验结果：

4

实验程序： 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、

实验结果：

5

实验程序： N1=8;

6

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、

实验结果：

7

实验程序： 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;

8

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、

实验结果：

9

实验程序： 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;

10

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) 实验结果：

11

实验程序： 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);

12

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) 实验结果：

13

实验程序： 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);

14

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、

实验结果：

15

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

16

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

17

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倍。如果信号的周期预先不知道，可以选择一个较大的采样频率，用试验法逐步增大采样频率，观察频谱特性。

18

一、实验目的

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

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

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

1hd(n)?2?

????Hd(ej?)ej?nd?

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

h(n)?hd(n)?(n)为

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

H(ej?)??h(n)e?j?nn?0

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

N?1h(n)??h(N?1?n)

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

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

设

19

H(k)?DFT[h(n)]H(k)?Hk(k)?jHI(k)2H(k)?HR(k)?HI2(k)

2?wk?kN 画图时，用20log|H(k)|打印幅度特性。第k点对应的频率为

为使曲线包络更接近H(ejw)的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33

时，可通过在h(n)末尾补零的方法，使长度变为64，再进行64点DFT，则可得到更精确的幅度衰减特性曲线。 四、上机实验内容

1、用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率Wc=

Π rad 。窗口N=15，433。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频曲线。观察3db带宽和20db带宽，总结窗口N对滤波特性的影响。

设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数

即

?j?a??e,???cj?Hd(e)???c??????0N?1a?21?1j?j?nhd(n)?Hd(e)e??2???2?sin?c(n?a)??(n?a)

20

????ce?j?ae?j?nd?c

五、实验结果 1、矩形窗，N=15

21

2、汉宁窗，N=15

22

23

3、海明窗，N=15

24

25

4、布莱克曼窗，N=15

26

27

5、矩形窗，N=33

28

29

elseif k==4

W=blackman(N); end

n=0:(N-1); wc=pi/4;

a=n-(N-1)/2+eps; hd=sin(wc*a)./(pi*a); h=hd.*(W)'; N2=1024;

H=fft(h,N2); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=unwrap(angle(H)); %得到相位 w=0:2*pi/N2:2*pi/N2*(N2-1); figure

subplot(2,2,1); n=0:(N-1); stem(n,h);

axis([0,N-1,-0.1,0.5]); xlabel('n'); ylabel('h(n)'); if N==15&&k==1

title('矩形窗N=15时的 h(n)'); elseif N==15&&k==2

title('汉宁窗N=15时的h(n)'); elseif N==15&&k==3

title('海明窗N=15时的h(n)'); elseif N==15&&k==4

title('布莱克曼窗N=15时的h(n)'); elseif N==33&&k==1

title('矩形窗N=33时的h(n)'); elseif N==33&&k==2

title('汉宁窗N=33时的h(n)'); elseif N==33&&k==3

title('海明窗N=33时的h(n)'); elseif N==33&&k==4

title('布莱克曼窗N=33时的h(n)'); end

subplot(2,2,2); plot(w,db);

axis([0,2,-100,5]); xlabel('w'); ylabel('db');

title('衰减特性（db）');

35

%grid on;

subplot(2,2,4); plot(w,abs(H)); axis([0,7,0,1.5]); title('幅频特性'); xlabel('w'); ylabel('幅值'); subplot(2,2,3); plot(w,pha); axis([0,4,-20,5]); title('相频特性'); xlabel('频率'); ylabel('相位'); %grid off;

N=input('请输入窗口长度N的值\\nN='); end

五、思考题

(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减， 如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。

① 通过阻带最小衰减确定使用哪一个窗函数

② 用阻带截止频率减去通带截止频率的绝对值还有对应窗函数的过渡带宽度确定阶数N

③ 构造希望逼近的频率响应函数 ④ 确定单位脉冲响应 ⑤ 加窗得到设计结果

(2) 如果要求用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd(n)。

① 由W1、W2求其理想线性相位带通滤波器频率响应 ； ② 用IDTFT求出hd(n)。

36

37

本文来源：https://www.bwwdw.com/article/edmx.html