elsenx2=0;len_x1=0;%若x1和x2同长，不补0
end%长者为补0长度的基础
lt=lmax;%先将x2补得与x1同长，再将两边补最大长度的0
u=[zeros（1,lt）,x2,zeros（1,nx2）,zeros（1,lt）];
t1=（-lt+1:
2*lt）;%先将x1补得与x2同长，再将左边补2倍最大长度的0
x1=[zeros（1,2*lt）,x1,zeros（1,nx1）];%将x1左右反折
hx1=fliplr（x1）;
N=length（hx1）;
y=zeros（1,3*lt）;%将y存储单元初始化
fork=0:
2*lt%动态演示绘图
p=[zeros（1,k）,hx1（1:
N-k）];%使hx1向右循环移位
y1=u.*p;%使输入和翻转移位的脉冲过渡函数逐项相乘
yk=sum（y1）;%相加
y（k+lt+1）=yk;%将结果放入数组y
subplot（4,1,1）;stem（t1,u）;
subplot（4,1,2）;stem（t1,p）;
subplot（4,1,3）;stem（t1,y1）;
subplot（4,1,4）;stem（k,yk）;%通过图形表示每一次卷积的结果
axis（[-10,50,0,5]）;holdon%在图形窗口上保留每一次运行的图形结果
pause（0.5）;%停滞0.5秒
end
第2章
例2-求以下序列的
变换：
，
，
，
程序如下：
symsaznw%把字符aznw定义为基本的符号对象
x1=a^n;X1=ztrans（x1）%求x1的z变换
x2=n-1;X2=ztrans（x2）
x3=（n*（n-1））/3;X3=ztrans（x3）
x4=exp（j*w*n）;X4=ztrans（x4）
运行结果如下：
X1=z/a/（z/a-1）
X2=z/（z-1）^2-z/（z-1）
X3=1/3*z*（z+1）/（z-1）^3-1/3*z/（z-1）^2
X4=z/exp（i*w）/（z/exp（i*w）-1）
例2-求下列函数的
的反变换：
，
，
程序如下：
symsnza%把anz定义为基本的符号对象
X1=z/（z-1）;x1=iztrans（X1）%求X1的z的反变换
X2=z/（1-z）^3;x2=iztrans（X2）
X3=z/（z-a）^2;x3=iztrans（X3）
运行结果如下：
x1=1
x2=1/2*n-1/2*n^2
x3=a^n*n/a
例2-已知
，试用部分分式法求
的反变换，并画
。
程序如下：
b=[1,0,0];
a=[1,-3,2];
[rpc]=residuez（b,a）;%把b（z）/a（z）展开成部分分式
N=15;
n=0:
N-1;
x=r
（1）*p
（1）.^n+r
（2）*p
（2）.^n;
stem（n,x）;
title（'x（n）'）;
在命令窗口中将显示：
r=
2
-1
p=
2
1
c=
0
所以多项式分解后表示为：
运行结果如下：
图2-部分分式法求
图
例2-研究
右半平面的实数极点对系统响应的影响。
已知系统的零-极点模型分别为：
，
，
求这些系统的零极点分布图以及系统的冲激响应，判断系统的稳定性。
程序如下：
z1=[0]';p1=[0.5]';k=1;%输入零极点增益模型零点向量z、极点向量p和增益系数k
[b1,a1]=zp2tf（z1,p1,k）;%求系数传递函数（tf）模型中分子、分母多项式的系数向量
subplot（3,2,1）,zplane（z1,p1）;%由零点和极点向量绘制零极点分布图
ylabel（'极点在单位圆内'）;
subplot（3,2,2）,impz（b1,a1,15）;%绘制数字系统的冲激响应，取样点数为15
z2=[0]';p2=[1]';
[b2,a2]=zp2tf（z2,p2,k）;
subplot（3,2,3）,zplane（z2,p2）;
ylabel（'极点在单位圆上'）;
subplot（3,2,4）,impz（b2,a2,15）;
z3=[0]';p3=[2]';
[b3,a3]=zp2tf（z3,p3,k）;
subplot（3,2,5）,zplane（z3,p3）;
ylabel（'极点在单位圆外'）;
subplot（3,2,6）,impz（b3,a3,15）;
运行结果如下：
图2-右半平面极点对系统响应影响图
由MATLAB结果图可见，这3个系统的极点均为实数且处于z平面的右半平面。
当极点处于单位圆内，系统的冲激响应曲线随着频率的增大而衰减；当极点处于单位圆上时，系统的冲击响应曲线为等幅振荡；当极点处于单位圆外时，系统的冲击响应曲线随着频率的增大而发散。
第3章
例3-已知两个周期序列分别为
，
，动态演示它们的周期卷积和
。
程序如下：
clear
n=0:
5;%建立时间向量序列n
xn1=[0,2,1,1,0,0];%建立xn1的序列主值
xn2=[0,1,1,0,0,0];%建立xn2的序列主值
N=length（xn1）;
nx=（-N:
3*N-1）;
hxn2=xn2（mod（nx,N）+1）;%将xn2序列周期严拓
u=[zeros（1,N）,xn2,zeros（1,2*N）];%按xn2周期严拓后的长度重建主值信号
xn12=fliplr（xn1）;%将xn1做左右反折
hxn1=xn12（mod（nx,N）+1）;%将xn1反折后的序列周期严拓
N1=length（hxn1）;
y=zeros（1,4*N）;%将y存储单元初始化
fork=0:
N-1%动态演示开始
p=[zeros（1,k+1）,hxn1（1:
N1-k-1）];%使hxn1向右循环移位
y1=u.*p;%使输入和翻转移位的脉冲过渡函数相乘
yk=sum（y1）;%相乘之后再相加
y（[k+1,k+N+1,k+2*N+1,k+3*N+1]）=yk;%将结果放入数组y
subplot（4,1,1）;stem（nx,hxn2）;
axis（[-1,3*N,0,1.1]）;ylabel（'x2（n）'）;
subplot（4,1,2）;stem（nx,p）;
axis（[-1,3*N,0,3.3]）;ylabel（'x1（n）'）;
subplot（4,1,3）;stem（k,yk）;%作图表示主值区每一次卷积的结果
axis（[-1,3*N,0,6.6]）;holdon%在图形窗口上保留每一次运行的图形结果
ylabel（'主值区'）;
subplot（4,1,4）;stem（nx,y）;
axis（[-1,3*N,0,6.6]）;ylabel（'卷积结果'）;
pause（3）;%停滞3秒
end
例3-求矩形序列的DTFT。
程序如下：
N=8;
xn=ones（1,N）;
n=0:
N-1;
w=linspace（0,2*pi,200）;%将[0,2π]频率区间分割成200份
X=xn*exp（-j*n'*w）;%离散时间的傅里叶变换
subplot（3,1,1）,stem（n,xn,'k'）;
ylabel（'x（n）'）;
subplot（3,1,2）,plot（w,abs（X）,'k'）;%显示幅度谱
axis（[-2*pi,2*pi,1.2*min（abs（X））,1.2*max（abs（X））]）;
ylabel（'幅度谱'）;
subplot（3,1,3）,plot（w,angle（X）,'k'）;%显示相位谱
axis（[-2*pi,2*pi,1.2*min（angle（X））,1.2*max（angle（X））]）;
ylabel（'相位谱'）;
运行结果如下：
图3-矩形序列的DTFT
例3-求矩形序列的DFT，并与上例进行比较。
程序如下：
N=8;
xn=ones（1,N）;
n=0:
N-1;
k=0:
N-1;
Xk=xn*exp（-j*2*pi/N）.^（n'*k）;%离散傅里叶变换
L=50;
yn=[xn,zeros（1,L-8）];%将x（n）有限长序列补足至50
m=0:
L-1;
t=0:
L-1;
Yk=yn*exp（-j*2*pi/L）.^（m'*t）;%离散傅里叶变换
subplot（2,2,1）,stem（n,xn）;%绘制原信号序列
title（'x（n）'）;
subplot（2,2,2）,stem（m,yn）;
title（'y（n）'）;
subplot（2,2,3）,stem（k,abs（Xk））;%傅里叶变换所对应的图形
title（'|X（k）|'）;
subplot（2,2,4）,stem（t,abs（Yk））;
title（'|Y（k）|'）;
运行结果如下：
图3-矩形序列的DFT
序列的DFT是DTFT周期内的等间隔采样，上图首先给出8点矩形序列的DTFT，然后给出通过补零加长为50点后的DTFT。
通过对比发现补零后，频谱的采样间隔变密，保留了更多的原信号的频谱特征。
第4章
例4-：
已知一个8点的时域非周期离散
信号和阶跃信号
，
，用
点进行FFT变换，作其时域信号图及信号频谱图。
程序如下：
n0=2;ns=0;ne=7;N=32;
n=ns:
ne;
W1=[（n-n0）==0];%建立时间信号
W2=[（n-n0）>=0];
i=0:
N-1;%频率样点自0开始
y1=fft（W1,N）;
Aw1=abs（y1）;%求幅度谱
y2=fft（W2,N）;
Aw2=abs（y2）;
subplot（2,2,1）;stem（n,W1）;
subplot（2,2,2）;stem（n,W2）;
subplot（2,2,3）;plot（i,Aw1）;
subplot（2,2,4）;plot（i,Aw2）;
运行结果如下：
图4-单位脉冲序列和单位阶跃信号及其频谱
例4-用FFT计算下列连续时间信号的频谱，并讨论
和
值对频谱特性的影响。
程序如下：
T0=[1,0.5,0.25,0.25];%输入不同的Ts值
N0=[256,256,256,1024];%输入不同的N值
form=1:
4;
Ts=T0（m）;N=N0（m）;%赋Ts和N值
n=0:
N-1;
D=2*pi/（Ts*N）;%计算模拟频率分辨率
xa=sin（n*Ts）+sin（1.1*n*Ts）+sin（1.18*n*Ts）;
k=floor（-（N-1）/2:
（N-1）/2）;
Xa=Ts*fftshift（fft（xa,N））;
[m,Xa
（1）]
subplot（2,2,m）;plot（k*D,abs（Xa）,'k'）;
axis（[0,3,1.1*min（abs（Xa））,1.1*max（abs（Xa））]）;
end
图4-FFT时不同
和
值对频谱特性的影响
和
，从图中可以看出，相同的
下，
越大，所取的时域信号长度
越长，则频率分辨率越高；相同的
下，
越大，则在频域抽样点数增加，则频率分辨率越高。
第5章
例5-已知系统的传递函数为
，将其从直接型转换为级联型和并联型。
从直接型转换为级联型，就是将系统传递函数模型（tf）转换为二次分式（sos）模型；从直接型转换为并联型，就是将系统的传递函数模型转换为极点留数（rpk）模型。
程序如下：
b=[4,-2,7,-2];%输入系统函数b参数
a=[1,-1.5,0.5,-0.75];%输入系统函数a参数
[sos,g]=tf2sos（b,a）%由直接型转换为级联型
[r,p,k]=residuez（b,a）%由直接型转换为并联型
运行结果如下：
sos=
1.0000-0.295901.0000-1.50000
1.0000-0.20411.68961.0000-0.00000.5000
g=
4
r=
4.2424
-1.4545+1.6713i
-1.4545-1.6713i
p=
1.5000
0.0000+0.7071i
0.0000-0.7071i
k=
2.6667
由sos和g的数据，可以写出级联型的表达式
由r、p、k的数据，可以列写出并联型的表达式
例5-已知一个FIR系统的传递函数为
，将其从级联型转换为横截型。
程序如下：
sos=[1,-2,3,1,0,0;1,-4,5,1,0,0];%输入sos参数
g=4;%输入g参数
[b,a]=sos2tf（sos,g）%由级联型转换为直接型
运行结果如下：
b=
4-2464-8860
a=
1000
由a和b的数据，可以写出级联型的表达式
第6章
例6-采用MATLAB直接法设计一个巴特沃斯型数字带通滤波器，要求：
;
。
描绘滤波器归一化的绝对和相对幅频特性、相频特性、零极点分布图，列出系统传递函数
程序如下：
ws1=0.2;ws2=0.8;%数字滤波器的阻带截止频率
ws=[ws1,ws2];
wp1=0.4;wp2=0.6;%数字滤波器的通带截止频率
wp=[wp1,wp2];
Rp=1;As=10;%输入滤波器的通阻带衰减指标
[n,wc]=buttord（wp,ws,Rp,As）%计算阶数n和截止频率
[b,a]=butter（n,Rp,wc）%直接求数字带通滤波器的系数
[H,w]=freqz（b,a）;%求数字系统的频率特性
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;%化为分贝值
subplot（2,2,1）,plot（w/pi,abs（H））;
title（'幅度相应'）;
subplot（2,2,2）,plot（w/pi,angle（H））;
title（'相位相应'）;
subplot（2,2,3）,plot（w/pi,dbH）;
title（'幅度相应（dB）'）;
subplot（2,2,4）,zplane（b,a）;
title（'零极点图'）;
程序执行结果：
n=2
wc=0.28630.7137
b=0.22920-0.458400.2292
a=1.0000-0.00000.2675-0.00000.1843
传递函数为
图形如下图所示
图6-巴特沃斯型数字带通滤波器设计
第7章
例7-用矩形窗设计一个FIR数字滤波器，要求：
N=16，截止频率为
，绘制理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。
程序如下：
wc=0.5pi;%输入设计指标
N=16
tao=（N-1）/2;
n=0:
N-1;
m=n-tao+eps;
hd=sin（wc*m）./（pi*m）;
windows=（boxcar（N））';%使用矩形窗，将列向量变为行向量
b=hd.*windows;%求FIR系统函数系数
[H,w]=freqz（b,1）;%求解频率特性
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;%化为分贝值
subplot（2,2,1）,stem（n,hd）;
axis（[0,N,1.1*min（hd）,1.1*max（hd）]）;title（'理想脉冲响应'）;
xlabel（'n'）;ylabel（'hd（n）'）;
subplot（2,2,2）,stem（n,windows）;
axis（[0,N,0,1.1]）;title（'窗函数特性'）;
xlabel（'n'）;ylabel（'wd（n）'）;
subplot（2,2,3）,stem（n,b）;
axis（[0,N,1.1*min（b）,1.1*max（b）]）;title（'实际脉冲响应'）;
xlabel（'n'）;ylabel（'h（n）'）;
subplot（2,2,4）,plot（w/pi,dbH）;
axis（[0,1,-80,10]）;title（'幅度频率响应'）;
xlabel（'频率（单位：
\pi）'）;ylabel（'H（e^{j\omega}）'）;
set（gca,'XTickMode','manual','XTick',[0,wc/pi,1]）;
set（gca,'YTickMode','manual','YTick',[-50,-20,-3,0]）;grid
运行结果如下：
图6-FIR数字带通滤波器设计
13．FDATool工具
FDATool是MATLAB中的一个图形化的滤波器设计与分析工具。
利用FDATool这一工具，我们可以进行FIR和IIR数字滤波器的设计，并且能够显示数字滤波器的幅频响应、相位响应以及零极点分布图等。
产生的数字滤波器系数在存储为文件后，可以直接提供给DSP程序代码调试工具—CCS或DSP存储器，以完成实际的数字滤波器的调试程序，从而实现实际的滤波器。
和在MATLAB命令窗口中输入命令设计滤波器相比，利用FDATool工具设计滤波器更加直观、方便。
在MATLAB命令窗口中输入命令“fdatool”将打开FDATool工作界面，如图1所示：
图1FDATool工作界面
下面通过一个实例来介绍FDATool工具的使用方法。
例子：
已知数据采样频率为1000Hz，现要设计一6阶的巴特沃斯低通滤波器，截止频率为200Hz，求其幅度响应、相位响应、脉冲响应、零极点图、滤波器系数等。
具体步骤如下
（1）根据任务首先确定滤波器种类、类型等指标。
本题应选用Lowpass、IIR、巴特沃思滤波器。
（2）根据指标给定的阶数【FilterOrder】一栏中应选择【SpecifyOrder】，并输入滤波器的阶数，输入采样频率、截止频率等指标。
（3）指标输入完后，按【DesignFilter】进行滤波器设计，将显示图2的结果
（4）观察滤波器的其它响应形式。
利用菜单按钮可以观察其幅度响应、相位响应、脉冲响应、零极点图、滤波器系数等。
图3给出了滤波器的幅度响应、相位响应、冲激响应、和零极点图。
显示当前滤波器的幅度响应。
显示当前滤波器的相位响应。
同时显示当前滤波器的幅度和相位响应。
显示当前滤波器的群延迟响应。
显示当前滤波器的相延迟。
显示当前滤波器的脉冲响应。
显示当前滤波器的阶跃响应。
显示当前滤波器的零极点的位置。
列出了当前滤波器的系数。
列出了当前滤波器的相关信息。
图2滤波器设计结果
a）幅度响应b）相位响应
c）冲激响应d）零极点图
图3滤波器的其它响应形式
a）幅度响应b）相位响应c）冲激响应d）零极点图
利用以上例子，我们可以在MATLAB命令窗口中输入命令来设计滤波器，并和使用fdatool工具设计的结果进行对比。
MATLAB中设计滤波器及查看相应特性的命令如下：
滤波器系数：
[b,a]=butter（6,200/1000*2）%Wn=fc/（Fs/2）
幅度响应、相位响应：
freqz（b,a,128,1000）
冲激响应：
impz（b,a）
零极点图：
zplane（b,a）
在MATLAB命令窗口中输入上述命令后，MATLAB将给出所要设计的滤波器的系数并绘制出该滤波器的幅度响应、相位响应、冲激响应和零极点图，图形结果如图4所示。
从图中可以看出使用Matlab命令设计的滤波器特性和使用fdatool工具设计出来的滤波器特性结果基本一致。
a）幅度响应b）相位响应
c）冲激响应d）零极点图
图4Matlab命令设计的滤波器特性
a）幅度响应b）相位响应c）冲激响应d）零极点图
14．SPTool工具
SPTool是MATLAB工具箱中自带的交互式图形信号处理工具，它包含了信号处理工具箱中的大部分函数，可以方便快捷地完成对信号、滤波器及频谱的分析、设计和浏览。
从工程角度看SPTool是一个非常实用的工具。
以下我们将结合例子来具体说明该工具的使用方法。
在MATLAB命令窗口中输入“sptool”，将打开SPTool工作界面，如图5所示：
下面通过一个实例说明SPATool工具的使用方法，
例子：
系统采样频率为1024Hz，产生100Hz和400Hz的合成正弦波，观察信号波形和频谱；设计滤波器去除400Hz分量，观察滤波后的信号波形和频谱。
（1）在MATLAB命令窗口中输入“sptool”，打开SPTool工
升级会员 