发布时间 : 星期二 文章数字信号处理上机指导 - 1dl更新完毕开始阅读a620bd4b852458fb770b56fc
7．y=filter(b , a , x)或[y , zf]=filter(b , a , x , xic)
系统以传递函数的形式表达，b、a为其分子、分母多项式的系数向量。信号x通过该系统后的零状态响应（即该系统对信号x的滤波），可用命令y = filter(b , a , x)或[y,zf]= filter(b , a , x)来完成，其中zf是系统的最后状态。
由于存储空间限制常常需要进行分布滤波，估算初始和最后条件对分布滤波来说是非常有用的。假设现有两部分数据，每部分有5000个点，如下：
x1=randn(5000 , 1);
x2=randn(5000 , 1);
可能第一个序列x1对应前10分钟采集的数据，第二个序列x2对应后10分钟采集的数据。那么整个序列为x=[x1 ; x2]。如果现在没有足够的空间存储连接起来的序列x,对x1、x2分别进行滤波，为了保证滤波后序列的连续性，采用滤波x1的最后条件作为滤波x2的初始条件。
[y1 , zf] = filter(b , a , x1);
y2 = filter(b , a , x2 , zf);
若计算系统在初始状态（Y , X）下的全响应，可用命令zf=filtic(b , a , Y , X) 将（Y , X）等效成输入xic，再用命令[y , zf] = filter(b , a , x , xic)来完成。
例 求解差分方程y(n)=[x(n)+x(n-1)+x(n-2)]/3+0.95y(n-1)-0.9025y(n-2)，n≥0。其中x (n) = cos (πn / 3)，且y (-1) = -2 ，y (-2)=-3。 用MATLAB语句实现如下： function y=chafen
b=[1,1,1]/3;
a=[1,-0.95,0.9025]; Y=[-2,-3];X=[0,0];
xic=filtic(b,a,Y,X); %初始状态等效成输入xic n=0:50;
x=cos(pi*n/3);
y=filter(b,a,x,xic); %全响应
plot(n,y) %或stem(n,y) 或者：
function chafen %用递推法求解 n=0:50; x0=cos(pi*n/3); x=[0,0,x0];
y=[-3,-2,zeros(1,51)]; for i=3:53
y(i)=(x(i)+x(i-1)+x(i-2))/3+0.95*y(i-1)-0.9025*y(i-2); end
z(1:51)=y(3:53); plot(n,z)
第6章 MATLAB在数字信号处理中的应用
§6-1圆周卷积的计算
对于无限长序列不能用MATLAB直接计算线性卷积，在MATLAB内部只提供了一个conv函数计算两个有限长序列的线性卷积。对于圆周卷积MATLAB内部没有提供现成的函数,我们可以按照定义式直接编程计算。 例6-1 已知两序列：
- 16 -
0.8n10?n?11 x(n)?? h(n)?? 0?n?5 ?0 ?其它其它??0求它们的线性卷积yl(n)=h(n)*x(n)和N点的圆周卷积y(n)=h(n)的关系。
x(n)，并研究两者之间
实现程序： （1）计算圆周卷积的函数
function yc=circonv(x1,x2,N) %实现两序列x1和x2的圆周卷积 if length(x1)>N
err0r('N must not be less than length of x1');
end
if length(x2)>N
err0r('N must not be less than length of x2');
end
x1=[x1,zeros(1,N-length(x1))]; %填充序列x1(n)使其长度为N x2=[x2,zeros(1,N-length(x2))]; %填充序列使x2(n)其长度为N m=[0:1:N-1];
x2=x2(mod(-m,N)+1); %生成x2的圆周反转序列 H=zeros(N,N);
for n=1:1:N %生成计算圆周卷积的矩阵 H(n,:)=cirshiftd(x2,n-1,N); % x2圆周右移n位 end
yc=x1*H '; %计算圆周卷积
function y=cirshiftd(x,m,N) %序列的圆周移位 if length(x)>N
error('The length of x must be less than N');
end
x=[x,zeros(1,N-length(x))]; %补零，长度变为N n=[0:1:N-1];
y=x(mod(n-m,N)+1); %得到输出
（2）研究两者之间的关系
function fun5 clear all;
[N,xn,hn]=fun50;
yln=conv(xn,hn); %直接用函数conv计算线性卷积
ycn=circonv(xn,hn,N); %用函数circonv计算N点的圆周卷积 ny1=[0:1:length(yln)-1]; ny2=[0:1:length(ycn)-1]; n=0:length(xn)-1; m=0:length(hn)-1; subplot(2,2,1);
stem(n,xn); xlabel('xn') subplot(2,2,2);
stem(m,hn); xlabel('hn') axis([0,16,0,4]); subplot(2,2,3);
stem(ny1,yln); xlabel('线性卷积') subplot(2,2,4);
stem(ny2,ycn);xlabel('圆周卷积')
function [N,xn,hn]=fun50 n=[0:1:11];
xn=0.8.^n; %生成x(n)
- 17 -
hn=ones(1,6); %生成h(n) N=12;
运行结果如右图所示。 §6-2 利用FFT实现线性卷积
若序列x1(n)、x2 (n)是长度分别为N1、N2的有限长序列，N点的圆周卷积yc(n)=x1(n)x2(n)，长度为N1+N2-1的线性卷积yl(n)=x1(n)*x2(n)。由DFT的性质可知：当N≥N1+N2-1时有yc(n)=yl (n)=IDFT{DFT [x1(n)]·DFT[x2 (n)]}。序列较长时DFT运算通常用快速算法FFT实现。
例6-2 用FFT实现上节例中的两序列的线性卷积。 实现程序如下：
%fun6.m
n=[0:1:11];m=[0:1:5];
N1=length(n); N2=length(m); xn=0.8.^n; %生成x(n) hn=ones(1,N2); %生成h(n) N=N1+N2-1; xk=fft(xn,N); hk=fft(hn,N); yk=xk.*hk; yn=ifft(yk,N);
if all(imag(xn)==0)&(all(imag(hn)==0)) yn=real(yn);
end
x=0:N-1;
stem(x,yn,'.'); 运行结果为如图。
§6-3 系统的单位冲激响应
线性移不变系统的特性可用它在输入信号是单位冲激序列时的输出即系统的单位抽样响应h(n)来表征。有了它，就可得到此线性移不变系统对任意输入x(n)的输出y(n)。 例6-3 一个因果线性移不变系统：y(n)=0.81y(n-2)+x(n)-x(n-2)。 求：（1）单位冲激响应；（2）单位阶跃响应；（3）绘制系统的幅频和相频响应。
由差分方程直接得出系统函数H (z) = (1-z – 2 ) / (1-0.81z – 2 )。从而可计算其它问题如下： function [h,y]=chongji
b=[1,0,-1]; a=[1,0,-0.81];
xh=[1,zeros(1,50)]; % 产生单位冲激序列 xn=ones(1,51); % 产生单位阶跃序列 h=filter(b,a,xh); % 单位冲激响应 y=filter(b,a,xn); % 单位阶跃响应 figure(1); n=0:50;
subplot(2,1,1); stem(n,h,'.r')
legend('冲激响应'); subplot(2,1,2); stem(n,y,'k');
legend('阶跃响应');
% 绘制系统的幅频和相频响应 figure(2);
freqz(b,a);
- 18 -
§6-4 IIR滤波器的设计与实现
Fig.2 Fig.1
无限长单位冲击响应（IIR）数字滤波器的设计方法有多种，对于一般条件下
使用的数字滤波器，其常用的设计方法是基于模拟滤波器变换原理的经典设计：首先，将数字滤波器的技术指标转换成相应的模拟滤波器技术指标；设计模拟低通滤波器；然后，将设计好的模拟滤波器转换成满足给定技术指标的数字滤波器（常用算法有脉冲响应不变法和双线性变换）。在MATLAB的数字信号处理工具箱中，提供的有关设计函数如下：
1．模拟低通滤波器设计
模拟低通滤波器的逼近有巴特沃思型、切比雪夫型、和考尔型，用如下的函数实现。 [Z,P,K]=buttap(n)；返回一个n阶巴特沃斯型归一化的模拟低通滤波器的零极点增益模型 [Z,P,K]=cheb1ap(n,Rp); n阶、通带内的最大衰减Rp、切比雪夫Ⅰ型
[Z,P,K]=cheb2ap(n,Rs); n阶、阻带内的最小衰减Rs、切比雪夫Ⅱ型
[Z,P,K]=elliap(n,Rp,Rs); n阶、通带内的最大衰减Rp、阻带内的最小衰减Rs、考尔型。 2．模拟低通滤波器阶数n的选择函数
滤波器阶数的选择在整个滤波器的设计中占有十分重要的地位和作用。根据需要选择合适的滤波器阶数，MATLAB工具箱中提供了对应于各类模拟低通滤波器的阶数选择函数，如巴特沃思型的buttord、切比雪夫型的cheb1ord、cheb2ord和考尔型的ellipord。这些函数的调用格式大同小异，仅以buttord为例加以说明。
[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)
输入参数：Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减。 输出参数：n为符合要求的滤波器最小阶数，Wn为巴特沃思型模拟低通滤波器3dB截止频率。’s’：表示模拟域。
3．零极点增益模型到传递函数模型的转换
[num,den]=zp2tf(Z,P,K)
输入参数：Z ，P，K分别表示系统的零极点增益模型的零点、极点和增益； 输出参数：num ，den分别为同一系统传递函数模型的分子和分母多项式系数。
4．模拟域的频率变换
将归一化的模拟低通滤波器转换成所需要类型（低通、高通、带通和带阻）的模拟滤波器，可分别用如下命令实现：
[b,a]=lp2lp(Bap,Aap,Wn)；把传递函数形式的归一化模拟低通滤波器原型转换成3dB截止频率为Wn的同型低通滤波器。
[b,a]=lphp(Bap,Aap,Wn)；转换成高通
[b,a]=lp2bp(Bap,Aap,W0,Bw)；转换成带通，W0：中心频率，Bw：带宽 [b,a]=lp2bs(Bap,Aap,W0,Bw)；转换成带阻
5．模拟滤波器数字化
[bz,az]=bilinear(b,a,Fs)：采用双线性变换法的映射关系。其中，Fs是采样频率。 [bz,az]=impinvar(b,a,Fs)：采用冲击响应不变法的映射关系。
例6-4 用双线性变换法设计一个Butterworth低通滤波器，要求其通带截止频率100Hz，阻带截止频率200Hz，通带衰减Rp小于2dB，阻带衰减大于15dB，采样频率Fs=500Hz。 MATLAB实现：
- 19 -