数字信号处理实验总结

实验一 离散信号及运算
一、 实验目的
1. 掌握MATLAB 语言的基本功能及实现方法；
2. 掌握MATLAB 中各种常用序列的表示和显示方法；
3. 熟练运用MATLAB 进行离散信号的各种运算。

二、 实验原理
我们所接触的信号大多为连续信号，而计算机及其他设备处理的大多为数字信号。

为了便于处理，往往要对信号进行处理使之变成离散数字信号。

对信号进行时间上的量化（即采样）是对信号作数字化处理的第一个环节，要求理解采样的原理和采样的性质，知道采样前后信号的变化及对离散信号和系统的影响。

三、 实验内容
1、用MATLAB 实现下列序列，并画出图形：
① 单位采样序列移位，100),3()(≤≤-=n n n x δ； 提示：实现单位采样序列：0001{
)(≠==n n n δ，可通过以下语句实现：x=zeros(1,N);x(1)=1; n=0:10;
x=[zeros(1,3),1,zeros(1,7)];
stem(n,x)； 01234567891000.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
② 单位阶跃序列移位，100),3()(≤≤-=n n u n x
提示：实现单位阶跃序列：0
001{)(≠==n n n u ，可通过以下语句实现：x=ones(1,N);
n=0:10; x=[zeros(1,3),1,ones(1,7)]; stem(n,x) 01234567891000.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
③正弦序列，
100),****2sin(*)(≤≤=n T n f A n x s π，其中A=2；f=10；
s T =0.005； A=2; f=10; Ts=0.005; n=0:10;
x=A*sin(2*pi*f*n*Ts); stem(n,x) 01234567891000.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
③ 指数序列，
100,9.0)(≤≤=n n x n n=0:10; x=0.9.^n; stem(n,x) 01234567891000.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
④ 复指数序列，0.05+j*pi/4*()e ,2020n x n n =-≤≤，画出该序列的实部、虚部，幅值和相位。

提示：可通过下列语句实现：
实部real(x)，虚部imag(x)，幅值abs(x)，相位angle(x)
n=-20:20; x=exp(0.05+j*pi/6*n);
xr=real(x); xi=imag(x); xm=abs(x); xa=angle(x);
figure;
subplot(411);stem(n,xr);title('实部');
subplot(412);stem(n,xi);title('虚部');
subplot(413);stem(n,xm);title('模');
subplot(414);stem(n,xa);title('相角'); -20
-15-10-505101520-20
2
实部-20
-15-10-505101520-20
2
虚部-20
-15-10-50510152001
2
模-20-15-10-505101520-505
相角
2、用MATLAB 实现两个序列相加：
序列1：x1=[1 0.5 0.3 0]，n1=1:4；
序列2：x2=[0.2 0.3 0.4 0.5 0.8 1]，n2=1:6；实现x=x1+x2，n=1:8，并画出x 的图形。

提示：MATLAB 中可用算术运算符“+”实现序列相加，但两个序列的长度必须相等。

如果序列长度不等，或者长度虽然相等但采样的位置不同，就不能运用“+”了。

当两序列的长度不等或位置不对应时，首先应使两者位置对齐，然后通过zeros 函数左右补零使其长度相等后再进行相加。

x1=[1,0.5,0.3,0];x2=[0.2,0.3,0.4,0.5,0.8,1];
n1=1:4; n2=1:6; n=1:8; x3=[x1,0,0]; x=x3+x2;
y=[x,0,0]; stem(n,y)
1234567800.2
0.4
0.6
0.8
1
1.2 1.4
3、用MATLAB 实现序列的反转：
实现)()(n x n y -=，序列x(n)采用100),3()(≤≤-=n n u n x ，并画出
y(n)的图形。

提示：可利用fliplr(x) 函数，例如：x=[1 2 3 4];y=fliplr(x);结果为y=[4 3 2 1]，要实
现)()(n x n y -=对fliplr(x) 函数进行合理运用。

n=0:10;
x=[zeros(1,3),1,ones(1,7)]; y=fliplr(x);
n1=-fliplr(n); stem(n1,y); -10-9-8-7-6-5-4-3-2-1000.1
0.20.3
0.4
0.5
0.6
0.7
0.8
0.9
1
4、序列的尺度变换，实现插值和抽取：
已知序列(){1,2,3,4,5,6,7,8},-43x n n =≤≤用MATLAB 分别实现下列尺度变换。

12()(2){()(/2)y n x n y n x n ==
提示：可对序列x(n)的下标进行取余计算，余数为零即为插值和抹去的点，函数如
下：mod(nx,m)，nx 为序列x(n)的下标，m 为插值或抽取的倍数。

clear all ;
x=[1,2,3,4,5,6,7,8];
n=-4:3;n1=n(1:2:length(n));
y1=x(1:2:length(x)); subplot(211); stem(n1,y1);
axis([-5 3 0 7]); y2=zeros(1,16);
for k=1:8
y2(2*k)=x(k);
End;
subplot(2,1,2);stem(-8:6,y2(2:end));
-5
-4-3-2-1012302
4
6-8-6-4-2024602
4
6
8
四、 思考题
1. 若用C 语言实现有限长序列的加法，编程如何实现？与MA TLAB 相比，优缺点，繁简度如何？
2. MATLAB 中可用算术运算符“+”实现序列相加，用算术运算符“*”实现序列相乘，试问用MATLAB 求任意序列相加、相乘时，分别应注意什么？
3. 分析stem(n,x)和plot(t,x)函数使用时的异同点。

实验二 离散时间系统的表达及计算
一、 实验目的
1. 掌握离散系统的性质、输入输出关系及卷积运算；
2. 掌握离散系统常系数线性差分方程的解法。

二、 实验原理
一个离散时间系统在数学上的定义是将输入序列x(n)映射成输出序列y(n)的唯一性变 换或运算。

它的输入是一个序列，输出也是一个序列，其本质是将输入序列转变成输出序列的一个运算。

在数字信号处理中，通常研究的离散系统大都为LSI 系统。

三、 实验内容
1.求两个序列的卷积（离散线性卷积）
对于一个LSI 系统，设其输入序列为)10()()(--=n u n u n x ，其中505≤≤-n 。

系统为单位冲
激响应
)(*9.0)(n u n n h ∧=，编程求出输出序列，并画出图形。

提示：输出序列应是输入序列和系统单位冲激相应序列的卷积，在Matlab 中可用conv 函数来实现两个序列的卷积。

但是使用conv 函数求卷积时，两序列卷积后的输出序列的下标不是任意的（具体是多少，查函数），所以需要自己构造下标。

n=-5:50;
u1=stepseq(0,-5,50);u2=stepseq(10,-5,50);
x=u1-u2; h=((0.9).^n).*u1;
[y,ny]=conv1(x,n,h,n); stem(ny,y);
function [x,n]=stepseq(n0,n1,n2);
if nargin~=3
disp('Usage:Y=stepseq(n0,n1,n2)');
elseif ((n0<n1)|(n0>n2)|(n1>n2))
error('arguments?must?satisfy?n1<=n0<=n2');
End;
n=[n1:n2];x=[(n-n0)>=0];
function [y,ny]=conv1(x,nx,h,nh)
if nargin~=4
disp('Usage:Y=conv_m(x,nx,h,nh)');
return ; end ;
nyb=nx(1)+nh(1); nye=nx(length(x))+nh(length(x));
ny=[nyb:nye]; y=conv(x,h); -2002040608010001
2
34
5
6
7
2.求两个序列的卷积（离散周期卷积）
已知序列~(){1,0,1,2,1}x n =，~
y(){1,1,0,1,2}n =，计算它们的周期卷积。

提示：将~()x n 和~y()n 的主值序列x(n)和y(n)作线性卷积，然后再将线性卷积序列以N 为周期进行周期延拓。

clear all ;
x=[1,0,1,2,1]; y=[1,1,0,1,2];
z=conv(x,y);
m=0:4;
h=[z(6:9),0]; h1=h+z(1:5); h2=[h1 h1 h1];
stem(h2)
05101501
2
3
4
56
3.已知LSI 系统的差分方程为：)()2(9.0)1()(n x n y n y n y =-+--，其中10020≤≤-n
求单位冲激响应)(n h ，并画出图形。

提示：可以使用filter 函数来完成，其调用格式为),,(x a b filter H =，注意b ，a ，x 参数都代表什么，如何设置。

clear all
b=[1]; a=[1 -1 0.9];
n=-20:100; N=120;
x=[zeros(1,20),1,zeros(1,100)]; y=filter(b,a,x);
stem(n,y); axis([-20 100 -2 2]); -20020406080100
-2-1.5
-1
-0.5
0.5
1
1.52
4.已知LSI 系统的差分方程为：)()2(9.0)1()(n x n y n y n y =-+--，其中10020≤≤-n
求单位阶跃响应)(n g ，并画出图形。

提示：可以使用filter 函数来完成，其调用格式为),,(x a b filter H =，注意b ，a ，x 参数都代表什么，如何设置。

b=[1]; a=[1 -1 0.9];
N=120; n=-20:100;
x=[zeros(1,20),1,ones(1,100)];
g=filter(b,a,x); stem(n,g); -2002040608010000.5
1
1.5
2
2.5
思考题
1. conv 函数来实现两个序列的卷积，其默认下标为？（10分）
2. 离散线性卷积与离散周期卷积的区别？（20分）
3. 解差分方程时，a 、b 分别代表什么？（10分）
当差分方程为)()1()(3/10)1(n x n y n y n y =++--时，a 、b 分别为多少？
4. 解释filter 函数？（10分）
实验三 z 变换及反变换
一．实验目的
1.通过Matlab 编程，熟悉z 变换定义及逆z 变换常用方法，加深对z 变换性质及其收敛域 的理解；
2.通过Matlab 编程实现离散信号及系统的z 域分析；
3.掌握利用Matlab 编程求解差分方程的方法。

二．实验原理
z 变换性质、求逆z 变换的方法、系统的z 域表示方法及差分方程的求解方法。

三．实验内容
1. z 正变换：已知序列x1=[1,0,5,6]，其中30≤≤n ；与x2=[1,3,5,2]，其中，30≤≤n ，求其卷积信号x=x1*x2的z 变换。

提示：例将序列x1进行Z 变换，则根据公式
∑=-⋅=Z 3
0)(1)](1[n n z n x n x 3210)3(1)2(1)1(1)0(1----⋅+⋅+⋅+⋅=z x z x z x z x
32106501----⋅+⋅+⋅+⋅=z z z z （与序列x1比较）
也就是说，没必要在Matlab 命令窗口中显示出完整序列z 变换的表达式，通过序列的位置信息和该位置信息上所对应的序列值，即可在报告中将该序列的z 变换写出来。

clear all
x1=[1,0,5,6];x2=[1,3,5,2];
syms z ;x=conv(x1,x2); s=0;
for n=1:length(x)
s=s+x(n)*z^(-n+1);
End;s
x = 1 3 10 23 43 40 12
s = 3/z + 10/z^2 + 23/z^3 + 43/z^4 + 40/z^5 + 12/z^6 + 1
2. z 反变换：已知1132)(-++=z z z X 与1225342)(-+++=z z z z X ，求)()()(213z X z X z X ⋅=的逆z
变换；
提示：思路与上题相同。

通过X 1(z)和X 2(z)找出序列x 1(n)和x 2(n)，直接求x 1(n)和x 2(n)的卷积，即是求X 1(z)·X 2(z)的逆z 变换。

注意序列x 1(n)和x 2(n)的n 的取值范围。

n1=-1:1; n2=-2:1;
x1=[1,2,3]; x2=[2,4,3,5];
y=conv(x1,x2); n=n1(1)+n2(1):n1(3)+n2(4); stem(n,y);
-3
-2.5-2-1.5-1-0.500.51 1.52
05
10
15
20
25
3.离散时间线性时不变系统的时域实现：已知因果LSI 系统的差分方程为：
)2()()2(81.0)(--+-=n x n x n y n y ，写出系统传递函数H(z)及其收敛域，编程实现以下内容：
① 系统冲激响应h(n)，并绘图； ② 系统的单位阶跃响应g(n)，并绘图； ③ 系统函数的幅频响应H(e jw )，并绘出幅频和相频特性。

提示：①冲激响应可用impz(b,a,n)函数，可查看Matlab Help 了解它的相关用法。

impz 函数：数字滤波器的冲激响应，用来求系统的冲激响应h(n)。

b 是差分方程中输入的系数，a 是差分方程中输出的系数。

注意缺项的系数前要添加系数0。

这个函数本身就有画图功能。

如果n 是一个整形变量，impz 计算的是在这些整数点位置上的冲激响应，且从0开始。

②单位阶跃响应可用stepz(b,a,n)函数，可查看Matlab Help 了解它的相关用法。

与impz 用法相似。

③幅频响应可用freqz （b,a,w ），其中w=[0:1:500]*pi/500和freqz(b,a,n)分别求， 可查看Matlab Help 了解它的相关用法。

①：b=[1 0 -1]; a=[1 0 -0.81]; impz(b,a,100);
0102030
405060708090
-0.2
0.2
0.4
0.6
0.8
1
1.2
n (samples)
A m p l i t u d e
Impulse Response
②：b=[1 0 -1]; a=[1 0 -0.81]; stepz(b,a,100);
10
20
30
40506070
80
90
00.10.20.30.40.50.60.7
0.80.9
1n (samples)
A m p l i t u d e
Step Response
③：a=[1,0,-0.81]; b=[1,0,-1]; freqz(b,a,50,'whole');
0.2
0.4
0.60.81 1.2 1.4 1.6
1.8
2
-40-2002040
Normalized Frequency (⨯π rad/sample)
P h a s e (d e g r e e s )
00.20.4
0.60.81 1.2 1.4 1.6 1.8
2
-2
-101
Normalized Frequency (⨯π rad/sample)
M a g n i t u d e (d B )
a=[1,0,-0.81]; b=[1,0,-1]; freqz(b,a,[0:1:500]*pi/500);
0.1
0.2
0.30.40.50.60.70.80.9
1
-100
-50050100
Normalized Frequency (⨯π rad/sample)
P h a s e (d e g r e e s )
0.1
0.2
0.30.40.50.60.70.80.9
1
-300
-200-1000100
Normalized Frequency (⨯π rad/sample)
M a g n i t u d e (d B )
四．思考题
1.画图语句后加grid，起什么作用？（10分）
2.画图时，不加axis语句，为什么也可以画出图形？此时横纵轴坐标怎样？（10分）
3.freqz（b,a,w）中，w=[0:1:500]*pi/500，若改为：
freqz(b,a,n)，n=[0:1:50]，会出现什么情况？为什么？（10分）
4.根据实验结果，写出系统函数的表达式H(z)，并写出收敛域。

（10分）
5.系统稳定性和因果性的判决条件？程序实现的主要思路？（10分）
实验四离散时间傅立叶变换DFT（离散时间信号的频域分析）
五、实验目的
4.深刻理解有限长序列的离散傅立叶变换的概念、计算及意义；
5.掌握从离散傅立叶级数到离散傅立叶变换的过程；
6.掌握离散傅立叶变换的性质，熟悉信号的频域转换。

六、实验原理
一个信号的时域转换可以通过傅立叶变换（DFT）来完成。

有限长序列的DFT是其Z 变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。

为了更好地理解DFT的概念，我们先了解周期序列及其离散傅立叶级数
（DFS ）表示。

七、 实验内容
1. 若)6
cos()(πn n x =是一个N=12的有限长序列，利用MATLAB 计算它的DFT 并画出图形。

提示：在计算机上实现信号的频谱分析及其他方面的处理工作时，对信号的要求是：在时域和频域都应该是离散的，且都应是有限长的，只有DFS （离散周期序列的傅立叶级数）在时域和频域都是离散的。

所以，DFT 实际上来是自于DFS 的，只不过仅在时域、频域各取一个周期而已。

本题要求自己编写一个),(N xn dft 函数，其中xn 代表序列，N 代表有限长序列的点数。

),(N xn dft 函数的算法采用序列求DFT 变换的基本公式，即：
∑
∑
-=-==-=1
01
)()2exp()()(N n N n nk N W n x nk N j n x k X π，其中，1,,1,0-=N k 。

由于X(k)为复数，所以画图时，应画X(k)的幅度特性abs()和相位特性angle()。

clc
N=12; n=0:N-1; k=0:N-1;
x=cos(n*pi/6); w=exp(-j*2*pi/N);
A=n'*k; X=x*w.^A; X1=abs(X); X2=angle(X); subplot(121);stem(n,X1);grid on ; subplot(122);stem(n,X2);grid on ;
0510150
1
2
3
4
5
6
7
051015
-4
-3-2-1
1
234
2. 对于实序列
200,)9.0()(<≤=n n x n
① 分解成偶部)(n x e 和奇部)(n x o ；
② 验证实序列的性质：)](Re[)]([k X n x DFT e = )](Im [)]([k X j n x DFT o = 提示：任意一个实序列都可以表示成偶对称序列与奇对称序列之和，即：)()()(n x n x n x o e +=。

令]))(()([21
)(N e n N x n x n x -+=，，这]))(()([2
1)(N o n N x n x n x --=样分别求)]([n x DFT e 和)]([n x DFT o ，再求)]([n x DFT ，将得到的X(k)取其)](Re[k X 和)](Im [k X j ，与上述结果相比较，验证实序列的上述性质。

本题要求自己编写一个circevod(x)函数，它的功能是将一个实序列分解成)(n x e 与)(n x o 之和，并求出)(n x e 与)(n x o 。

其中需要有判断语句，判断这个序列是否是实序列；
N n N x ))((-的求法可参考mod(-n,N)函数。

clear all ;
n=0:19; N=length(n); m=mod(-n,N)+1;
x1=0.9.^n; x2=x1(m); xe=1/2*(x1+x2); xo=1/2*(x1-x2); subplot(321);stem(n,xe);title('x1(n)的偶部'); subplot(322);stem(n,xo);title('x1(n)的奇部'); Xe=dft(xe,N); Xo=dft(xo,N);
X=dft(x1,N); ReX=real(X); ImX=imag(X); subplot(323);stem(n,Xe);title('xe(n)的DFT'); subplot(324);stem(n,ReX);title('X 的实部'); subplot(325);stem(n,imag(Xo));title('xo 的DFT'); subplot(326);stem(n,ImX);title('X 的虚部');
0510152000.51x1(n)的偶部
05101520
-0.500.5x1(n)的奇部
05
1015200510xe(n)的DFT
05
101520
0510X 的实部
05101520-5
05xo 的DFT
05101520
-5
05X 的虚部
离散傅里叶变换的应用
1用DFT 计算线性卷积：（了解） 算法流程如下：
补零到N 点
补零到N 点N 点DFT
N 点DFT
11()(01)
x n n N ≤≤-22()(01)
x n n N ≤≤-121
N N N =+-N 点IDFT
1X (k)
2X (k)
T(k)
()()
t n g n =
数字信号处理的优势是“实时实现”，即信号进来经处理后马上输出去。

然而：
()()()()()k y n x n h n x k h n k ∞
=-∞
=*=
-∑
引入两个问题：
● x(n)没有全部进入，如何实现卷积？ ● 全部进入再卷积，又如何保证实时实现？
解决方法通常是将较长序列进行分段，然后计算每段子序列和较短序列的线性卷积，最后再将各段线性卷积结果序列进行相加得到结果。

这类方法包括：重叠相加法和重叠
保留法。

2用DFT 对模拟信号进行谱分析
编程实现：给定模拟信号)8cos(5)4sin(2)(t t t x ππ+=，以)1:0(01.0-==N n n t 毫秒进行取样， 用DFT 进行信号频谱分析：
① 若要能区别该信号中的两个频率分量，试问取样信号的长度至少为多少？共取多少采样点？②所用DFT 的点数N 分别等于128、256、1024，画出信号的N 点DFT 的幅度谱。

③讨论幅度谱结果，N 为多少时能分辨出信号中的所有频率分量？ clear;
N1=64;n1=0:N1-1; N2=256;n2=0:N2-1; N3=1024;n3=0:N3-1; t=0.01*n1;
x=2*sin(4*pi*t)+5*cos(8*pi*t); X=dft(x,N1); subplot(321);plot(t,x);title(' x(t)---N=64');
subplot(322);plot(t,abs(X));title(' abs(X)---N=64'); t=0.01*n2;
x=2*sin(4*pi*t)+5*cos(8*pi*t); X=dft(x,N2); subplot(323);plot(t,x);title('x(t)---N=256');
subplot(324);plot(t,abs(X));title(' abs(X)---N=256'); t=0.01*n3;
x=2*sin(4*pi*t)+5*cos(8*pi*t); X=dft(x,N3); subplot(325);plot(t,x);title('x(t)---N=1024');
subplot(326);plot(t,abs(X));title(' abs(X)---N=1024');
00.20.4
0.60.8-10010 x(t)---N=64
00.20.40.60.8
100
200 abs(X )---N=64
12
3-10010x(t)---N=256
123
500
1000 abs(X )---N=256
051015-10
010x(t)---N=1024
051015
2000
4000 abs(X )---N=1024
八、 思考题
1. 圆周卷积与线性卷积的关系，二者是否相同？在什么条件下相同？
2. 离散信号的卷积运算有何用途？
3. DFT 变换中补零对信号频谱的分辨率有影响吗？补零后频谱图有何变化？
4. 实验内容3中，若要能区别该信号中的两个频率分量，取样信号的长度至少为多少？ N 分别为60、500时能否分辨出信号中所有频率分量？与（2）的结果比较。

实验五 离散傅立叶变换和快速FFT
一．实验目的
1.在理论学习的基础上，通过本实验，进一步加深对DFT的算法原理及性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的结果必然满足DFT的基本性质）；
2.熟悉并掌握按时间/频率抽取FFT算法原理和子程序的应用；
3.学习用FFT对连续信号和时域离散信号进行谱分析的方法；
4.了解应用FFT进行信号频谱分析过程中可能出现的问题（如：混叠、泄露、栅栏效应等），
分析其原因，以便在实际中正确应用FFT。

二．实验原理
1.FFT 原理：DFT如果直接计算的话，计算量非常大，而且不利于信号的实时处理，在
实际应用中遇到很大的困难，由此出现了很多快速的计算DFT的方法，在此我们以基2的时间抽取快速傅立叶算法为例。

2.混叠：序列的频谱是被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会
发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱，避免混叠现象的唯一方法是保证采样速率足够高。

3.泄露：用截短的序列来近似很长甚至是无限长的序列，这样可以使用较短的DFT对信
号进行频谱分析，所得的频谱是原序列频谱的扩展。

为了减少泄露的影响，可以选择适当的窗函数使频谱的扩散减至最小。

4.栅栏效应：DFT是对单位圆上Z变化的均匀采样，所以它不可能将频谱视为一个连续
函数，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样一些频谱的峰点或谷点被“栅栏”所拦住，不能被观察到。

减小栅栏效应的方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“栅栏”
的位置，使频谱的峰点或谷点暴露出来。

三．实验内容
1.求序列x=[5 2 7 4 1 1 3 9]的快速傅立叶变换，并绘出幅频特性曲线； 提示：采用fft(x,N)函数，N 为FFT 的点数，当序列长度小于N 时，系统自动在序列末尾补零；当序列长度大于N 时，系统自动截断序列多余的部分。

x=[5 2 7 4 1 1 3 9];
N=length(x); xk=fft(x,N); n=0:N-1; stem(n,abs(xk));
01234567
5
10
15
20
25
30
35
2.设一序列中含有两种频率成分，f 1=2 Hz ，f 2=2.05Hz ，采样率取为fs=10Hz ，即 x(n)=sin(2πf 1n/f s )+sin(2πf 2n/f s )，根据公式N
f s 2<
2
1f f -，要区分出这两种频率成分，N 必
须满足多少？
⑴取x(n)(0≤n<128)时，计算128点FFT ，并绘出幅频特性曲线；
⑵取x(n)(0≤n<128)时，补384个0，计算512点FFT ，并绘出幅频特性曲线； ⑶取x(n)(0≤n<512)时，计算512点FFT ，并绘出幅频特性曲线；
⑷分别改变FFT 变换的点数和采样时间，对结果进行分析，对比看有何不同。

提示：补零时可采用zeros 函数；(4)中，N ，t 的取值为N=512，t=0.1n ；N=1024，t=0.1n ；N=512，t=0.05n ；N=1024，t=0.05n 四种情况。

f1=2; f2=2.05; fs=10; N1=128; N2=512; n1=0:N1-1; n2=0:N2-1;
n3=n1*fs/N1; n4=n2*fs/N2; n5=n2*fs/N2; xn1=sin(2*pi*f1.*n1/fs)+sin(2*pi*f2.*n1/fs);
xn2=sin(2*pi*f1.*n2/fs)+sin(2*pi*f2.*n2/fs); xk1=fft(xn1,N1);xk2=fft(xn1,N2);xk3=fft(xn2,N2); subplot(311);plot(n3,abs(xk1));title('128点'); subplot(312);plot(n4,abs(xk2));title('512点补零'); subplot(313);plot(n5,abs(xk3));title('512点');
012345
678910
050100128点
01234
5678910
050100512点补零
012345678910
200400512点
f1=2; f2=2.05; fs1=10; fs2=20; N1=512; N2=1024;
n1=0:N1-1; n2=0:N2-1; n3=n1*fs1/N1; n4=n1*fs2/N1; n5=n2*fs1/N2; n6=n2*fs2/N2; xn1=sin(2*pi*f1.*n1/fs1)+sin(2*pi*f2.*n1/fs1); xn2=sin(2*pi*f1.*n2/fs1)+sin(2*pi*f2.*n2/fs1); xn3=sin(2*pi*f1.*n1/fs2)+sin(2*pi*f2.*n1/fs2); xn4=sin(2*pi*f1.*n2/fs2)+sin(2*pi*f2.*n2/fs2); xk1=fft(xn1,N1);xk2=fft(xn2,N1); xk3=fft(xn3,N2);xk4=fft(xn4,N2);
subplot(411);plot(n3,abs(xk1));title('512点 t=0.1n'); subplot(412);plot(n4,abs(xk2));title('512点 t=0.05n'); subplot(413);plot(n5,abs(xk3));title('1024点 t=0.1n');
subplot(414);plot(n6,abs(xk4));title('1024点 t=0.05n');
1
2
3
4567
8
9
10
020*******点 t=0.1n
2
4
6
8101214
16
18
20
020*******点 t=0.05n
1
2
3
4567
8
9
10
020********点 t=0.1n
2
4
6
8
10
12
14
16
18
20
050010001024点 t=0.05n
3.计算序列y(n)的自相关函数，其中：
x(n)=sin(2π2n/16)，y(n)=e(n)+x(n)，式中e(n)为一随机噪声。

提示：随机噪声e(n)可通过rand(1,N)函数构造，rand 函数产生[0,1]分布的随机数，均值为0.5，(1,N)表示的是1行N 列；本程序要求编写产生均值为0，分布在[-0.5,0.5] 的随机数；求y(n)的自相关函数，可采用xcorr(yn,yn)函数实现。

本实验要求取180≤≤n ，画出x(n)，y(n)和y(n)的自相关函数ry 的波形。

N=500; n=0:N-1;
xn=sin(2*pi*2.*n/16); en=rand(1,N)-0.5; yn=en+xn; zn=xcorr(yn,yn,'unbiased'); N1=length(xn)+length(yn); n1=0:N1-2;
subplot(3,1,1);plot(n,xn);title('xn') subplot(3,1,2);plot(n,yn);title('yn') subplot(3,1,3);plot(n1,zn);title('自相关')
050100150200250300350400450500
-101xn
050100150200
250300350400450500
-202yn
01002003004005006007008009001000
-1
01自相关
四．思考题
4. 怎样知道FFT 比DFT 快？（从理论上、程序上分析写）
5. 简述FFT 变换的点数和采样时间对频谱分析的影响？
6. 抽样定理对正弦信号的适用条件？以及点数N 的选取？
实验六 数字滤波器组的设计——IIR 数字滤波器的设计
一．实验目的
1. 通过实验熟悉各种数字滤波器的计算机实验方法；
2. 通过分析实验结果加深对各类数字滤波器频率特性的知性认识。

二．实验原理
数字滤波器组的设计实际上是确定其系统函数并实现的过程。

我们在进行滤波器设计时，需要确定其性能指标。

理想滤波器是物理不可实现的（由于从一个频带到另一个频带之间的突变）；要物理可实现，用从一个带到另一个带之间设置一个过渡带且在通带和止带内也不应该严格为1或0，应给以较小容限。

三．实验内容
设计一个界面，完成滤波器组的设计。

具体要求：
（1）归一化模拟低通滤波器的设计
模拟低通滤波器的技术指标有Ωp、Ωs、Ωc、αp、αs。

其中，Ωp为通带截止频率、Ωs为阻带起始频率、Ωc为3dB通带截止频率，αp为通带内的最大衰减，αs为阻带内的最小衰减，αp、αs一般用dB来表示。

归一化模拟低通滤波器的所有频率都用通带截止频率Ωp作归一化，即归一化后通带截止频率为1。

用λ表示归一化后模拟滤波器的频率，即：λ=Ω/Ωp。

a)巴特沃斯滤波器
3dB通带截止频率和阶数N可以完全描述一个巴特沃斯归一化模拟低通滤波器。

编程实现：滤波阶次分别取N=2、5、10、20，在一张图中绘出巴特沃斯低通模拟原型滤波器的平方幅频响应曲线。

提示：采用[z0,p0,k0]=buttap(N)：N为巴特沃斯归一化模拟低通滤波器的阶数，返回值为滤波器组的零点数z0、极点数组p0和增益k0。

[B,A]=zp2tf(z,p,k)该函数从系统的零极点及增益得到滤波器系统函数的分析多项式系数向量B和分母多项式系数向量A。

[H,w]=freqs(b,a,w)函数求模拟系统的频率响应，用法相似于freqz。

其中，w=0:2*pi*100:2*pi*30000，画幅频特性曲线时注意：①H包含幅频和相频，幅频要取abs(H)；②横轴若表示频率，应取w/(2*pi)。

clear all;
N=[2,5,10,20]; w=0:2*pi*100:2*pi*30000;
for i=1:4
[z,p,k]=buttap(N(i)); [b,a]=zp2tf(z,p,k); [h,w]=freqs(b,a,length(w)); plot(w,abs(h));axis([0 2 0 1]); hold on ; end ;
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
00.10.20.30.40.50.60.70.80.91
b) 切比雪夫I 型
编程实现：滤波阶次分别取N=2、4、6、8，在一张图中绘出切比雪夫I 型低通模拟原型滤波器的平方幅频响应曲线。

吃
提示：采用[z0,p0,k0]=cheb1ap(N,Rp)，N 为切比雪夫I 型归一化模拟低通滤波器的阶数，Rp 为通带最大衰减。

[B,A]=zp2tf(z,p,k)和[H,w]=freqs(b,a,w)函数；其中Rp=1dB 。

clear all ；
N=[2,4,6,8]; Rp=1; w=0:2*pi*100:2*pi*30000; for i=1:4
[z,p,k]=cheb1ap(N(i),Rp); [b,a]=zp2tf(z,p,k); [h,w]=freqs(b,a,length(w)); plot(w,abs(h));axis([0 2 0 1]) hold on ; end ;
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
00.10.20.30.40.50.60.70.80.91
（2）IIR 数字滤波器设计
实际数字滤波器指标
归一化低通滤波器设计指标
H a(p)
H a(s)
H (z)
MATLAB 或查表
频率变换
脉冲响应不变或双线性变换
IIR 数字滤波器的设计过程
①冲激响应不变法
编程实现：设计巴特沃斯低通滤波器，通带纹波为Rp=1dB ，通带上限角频率p ω=0.2π，阻带下限角频率s ω=0.3π，阻带最小衰减s a =15dB ，画出设计后的数字低通滤波器的幅频特性曲线。

提示：采用[N, OmgC]=buttord(OmgP,OmgS,Rp,As,‘s ’)计算巴特沃斯低通滤波器的阶数N 和3dB 通带截止频率、[bs,as]=butter(N, OmgC,‘s ’)计算模拟滤波器系统函数的多项式，[Bz,Az]=impinvar(bs,as,1/T) 计算数字滤波器系统函数的多项式，其中，采样周期T=1。

[H,w]=freqz(Bz,Az,w)函数，其中，w=0:0.05*pi*pi 。

冲激响应不变法设计IIR 滤波器时，如果设计指标以数字滤波器的指标给出时，设计结果与采样间隔无关。

即：不管T 的取值如何，Bz 、Az 的结果相同。

clear all ;
Rp=1; As=15; w=0:0.005*pi:pi; T=1; wp=0.2*pi; ws=0.3*pi; OmgP=wp/T; OmgS=ws/T;
[N,OmgC]=buttord(OmgP,OmgS,Rp,As,'s'); [bs,as]=butter(N,OmgC,'s');
[bz,az]=impinvar(bs,as,1/T); [H,w]=freqz(bz,az,w); plot(w/pi,abs(H));
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
00.10.20.30.40.50.60.70.80.91
②双线性z 变换法
编程实现：设计巴特沃斯低通滤波器，p f =100Hz ，s f =300Hz ，p a =3dB ，s a =20dB ，抽样频率Fs=1000Hz ，画出设计后的数字低通滤波器的幅频特性曲线。

提示：采用[N, OmgC]=buttord(wap,was,Rp,As,‘s ’)、[bs,as]=butter(N, OmgC,‘s ’)，说明：若Ω和w 在映射的过程中乘了系数2/T ，调用[Bz,Az]=bilinear(bs,as,Fs)，否则调用 [Bz,Az]=bilinear(bs,as,Fs/2)和[H,w]=freqz(Bz,Az,w)函数。

fs=1000; fp=100; fst=300; Rp=3; As=20; w=0:0.005*pi:pi;
wap=tan(fp*2*pi/fs/2)*fs*2; was=tan(fst*2*pi/fs/2)*fs*2; [N,OmgC]=buttord(wap,was,Rp,As,'s');
[bs,as]=butter(N,OmgC,'s'); [bz,az]=bilinear(bs,as,fs); [H,w]=freqz(bz,az,w); plot(w,abs(H));
0.5
1
1.5
2
2.5
3
3.5
00.10.20.30.40.50.60.70.80.91
四．思考题
1.冲激响应不变法设计滤波器的步骤？优缺点？
2.双线性z 变换法设计滤波器的步骤？优缺点？
实验六 数字滤波器组的设计——FIR 数字滤波器的设计
一、实验目的
1、掌握用窗函数法和频率采样法设计FIR 数字滤波器的原理和方法；
2、熟悉线性相位FIR 滤波器的幅频特性和相频特性；
3、了解不同窗函数对滤波器性能的影响。

4、了解过渡带采样点对滤波器性能的影响。

二、实验原理
FIR 数字滤波器的系统函数为1
0()()N n
n H z h n z --==∑，该滤波器总是因果稳定的系统，且
当满足()(1)h n h N n =--或()(1)h n h N n =---的对称条件时，FIR 滤波器具有线性相位。

其在数据通信、图像处理、语音信号处理等实际领域得到广泛应用。

FIR 数字滤波器的设计方法主要有窗函数法和频率采样法。

（1）窗函数法
滤波器设计是以理想滤波器的特性为基准，使待设计的实际滤波器的特性逼近理想滤波器的特性。

以低通滤波器为例，一个截止频率为c ω、相位响应为零的理想低通滤波器的频率响应为：
j d 1
(e )0
c
c H ω
ωωωωπ≤⎧⎪=⎨
<≤⎪⎩
求离散时间傅里叶反变换可得理想冲激响应为：
--sin 11()()22c c j j n
j n c d d n h n H e e d e d n n
πωωωωπωωωωπππ===-∞<<+∞⎰⎰
理想冲激响应()d h n 是时间上的无限长序列，信号的绝大部分能量集中在零时刻附近，随着时间的增加，能量越来越小，最后趋于零。

因此，可以用一个窗函数对()d h n 进行截短，原则上尽量保留()d h n 的最高能量部分，以达到用时间上的有限长序列逼近理想冲激响应的目的，这就是窗函数法设计FIR 数字滤波器的基本思想。

若窗函数用()w n 表示，则有：
()()()w d h n h n w n =
从物理可实现的角度考虑，此时得到的有限长序列()w h n 并非因果序列，需要将其进行移位，最终得到FIR 数字低通滤波器的冲击响应序列()h n ，即：
()()w h n h n N =-
若相位响应不为零，而是令()ϕωωτ=-，即()ϕω具有线性相位，则：
-j j d e (e )0ωτωωωωωπ⎧≤⎪=⎨<≤⎪⎩
c c H ，其中1
2τ-=N ，N 为滤波器阶数。

-sin[()]
1()2()ωωτ
ωωωτωπ
πτ--=
=
-⎰c
c
j j n
c d n h n e
e d n
窗函数法的设计思想是从时域出发，截取有限长的一段冲激响应作为()H z 系数，冲激响应长度N 就是系统函数()H z 的阶数。

只要N 足够长，并且截取方法合理，总能够满足频域的要求。

窗函数法设计FIR 数字滤波器的基本步骤为：
Step1：确定数字滤波器理想特性()j d H e ω
；
Step2：由()j d H e ω
求出单位冲激响应()d h n ；。

Step3：由过渡带宽及阻带最小衰减的要求选择适当的窗函数，并根据线性相位条件确定窗函数的长度N ；
如()h n 偶对称，N 为奇数，或()h n 奇对称，N 为偶数。

ωωω∆=-p s ，对于海明窗，2 3.3
N0πω⨯=∆，若选()h n 偶对称，N 为奇数，在MATLAB 中，可用下面语句实现：
deltaw=ws-wp;%求过渡带宽度
N0=ceil(6.6*pi/deltaw);%/窗宽由过渡带决定，课本342页表7-3， N=N+mod(N0+1,2);
wdham=(hamming(N))';%求窗函数
在MATLAB 中，可由w=boxcar(N)（矩形窗，返回值w 是一个长度为N 的矩形窗序列）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗，最有用的窗结构之一，对于给定的阻带衰减，它能够提供最陡峭的过渡带，beta 是窗函数的形状参数，通过改变beta 可以对主瓣宽度）等函数来实现窗函数设计法中所需的窗函数。

Step4：求出单位冲激响应h (n )：h (n )=h d (n )w (n ) 0≤n ≤N -1
Step5：对h (n )作离散时间傅立叶变换，得到()j H e ω，分析幅频特性是否满足要求。

（2）频率采样法
频率采样法是先对理想频响
()
j d H e ω进行采样，得到采样值
()
H k ，即
2()()(
),
0,1,,1j d d k N
H k H k H e k N ωπ
ω====-，再对频域序列()H k 做IDFT 得到时域的滤波器
序列：()[()]h n IDFT H k =，这种设计过程是在频域进行抽样得到离散的频域点，称之为频率抽样法。

频率抽样法的设计思路：
()
j d H e ω()(2/)
d H k H k N π=()
h n 频率取样
IDFT
FIR 数字滤波器的频率抽样法设计思路图
频率采样法设计FIR 数字滤波器的基本步骤为：
Step1：根据设计要求选择滤波器的种类，确定合适的滤波器阶数N [0,2)π； Step2：然后对其在上进行N 点等间隔采样得到()H k ；
Step3：将()H k 带入内插公式得到所设计滤波器的频率响应()j H e ω；
在 MATLAB 中可由函数h=real(ifft(H,N))和实现。

频率采样法得到滤波器，在采样点上滤波器的实际频率响应是严格地和理想频率响应数值相等的。

但是在采样点之间的频响则是由各采样点的加权内插函数的延伸叠加而成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状。

理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小。

为解决上述问题，可在频率响应的过渡带内插入一个、两个或三个采样点，增加过渡带，减小频带边缘的突变，减小通带和阻带的波动，因而可增大阻带最小衰减。

三、实验内容
1、数字滤波器的技术指标如下：0.4P ωπ=、0.5P dB α=、0.6s ωπ=、50S dB α=，采用窗函数法来设计一个FIR 数字滤波器。

画出各种滤波器的幅频和相频特性曲线。

（1）选择一个合适的窗函数，取 N ＝15，观察所设计滤波器的幅频特性，分析是否满足设计要求；
（2）取N ＝35，重复上述设计，观察幅频和相频特性的变化，分析长度N 变化的影响；
clear all ;
N=35; wp=0.4*pi; ws=0.6*pi; wc=1/2*(wp+ws); n=[0:(N-1)]; tao=(N-1)/2; m=n-tao+eps; hd=sin(wc*m)./(pi*m); deltaw=ws-wp; N0=ceil(6.6*pi/deltaw); N=N+mod(N0+1,2); wdham=(hamming(N))'; h=hd.*wdham; [H,w]=freqz(h,[1],1000,'whole');
H=(H(1:1:501))'; w=(w(1:1:501))'; mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); subplot(221);stem(n,hd); subplot(222);stem(n,wdham); subplot(223);stem(n,h); subplot(224);plot(w,db);
10
20
30
40
-0.2
00.20.40.60
10
20
30
40
00.20.40.60.81
10
20
30
40
-0.2
00.20.40.60
1
2
3
4
-150
-100-50050
（3）保持 N ＝35不变，改变窗函数为blackman 窗，观察并记录窗函数对滤波器幅频特性的影响，比较两种窗的特点。

提示：m=n-tao+eps; %eps 为一个很小的数，避免下一步中分母为0；
hd=sin(wc*m)./(pi*m); %求理想脉冲响应。

clear all ;
N=35; wp=0.4*pi; ws=0.6*pi; wc=1/2*(wp+ws); n=[0:(N-1)]; tao=(N-1)/2; m=n-tao+eps;
hd=sin(wc*m)./(pi*m); deltaw=ws-wp; N0=ceil(11*pi/deltaw); N=N+mod(N0+1,2); wdham=(blackman(N))'; h=hd.*wdham; [H,w]=freqz(h,[1],1000,'whole');
H=(H(1:1:501))'; w=(w(1:1:501))'; mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); subplot(221);stem(n,hd); subplot(222);stem(n,wdham); subplot(223);stem(n,h); subplot(224);plot(w/pi,db);
10
20
30
40
-0.2
00.20.40.60
10
20
30
40
00.20.40.60.81
10
20
30
40
-0.2
00.20.40.60
0.5
1
-200
-150-100-50050
2、数字滤波器的技术指标如下：0.4P
ωπ=、0.5P dB α=、0.6s ωπ=、50S dB α=，
采用频率采样法来设计一个FIR 数字滤波器。

画出各种滤波器的频率样本、设计滤波器的脉冲响应和幅频特性。

（1）取 N ＝35，观察所设计滤波器的幅频特性，最小阻带衰减为多少，是否满足设计要求？
N=35; N=N+mod(N+1,2);
wc=0.5*pi; k=0:N-1; w=2*pi*k/N; N1=fix(wc/(2*pi/N)); N2=N-2*N1-1; A=[ones(1,N1+1),zeros(1,N2),ones(1,N1)]; theta=-pi*[0:N-1]*(N-1)/N;
H=A.*exp(j*theta); h=real(ifft(H));。