3．单位抽样序列
可以通过借助MATLAB中的零矩阵函数zeros表示。
全零矩阵zeros（1,N）产生一个由N个零组成的列向量，对于有限区间的
可以通过以下MATLAB程序表示
%单位抽样序列实现程序
k=-30:
30;
delta=[zeros（1,30）,1,zeros（1,30）];
stem（k,delta）
4．单位阶跃序列
可以通过借助MATLAB中的单位矩阵函数ones表示。
单位矩阵ones（1,N）产生一个由N个1组成的列向量，对于有限区间的
可以通过以下MATLAB程序表示
%单位阶跃序列实现程序
k=-30:
30;
uk=[zeros（1,30）,ones（1,31）];
stem（k,uk）
程序运行结果如图1.24所示。
例1.7编写程序来产生下列基本脉冲序列。
（1）单位抽样序列，起点ns=0，终点nf=10，在n0=3处有一单位脉冲。
（2）单位阶跃序列，起点ns=0，终点nf=10，在n0=3前为0，在ns=3后为1。
（3）复指数序列，σ=-0.2，ω0=0.5。
解：
程序清单如下。
clear,n0=3;ns=0;nf=10;
n1=[ns:
nf];x1=[zeros（1,n0-ns）,1,zeros（1,nf-n0）];%单位抽样序列的产生
n2=[ns:
nf];x2=[zeros（1,n0-ns）,ones（1,nf-n0+1）];%单位阶跃序列的产生
n3=[ns:
nf];x3=exp（（-0.2+0.5j）*n3）;%复指数序列的产生
subplot（2,2,1）,stem（n1,x1）;title（'单位脉冲序列δ（n-3）'）;%画图表示单位抽样序列
subplot（2,2,2）,stem（n2,x2）;title（'单位阶跃序列u（n-3）'）;%画图表示单位阶跃序列
subplot（2,2,3）,stem（n3,real（x3））;line（[0,10],[0,0]）
title（'复指数序列'）;ylabel（'实部'）;%画图表示复指数序列的实部
subplot（2,2,4）,stem（n3,imag（x3））;line（[0,10],[0,0]）
title（'复指数序列'）;ylabel（'虚部'）;
实验二利用MATLAB实现频谱分析
例1已知有限长序列
的长度
，且
，用FFT求
，再用IFFT求
。
解：
利用快速傅里叶变换函数求解的MATLAB实现程序清单如下
clear
xn=[1,2,-1,3];
X=fft（xn）
x=ifft（X）
程序运行结果如下
X=
5.00002.0000+1.0000i-5.00002.0000-1.0000i
x=
12-13
例2设
是由两个正弦信号及白噪声的叠加，试用FFT文件对其作频谱分析。
图2.1运行结果
解：
程序清单如下
图4.14
应用实例
%产生两个正弦加白噪声；
N=256;
f1=.1;f2=.2;fs=1;
a1=5;a2=3;
w=2*pi/fs;
x=a1*sin（w*f1*（0:
N-1））+a2*sin（w*f2*（0:
N-1））+randn（1,N）;
%应用FFT求频谱；
subplot（2,2,1）;
plot（x（1:
N/4））;
title（'原始信号'）;
f=-0.5:
1/N:
0.5-1/N;
X=fft（x）;
y=ifft（X）;
subplot（2,2,2）;
plot（f,fftshift（abs（X）））;
title（'频域信号'）;
subplot（2,2,3）;
plot（real（x（1:
N/4）））;
title（'时域信号'）;
运行结果如图4.14所示，该程序同时完成了傅里叶变换与傅里叶反变换。
例3设
为长度
的矩形序列，用MATLAB程序分析FFT取不同长度时
的频谱变化。
解：
时
的FFTMATLAB实现程序如下
x=[1,1,1,1,1,1];
N=8;
y1=fft（x,N）;
n=0:
N-1;
subplot（3,1,1）;stem（n,abs（y1）,'.k'）;axis（[0,9,0,6]）;
N=32;
y2=fft（x,N）;
n=0:
N-1;
subplot（3,1,2）;stem（n,abs（y2）,'.k'）;axis（[0,40,0,6]）;
N=64;
y3=fft（x,N）;
n=0:
N-1;
subplot（3,1,3）;stem（n,abs（y3）,'.k'）;axis（[0,80,0,6]）;
运行结果
的频谱如图2.2所示
图2.2不同长度时
的频谱图
例3
，求其傅立叶变换
。
解：
严格说，在MATLAB中不使用symbolic工具箱是不能分析模拟信号的，但是当采样时间间隔充分小的时候，可产生平滑的图形，当时间足够长，可显示出所有的模型，也就是可以近似的分析。
此例中，
为fh=2000Hz的带限信号，因此取
程序清单如下。
Dt=0.00005;t=-0.005:
Dt:
0.005;%模拟信号
xa=exp（-1000*abs（t））;
Wmax=2*pi*2000;K=500;k=0:
1:
K;W=k*Wmax/K;%连续时间傅立叶变换
Xa=xa*exp（-j*t'*W）*Dt;Xa=real（Xa）;
W=[-fliplr（W）,W（2:
501）];
Xa=[fliplr（Xa）,Xa（2:
501）];
figure
（1）
subplot（2,1,1）;plot（t*1000,xa,'.'）;
xlabel（'tinmsec'）;ylabel（'xa（t）'）;gtext（'模拟信号'）;
subplot（2,1,2）;plot（W/（2*pi*1000）,Xa*1000,'.'）;
xlabel（'FrequenceinKHz'）;ylabel（'Xa（jw）*1000'）;gtext（'连续时间傅立叶变换'）;
程序运行结果如图2.3所示。
图2.3模拟信号及傅立叶变换曲线
例4以上例中的
说明采样频率对频域特性的影响，分别取采样频率为fs=5000Hz和fs=1000Hz，绘出
曲线。
解：
程序清单如下。
Dt=0.00005;t=-0.005:
Dt:
0.005;%模拟信号
xa=exp（-1000*abs（t））;
Ts=0.0002;n=-25:
1:
25;%离散时间信号
x=exp（-1000*abs（n*Ts））;
K=500;k=0:
1:
K;w=pi*k/K;%离散时间傅立叶变换
X=x*exp（-j*n'*w）;X=real（X）;
w=[-fliplr（w）,w（2:
501）];
X=[fliplr（X）,X（2:
501）];
figure
subplot（2,1,1）;plot（t*1000,xa,'.'）;ylabel（'xa1（t）'）;
gtext（'离散信号'）;holdon
stem（n*Ts*1000,x）;holdoff
subplot（2,1,2）;plot（w/pi,X,'.'）;ylabel（'X1（jw）'）;
gtext（'离散时间傅立叶变换'）;gtext（'Ts=0.2msec'）
程序运行结果如图2.4所示。
图2.4fs=5000Hz时的离散信号及傅立叶变换曲线
对上面的MATLAB程序稍加修改，就可以得到fs=1000Hz时的
，程序运行结果如图2.5所示。
可以看出，当fs=5000Hz时，
的曲线和例1.14中的
曲线完全一致，没有发生混叠现象，而当fs=1000Hz时，
的曲线与例1.14中的
曲线不同，产生了混叠的现象。
图2.5fs=1000Hz时的离散信号及傅立叶变换曲线
例5对上列例产生的离散序列x1（n）和x2（n），采用sinc函数进行内插重构。
解：
程序清单如下
Ts1=0.0002;Fs1=1/Ts1;n1=-25:
1:
25;nTs1=n1*Ts1;%离散时间信号
x1=exp（-1000*abs（nTs1））;
Ts2=0.001;Fs2=1/Ts2;n2=-5:
1:
5;nTs2=n2*Ts2;
x2=exp（-1000*abs（nTs2））;
Dt=0.00005;t=-0.005:
Dt:
0.005;%模拟信号重构
xa1=x1*sinc（Fs1*（ones（length（nTs1）,1）*t-nTs1'*ones（1,length（t））））;
xa2=x2*sinc（Fs2*（ones（length（nTs2）,1）*t-nTs2'*ones（1,length（t））））;
subplot（2,1,1）;plot（t*1000,xa1,'.'）;ylabel（'xa1（t）'）;
title（'从x1（n）重构模拟信号xa（t）'）;holdon
stem（n1*Ts1*1000,x1）;holdoff
subplot（2,1,2）;plot（t*1000,xa2,'.'）;ylabel（'xa2（t）'）;
title（'从x2（n）重构模拟信号xa（t）'）;holdon
stem（n2*Ts2*1000,x2）;holdoff
程序运行结果如图1.37所示。
将图2.6中的重构后的模拟信号曲线与图1.34中的原始模拟信号曲线进行比较，可以看出用离散信号x1（n）重构出的模拟信号与原始信号误差很小，而用离散信号x2（n）重构出的模拟信号误差较大，这是因为x2（n）的采样频率为fs=1000Hz，不满足采样定理fs>2fh，因此采样信号的频谱发生混叠造成取样信号不能无失真的恢复出原始信号。
图2.6从x1（n）和x2（n）分别重构模拟信号
实验三IIR数字滤波器的设计MATLAB实现
目前，设计IIR数字滤波器的通用方法是先设计相应的低通滤波器，然后再通过双线性变换法和频率变换得到所需要的数字滤波器。
模拟滤波器从功能上分有低通、高通、带通及带阻四种，从类型上分有巴特沃兹（Butterworth）滤波器、切比雪夫（Chebyshev）I型滤波器、切比雪夫II型滤波器、椭圆（Elliptic）滤波器以及贝塞尔（Bessel）滤波器等。
1应用MATLAB设计IIR数字滤波器相关文件
下面给出与IIR数字滤波器设计有关的MATLAB文件。
1．buttord.m
用来确定数字低通或模拟低通滤波器的阶次，其调用格式分别是
（1）[N,Wn]=buttord（Wp,Ws,Rp,Rs）
（2）[N,Wn]=buttord（Wp,Ws,Rp,Rs,‘s’）
格式
（1）对应数字滤波器，式中Wp,Ws分别上通带和阻带的截止频率，实际上它们是归一化频率，其值在0~1之间，1对应抽样频率的一半。
对低通和高通滤波器，Wp,Ws都是标量，对带通和带阻滤波器，Wp,Ws都是
的向量。
Rp,Rs分别是通带和阻带的衰减，单位为dB。
N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。
格式
（2）对应模拟滤波器，式中各个变量的含义与格式
（1）相同，但Wp,Ws及Wn的单位为rad/s，因此，它们实际上是频率。
2．buttap.m
用来设计模拟低通原型滤波器
，其调用格式是
[z,p,k]=buttap（N）
N是欲设计的低通原型滤波器的阶次，z、p和k分别是设计出的
的极点、零点及增益。
3．lp2lp.m
4．lp2hp.m
5．lp2bp.m
6．lp2bs.m
以上4个文件用来将模拟低通原型滤波器
分别转换为低通、高通、带通、及带阻滤波器。
其调用格式分别是
（1）[B,A]=lp2lp（b,a,Wo）或[B,A]=lp2hp（b,a,Wo）
（2）[B,A]=lp2bp（b,a,Wo,Bw）或[B,A]=lp2bs（b,a,Wo,Bw）
式中b，a分别是模拟低通原型滤波器
有分子、分母多项式的系数向量，B，A分别是转换后的
有分子、分母多项式的系数向量；在格式
（1）中，Wo是低通或高通滤波器的截止频率；在格式
（2）中Wo是带通或带阻滤波器的中心频率，Bw是其带宽。
7．bilinear.m
实现双线性变换，即由模拟滤波器
得到数字滤波器
。
其调用格式是
[Bz,Az]=bilinear（B,A,Fs）
式中B、A分别是
的分子、分母多项式的系数向量；Bz、Az分别是
的分子、分母多项式的系数向量，Fs是抽样频率。
8．butter.m
用来直接设计巴特沃兹数字滤波器，实际上它把buttord.m，buttap.m，lp2lp.m及bilinear.m等文件都包含进去，从而使设计过程更简捷。
其调用格式是
（1）[B,A]=butter（N,Wn）
（2）[B,A]=butter（N,Wn,’high’）
（3）[B,A]=butter（N,Wn,’stop’）
（4）[B,A]=butter（N,Wn,’s’）
格式
（1）~（3）用来设计数字滤波器，B、A分别是
的分子、分母多项式的系数向量，Wn是通带截止频率，范围在0~1之间，1对应抽样频率的一半。
若Wn是标量，则格式
（1）用来设计低通数字滤波器。
若Wn是
的向量，则格式
（1）用来设计数字带通滤波器；格式
（2）用来设计数字高通滤波器；格式（3）用来设计数字带阻滤波器。
格式（4）用来设计模拟滤波器。
9．cheb1ord.m
求切比雪夫I型滤波器的阶次。
10．cheb1ap.m
用来设计原型切比雪夫I型模拟滤波器。
11．cheby1.m
直接设计切比雪夫I型滤波器。
以上3个文件的调用格式与对应的巴特沃兹滤波器的文件类似。
12．cheb2ord.m
13．ellipord.m
14．cheb2ap.m
15．ellipap.m
16．besselap.m
17．cheby2.m
18．ellip.m
19．besself.m
以上分别为切比雪夫II型滤波器、椭圆及贝塞尔滤波器设计函数。
其格式与切比雪夫I型滤波器和巴特沃兹滤波器类同。
此外，与本章内容有关的MATLAB文件还有：
20．impinvar.m
用冲激响应不变法实现
到
及
到z的转换。
21．maxflat.m
设计广义巴特沃兹低通滤波器。
22．yulewalk.m
利用最小平方方法设计Yule-Walker滤波器。
例1利用巴特沃斯模拟滤波器设计数字高通滤波器，要求通带截止频率为0.6π，通带内衰减不大于1dB，阻带起始频率为0.4π，阻带内衰减不小于15dB,，采样周期Ts=1s。
%usingbutterworthtodesignhighpassDF
clearall;
Wp=0.6*pi;
Ws=0.4*pi;
Ap=1;
As=15;
[N,wn]=buttord（Wp/pi,Ws/pi,Ap,As）;%计算巴特沃斯滤波器阶次和截止频率
[b,a]=butter（N,wn,'high'）;%频率变换法设计巴特沃斯高通滤波器、
[C,B,A]=sdir2cas（b,a）%把滤波器的系统函数转换为成级联型
[db,mag,pha,w]=freqz_m2（b,a）;%数字滤波器的频率响应
subplot（211）;
plot（w/pi,mag）;grid;
axis（[0,1,0,1.2]）;
title（'数字滤波器的幅频响应'）;
subplot（212）;
plot（w/pi,db）;grid;
axis（[0,1,10,-200]）;
title（'数字滤波器的幅频响应dB'）;
function[C,B,A]=sdir2cas（b,a）;
%把滤波器的系统函数转换为成级联型
%a---数字滤波器系统函数分母系数
%b---数字滤波器系统函数分子系数
%C---数字滤波器阶联结构的常数因子
%B---数字滤波器阶联结构的分母系数
%A---数字滤波器阶联结构的分子系数
Na=length（a）-1;Nb=length（b）-1;
b0=b
（1）;b=b/b0;
a0=a
（1）;a=a/a0;
C=b0/a0;
p=cplxpair（roots（a））;K=floor（Na/2）;
ifK*2==Na
A=zeros（K,3）;
forn=1:
2:
Na
Arow=p（n:
1:
n+1,:
）;Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
elseifNa==1
A=[0real（poly（p））];
else
A=zeros（K+1,3）;
forn=1:
2:
2*K
Arow=p（n:
1:
n+1,:
）;Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
A（K+1,:
）=[0real（poly（p（Na）））];
end
z=cplxpair（roots（b））;K=floor（Nb/2）;
ifNb==0
B=[00poly（z）];
elseifK*2==Nb
B=zeros（K,3）;
forn=1:
2:
Nb
Brow=z（n:
1:
n+1,:
）;Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
elseifNb==1
B=[0real（poly（z））];
else
B=zeros（K+1,3）;
forn=1:
2:
2*K
Brow=z（n:
1:
n+1,:
）;Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
B（K+1,:
）=[0real（poly（z（Nb）））];
end
例
升级会员 