数字信号处理实验报告汇总

目 录

实验1 离散时间信号的频域分析-----------------------2

实验2 FFT

实验3 IIR

实验4 FIR

算法与应用-------------------------------7 数字滤波器的设计------------------------12 数字滤波器的设计------------------------17 - 1 -

实验1 离散时间信号的频域分析

一．实验目的

信号的频域分析是信号处理中一种有效的工具。在离散信号的时域分析中，通常将信号表示成单位采样序列（n）的线性组合，而在频域中，将信号表示成复变量e

jn或 e

j2nN的线性组合。通过这样的表示，可以将时域的离散序

列映射到频域以便于进一步的处理。

在本实验中，将学习利用MATLAB计算离散时间信号的DTFT和DFT,并加深对其相互关系的理解。

二、实验原理

(1)DTFT和DFT的定义及其相互关系。序列x(n)DTFT定义为X(ejw)=

x(n)enjnX(ejw)它是关于自变量的复函数，且是以2为周期的连续函数。

可以表示为X(ejw)Xre(ejw)jXim(ejw),其中，Xre(ejw)和Xim(ejw)分别是

X(ejw)实部和虚部；|X(ejw)| 和还可以表示为 X(ejw)=|X(ejw)|ej(w)，其中，

(w)argX(ej)分别是X(ejw)的幅度函数和相位函数；它们都是的实函数，也是以2为周期的周期函数。

序列x(n)的N点DFT定义为X(k)X(ej2NkN1N0)x(n)ej2NknN1x(n)WNkn，n0X(k)是周期为N的序列。X(ej)与X(k)的关系：X(k)是对X(ej)）在一个周期

jw|X(k)X(e)|中的谱的等间隔N点采样，即

jX(e)可以 ，而kw2N通过对X(k)内插获得，即

- 2 -

X(e)jw1N1k0NX(k)sin(N2k)j(2k/N)2(N1)/2•e

N2ksin()2N(2)使用到的MATLAB命令有基于DTFT离散时间信号分析函数以及求解序列的DFT函数。

1）基于 DTFT离散时间信号分析函数有：freqz,real,imag,abs,angle。函数freqz可以用来计算一个以ej的有理分式形式给出的序列的DTFT值。freqz的形式多样，常见的有H=freqz(num,den,w)，其中num表示序列有理分式DTFT的分子多项式系数，den表示分母多项式系数（均按z的降幂排列），矢量w表示在0~2中给定的一系列频率点集合。freqz函数的其他形式参见帮助文件。在求出DTFT值后，可以使用函数real,imag,abs和angle分别求出并绘出起实

10.96ej0.9028ej2)部、虚部和相位谱。例如X(e)）=利用函数freqz计jj211.56e0.8109e)j算出H(ej)，然后利用函数abs和angle分别求出幅频特性与相位特性最后利用plot命令绘出曲线。

2）求解序列DFT的函数有：fft,ifft。函数fft（x）可以计算R点序列的R点DFT值；而fft（x，N）则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R三、实验设备

计算机、MATLAB软件

四、实验内容

（1)编程计算并画出下面DTFT的实部，虚部、幅度和相位谱。

0.13130.1553ejw0.1313ej2w0.0518ej3wX(e) jwj2wj3w11.2828e1.0388e0.3418ejw程序如下：num=[0.1313 -0.1553 0.1313 0.0518];

den=[1 1.2828 1.0388 0.3418];

- 3 -

w=0:0.001:2*pi; H=freqz(num,den,w); figure subplot(221) plot(w,real(H)); title('实部') grid on; subplot(222) plot(w,imag(H)) title('虚部') grid subplot(223) plot(w,abs(H)) title('幅度') grid subplot(224) plot(w,abs(H)) title('相位谱') grid on;

波形如下：

图1-1 DTFT的实部，虚部、幅度和相位谱

- 4 -

(2)计算32点序列x(n)=cos

13n,0n31的32点和64点DFT，分别绘出32幅度谱图形，并绘出该序列的DTFT图形。

程序如下：n1=0:15;

n2=0:31;

x1=cos((13*pi*n1)/32) x2=cos((13*pi*n2)/32) X1=fft(x1); X2=fft(x2); subplot(211) plot(n1,abs(X1)) title('32点DFT') grid on; subplot(212) plot(n2,abs(X2)) title('64点DFT') grid on;

波形如下：

图1-2 序列DFT幅度谱图形和DTFT图形

DTFT程序如下：

- 5 -

A=1; n=0:31;

x=cos((13*pi*n)/32); B=x;

w=0:0.01:2*pi*2; [H]=freqz(B,A,w); magH=abs(H); subplot(2,1,1); plot(w,magH);grid; ylable('Magnitude'); subplot(2,1,2); plot(w,phaH);grid; xlable('w'); ylable('phase'); grid on;

波形如下：

图1-3 DTFT图形

- 6 -

实验2 FFT算法与应用

一、实验目的

在理论学习的基础上，通过本次实验。加深对快速傅里叶变换的理解，熟悉FFT算法。熟悉应用FFT对典型信号进行频谱分析的方法。了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。

二．实验原理

在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。无限长序列往往也可以用有限长序列来逼近。对于有限长序列可以使用离散傅里叶变换（DFT），这一变换可以很好地反映序列的频域特性，并且容易利用快速算法在计算机的上的实现。当序列的长度是N时，定义离散傅里叶变换为

N1n02X（k）=DFT[x(n)]=x(n)eN1n0jNkn

WNkn =x(n)式中，WNej2N的反变换定义为

1N10 x(n)IDFT[X(k)]＝根据上式令ZWNk，则有 X(z)N1zWNkNX(k)WNkn kx(n)WNknn0DFT[x(n)]

所以，X（k）是在z变换在单位圆上的等间隔采样，或者说是序列傅里叶变换的等间隔采样。时域采样在满足乃奎斯特定理时，就是不会发生频谱混叠；同样，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列的混叠。

DFT是对序列傅里叶变换的等间距采样，因此可以利用于序列的频谱分析。在运用于DFT进行频谱分析的时候可能有三种误差，分析如下：

（１）混叠现象

序列的频谱是采样信号频谱的周期延拓，周期是

2，因此当采样速率不满T足乃奎斯特定理时，即采样频率小于两倍的信号频率时，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实地反映原信号的频谱。

（２）泄露现象

实际中的信号序列往往很长，甚至是无限长序列，为了方便，往往用截断的

- 7 -

序列来近似它们。这样可以使用较短的DFT来对信号进行频谱分析。这种截断等价于给原信号序列乘以一个矩形窗函数，而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。

（3）栅栏响应

因为DFT是对单位圆上z变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度来看，用DFT来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真实的频谱。减小栅栏效应的一个方法是在源程序的末端补一些零值，从而变动DFT的点数。

函数fft（x）可以计算R带你序列的R点的DFT值；而fft（x，N）则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R三． 实验设备

计算机，MATLAB软件。

四． 实验内容

编制信号产生子程序及本实验的频谱分析主程序。实验中需要用到的基本信号包括：

（1）三角波序列：

x(n)=n+13,0n3 13-n ,4n7 0,其它

程序如下：n2=0:7;

x2=[13,14,15,16,9,8,7,6]; subplot(221); stem(n2,x2); title('x2序列'); grid on; k1=0:7; y21=fft(x2,8); magy21=abs(y21); subplot(222); stem(k1,magy21); title('x2的8点FFT'); grid on; k2=0:15; y22=fft(x2,16); magy22=abs(y22); subplot(224);

- 8 -

stem(k2,magy22); title('x2的16点FFT'); grid on;

波形如下：

图2-1 三角波序列以及它的8点和16点FFT

（2）利用DFT计算下面两序列的线性卷积x(n)={3,-4,6,0,6,-4},h(n)={1,-4,

3,4}

程序如下：clear all; close all; clc; N=9;

g=[3 -4 6 0 6 -4]; h=[1 -4 3 4]; x=conv(g,h) gf=fft(g,N); hf=fft(h,N); xx=gf.*hf X=ifft(xx,N); grid on; 波形如下：

- 9 -

图2-2 用DFT计算的序列x(n)与h(n)线性卷积图

（3）已知某序列x（n）在单位圆上的N=64等分样点的z变换为X（Zk）=X（k）

=

1，k=0,1,2…,63。试分别用N点IFFT程序和N点FFTj2k/N10.13e程序计算x(n)IDFT[X(k)],绘出x(n)及|X(k)|。 程序如下：

k=0:63; N=64;

X=1./(1-0.13*exp(-j*2*pi*k/N)); x=ifft(X,64) stem(abs(x));

xlabel('n');ylabel('x(n)');title('25FF.单位圆上等分样点的Z变换图') grid on;

波形如下：

- 10 -

图2-3 单位圆等分样点Z变换图

- 11 -

实验3 IIR数字滤波器的设计

一． 实验目的

从理论上讲，任何的线性时不变（LTI）离散时间系统都可以看做一个数字滤波器，因此设计数字滤波器实际就是设计离散时间系统。数字滤波器包括IIR（无限冲击响应）和FIR（有限冲击响应）型，在设计时通常采用不同的方法。

本实验通过使用MATLAB函数对数字滤波器进行设计和实现，要求掌握IIR数字巴特沃斯滤波器,数字切比雪夫滤波器的设计原理，设计方法和设计步骤；能根据给定的滤波器指标进行滤波器设计；同时也加深学生对数字滤波器的常用指标和设计过程的理解。

二． 实验原理

在IIR滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器Ha(s),然后根据一定的准则将比模拟原型滤波器转换为数字滤波器Hd(z),即我们需要设计的数字滤波器。

IIR滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定主要是在模拟原型滤波器设计中进行的。

IIR数字滤波器的设计方法如下：

（1）冲击响应不变法。对满足设计要求的模拟原型滤波器Ha(s)进行部分分式展开为

Ha(s)sk1NAkks ;基于

Ah(n)=h(nT)，可得H(z)1eNdadk1kskTz1

(2)双线性变换法。对设计要求中给出的边界频率进行预畸处理，然后用得到的频率进行模拟滤波器设计，得到模拟原型滤波器Ha(s);用双线性变换法求出数字滤波器Hd(z)Ha(s)z1z1。

1z1一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。冲击响应不变法一个重要特点是频率坐标的变换是线性的，因此如果模拟滤波器的频率响应带限于折叠频率的话，则通过变换后滤波器的频率响应可不失真地反映原响应与频率的关系。

与冲激响应不变法比较，双线性变换的主要优点是靠频率的非线性关系得到s平面与z平面的单值一一对应关系，整个值对应于单位圆一周。所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。

MATLAB提供了一组标准的数字滤波器设计函数，大大简化了滤波器的设计过程。

- 12 -

（1） butter。功能：Butterworth模拟/数字滤波器设计如图；格式：[b,a]=butter(n,wn,‘ftype’,’s’),[b,a]=butter(n,wn,‘ftype’)。

（2） Cheby1、cheby2。功能：chebyshev、chebyshev型模拟/数字滤波器设计如图所示。格式：

[b,a]=cheby1(n,Rp,wn,ftype,),[b,a]=cheby2(n,Rs,wn,ftype)

IIR滤波器设计函数还包括：

buttord,chebwin,cheb1ord,cheb2ord,ellip,ellipord等。还可以用

[N,Wn]=buttord(Wp,Ws,Rp,Rs)的MATLAB命名估算一个Butterworth滤波器的阶数。

三．实验设备

计算机，MATLAB软件。 四．实验内容

利用MATLAB编程方法或利用MATLAB中的fdatool工具设计不同功能的IIR数字滤波器。

（1）基于chebyshev型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求参数为通带截止频率Wp=0.213；通带最大衰减Ap=1dB;阻带截止频率Ws=0.413；阻带最小衰减As=35dB。

程序如下：

clear; wp=0.213*pi; ws=0.413*pi; Ap=1; As=35; T=1;

Rip=10^(-Ap/20); Atn=10^(-As/20); OmegaP=wp*T; OmegaS=ws*T;

[n,Wn]=cheb1ord(OmegaP,OmegaS,1,35,'s'); [b,a]=cheby1(n,1,Wn,'low','s'); freqs(b,a);

[bz,az]=impinvar(b,a,T); [H,W]=freqz(bz,az,512,T);

- 13 -

ma=20*log10(abs(H)),pha=20*log10(unwrap(angle(H))),hi=impz(bz,az)

ni=step(bz,az);

subplot(2,2,1),plot(W,ma); title('幅频特性(dB)'); xlabel('w(/pi)'); ylabel('dB');

subplot(2,2,2),plot(W,pha); title('相频特性'); xlabel('w(/pi)'); ylabel('pha(/pi)'); subplot(2,2,3),plot(hi); title('单位冲击响应'); xlabel('n'); ylabel('h(n)');

subplot(2,2,4),plot(ni); title('单位阶跃响应'); xlabel('n'); ylabel('h(n)'); bz,az;

波形如下：

图3-1 基于chebyshevI模型的数字滤波器

- 14 -

;

(2)基于Butterworth型模拟滤波器原型使用双线性变换方法设计数字滤波器的，要求参数为通带截止频率Wp=0.24；通带最大衰减Ap=1.3dB;阻带截止频率Ws=0.13；阻带最小衰减As=40dB。

程序如下：

clear;

wp=0.4*pi; ws=0.13*pi; Ap=1.3; As=40; fs=1; ts=1/fs; wp2=2*fs*tan(wp/2*ts); ws2=2*fs*tan(ws/2*ts); [n,wn]=buttord(wp2,ws2,Ap,As,'s') [z,p,k]=buttap(n) [bap,aap]=zp2tf(z,p,k) [b,a]=lp2lp(bap,aap,wn) [bz,az]=bilinear(b,a,fs); [h,w]=freqz(bz,az); subplot(2,1,1); plot(w/pi,abs(h)); grid on xlabel('频率'); ylabel('幅度'); subplot(2,1,2);

plot(w/pi,20*log10(abs(h))); grid on xlabel('频率');

ylabel('幅度');

波形如下:

- 15 -

图3-2 基于Butterworth模型的数字滤波器

- 16 -

实验4 FIR数字滤波器的设计

一．实验目的

掌握用窗函数法设计FIR数字滤波器的原理与其设计步骤；熟悉线性相位数字滤波器的特性。学习编写数字滤波器的设计程序的方法，并能进行正确编程；根据给定的滤波器指标，给出设计步骤。

二．实验原理

如果系统的冲击响应hd(n)为已知，则系统的输入-输出关系为: y(n)x(n)*hd(n)

对于低通滤波器，只要设计出低通滤波器的冲击响应函数，就可以由式得到系统的输出了。假设所希望的数字滤波器的频率响应为Hd(ejw)，它是频域的周期函数，周期为2，那么它与Hd(ejw)相对应的傅里叶系数为:

1 hd(n)2Hd(ejw)ejnwdw

但是将hd（n）作为滤波器脉冲响应有两个问题：一是它是无限长的，与FIR数字滤脉冲响应有限长这一前提不一致；二是它是非因果的，

hd(n)0,n0。对此，要采取两项措施：一是将hd(n)截短；二是将其往右移。由此得到h2(n)的实际频域响应Hd(e)jwN10h(n)ejnw，与理想频域响应Hd(ejw)n2相近，但不完全一致。理论证明上述现象是对hd(n)进行简单截短处理的必然结果，一般称为吉布斯现象，为尽可能的减少吉布斯现象，应对hd(n)进行加窗截

WN(n)作为FIR滤波器的系数。其中N为窗口WN(n)的长度。取，即以h(n)hd(n)窗口函数的形状和窗口长度N决定了窗函数法设计出的FIR滤波器的性能。

设计时，要根据阻带的最小衰减和过渡带宽来选择恰当的窗函数类型和窗长度N。常用窗函数有矩形窗、海明窗、和布莱克曼窗等。 窗函数设计FIR滤波器步骤如下：

（1）给定理想频率响应Hd(ejw)的幅频特性和相频特性；

（2）求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ejw)采样，并对其求IDFT的hM(n)，用hM(n)代替hd(n)；

（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗长度N；

WN(n),求所需设计的FIR滤波器单位脉冲响应； （4）根据h(n)hd(n)（5）求Hd(ejw)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。

- 17 -

三．实验设备

计算机，MATLAB软件。

四．实验内容

(1)用海明窗设计一个13阶的FIR的低通滤波器，wc=0.13

程序如下： clear all;%海明窗 N=13; wc=0.13*pi;

h=fir1(N-1,wc/pi,'bandpass',hamming(N)); [h1,wc]=freqz(h,1);

subplot(2,1,1);plot(wc/pi,20*log10(abs(h1))); axis([0,1,-80,10]); grid;

xlabel('归一化频率/\\pi'); ylabel('幅度/dB'); title('N=13,海明窗'); 波形如下：

图4-1 用海明窗设计的低通滤波器

(2)用布莱克曼窗设计一个13阶的FIR的低通滤波器，wc=0.13

程序如下：

clear all;%布莱克曼窗 N=13; wc=0.13*pi;

h=fir1(N-1,wc/pi,'bandpass',blackman(N)); [h1,wc]=freqz(h,1);

subplot(2,1,1);plot(wc/pi,20*log10(abs(h1))); axis([0,1,-80,10]);

- 18 -

grid;

xlabel('归一化频率/\\pi'); ylabel('幅度/dB'); title('N=13,布莱克曼窗'); 波形如下：

图4-2 用布莱克曼窗设计的低通滤波器

(3)用矩形窗设计线性相位低通滤波器：0.13w,Hd(ejw)e-j（w-），其余为0。

程序如下： clear all;%矩形窗 N=13; wc=0.13*pi;

h=fir1(N-1,wc/pi,'bandpass',boxcar(N)); [h1,wc]=freqz(h,1);

subplot(2,1,1);plot(wc/pi,20*log10(abs(h1))); axis([0,1,-80,10]); grid;

xlabel('归一化频率/\\pi'); ylabel('幅度/dB'); title('N=13,矩形窗'); 波形如下：

- 19 -

图4-3 用矩形窗设计的线性相位低通滤波器

- 20 -