数字信号处理讲义
MATLAB系统分析
MATLAB信号处理工具箱提供了许多工具函数来分析滤波器（系统）的特性，包括时域分析，对任意输入响应、脉冲响应等；频域分析，幅值响应、相位响应、零极点位置；其他特性，群延迟等。滤波器时域和频域分析是设计各类滤波器、评价滤波器性能和滤波器应用基础。这里简单介绍一下有关函数。
一．时间响应分析工具filter,filtic,fftfilt和impulse
1．Filter（）
函数filter用于实现IIR和FIR滤波器对数据滤波，常用来计算滤波器对输入的响应，函数调用格式为：
y=filter(b,a,x)
其中，b，a分别是滤波器系统函数H(z)的分子和分母多项式系数向量；x
为滤波器输入向量或矩阵；
函数filter可用于全极点、全零点、零极点系统。 2．Filtic（）
函数filtic用于从系统过去值y和x求系统状态初始值，调用格式：
z=filtic(b,a,y,x) z=filtic(b,a,y)
其中，y=[y(n-1),y(n-2),..,y(-na)]为输出y的过去向量；x=[x(n-1),x(n-2),..,x(-nb)]为输出x的过去向量;b,a分别为分子、分母多项式系数向量；z为系统初始状态。如果，n<=-1时，x(n)=0,则filtic函数中不用指定x。
22
数字信号处理讲义
【例3－2】求解 y(n)?n31y(n?1)?y(n?2)?x(n),n?0 22其中 x(n)??14?u(n)，初始条件为y(-1)=4和y(-2)=10。
解：用MATLAB计算全响应：
n=[0:7]; x=(1/4).^n; xic=[1,2]; b=[1];a=[1,1.5,0.5]; format long y1=filter(b,a,x,xic);
Matlab提供了filtic函数计算 对于上例： Y=[4,10]; xic=filter(b,a,Y)
3．Fftfilt（）
函数fftfilt利用效率高的基于FFT重叠相加法实现对数据滤波，该函数只适用于FIR滤波器，函数调用格式为：
y=fftfilt(b,x) y=fftfilt(b,x,n)
式中，b为FIR滤波器系数向量；x为输入数据；n为FFT长度， 省时，
函数选择最佳的FFT长度；y为滤波器输出。
Y=fftfilt(b,x)等价于y=filter(b,1,x)。
【例3－3】：FIR低通滤波器截止频率为200Hz，采样频率Fs=1000Hz，对信号
x(t)=sin(2πf1t)+sin(2πf2t)滤波，f1=50Hz,f2=250Hz,求滤波器输出。
%Matlab Program 3.2
23
数字信号处理讲义
%Create Orignal signal data N=1000;Fs=1000; n=[0：N-1];t=n/.Fs;
x=sin(2*pi*50*t)+sin(2*pi*240*t); %Create a FIR filter with lawpass b=fir1(40,200/500); %Filtering the data yfft=fftfilt(b,x,256); ％Output n1=81:241;
t1=t(n1);x1=x(n1);y1=yfft(n1); subplot(221);plot(t1,x1); tilte(‘Original Signal’); subplot(222);plot(t1,y1); title(‘Signal after the filter’); grid;
执行上面这个程序，分析程序结果。 4．Impz（）
impz用于产生数字系统的脉冲响应。基本调用格式为：
[h,t]=impz(b,a,n,Fs)
b,a分别为系统传递函数分子和分母多项式系数向量；n为采样点数；Fs为采样周期，缺省值＝1；h为单位脉冲响应向量；t为和h对应的时间向量。当函数输出缺省时，绘制系统脉冲响应图；当n缺省时，函数自动选择n值。
在Matlab中，h=impz(b,a,n)等价于下面语句：
24
数字信号处理讲义
x=[1 zeros(1,n)]; h=filter(b,a,x);
因此，不少场合下，用后一种方法求系统的脉冲响应更为灵活。
5． Conv（）
线性时不变系统的输入输出关系还可通过单位脉冲响应h(n)表示：
y(n)=x(n)*h(n)=?x(k)h(n?k)
k?0N所以可直接采用MATLAB中的函数conv或者用可求出带下标的序列卷积的函数conv_m来求系统对输入的响应。
二．
1． Freqs()
函数freqs用于求模拟滤波器的频率响应。调用格式为
h=freqs(b,a,ω)； [h,ω]=freqs(b,a,n); freqs(b,a)
其中，b,a分别为模拟滤波器传递函数分子和分母多项式向量；n 为频率点数(整数),缺省＝200；h为频率响应,复数；ω频率向量，实数。
函数缺省时，绘制模拟滤波器的幅频响应图和相频响应图。 模拟滤波器的系统函数形式为
频率响应
b(1)snb?b(2)snb?1?...?b(nb?1)H(z)?????(式3?7) nana?1a(1)s?a(2)s?...?a(na?1)输出h为模拟滤波器的频率响应H(ejω)。 2．Freqz（）
函数freqz用于求数字滤波器的频率响应。调用格式为
[h,ω]=freqz(b,a,n); [h,f]=freqz(b,a,n,Fs);
[h,ω]=freqz(b,a,n,’whole’); [h,f]=freqz(b,a,n,’whole’,Fs);
25
数字信号处理讲义
[h]=freqz(b,a,f,Fs); freqz(b,a);
说明：
（1）函数输入：b,a分别为数字滤波器传递函数分子和分母多项式向
量；n为频率点数(整数),最好为2的幂，缺省＝512；Fs为采样频率，Hz；f为给定的频率矢量；’whole 表示返回的ω值为包含z整个单位圆频率矢量，即0～2π；缺省时，ω仅为包含z平面上半单位圆（0~π）之间等间距N个点频率矢量。
函数返回：h为频率响应,复数；ω为n点频率向量（弧度），返
回范围与输入参量‘whole’有关；f为n点频率向量(Hz)。
函数缺省时，绘制模拟滤波器的幅频响应图和相频响应图。 （2）该函数适用于和式3－1形式的数字滤波器。 （3）该函数执行补零的FFT算法，即h=fft(b,n)./fft(a,n)
（4）Matlab提供两个函数ABS和ANGLE由频率响应H(ejω)求幅频响应| H(ejω)|和相频响应arg[H(ejω)].
（5）函数freqz输出的频率向量ω在0～2π之间，为了获得一个滤波器真正的相频特性图，要对相位角ω进行修正。为此Matlab信号处理工具箱提供一个函数UNWRAP. 3．UNWRAP.
函数UNWARP展开由函数freqz产生的频率ω.调用格式为：
P=unwrap(ω) 三．
零极点图
系统零极点位置决定了系统稳定性和性能，因此，考察系统零极点位置是分析系统特性重要方面之一。
26