数字信号处理 实验报告
实验一 信号（模拟、数字）的输入输出实验
（常见离散信号产生和实现）
一、实验目的
1．加深对常用离散信号的理解；
2．掌握matlab 中一些基本函数的建立方法。
二、实验原理 1. 单位抽样序列
δ(n ) =⎨
⎧1⎩0
n =0n ≠0
在MATLAB 中可以利用zeros()函数实现。
x =zeros (1, N );
x (1) =1;
如果δ(n ) 在时间轴上延迟了k 个单位，得到δ(n -k ) 即：
δ(n -k ) =⎨
2．单位阶跃序列
⎧1⎩0
n =k n ≠0
n ≥0⎧1
u (n ) =⎨
n
在MATLAB 中可以利用ones()函数实现。
x=ones(1,N)
3．正弦序列
x (n ) =A sin(2πfn /Fs +ϕ)
在MATLAB 中,
n=0:N-1;
x=A*sin(2*pi*f*n/Fs+fai)
4．复指数序列
x (n ) =r ⋅e j ϖn
在MATLAB 中,
n=0:N-1;
x=r*exp(j*w*n) 5．指数序列
x (n ) =a n
在MATLAB 中,
n=0:N-1;
x=a.^n
三、实验内容实现和图形生成 1．五种基本函数的生成 程序如下： (1)单位抽样序列
% 单位抽样序列和延时的单位抽样序列 n=0:10;
x1=[1 zeros(1,10)];x2=[zeros(1,5) 1 zeros(1,5)]; subplot(1,2,1);
stem(n,x1);xlabel ('时间序列n');ylabel('振幅');title('单位抽样序列x1');
subplot(1,2,2);
stem(n,x2); xlabel('时间序列n');ylabel('振幅');title('延时了5的单位抽样序列');
单位抽样序列x1
2
2
延时了5的单位抽样序列
1.5
1.5
11
振幅
0.5
振幅
5
时间序列n
10
0.5
00
-0.5-0.5
-1-1
5
时间序列n
10
（2）单位阶跃序列 n=0:10;
u=[ones(1,11)];
stem(n,u);xlabel ('时间序列n');ylabel('振幅');title('单位阶跃序列'); 所得的图形如下所示：
振幅
1
2
3
4
56时间序列n
7
8
9
10
（3）正弦函数 n=1:30;
x=2*sin(pi*n/6+pi/3);
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('正弦函数序列x=2*sin(pi*n/6+pi/3)');
21.510.5
振幅
0-0.5-1-1.5-2
时间序列n
(4)复指数序列 n=1:30; x=2*exp(j*3*n);
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('复指数序列x=2*exp(j*3*n)'); 图形如下：
复指数序列x=2*exp(j*3*n)
21.5
10.5
振幅
0-0.5
-1-1.5
-2
时间序列n
（5）指数序列 n=1:30;
x=1.2.^n;
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('指数序列x=1.2.^n');
指数序列x=1.2.n
250
200
150
振幅
100500时间序列n
2．绘出信号x (n ) =1. 5sin(2π*0. 1n ) 的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期? 程序如下： n=0:40;
x1=1.5*sin(2*pi*0.1*n);x2=sin(0.9*n); subplot(1,2,1);
stem(n,x1); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x1=1.5*sin(2*pi*0.1*n)'); subplot(1,2,2);
stem(n,x2); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x2=sin(0.9*n)'); 运行结果如下：
正弦序列x1=1.5*sin(2*pi*0.1*n)
正弦序列x2=sin(0.9*n)
振幅
振幅
10
2030
时间序列n
40
时间序列n
由上图看出：x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。理论计算中对第一个，N ＝2*ｐｉ/（0.1*ｐｉ）＝10，第二个0.9不是pi 的倍数，所以不是周期的。因此可以看出，实验结果和理论相符。
3．x(n)=[2,3,1,1,2,-1,0,3]，-2≤n ≤5；h(n)=[2,4,1,-2,0,-1]，-3≤n ≤2，手工计算和MATLAB 计算卷积y(n)=x(n)*h(n)。
4. 如果x (n ) 、h (n ) 的起点不为0，则采用conv_m计算卷积； 编写conv_m函数：
function[y,ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];
y=conv(x,h); %MATLAB自带的函数 在命令窗口输入：
x=[2,3,1,1,2,-1,0,3];nx=[-2:5]; h=[2,4,1,-2,0,-1];nh=[-3:2]; [y,ny]=conv_m(x,nx,h,nh)
stem(ny,y,'.');xlabel('时间序号n');title('卷积和y(n)=x(n)*h(n)');
卷积和y(n)=x(n)*h(n)
20
15
10
5
-5
-10-6
-4-2
02时间序号n
468
附：手工计算
根据不进位卷积运算得:y(n)=[4,14,16,5,3,3,-7,0,13,1,-5,0,-3], -5≤n ≤7，图形如下
四、心得体会
通过此次试验，我深切体会到matlab 在数字信号处理中应用广泛，试验中也掌握了matlab 中一些基本函数的建立方法 ，加深了对常用离散信号的理解 ，同时课本上抽象的知识也变得更加形象化，枯燥的理论知识瞬间变得生动有趣起来。
实验二 FFT频谱分析及应用
一、实验目的
1．通过实验加深对FFT 的理解；
2．熟悉应用FFT 对典型信号进行频谱分析的方法。 三、实验原理与方法
在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT 不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。
有限长序列的DFT 是其z 变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT 是DFT 的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。
在MATLAB 信号处理工具箱中的函数fft （x,N ）, 可以用来实现序列的N 点快速傅立叶变换。经函数fft 求得的序列一般是复序列，通常要求出其幅值和相位。MATLAB 中提供了求复数的幅值和相位的函数：abs 、angle ，这些函数一般和fft 同时使用。
1．模拟信号x (t ) =2sin(4πt ) +5cos(8πt ) ，以t =0. 01n (n =0:N -1) 进行采样，求： （1）N ＝40点FFT 的幅度频谱，从图中能否观察出信号的2个频谱分量? （2）提高采样点数，如N ＝128，再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT 频谱分析结果与理论上是否一致？
2．有限长序列x(n)={2,1,0,1,3};h(n)={1,3,2,1}，试利用FFT 实现由DFT 计算线性卷积，并与线性卷积直接计算（conv ）的结果进行比较。 1．程序1： %xh11 N=40;n=0:N-1; t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t); k=0:N/2;w=2*pi/N*k; X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');title('signal x(n)'); subplot(2,1,2);plot(w/pi,magX);title('FFT N=40'); xlabel('w (unit :pi)');ylabel('|X|');grid
signal x(n)
5
10
15
20FFT N=40
100
25
30
35
40
|X |
50
00.10.20.30.4
0.50.6w (unit :pi)
0.70.80.91
（2）clear clc %xh11
%N=40;%取40个点，时间太短，现象不够精确 N=1000;%取1000个点 n=0:N-1; T=0.01; Fs=1/T; t=T*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t); k=0:N/2;f=Fs/N*k;%w=2*pi/N*k; X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');title('signal x(n)'); subplot(2,1,2);plot(f,magX);%plot(w/pi,magX); title('FFT N=40');
xlabel('w (unit :pi)');ylabel('|X|'); grid
axis([0 5 0 3000])
signal x(n)
[**************]0FFT N=40
[**************]0
30002000
|X |
1000
00.511.52
2.5
w (unit :pi)
33.544.55
2．程序2： %xh12
% 用FFT 实现由DFT 计算线性卷积 x=[2 1 0 1 3];h=[1 3 2 1]; L=length(x)+length(h)-1; XE=fft(x,L);HE=fft(h,L); y1=ifft(XE.*HE);
%画出由圆周卷积计算线性卷积结果及误差 k=0:L-1;
subplot(211);stem(k,real(y1));axis([0 8 0 15]); title('Result of Linear Convolution'); xlabel('Time index k');ylabel('Amplitude');
y2=conv(x,h); error=y1-y2; subplot(212); stem(k,abs(error));
xlabel('Time index k');ylabel('Amplitude'); title('Error Magnitude');
Result of Linear Convolution
A m p l i t u d e
0x 10
-16
123
45Time index kError Magnitude
678
A m p l i t u d e
Time index k
线性卷积直接计算：
function[y,ny]=conv_m(x,nx,h,nh) %改进卷积程序 nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];
y=conv(x,h); %MATLAB自带的函数 在命令窗口输入： x=[2,1,0,1,3];nx=[0:5];
h=[1,3,2,1];nh=[0:3]; [y,ny]=conv_m(x,nx,h,nh)
stem(ny,y,'.');xlabel('时间序号n');title('卷积和y(n)=x(n)*h(n)');
卷积和y(n)=x(n)*h(n)
012
34时间序号n
567
比较结果：这两个有限长序列的圆周卷积等于对两个序列的线性卷积。但当满足一定条件时才相等，若不满足，则这两个有限长序列的圆周卷积等于对两个序列的线性卷积在时间上的混叠。
实验三 IIR数字滤波器的设计
一、实验目的
1．掌握脉冲响应不变法和双线性变换法设计IIR 数字滤波器的原理和方法； 2．观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别。 三、实验原理与方法和手段 1．脉冲响应不变法
所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和ha(t)的采样值，即：h (n ) =h a (t ) t =nT =h a (nt ) ，其中，T 为采样周期。
在脉冲响应不变法中，模拟角频率和数字角频率的变换关系为：ω=ΩT ，
可见，Ω和ω之间的变换关系为线性的。
在MATLAB 中，可用函数impinvar 实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：
[B,A]=impinvar(b,a,fs1) [B,A]=impinvar(b,a)
其中，b 、a 分别为模拟滤波器的分子和分母多项式系数向量；fs1为采样频率（Hz ），缺省值fs=1Hz；B 、A 分别为数字滤波器分子和分母多项式系数向量。 2．双线性变换法：
21-z -1
由于s 平面和z 平面的单值双线性映射关系为s ＝，其中T 为采样
T 1+z -1
周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。
在双线性变换中，模拟角频率和数字角频率的变换关系为：Ω=见，Ω和ω之间的变换关系为非线性的。
在MATLAB 中，可用函数bilinear 实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：[B,A]=bilinear(b,a,fs1) 3．数字滤波器设计
（1）定技术指标转换为模拟滤波器设计性能指标。
（2）估计满足性能指标的模拟相应滤波器性能阶数和截止频率。 利用MATLAB 中buttord 、cheb1ord 、cheb2ord 、ellipord 等函数， 调用格式如：[N , Ωc ]=buttord (Ωp , Ωs , αp , αs , ' s ' )
其中，Ωp 为通带边界频率，rad/s；Ωs 为阻带边界频率，rad/s;αp 为带通波动，dB; αs 为阻带衰减，dB ；‘s ’表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωc 为模拟滤波器的截止频率(-3dB频率) ，rad/s。函数适用低通、高通、带通、带阻滤波器。 （3）设计模拟滤波器。
MATLAB 信号处理工具箱提供了模拟滤波器设计的完全工具函数：butter 、cheby1，cheby2、ellip 、besself 。用户只需一次调用就可完成低通、高通、带通、带阻滤波器设计。 调用格式如：[b,a]=butter(N,Ωc, ’ftype ’, ‘s ’) ，其中，’ftype ’
2ω
tg ，可T 2
为滤波器类型：
‘high ’表示高通滤波器，截止频率为Ωc ；
‘stop ’表示带阻滤波器，Ωc=[Ω1 Ω2] (Ω1
（４）利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。
五、实验步骤 具体步骤如下：
1、查看帮助文件，了解相关函数的调用格式。 2、实验题目
1. 设计模拟巴特沃斯低通滤波器，fp=300Hz，αp=1dB，fs=800Hz，αs=20dB。 %xh13模拟低通滤波器技术指标 αp=1; αs=20; fp=300; fs=800; Ωap=2*pi*fp; Ωas=2*pi*fs;
%设计模拟巴特沃斯低通滤波器
[N,Ωac]=buttord(Ωap, Ωas, αp, αs,'s')
[b,a]=butter(N,Ωac,'s') %设计模拟巴特沃斯低通滤波器，Ωap 为通带边界频率，rad/s；Ωas 为阻带边界频率，rad/s;αp 为通带最大衰减，dB ；αs 为阻带最小衰减，dB ；'s' 表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωac 为模拟滤波器的截止频率(-3dB频率) ，rad/s；b 、a 分别为模拟滤波器的系统函数分子和分母多项式系数向量； [H,Ω]=freqs(b,a); %求模拟滤波器的频率响应 %绘制频响幅度谱
plot(Ω/2/pi,20*log10(abs(H))); %横轴为频率，单位：HZ; 纵轴频响幅度，
单位：dB
axis([0,1500,-50,0]);
xlabel('频率Hz');ylabel(' H幅值dB');
H 幅值d B
[1**********]00
2：设计模拟巴特沃斯高通滤波器，fp=800Hz，αp=1dB，fs=300Hz，αs=20dB。 %xh14模拟高通滤波器技术指标 αp=1; αs=20; fp=800; fs=300; Ωp=2*pi*fp; Ωs=2*pi*fs;
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s')
[b,a]=butter(N,Ωc,'high','s') %设计模拟巴特沃斯高通滤波器，Ωp 为通带边界频率，rad/s；Ωs 为阻带边界频率，rad/s；αp 为通带最大衰减，dB ；αs 为阻带最小衰减，dB ；'s' 表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωc 为模拟滤波器的截止频率(-3dB频率) ，rad/s；b 、a 分别为模拟高通滤波器的系统函数分子和分母多项式系数向量； [H,Ω]=freqs(b,a);
plot(Ω/2/pi,20*log10(abs(H))); axis([0,1000,-50,0]);
10001200
频率Hz
[**************]0
xlabel('频率Hz');ylabel(' H幅值dB');
-10
-20
H 幅值d B
-30
-40
-50
-60
-70
[1**********]00
500600频率Hz
[1**********]00
3：fp=0.１kHZ, αp=1dB，fs=0.3kHZ，αs=25dB，T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth 数字低通滤波器。 %xh15采用冲击响应不变法 αp=1; αs=25; fp=100; fs=300; Ωp=2*pi*fp;
Ωs=2*pi*fs; %数字滤波器技术指标要求转化成模拟滤波器技术指标要求 fs1=1000; %采样频率 %设计模拟滤波器
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s'); [b,a]=butter(N,Ωc,'s');
[B,A]=impinvar(b,a,fs1) %用冲击响应不变法将模拟滤波器变换成数字滤波器。B 、A 分别为数字滤波器的系统函数分子和分母多项式系数向量； [H1,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应
%绘制数字滤波器频响幅度谱
subplot(211);plot(w*fs1/2/pi,20*log10(abs(H1))); axis([0,1000,-100,0]);
xlabel('频率Hz');ylabel(' H1幅值dB'); %采用双线性变换法
%数字滤波器的技术指标要求 αp=1; αs=25; fp=100; fs=300;
fs1=1000; %采样频率
wp=2*pi*fp/fs1; %数字滤波器通带边界频率 ws=2*pi*fs/fs1; %数字滤波器阻带边界频率 %变换为同类型模拟滤波器的技术指标要求
Ωp=2*fs1*tan(wp/2); %同类模拟滤波器通带边界频率 Ωs=2*fs1*tan(ws/2); %同类模拟滤波器阻带边界频率 %设计同类型模拟滤波器
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s'); [b,a]=butter(N,Ωc,'s');
%用双线性变换法将模拟滤波器变换成数字滤波器 [B,A]=bilinear(b,a,fs1)
[H2,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应 %绘制数字滤波器频响幅度谱
subplot(212);plot(w*fs1/2/pi,20*log10(abs(H2))); axis([0,1000,-100,0]);
xlabel('频率Hz');ylabel('H2幅值dB');
H 1幅值d B
-50
-100
[1**********]00
500600频率Hz
[1**********]00
H 2幅值d B
-50
-100
100
200
300
400500频率Hz
600
700
800
900
双线性变换法和脉冲不变法的特点和区别：
脉冲响应不变法：1，模拟频率到数字频率的转换时线性的。
2，数字滤波器单位脉冲响应的数字表示近似原型的模拟滤波器单位脉冲响应，因此时域特性逼近好，但会产生频谱混叠现象，只适合带限滤波器
双线性变换法：克服多值映射得关系，可以消除频率的混叠，但是因为是非线性的，所以在高频处有较大的失真。
实验四 FIR数字滤波器的设计
一、实验目的
1．掌握用窗函数法设计FIR 数字滤波器的原理和方法； 2．熟悉线性相位FIR 滤波器的幅频特性和相频特性； 3．了解不同窗函数对滤波器性能的影响。 三、实验原理与方法和手段
窗函数设计法是一种把一个长序列变成有限长的短序列的设计方法，是在时域进行的。用窗函数法设计FIR 数字滤波器时，先根据W c和N 求出相应的的理想滤波器的单位脉冲响应hd(n)。
1
hd (n ) =
2π
⎰
w c
-w c
H d (e ) e
jw jwn
sin(w c (n -a )) dw =
π(n -a )
因为hd(n)一般是非因果的，且无限长，物理上是不可实现的。为此可选择适当的窗函数w(n)截取有限长的hd(n)， 即h(n)= hd(n)w(n)，只要阶数足够长，截取的方法合理，总能够满足频域的要求。
实际中常用的窗函数有矩形（Boxcar ）窗、三角（Bartlett ）窗、汉宁(Hanning)窗、汉明(Hamming) 窗和布莱克曼（Blackman ）窗。这些窗函数各有优缺点, 所以要根据实际情况合理选择窗函数类型。
1．窗函数法设计线性相位FIR 滤波器的一般步骤为：
jw
H (e ) 的特性； （1）确定理想滤波器d
（2）由Hd （e jw ）求出h d (n ) ；
（3）根据过渡带宽度和阻带最小衰减，借助窗函数确定窗的形式及N 的大小，即选择适当的窗函数，并根据线性相位条件确定窗函数的长度N ；在MATLAB 中，可由w=boxcar(N)（矩形窗）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗）等函数来实现窗函数设计法中所需的窗函数。
（4）由h(n)= h d (n ) .w(n), 0≤n ≤ N-1，得出单位脉冲响应h(n); （5）对h(n)作离散时间傅立叶变换，得到H(e ) 。
2．在MATLAB 中，可以用b=fir1(N,Wn,’ftype’,taper) 等函数辅助设计FIR 数字滤波器。N 代表滤波器阶数；Wn 代表滤波器的截止频率(归一化频率) ，当设计带通和带阻滤波器时，Wn 为双元素相量；ftype 代表滤波器类型，如’high ’高通，’stop ’带阻等；taper 为窗函数，默认为海明窗，窗函数实现需要用窗函数blackman, hamming,hanning chebwin, kaiser产生。 五、实验步骤
在“开始--程序”菜单中，找到MATLAB 程序，运行启动；
jw
进入MATLAB 后 ，在Command Window中输入自己编写的主程序，并执行； 记录运行结果图形，作分析对比。
具体步骤如下：
1. 用窗函数法设计一线性相位FIR 低通滤波器，要求通带截止频率wc=冲击响应h （n ）的长度N ＝21。绘出h(n)及其幅频响应曲线.
sin((n -a ))
3 设计分析:理想滤波器单位冲击响应为hd (n ) =；为了满足线性相π(n -a ) N -1
位FIR 滤波器条件h(n)= h(N-1-n)，要求a==10。
2
π
，单位3
π
%xh19
N=21;wc=pi/3; %理想低通滤波器参数 n=0:N-1; a=(N-1)/2;
hdn=sin(wc*(n-a))/pi./(n-a); %计算理想低通滤波器单位冲击响应hd(n) if rem(N,2)~=0 hdn(a+1)=wc/pi; end %N为奇数时, 处理n=a点的0/0型 wn=hamming(N); %hamming窗 hn=hdn.*wn';
subplot(121); %绘出h(n)及幅频特性曲线 stem(n,hn,'.');
xlabel('n');ylabel('h(n)'); title('hamming窗设计的h(n)'); grid on; subplot(122); hw=fft(hn,512); w=2*[0:511]/512;
plot(w,20*log10(abs(hw)));
xlabel('w/pi');ylabel('Magnitude(dB)'); title('hamming窗设计的幅频特性'); grid on;
21
hamming 窗设
计的h(n)
20
hamming 窗设计的幅频特性
-20
M a g n i t u d e (d B )
5
10n
15
20
-40
h (n )
-60
-80
-100
-120
00.5
2．针对一个含有5Hz 、15Hz 和30Hz 的混和正弦波信号，设计一个FIR 带通滤波器。
参数要求:采样频率fs=100Hz，通带下限截止频率fc1=10Hz，通带上限截止频率fc2=20Hz，过渡带宽6Hz ，通阻带波动0.01，采用凯塞窗设计。 %xh20
fc1=10; fc2=20; fs=100;
[n,Wn,beta,ftype]=kaiserord([7 13 17 23],[0 1 0], [0.01 0.01 0.01],100) window=kaiser(n+1,beta); %使用kaiser 窗函数
b=fir1(n, Wn,window); %使用标准频率响应的加窗设计函数fir1 freqz(b,1,512); %数字滤波器频率响应 t=(0:100)/fs;
s=sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30); sf=filter(b,1,s); %对信号s 进行滤波 figure
subplot(2,1,1); plot(t,s) subplot(2,1,2); plot(t,sf)
22
1w/pi
1.52
50
M a g n i t u d e (d B )
-50-100-150
0.1
0.2
0.30.40.50.60.70.8Normalized Frequency (⨯π rad/sample)
0.9
1
500
P h a s e (d e g r e e s )
-500
-1000
00.10.2
0.30.40.50.60.70.8Normalized Frequency (⨯π rad/sample)
0.91
42
0-2-4
00.10.20.30.40.50.60.70.80.91
21
0-1-2
00.10.20.30.40.50.60.70.80.91
23
数字信号处理 实验报告
实验一 信号（模拟、数字）的输入输出实验
（常见离散信号产生和实现）
一、实验目的
1．加深对常用离散信号的理解；
2．掌握matlab 中一些基本函数的建立方法。
二、实验原理 1. 单位抽样序列
δ(n ) =⎨
⎧1⎩0
n =0n ≠0
在MATLAB 中可以利用zeros()函数实现。
x =zeros (1, N );
x (1) =1;
如果δ(n ) 在时间轴上延迟了k 个单位，得到δ(n -k ) 即：
δ(n -k ) =⎨
2．单位阶跃序列
⎧1⎩0
n =k n ≠0
n ≥0⎧1
u (n ) =⎨
n
在MATLAB 中可以利用ones()函数实现。
x=ones(1,N)
3．正弦序列
x (n ) =A sin(2πfn /Fs +ϕ)
在MATLAB 中,
n=0:N-1;
x=A*sin(2*pi*f*n/Fs+fai)
4．复指数序列
x (n ) =r ⋅e j ϖn
在MATLAB 中,
n=0:N-1;
x=r*exp(j*w*n) 5．指数序列
x (n ) =a n
在MATLAB 中,
n=0:N-1;
x=a.^n
三、实验内容实现和图形生成 1．五种基本函数的生成 程序如下： (1)单位抽样序列
% 单位抽样序列和延时的单位抽样序列 n=0:10;
x1=[1 zeros(1,10)];x2=[zeros(1,5) 1 zeros(1,5)]; subplot(1,2,1);
stem(n,x1);xlabel ('时间序列n');ylabel('振幅');title('单位抽样序列x1');
subplot(1,2,2);
stem(n,x2); xlabel('时间序列n');ylabel('振幅');title('延时了5的单位抽样序列');
单位抽样序列x1
2
2
延时了5的单位抽样序列
1.5
1.5
11
振幅
0.5
振幅
5
时间序列n
10
0.5
00
-0.5-0.5
-1-1
5
时间序列n
10
（2）单位阶跃序列 n=0:10;
u=[ones(1,11)];
stem(n,u);xlabel ('时间序列n');ylabel('振幅');title('单位阶跃序列'); 所得的图形如下所示：
振幅
1
2
3
4
56时间序列n
7
8
9
10
（3）正弦函数 n=1:30;
x=2*sin(pi*n/6+pi/3);
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('正弦函数序列x=2*sin(pi*n/6+pi/3)');
21.510.5
振幅
0-0.5-1-1.5-2
时间序列n
(4)复指数序列 n=1:30; x=2*exp(j*3*n);
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('复指数序列x=2*exp(j*3*n)'); 图形如下：
复指数序列x=2*exp(j*3*n)
21.5
10.5
振幅
0-0.5
-1-1.5
-2
时间序列n
（5）指数序列 n=1:30;
x=1.2.^n;
stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('指数序列x=1.2.^n');
指数序列x=1.2.n
250
200
150
振幅
100500时间序列n
2．绘出信号x (n ) =1. 5sin(2π*0. 1n ) 的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期? 程序如下： n=0:40;
x1=1.5*sin(2*pi*0.1*n);x2=sin(0.9*n); subplot(1,2,1);
stem(n,x1); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x1=1.5*sin(2*pi*0.1*n)'); subplot(1,2,2);
stem(n,x2); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x2=sin(0.9*n)'); 运行结果如下：
正弦序列x1=1.5*sin(2*pi*0.1*n)
正弦序列x2=sin(0.9*n)
振幅
振幅
10
2030
时间序列n
40
时间序列n
由上图看出：x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。理论计算中对第一个，N ＝2*ｐｉ/（0.1*ｐｉ）＝10，第二个0.9不是pi 的倍数，所以不是周期的。因此可以看出，实验结果和理论相符。
3．x(n)=[2,3,1,1,2,-1,0,3]，-2≤n ≤5；h(n)=[2,4,1,-2,0,-1]，-3≤n ≤2，手工计算和MATLAB 计算卷积y(n)=x(n)*h(n)。
4. 如果x (n ) 、h (n ) 的起点不为0，则采用conv_m计算卷积； 编写conv_m函数：
function[y,ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];
y=conv(x,h); %MATLAB自带的函数 在命令窗口输入：
x=[2,3,1,1,2,-1,0,3];nx=[-2:5]; h=[2,4,1,-2,0,-1];nh=[-3:2]; [y,ny]=conv_m(x,nx,h,nh)
stem(ny,y,'.');xlabel('时间序号n');title('卷积和y(n)=x(n)*h(n)');
卷积和y(n)=x(n)*h(n)
20
15
10
5
-5
-10-6
-4-2
02时间序号n
468
附：手工计算
根据不进位卷积运算得:y(n)=[4,14,16,5,3,3,-7,0,13,1,-5,0,-3], -5≤n ≤7，图形如下
四、心得体会
通过此次试验，我深切体会到matlab 在数字信号处理中应用广泛，试验中也掌握了matlab 中一些基本函数的建立方法 ，加深了对常用离散信号的理解 ，同时课本上抽象的知识也变得更加形象化，枯燥的理论知识瞬间变得生动有趣起来。
实验二 FFT频谱分析及应用
一、实验目的
1．通过实验加深对FFT 的理解；
2．熟悉应用FFT 对典型信号进行频谱分析的方法。 三、实验原理与方法
在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT 不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。
有限长序列的DFT 是其z 变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT 是DFT 的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。
在MATLAB 信号处理工具箱中的函数fft （x,N ）, 可以用来实现序列的N 点快速傅立叶变换。经函数fft 求得的序列一般是复序列，通常要求出其幅值和相位。MATLAB 中提供了求复数的幅值和相位的函数：abs 、angle ，这些函数一般和fft 同时使用。
1．模拟信号x (t ) =2sin(4πt ) +5cos(8πt ) ，以t =0. 01n (n =0:N -1) 进行采样，求： （1）N ＝40点FFT 的幅度频谱，从图中能否观察出信号的2个频谱分量? （2）提高采样点数，如N ＝128，再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT 频谱分析结果与理论上是否一致？
2．有限长序列x(n)={2,1,0,1,3};h(n)={1,3,2,1}，试利用FFT 实现由DFT 计算线性卷积，并与线性卷积直接计算（conv ）的结果进行比较。 1．程序1： %xh11 N=40;n=0:N-1; t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t); k=0:N/2;w=2*pi/N*k; X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');title('signal x(n)'); subplot(2,1,2);plot(w/pi,magX);title('FFT N=40'); xlabel('w (unit :pi)');ylabel('|X|');grid
signal x(n)
5
10
15
20FFT N=40
100
25
30
35
40
|X |
50
00.10.20.30.4
0.50.6w (unit :pi)
0.70.80.91
（2）clear clc %xh11
%N=40;%取40个点，时间太短，现象不够精确 N=1000;%取1000个点 n=0:N-1; T=0.01; Fs=1/T; t=T*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t); k=0:N/2;f=Fs/N*k;%w=2*pi/N*k; X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');title('signal x(n)'); subplot(2,1,2);plot(f,magX);%plot(w/pi,magX); title('FFT N=40');
xlabel('w (unit :pi)');ylabel('|X|'); grid
axis([0 5 0 3000])
signal x(n)
[**************]0FFT N=40
[**************]0
30002000
|X |
1000
00.511.52
2.5
w (unit :pi)
33.544.55
2．程序2： %xh12
% 用FFT 实现由DFT 计算线性卷积 x=[2 1 0 1 3];h=[1 3 2 1]; L=length(x)+length(h)-1; XE=fft(x,L);HE=fft(h,L); y1=ifft(XE.*HE);
%画出由圆周卷积计算线性卷积结果及误差 k=0:L-1;
subplot(211);stem(k,real(y1));axis([0 8 0 15]); title('Result of Linear Convolution'); xlabel('Time index k');ylabel('Amplitude');
y2=conv(x,h); error=y1-y2; subplot(212); stem(k,abs(error));
xlabel('Time index k');ylabel('Amplitude'); title('Error Magnitude');
Result of Linear Convolution
A m p l i t u d e
0x 10
-16
123
45Time index kError Magnitude
678
A m p l i t u d e
Time index k
线性卷积直接计算：
function[y,ny]=conv_m(x,nx,h,nh) %改进卷积程序 nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];
y=conv(x,h); %MATLAB自带的函数 在命令窗口输入： x=[2,1,0,1,3];nx=[0:5];
h=[1,3,2,1];nh=[0:3]; [y,ny]=conv_m(x,nx,h,nh)
stem(ny,y,'.');xlabel('时间序号n');title('卷积和y(n)=x(n)*h(n)');
卷积和y(n)=x(n)*h(n)
012
34时间序号n
567
比较结果：这两个有限长序列的圆周卷积等于对两个序列的线性卷积。但当满足一定条件时才相等，若不满足，则这两个有限长序列的圆周卷积等于对两个序列的线性卷积在时间上的混叠。
实验三 IIR数字滤波器的设计
一、实验目的
1．掌握脉冲响应不变法和双线性变换法设计IIR 数字滤波器的原理和方法； 2．观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别。 三、实验原理与方法和手段 1．脉冲响应不变法
所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和ha(t)的采样值，即：h (n ) =h a (t ) t =nT =h a (nt ) ，其中，T 为采样周期。
在脉冲响应不变法中，模拟角频率和数字角频率的变换关系为：ω=ΩT ，
可见，Ω和ω之间的变换关系为线性的。
在MATLAB 中，可用函数impinvar 实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：
[B,A]=impinvar(b,a,fs1) [B,A]=impinvar(b,a)
其中，b 、a 分别为模拟滤波器的分子和分母多项式系数向量；fs1为采样频率（Hz ），缺省值fs=1Hz；B 、A 分别为数字滤波器分子和分母多项式系数向量。 2．双线性变换法：
21-z -1
由于s 平面和z 平面的单值双线性映射关系为s ＝，其中T 为采样
T 1+z -1
周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。
在双线性变换中，模拟角频率和数字角频率的变换关系为：Ω=见，Ω和ω之间的变换关系为非线性的。
在MATLAB 中，可用函数bilinear 实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：[B,A]=bilinear(b,a,fs1) 3．数字滤波器设计
（1）定技术指标转换为模拟滤波器设计性能指标。
（2）估计满足性能指标的模拟相应滤波器性能阶数和截止频率。 利用MATLAB 中buttord 、cheb1ord 、cheb2ord 、ellipord 等函数， 调用格式如：[N , Ωc ]=buttord (Ωp , Ωs , αp , αs , ' s ' )
其中，Ωp 为通带边界频率，rad/s；Ωs 为阻带边界频率，rad/s;αp 为带通波动，dB; αs 为阻带衰减，dB ；‘s ’表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωc 为模拟滤波器的截止频率(-3dB频率) ，rad/s。函数适用低通、高通、带通、带阻滤波器。 （3）设计模拟滤波器。
MATLAB 信号处理工具箱提供了模拟滤波器设计的完全工具函数：butter 、cheby1，cheby2、ellip 、besself 。用户只需一次调用就可完成低通、高通、带通、带阻滤波器设计。 调用格式如：[b,a]=butter(N,Ωc, ’ftype ’, ‘s ’) ，其中，’ftype ’
2ω
tg ，可T 2
为滤波器类型：
‘high ’表示高通滤波器，截止频率为Ωc ；
‘stop ’表示带阻滤波器，Ωc=[Ω1 Ω2] (Ω1
（４）利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。
五、实验步骤 具体步骤如下：
1、查看帮助文件，了解相关函数的调用格式。 2、实验题目
1. 设计模拟巴特沃斯低通滤波器，fp=300Hz，αp=1dB，fs=800Hz，αs=20dB。 %xh13模拟低通滤波器技术指标 αp=1; αs=20; fp=300; fs=800; Ωap=2*pi*fp; Ωas=2*pi*fs;
%设计模拟巴特沃斯低通滤波器
[N,Ωac]=buttord(Ωap, Ωas, αp, αs,'s')
[b,a]=butter(N,Ωac,'s') %设计模拟巴特沃斯低通滤波器，Ωap 为通带边界频率，rad/s；Ωas 为阻带边界频率，rad/s;αp 为通带最大衰减，dB ；αs 为阻带最小衰减，dB ；'s' 表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωac 为模拟滤波器的截止频率(-3dB频率) ，rad/s；b 、a 分别为模拟滤波器的系统函数分子和分母多项式系数向量； [H,Ω]=freqs(b,a); %求模拟滤波器的频率响应 %绘制频响幅度谱
plot(Ω/2/pi,20*log10(abs(H))); %横轴为频率，单位：HZ; 纵轴频响幅度，
单位：dB
axis([0,1500,-50,0]);
xlabel('频率Hz');ylabel(' H幅值dB');
H 幅值d B
[1**********]00
2：设计模拟巴特沃斯高通滤波器，fp=800Hz，αp=1dB，fs=300Hz，αs=20dB。 %xh14模拟高通滤波器技术指标 αp=1; αs=20; fp=800; fs=300; Ωp=2*pi*fp; Ωs=2*pi*fs;
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s')
[b,a]=butter(N,Ωc,'high','s') %设计模拟巴特沃斯高通滤波器，Ωp 为通带边界频率，rad/s；Ωs 为阻带边界频率，rad/s；αp 为通带最大衰减，dB ；αs 为阻带最小衰减，dB ；'s' 表示为模拟滤波器；函数返回值N 为模拟滤波器的最小阶数；Ωc 为模拟滤波器的截止频率(-3dB频率) ，rad/s；b 、a 分别为模拟高通滤波器的系统函数分子和分母多项式系数向量； [H,Ω]=freqs(b,a);
plot(Ω/2/pi,20*log10(abs(H))); axis([0,1000,-50,0]);
10001200
频率Hz
[**************]0
xlabel('频率Hz');ylabel(' H幅值dB');
-10
-20
H 幅值d B
-30
-40
-50
-60
-70
[1**********]00
500600频率Hz
[1**********]00
3：fp=0.１kHZ, αp=1dB，fs=0.3kHZ，αs=25dB，T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth 数字低通滤波器。 %xh15采用冲击响应不变法 αp=1; αs=25; fp=100; fs=300; Ωp=2*pi*fp;
Ωs=2*pi*fs; %数字滤波器技术指标要求转化成模拟滤波器技术指标要求 fs1=1000; %采样频率 %设计模拟滤波器
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s'); [b,a]=butter(N,Ωc,'s');
[B,A]=impinvar(b,a,fs1) %用冲击响应不变法将模拟滤波器变换成数字滤波器。B 、A 分别为数字滤波器的系统函数分子和分母多项式系数向量； [H1,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应
%绘制数字滤波器频响幅度谱
subplot(211);plot(w*fs1/2/pi,20*log10(abs(H1))); axis([0,1000,-100,0]);
xlabel('频率Hz');ylabel(' H1幅值dB'); %采用双线性变换法
%数字滤波器的技术指标要求 αp=1; αs=25; fp=100; fs=300;
fs1=1000; %采样频率
wp=2*pi*fp/fs1; %数字滤波器通带边界频率 ws=2*pi*fs/fs1; %数字滤波器阻带边界频率 %变换为同类型模拟滤波器的技术指标要求
Ωp=2*fs1*tan(wp/2); %同类模拟滤波器通带边界频率 Ωs=2*fs1*tan(ws/2); %同类模拟滤波器阻带边界频率 %设计同类型模拟滤波器
[N,Ωc]=buttord(Ωp, Ωs, αp, αs,'s'); [b,a]=butter(N,Ωc,'s');
%用双线性变换法将模拟滤波器变换成数字滤波器 [B,A]=bilinear(b,a,fs1)
[H2,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应 %绘制数字滤波器频响幅度谱
subplot(212);plot(w*fs1/2/pi,20*log10(abs(H2))); axis([0,1000,-100,0]);
xlabel('频率Hz');ylabel('H2幅值dB');
H 1幅值d B
-50
-100
[1**********]00
500600频率Hz
[1**********]00
H 2幅值d B
-50
-100
100
200
300
400500频率Hz
600
700
800
900
双线性变换法和脉冲不变法的特点和区别：
脉冲响应不变法：1，模拟频率到数字频率的转换时线性的。
2，数字滤波器单位脉冲响应的数字表示近似原型的模拟滤波器单位脉冲响应，因此时域特性逼近好，但会产生频谱混叠现象，只适合带限滤波器
双线性变换法：克服多值映射得关系，可以消除频率的混叠，但是因为是非线性的，所以在高频处有较大的失真。
实验四 FIR数字滤波器的设计
一、实验目的
1．掌握用窗函数法设计FIR 数字滤波器的原理和方法； 2．熟悉线性相位FIR 滤波器的幅频特性和相频特性； 3．了解不同窗函数对滤波器性能的影响。 三、实验原理与方法和手段
窗函数设计法是一种把一个长序列变成有限长的短序列的设计方法，是在时域进行的。用窗函数法设计FIR 数字滤波器时，先根据W c和N 求出相应的的理想滤波器的单位脉冲响应hd(n)。
1
hd (n ) =
2π
⎰
w c
-w c
H d (e ) e
jw jwn
sin(w c (n -a )) dw =
π(n -a )
因为hd(n)一般是非因果的，且无限长，物理上是不可实现的。为此可选择适当的窗函数w(n)截取有限长的hd(n)， 即h(n)= hd(n)w(n)，只要阶数足够长，截取的方法合理，总能够满足频域的要求。
实际中常用的窗函数有矩形（Boxcar ）窗、三角（Bartlett ）窗、汉宁(Hanning)窗、汉明(Hamming) 窗和布莱克曼（Blackman ）窗。这些窗函数各有优缺点, 所以要根据实际情况合理选择窗函数类型。
1．窗函数法设计线性相位FIR 滤波器的一般步骤为：
jw
H (e ) 的特性； （1）确定理想滤波器d
（2）由Hd （e jw ）求出h d (n ) ；
（3）根据过渡带宽度和阻带最小衰减，借助窗函数确定窗的形式及N 的大小，即选择适当的窗函数，并根据线性相位条件确定窗函数的长度N ；在MATLAB 中，可由w=boxcar(N)（矩形窗）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗）等函数来实现窗函数设计法中所需的窗函数。
（4）由h(n)= h d (n ) .w(n), 0≤n ≤ N-1，得出单位脉冲响应h(n); （5）对h(n)作离散时间傅立叶变换，得到H(e ) 。
2．在MATLAB 中，可以用b=fir1(N,Wn,’ftype’,taper) 等函数辅助设计FIR 数字滤波器。N 代表滤波器阶数；Wn 代表滤波器的截止频率(归一化频率) ，当设计带通和带阻滤波器时，Wn 为双元素相量；ftype 代表滤波器类型，如’high ’高通，’stop ’带阻等；taper 为窗函数，默认为海明窗，窗函数实现需要用窗函数blackman, hamming,hanning chebwin, kaiser产生。 五、实验步骤
在“开始--程序”菜单中，找到MATLAB 程序，运行启动；
jw
进入MATLAB 后 ，在Command Window中输入自己编写的主程序，并执行； 记录运行结果图形，作分析对比。
具体步骤如下：
1. 用窗函数法设计一线性相位FIR 低通滤波器，要求通带截止频率wc=冲击响应h （n ）的长度N ＝21。绘出h(n)及其幅频响应曲线.
sin((n -a ))
3 设计分析:理想滤波器单位冲击响应为hd (n ) =；为了满足线性相π(n -a ) N -1
位FIR 滤波器条件h(n)= h(N-1-n)，要求a==10。
2
π
，单位3
π
%xh19
N=21;wc=pi/3; %理想低通滤波器参数 n=0:N-1; a=(N-1)/2;
hdn=sin(wc*(n-a))/pi./(n-a); %计算理想低通滤波器单位冲击响应hd(n) if rem(N,2)~=0 hdn(a+1)=wc/pi; end %N为奇数时, 处理n=a点的0/0型 wn=hamming(N); %hamming窗 hn=hdn.*wn';
subplot(121); %绘出h(n)及幅频特性曲线 stem(n,hn,'.');
xlabel('n');ylabel('h(n)'); title('hamming窗设计的h(n)'); grid on; subplot(122); hw=fft(hn,512); w=2*[0:511]/512;
plot(w,20*log10(abs(hw)));
xlabel('w/pi');ylabel('Magnitude(dB)'); title('hamming窗设计的幅频特性'); grid on;
21
hamming 窗设
计的h(n)
20
hamming 窗设计的幅频特性
-20
M a g n i t u d e (d B )
5
10n
15
20
-40
h (n )
-60
-80
-100
-120
00.5
2．针对一个含有5Hz 、15Hz 和30Hz 的混和正弦波信号，设计一个FIR 带通滤波器。
参数要求:采样频率fs=100Hz，通带下限截止频率fc1=10Hz，通带上限截止频率fc2=20Hz，过渡带宽6Hz ，通阻带波动0.01，采用凯塞窗设计。 %xh20
fc1=10; fc2=20; fs=100;
[n,Wn,beta,ftype]=kaiserord([7 13 17 23],[0 1 0], [0.01 0.01 0.01],100) window=kaiser(n+1,beta); %使用kaiser 窗函数
b=fir1(n, Wn,window); %使用标准频率响应的加窗设计函数fir1 freqz(b,1,512); %数字滤波器频率响应 t=(0:100)/fs;
s=sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30); sf=filter(b,1,s); %对信号s 进行滤波 figure
subplot(2,1,1); plot(t,s) subplot(2,1,2); plot(t,sf)
22
1w/pi
1.52
50
M a g n i t u d e (d B )
-50-100-150
0.1
0.2
0.30.40.50.60.70.8Normalized Frequency (⨯π rad/sample)
0.9
1
500
P h a s e (d e g r e e s )
-500
-1000
00.10.2
0.30.40.50.60.70.8Normalized Frequency (⨯π rad/sample)
0.91
42
0-2-4
00.10.20.30.40.50.60.70.80.91
21
0-1-2
00.10.20.30.40.50.60.70.80.91
23