数字信号处理实验报告

姓名： 学号： 0905130129 专业班级：通信工程1301班 学院： 信息科学与工程学院 指导老师： 李宏

目录

实验一 常见离散时间信号的产生和频谱分析 ................. 3 一、实验目的 .......................................... 3 二、实验原理 .......................................... 3 三、实验内容及要求 .................................... 6 四、实验用MATLAB函数介绍 ............................. 7 五、 实验代码、结果及分析 ............................. 7 实验二 数字滤波器的设计 ............................... 20 一、实验目的 ......................................... 20 二、实验原理 ......................................... 20 三、实验内容及要求 ................................... 24 四、实验用MATLAB函数介绍 ............................ 24 五、 实验代码、结果及分析 ............................ 25

实验一 常见离散时间信号的产生和频谱分析

一、实验目的

（1） 熟悉MATLAB应用环境，常用窗口的功能和使用方法； （2） 加深对常用离散时间信号的理解； （3） 掌握简单的绘图命令； （4） 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶

变换对离散信号进行频域分析。

二、实验原理

（1） 常用离散时间信号

a）单位抽样序列

n01 (n)

0n0如果(n)在时间轴上延迟了k个单位，得到(nk)即：

(nk)

10nk n0b）单位阶跃序列

n01 u(n) n0010nN1c）矩形序列 RN(n)

其他0d）正弦序列

x(n)Asin(wn)

e）实指数序列

nx(n)au(n)

f）复指数序列

x(n)e(jw)n

（2）离散傅里叶变换：

设连续正弦信号x(t)为

x(t)Asin(0t)

这一信号的频率为f0，角频率为02f0，信号的周期为T012f00。

如果对此连续周期信号x(t)进行抽样，其抽样时间间隔为T，抽样后信号以x(n)表示，则有x(n)x(t)tnTAsin(0nT)，如果令w为数字频率，满足w00T012f0，其中fs是抽样重复频率，简称抽样频率。

fsfs为了在数字计算机上观察分析各种序列的频域特性，通常对

X(ejw)在0,2上进行M点采样来观察分析。 对长度为N的有限长序

列x(n), 有

X(ejwk)x(n)ejwkn

n0N1其中 wk2k，k0,1,,M1 Mk通常M应取得大一些，以便观察谱的细节变化。取模|X(ejw)|可绘出幅频特性曲线。

（3）用DFT进行普分析的三种误差

三种误差：混叠现象、泄露现象、栅栏效应

a) 混叠现象

当采样频率小于两倍信号（这里指是信号）最大频率时，经

过采样就会发生频谱混叠，这使得采样后的信号序列频谱不能真实地反映原信号的频谱。所以在利用DFT分析连续信号的频谱时，必须注意这一问题。避免混叠现象的唯一方法是保证采样速率足够高，使频谱交叠现象不致出现。也就是说，在确定采样频率之前，必须对信号的性质有所了解，一般在采样前，信号通过一个防混叠低通滤波器。 b) 泄漏现象

实际中的信号序列往往很长,为了方便我们往往用截短的

序列来近似它们,这样可以使用较短的DFT来对信号进行频谱分析,这种截短等价于给原信号序列乘以一个矩形窗函数。

泄漏是不能与混叠完全分离开的,因为泄漏导致频谱的扩散，从而造成混叠。为了减小泄漏的影响，可以选择适当的窗函数，使频谱的扩散减到最小。 c) 栅栏效应

因为DFT是对单位圆上Z变换的均匀采样，所以他不可能

将频谱视为一个连续函数。这样就产生了栅栏效应，就一定意义上看，DFT来观看频谱就好像通过一个尖桩的栅栏来观看一个图景一样，只能在离散点上看到真实频谱，这样就可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所挡住，不能被我们观察到。减小栅栏效应的一个方法就是借助在原序列的末端添补一些零值，从而变动DFT的点数。这一方法实际上是人为地改变了对真实谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位

置，从而使得频谱的峰点或者谷点暴露出来。当然，这是每根谱线所对应的频率和原来的不同了。

综上所述，DFT可以用于信号的频谱分析，但必须注意可

能产生的误差，在应用过程中要尽可能减少和消除这些误差的影响。

三、实验内容及要求

（1）复习常用离散时间信号的有关内容；

（2） 用MATLAB编程产生上述任意3种序列（长度可输入确定，对(d) (e)(f)中的参数可自行选择），并绘出其图形；

（3） 混叠现象 对连续信号x(t)sin(2*pi*f01*t)其中，f01500Hz进行采样，分别取采样频率fs2000Hz,1200Hz,800Hz，观察|X(ejw)|的变化，并做记录(打印曲线),观察随着采样频率降低频谱混叠是否明显存在，说明原因。 （4）截断效应

给定x(n)cos(n)，截取一定长度的信号y(n)x(n)w(n)，w(n)为

4窗函数，长度为N，w(n)RN(n)。分别取N=6，8，12，计算y(n)的N点DFT变换，画出其幅频特性曲线；做2N点DFT变换，分析当N逐渐增大时，分析是否有频谱泄露现象、主瓣的宽度变化？如何减小泄露？

（5）栅栏效应

给定x(n)R4(n)，分别计算X(ejw)在频率区间0,2上的16点、

32点、64点等间隔采样，绘制X(ejw)采样的幅频特性图，分析栅栏效应，如何减小栅栏效应？

四、实验用MATLAB函数介绍

（1）数字信号处理中常用到的绘图指令(只给出函数名，具体调用格式参看help)

figure(); plot(); stem(); axis(); grid on; title(); xlabel(); ylabel(); text(); hold on; subplot() （2）离散时间信号产生可能涉及的函数

zeros(); ones(); exp(); sin(); cos(); abs(); angle(); real(); imag();

五、实验代码、结果及分析

（1）复习常用离散时间信号的有关内容；

在时间上依次出现的数值序列,例如，{…，0.5，1，2，-1，0，5，…}。

相邻两个数之间的时间间隔可以是相等的，也可以是不等的。在前一情况下，设时间间隔为T秒,则离散信号可用符号x(nT)来表示。在间隔T归一化为1的条件下,T可以省略,即将x(nT)表示为x(n)。x(n)既可表示整个序列, 也可表示离散信号在nT瞬间的值。大多数离散时间信号是由对连续时间信号采样得到的，取值上可以仍然取连续值。

（2）用MATLAB编程产生上述任意3种序列（长度可输入确定，对(d)

(e)(f)中的参数可自行选择），并绘出其图形； 1、单位阶跃序列

n=-10:10;

xn=heaviside(n); xn(n==0)=1; plot(n,xn); stem(n,xn);

axis([-10 10 0 1.1]); title('单位阶跃序列'); xlabel('n'); ylabel('u(n)')

2、正弦序列

n=-15:0.3:15; x=sin(pi*n/5); stem(n,x);

line([-16,16],[0,0]) axis([-16,16,-1.2,1.2]); title('正弦序列') xlabel('n'); ylabel('x(n)')

3、实指数序列

n=0:10;

xn=(1.5).^n; plot(n,xn); stem(n,xn);

axis([0 10 0 65]); title('实指数序列'); xlabel('n'); ylabel('u(n)')

（3） 混叠现象 对连续信号x(t)sin(2*pi*f01*t)其中，f01500Hz进行采样，分别取采样频率fs2000Hz,1200Hz,800Hz，观察|X(ejw)|的变化，并做记录(打印曲线),观察随着采样频率降低频谱混叠是否明显存在，说明原因。 ①采样频率为2000Hz

N=200; T=0.1;

t=linspace(0,T,N); x=sin(2*pi*500*t); dt=t(2)-t(1); f=1/dt; X=fft(x); F=X(1:N/2+1); f=f*(0:N/2)/N; subplot(2,1,1)

plot(t,x)

title('x=sin(2*pi*500*t)') xlabel('t') ylabel('振幅')

axis([0,0.1,-1,1]); subplot(2,1,2) plot(f,abs(F))

title('采样频率为2000Hz') xlabel('频率');

ylabel('|X(e^{jw})|')

②采样频率为1200Hz

N=120; T=0.1;

t=linspace(0,T,N); x=sin(2*pi*500*t); dt=t(2)-t(1); f=1/dt; X=fft(x);

F=X(1:N/2+1); f=f*(0:N/2)/N; subplot(2,1,1)

plot(t,x)

title('x=sin(2*pi*500*t)') xlabel('t') ylabel('振幅') axis([0,0.1,-1,1]); subplot(2,1,2) plot(f,abs(F))

title('采样频率为1200Hz') xlabel('频率');

ylabel('|X(e^{jw})|')

③采样频率为800Hz

N=80; T=0.1;

t=linspace(0,T,N); x=sin(2*pi*500*t); dt=t(2)-t(1); f=1/dt; X=fft(x);

F=X(1:N/2+1); f=f*(0:N/2)/N; subplot(2,1,1) plot(t,x)

title('x=sin(2*pi*500*t)') xlabel('t') ylabel('振幅') axis([0,0.1,-1,1]); subplot(2,1,2) plot(f,abs(F))

title('采样频率为800Hz') xlabel('频率');

ylabel('|X(e^{jw})|')

结果分析：因为连续信号的最高频率为f01500Hz，根据取样定理，采样频率必须满足Fs>=2f01，否则会在折叠频率Fs/2处出现频谱混叠。通过实验，当采样频率Fs为200HZ，1200HZ时，峰值在500HZ处，没有发生混叠。但当取样频率为800HZ时，峰值在300HZ处，产生了较大的误差。 （4）截断效应

给定x(n)cos(n)，截取一定长度的信号y(n)x(n)w(n)，w(n)为窗函数，

4长度为N，w(n)RN(n)。分别取N=6，8，12，计算y(n)的N点DFT变换，画出其幅频特性曲线；做2N点DFT变换，分析当N逐渐增大时，分析是否有频谱泄露现象、主瓣的宽度变化？如何减小泄露？

n=0:1:99; N=6;

u = [ones(1,N-1) zeros(1,100-N+1)]; xn=cos(pi/4*n).*u; w0=pi/4;

Xk=fft(xn,2*N);

kx=[-N:1:N-1]/(2*N)*2; figure subplot(2,1,1) stem(xn)

title('N=6截断信号') xlabel('Time index n'); ylabel('Amplitude'); subplot(2,1,2) plot(kx,abs(Xk)) title('频谱图') xlabel('w / pi'); ylabel('|X(e^jw)|');

结果分析：由以上多图可知，随着矩形窗的N增大，主瓣宽度变

窄，分辨率提高，泄露也相继减少；随N减少，即取样长度比较小时，波形泄露比较严重，无法反映实际波形。另外，为了减小谱间干扰，应用其他形状的窗函数w（n）代替矩形窗会降低泄露程度。 （5）栅栏效应

给定x(n)R4(n)，分别计算X(ejw)在频率区间0,2上的16点、32点、64点等间隔采样，绘制X(ejw)采样的幅频特性图，分析栅栏效应，如何减小栅栏效应？

n=0:1:10;

xn=[ones(1,4),zeros(1,7)]; Xk16=fft(xn,16); Xk32=fft(xn,32); Xk64=fft(xn,64); subplot(2,2,1); stem(n,xn);

title('(a) x_1 (n)'); xlabel('n');

ylabel('x_1 (n)'); k=0:15; wk=2*k/16;

subplot(2,2,2); stem(wk,abs(Xk16));

title('(c) 16点DFT的幅频特性图'); xlabel('\\omega/\\pi'); ylabel('幅度'); k=0:31; wk=2*k/32;

subplot(2,2,3); stem(wk,abs(Xk32));

title('(d) 32点DFT的幅频特性图'); xlabel('\\omega/\\pi'); ylabel('幅度'); k=0:63; wk=2*k/64;

subplot(2,2,4); stem(wk,abs(Xk64));

title('(d) 64点DFT的幅频特性图'); xlabel('\\omega/\\pi'); ylabel('幅度');

结果分析：由栅栏效应可知，采样点的间隔内的情况是看不到的。所以对于有限长序列，我们在原序列尾部补零，从而增大频域采样点数和采样点位置。只要采样频率Fs足够高，且采样点数足够多，分辨率就会提高到认为DFT的离散谱近似代表原信号频谱。

实验二 数字滤波器的设计

一、实验目的

（1）熟悉用双线性变换法和脉冲响应不变法设计IIR数字滤波器的原理与方法； （2）学会调用MATLAB信号处理工具中滤波器设计函数，设计各种IIR滤波 器，学会根据滤波需求确定滤波器指标参数； （3）掌握用窗函数法设计FIR数字滤波器的原理和方法；

二、实验原理

（一）IIR滤波器 模拟滤波器Ha(s)设计

巴特沃兹滤波器的振幅平方函数A为

2A2()Ha(j)211()2Nc （1）

其传输函数为

Ha(s)NCN1p0

SpCe(ss12p1j()22Np) （2）

（3）

p0,1,(2N1)

首先确定技术指标： i. ii.

通带中允许的最大衰减ap和通带截止频率p； 阻带允许的最小衰减as和阻带起始频率s。

由式（8-11）可得：

pap10log10log12Ha(jp)C12N

（4）

sas10log10log12Ha(js)C12N

（5）

得到

p1C2N10ap/10 （6） （7）

s1C2N10as/10再利用上面两式得到

psa/1010p110as/101 N令

sp，

10p1k10as/101

a/10则

NlogKlog （8）

已知p,ap,s,as，可由式(8)求出滤波器的阶数N。求出的N可能有小数部分一般取大于等于N的最小整数。关于3dB截止频率c，有时在技术指标中给出，如果没有给出可以按照式（6）或式（7）求出。

根据以上所述，巴特沃兹滤波器的设计步骤为： i. ii. iii. iv.

根据要求p,ap,s,as，由式(8)求出阶数N； 由式（6）或式(7)求出3dB截止频率c； 由式（3）求出N个极点； 由式（2）写出传递函数Ha(s)。

双线性变换法和脉冲响应不变法实现将数字频率和模拟频率相互映射，具体见教材相关介绍。 （二）FIR滤波器

设所希望得到的滤波器的理想频率响应为Hd(ejw)。那么FIR滤波

器的设计就在于寻找一个传递函数Hd(ejwjw在)h(n)ejw去逼近Hd(e)。

n0这种逼近中最直接的一种方法是从单位取样响应序列h(n)着手，使我们知道h(n)可以从理想频率响应h(n)逼近理想的单位取样响应hd(n)。

Hd(ejw)通过傅里叶反变换来得到，即：

Hd(e)hd(n)ejwnjwnhd(n)1220Hd(ejw)ejwndw （9）

但是一般来说，这样得到的单位取样响应hd(n)往往都是无限长序列；而且是非因果的。以一个截止频率为wc的线性相应位理想低通为例来说明。设低通滤波器的时延为，即：

jweHd(e)0jwwwcwcw （10）

则

hd(n)12wcwcejwejwndw

sin[wc(n)](n)

这是一个以为中心的偶对称的无限长非因果序列。这样一个无限长的序列怎样用一个有限长序列去近似呢？最简单的办法就是直接截取它的一段来代替它。例如把n0到nN1的一段截取来作为

h(n)，但是为要保证所得到的是线性相位滤波器。必须满足h(n)的对

称性，所以时延应该取h(n)长度的一半，即(N1)/2

h(n)0nN1h(n)d其它n0

这种直接截取的办法可以形象地想象为，h(n)好比是通过一个“窗口”所看到的一段hd(n)。h(n)中表达为hd(n)和一个“窗口函数”的乘积。在这里，窗口函数就是矩形脉冲函数RN(n)，即

h(n)hd(n)RN(n)

但是一般来说，窗口函数并不一定是矩形函数，可以在矩形以内还对hd(n)作一定的加权处理，因此，一般可以表示为

h(n)hd(n)w(n)

这里w(n)就是窗口函数。这种对理想单位取样响应加窗的处理对频率响应会产生以下三点影响：

a) 使理想特性不连续的边沿加宽，形成一过渡带，过渡带的宽度取决于窗口频谱的主瓣宽度。

b) 在过渡带两旁产生肩峰和余振，它们取决于窗口频谱的旁瓣；旁瓣越多，余振也越多；旁瓣相对值越大，肩峰则越强。 c) 增加截取长度N，只能缩小窗口频谱的主瓣宽度而不能改变旁瓣的相对值；旁瓣与主瓣的相对关系只决定于窗口函灵敏的形状。因此增加N，只能相对应减小过渡带宽。而不能改变肩峰值。肩峰值的大小直接决定通带内的平稳和阻带的衰减，对滤波器性能有很大关系。例如矩形窗的情况下，肩峰达8.95%，致使阻带最小衰减只有21分贝，这在工程上往往是不够的。怎样才能改善阻带的衰减特性呢？只能从改善窗口函数的形状上找出路，所以希

望的窗口频谱中应该减少旁瓣，使能量集中在主瓣，这样可以减少肩峰和余振，提高阻带的衰减。而且要求主瓣宽度尽量窄，以获得较陡的过渡带，然而这两个要求总不能兼得，往往需要用增加主瓣宽度带换取决瓣的抑制，于是提出了海明窗、凯塞-贝塞尔窗、切比雪夫窗等窗口函数。

三、实验内容及要求

(1) 设计IIR数字滤波器，绘制幅频响应特性曲线H(ejw)。

基于脉冲响应和双线性变换设计低通数字滤波器：要求其通带截至频率100Hz，阻带截至频率200Hz，通带衰减Rp小于2dB，阻带衰减Rs大于15dB，采样频率Fs=500HZ。

(2) 设计FIR数字滤波器，绘制幅频响应特性曲线H(ejw)。

(a)用矩形窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率wp0.2rad，阻带截止频率ws0.3rad，阻带衰减为20dB，通带衰减为3dB，

(b) 当阻带衰减不小于40dB，通带衰减不大于3dB，应选择何种窗函数，画出其幅度响应H(ejw)的曲线。

四、实验用MATLAB函数介绍

（1）数字信号处理中常用到的绘图指令(只给出函数名，具体调用格式参看help) figure(); plot();grid on; title(); xlabel(); ylabel();hold on; subplot()

（2）滤波器设计可能用到的函数

buttord() ; buttap(); zp2tf(); lp2lp(); lp2bp(); lp2bs(); lp2hp(); bilinear(); freqz(); butter(); impinvar(); fir1(); ceil(); load()

五、实验代码、结果及分析

（1）设计IIR数字滤波器，绘制幅频响应特性曲线H(ejw)。 基于脉冲响应和双线性变换设计低通数字滤波器：要求其通带截至频率100Hz，阻带截至频率200Hz，通带衰减Rp小于2dB，阻带衰减Rs大于15dB，采样频率Fs=500HZ。 ①脉冲响应不变法

Fs=500; fp=100; fs=200;

Wp=fp*2*pi/Fs; Ws=fs*2*pi/Fs; ap=2; as=15;

[n,wc]=buttord(Wp,Ws,ap,as,'s'); [b,a]=butter(n,wc,'s'); [bz,az]=impinvar(b,a); [h,w]=freqz(bz,az,Fs); plot(w/pi,abs(h)) xlabel('频率w/pi'); ylabel('幅度/dB');

title('脉冲响应不变法');

②双线性变换法

Rp=2; Rs=15; Fs=500; Ts=1/Fs;

wp=100*2*pi*Ts; ws=200*2*pi*Ts; wp1=2*tan(wp/2)/Ts; ws1=2*tan(ws/2)/Ts;

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); [Z,P,K]=buttap(N); [B,A]=zp2tf(Z,P,K); [b,a]=lp2lp(B,A,Wn);

[bz,az]=bilinear(b,a,Fs); [H,W]=freqz(bz,az); plot(W/pi,abs(H)); xlabel('频率w/pi'); ylabel('幅度/dB'); title('双线性变换法');

（2）设计FIR数字滤波器，绘制幅频响应特性曲线H(ejw)。 (a)用矩形窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率wp0.2rad，阻带截止频率ws0.3rad，阻带衰减为20dB，通带衰减为3dB，

wp=0.2*pi; ws=0.3*pi; Bt=ws-wp;

N=ceil(4*pi/Bt); Wn=(0.2+0.3)*pi/2;

hn=fir1(N,Wn/pi,boxcar(N+1)) subplot(2,1,1) stem(hn,'.') xlabel('n') ylabel('h(n)') subplot(2,1,2)

[h,w]=freqz(hn,1,512);

plot(w/pi,20*log10(abs(h))) xlabel('w/pi') ylabel('幅度/dB')

(b) 当阻带衰减不小于40dB，通带衰减不大于3dB，应选择何种窗函数，画出其幅度响应H(ejw)的曲线。

wp=0.2*pi; ws=0.3*pi; Bt=ws-wp;

N=ceil(8*pi/Bt); Wn=(0.2+0.3)*pi/2;

hn=fir1(N,Wn/pi,hanning(N+1)) subplot(2,1,1) stem(hn,'.') xlabel('n') ylabel('h(n)') subplot(2,1,2)

[h,w]=freqz(hn,1,512);

plot(w/pi,20*log10(abs(h))) xlabel('w/pi') ylabel('幅度/dB')

结果分析：由阻带衰减不小于40dB可知，应该选择汉宁窗。