Matlab实验指导书ok - 图文

数字信号处理实验指导书

（基于MATLAB）

尚秋峰 宋文妙

华北电力大学 二○○六年十一月

前 言

MATLAB 是一套功能强大的工程计算及数据处理软件，广泛应用于工业，电子，医疗和建筑等众多领域。它是一种面向对象的，交互式程序设计语言，其结构完整，又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，MATLAB提供了方便的绘图功能，便于用户直观地输出处理结果。

本课程实验要求学生掌握用MATLAB语言仿真数字信号处理算法的基本方法，加深对教学内容的理解。

说明：本书给出的为本课程基本实验，教师可根据学生情况安排内容延伸实验。

目 录

实验一 序列的产生和运算 ........................................... 1 实验二 因果性数字系统的时域实现 ................................... 5 实验三 离散傅里叶变换（DFT）及其快速算法（FFT） ................... 10 实验四 FFT的典型应用 ............................................. 17 实验五 FIR滤波器设计 ............................................. 23 实验六 IIR滤波器设计 ............................................. 29 附 录 ............................................................ 32

实验一 序列的产生和运算

一、实验目的：

1、熟悉MATLAB环境，掌握基本编程方法。 2、熟悉MATLAB中序列产生和运算的基本函数。 二、实验内容： 1、MATLAB入门

（1）了解MATLAB 的工作窗口

进入MATLAB环境，选择View中的各个选项，可以打开命令窗口Command Window、变量浏览器Workspace、路径浏览器Current Directory，如下图示：

图1－1 Matlab 工作窗口

MATLAB命令窗口，是键入指令的地方也是MATLAB显示计算结果的地方。在符号>>之后键入“1+2+3”，或者“x=1+2+3”，显示结果有什么不同？如果在上述的例子结尾加上“;”，则计算结果不会在窗口自动显示，要显示计算结果须键入该变量x。

进行以上操作后打开变量浏览器，观察里面的变量名和变量的值。

在MATLAB环境选择File 再选择New，即进入程序编辑/调试窗口，如下图示。存档时必须以.m名称储存。要执行 M-file 可以在命令窗口下直接键入该文件名；或是选择Debug

1

下的Run M-file来执行M-file。如果要修改 M-file 可以选择功能表上的Open M-file ，即可搜寻要修改的 M-file，修改后再存档。尝试编写程序文件，完成以上操作。

图1－2 程序编辑/调试窗口

练习使用路径浏览器打开特定目录下的M-file。

M-file还可以用来定义函数，然后储存起来，就可以和那些内建的函数（如sin, cos,log等）一样的自由使用。举例来说，我们可以定义一函数cirarea来计算圆的面积，以下的 M-file: cirarea.m就是定义这个函数的：

% M-file function, cirarea.m

% Calculate the area of a circle with raduis r % r can be a scalar or an array

function c=cirarea(r)

c=pi*r.^2;

注意，M-file定义的函数语法上有一些规定：

? 第一行指令以function这个字做为起头，接着是输出的变量，等号，函数名称，输入的

变量是放在括号之内。function out1=userfun(in1)，这行的out1是输出的变量，userfun是函数名称，in1是输入的变量。function [out1, out2]=serfun(in1, in2) 如果输出变量 [out1,out2]和输入变数(in1, in2)不只一个时，则在输出变量部份须加上[ ]。 ? 上述的输入变量是在调用函数时输入的，而输出的变量即是函数的返回值。 ? 函数名称的取法的规定与一般变量相同。

? 在定义函数之前，最好加上注解行来说明这个函数的特色及如何使用，这样，使用指令

如help cirarea，该函数的注解行会出现在指令视窗。

2

尝试编写函数文件，并在你编写的程序文件中应用此函数。

MATLAB会将绘图结果展示在另一个视窗，称为MATLAB Figure Windows，如下图示。如果你看不到此视窗，可以进入Windows再选择Figure。在图形窗口中，可以利用Edit菜单中的选项来改变显示效果以及拷贝图形，尝试这些操作。

图1－3 MATLAB 图形视窗

（2）寻求帮助

在MATLAB系统中相关的线上(on-line)求助方式有三种：

? 利用help指令，如果你已知要找的主题为何的话，直接键入help 。所以即使身

旁没有使用手册，也可以使用help指令查询不熟悉的指令或是题材之用法，例如help sqrt, help topic。

? 利用lookfor指令，它可以从你键入的关键字(key-word)（即使这个关键字并不是

MATLAB的指令）列出所有相关的主题，例如lookfor cosine, lookfor sine。

? 利用指令视窗的功能选单中的Help，从中选取Table of Contents（目录）或是Index（索

引）。

? 练习以上寻求帮助的方法。

2、波形产生

学习附录中的有关知识，完成下列练习。

练习一：建立一个MATLAB函数impseq.m，用来生成迟延的单位脉冲序列。其输入参数为：序列起始位置ns，序列终止位置nf，脉冲位置np。

练习二：编写MATlAB程序，以产生下图所示的波形。并将序列绘制出来。

3

练习三：设x＝[2 3 4 5]，求将它延拓6个周期所得的序列。 练习四：设x＝[2 3 4 5]，实现以下运算并绘制波形

练习五： 编写程序实现下列函数，并绘制波形

x(t)?2sin(4?t)?5sin(8?t)?0.8w(t)t以t=0.01（n=0:N-1）进行取样，N=64 w(n)?randn(1,N)三、思考题：

算术运算符*和．*之间的区别是什么？

4

以

实验二 因果性数字系统的时域实现

一、 实验目的

编制实现下列差分方程的程序:

y(n)=∑bkx(n-k)+∑aky(n-k)

二、 实验内容

1、编制nonrec.m函数文件,实现FIR滤波y(n)=h(n)*x(n)。这里给定h(n)=R8(n)，x(n)= nR16(n)，求y(n)。

nonrec.m函数文件: function y=nonrec(x,h)

x=[x,zeros(1,length(h)-1)]; %补零 w=zeros(1,length(h)); for i=1:length(x) for j=length(h):-1:2

w(j)=w(j-1); %得到每一延时单元输出变量 end w(1)=x(i);

y(i)=w*h.’; %h与w对应相乘 end

主程序文件: x=0:15; h=ones(1,8); y=nonrec(x,h); n=0:22; stem(n,y);

图2－1 卷积运算实现FIR滤波

5

分析：线性卷积y(n)=x(n)*h(n)的长度为16+8-1=23，可利用y(n)=∑h(m)x(n-m)直接计算得

n(n+1)/2, n≤7 y(n)= 4(2n-7), 8≤n≤15 (n+8)(23-n)/2, 16≤n≤22

即 y=[ 0 1 3 6 10 15 21 28 36 44 52 60 68 76 84 92 84 75 54 42 29 15] ,与曲线相符。

2、编制rec.m函数文件,实现纯递归IIR滤波 y(n)=x(n)+∑aky(n-k)。 这里给定a1=2rcosw0, a2=-r2

, r=0.95, w0=π/8, 求单位抽样响应h(n).

rec.m函数文件:

function y=rec(x,a,n)

x=[x,zeros(1,n-length(x))]; %补零到所需长度 sum=0;

w=zeros(1,length(a)); for i=1:n

y(i)=sum+x(i); % 递归 for j=length(a):-1:2

w(j)=w(j-1); %延时 end w(1)=y(i); sum=w*a'; end 主程序文件:

x=[1];

a=[2*0.95*cos(pi/8),-0.95^2];

h=rec(x,a,75); %取h(n)的长度为75点 n=0:74; stem(n,h);

图2－2 纯递规IIR滤波

分析计算:

6

由题意, a1=2*0.95*cos(π/8), a2=-0.95, 所以, 得到系统函数 H(z)=1/[1-1.9cos(π/8)z+0.95z],

做逆Z变换得 h(n)=0.95ncos(πn/8)+ctg(π/8)*0.95nsin(πn/8), 利用MATLAB直接画h(n), 即,使用下列语句 n=0:74;

h=0.95.^n.*cos(pi.*n./8)+cot(pi/8).*(0.95.^n).*sin(pi.*n./8); stem(n,h);

得到如图2－3的曲线

-1

2

-2

2

图2－3

此曲线与用rec函数所画曲线完全一致, 可见, 用MATLAB编制的FIR滤波程序与理论计算所得结果是相同的.

3、用nonrec.m和rec.m函数编制dfilter.m 函数文件, 构成完整的一般IIR滤波程序, 并完成下列信号滤波: x(n)=cos(2πn/5)R64(n)

这里给定系统函数 H(z)=(1-2z-1+z-2)/(1-0.5z-1+0.5z-2) 求y(n).

dfilter.m函数文件:

function y=dfilter(x,b,a,n) y1=nonrec(x,b); y=rec(y1,a,n); 主程序文件: n=0:63; x=cos(2*pi/5*n);

b=[1,-2,1]; %由H(z)得到系数a,b

7

a=[0.5,-0.5];

y=dfilter(x,b,a,64); %取y(n)的长度为64点 stem(n,y);

图2－4

4．用help查看内部函数filter.m, 了解调用格式, 重做3, 并与你编写的dfilter.m结果进行比较.

用help可以看到内部函数为Y = FILTER(B,A,X), 且有 “The filter is a Direct Form II Transposed implementation of the standard difference equation”:

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

-a(2)*y(n-1) - ... - a(na+1)*y(n-na) 因此, 调用内部函数filter时, 要对原系数a做适当变化:

a1(1)=1, a1(i)=-a(i-1), i>1. n=0:63; x=cos(2*pi/5*n); b=[1,-2,1];

a=[0.5,-0.5];

y=dfilter(x,b,a,64); %调用自己编的dfilter函数

a1=[1,-0.5,0.5]; %a变为a1, 用于调用内部函数filter y1=filter(b,a1,x); subplot(2,1,1); stem(n,y); subplot(2,1,2); stem(n,y1);

8

图2－5

三、 思考题

1、 内部函数filter.m的调用格式是什么?与你编写的dfilter.m调用格式是否一致?差异在何处?

2、 实验中如何合理地选取单位抽样响应h(n)的点数?实验内容及要求中第三项的物理概念是什么?给定的H(z)是具有什么性质的滤波器?

3、 为实现各实验内容中的滤波器，你编写程序时采用的什么结构？

4、 试利用Matlab中的卷积函数conv完成实验内容1,并与本书的例程进行比较。

9

实验三 离散傅里叶变换（DFT）及其快速算法（FFT）

一、实验目的

1、学习DFT和FFT的原理及编程实现方法。 2、熟悉MATLAB编程的特点。 二、实验内容

1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X (ej?)，并比较三种程序计算机运行时间。

（1）用for loop 语句的M函数文件dft1.m，用循环变量逐点计算X(k)； （2）编写用MATLAB矩阵运算的M函数文件dft2.m, 完成下列矩阵运算； （3）调用FFT库函数，直接计算X(k)；

（4）分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换X(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。

（一）实验例程 M函数文件如下： dft1.m:

function[Am,pha]=dft1(x)

N=length(x);

w=exp(-j*2*pi/N); for k=1:N sum=0; for n=1:N

sum=sum+x(n)*w^((k-1)*(n-1)); end

Am(k)=abs(sum); pha(k)=angle(sum); end

dft2.m:

function[Am,pha]=dft2(x)

N=length(x);

n=[0:N-1]; k=[0:N-1];

w=exp(-j*2*pi/N); nk=n'*k;

10

wnk=w.^(nk); Xk=x*wnk; Am=abs(Xk); pha=angle(Xk);

dft3.m:

function[Am,pha]=dft3(x) Xk=fft(x); Am=abs(Xk); pha=angle(Xk);

源程序及运行结果：

（1）x=[ones(1,8),zeros(1,248)];

t=cputime;

[Am1,pha1]=dft1(x); t1=cputime-t n=[0:(length(x)-1)]; w=(2*pi/length(x))*n; figure(1)

subplot(2,1,1), plot(w,Am1,'b'); grid; title('Magnitude part'); xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

subplot(2,1,2), plot(w,pha1,'r'); grid; title('Phase Part');

xlabel('frequency in radians'); ylabel('argX[exp(jw)]/radians');

11

图3－1 DFT1

记录运行时间：t1=

2）x=[ones(1,8),zeros(1,248)];

t=cputime;

[Am2,pha2]=dft2(x); t2=cputime-t n=[0:(length(x)-1)]; w=(2*pi/length(x))*n; figure(2)

subplot(2,1,1), plot(w,Am2,'b'); grid; title('Magnitude part'); xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

subplot(2,1,2), plot(w,pha2,'r'); grid; title('Phase Part');

xlabel('frequency in radians'); ylabel('argX[exp(jw)]/radians');

（ 12

图3－2 DFT2

记录运行时间：t2=

（3）x=[ones(1,8),zeros(1,248)]; t=cputime;

[Am3,pha3]=dft3(x); t3=cputime-t; n=[0:(length(x)-1)]; w=(2*pi/length(x))*n; figure(3)

subplot(2,1,1), plot(w,Am3,'b'); grid; title('Magnitude part'); xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

subplot(2,1,2), plot(w,pha3,'r'); grid; title('Phase Part');

xlabel('frequency in radians'); ylabel('argX[exp(jw)]/radians');

13

图3－2 DFT2

记录运行时间：t3=

从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。

（二）练习

改变输入信号类型和点数，重做例程中的内容，并整理总结实验结果，给出结论。 2、实现序列的内插和抽取所对应的傅里叶变换。

给定序列x(n)=[cos(π/36×n)+cos(1.5π/36×n)]R128(n)，做N=128点的傅里叶变换，并求

x1(n)=x(4n) 和 x2(n)??对应的傅里叶变换(N=128点)。

（一）实验例程 源程序如下：

N=32; %抽取32点 t=0:N-1;

x2=cos(pi/9*t)+cos(1.5*pi/9*t);

Xk2=fft(x2,128); %做128点傅里叶变换 Am2=abs(Xk2); n=[0:127];

14

?x(n)?0n?4kn?4k

w=(2*pi/128)*n; figure(2) plot(w,Am2,'b'); title('Magnitude part'); xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

x3=zeros(1,128); %初始化数组为零 N=32;

for i=0:N-1 %n=4i时赋值 x3(1,4*i+1)=cos(pi/36*i)+cos(1.5*pi/36*i); end

Xk3=fft(x3); %做128点fft Am3=abs(Xk3); n=[0:127]; w=(2*pi/128)*n; figure(3) plot(w,Am3,'b'); title('Magnitude part'); xlabel('frequency in radians'); ylabel('|X(exp(jw))|'); 抽取后的离散傅里叶变换：

图3－4 抽取的离散傅里叶变换

内插后的离散傅里叶变换：

15

图3－3 内插的离散傅里叶变换

（二）练习

仿照例程编写给定序列x(n)=[cos(π/36×n)+cos(1.5π/36×n)]R128(n)的128点傅里叶变换程序，比较这三个计算结果得到的幅频特性图，分析其差别产生的原因。

三、思考题

实验内容4中，分别求X1(e

j?w

)和X2(e

j?w

)与原X(e)的关系，并说明为什么抽取后能分

j?w

辨出两个频率分量？有无任何条件的限制？

16

实验四 FFT的典型应用

一、实验目的

1、学习FFT的应用于卷积运算和谱分析时的编程实现方法。 2、熟悉MATLAB编程的特点。 二、实验内容

1、利用FFT实现两序列的卷积运算，并研究FFT点数与混叠的关系。

给定x(n)=R16(n)，h(n)=R8(n)，用FFT和IFFT分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。

（一）实验例程 源程序如下：

N=16; x=[0:N-1]; h=ones(1,8);

线性卷积：Xk1=fft(x,23); %做23点fft

Hk1=fft(h,23); Yk1=Xk1.*Hk1; y1=ifft(Yk1);

n=0:22; figure(1) stem(n,y1)

运行结果： 线性卷积：

图3－1 线性卷积

17

（二）练习

改写上面例程做18点fft,由程序运行结果总结前多少个点为混叠结果，后多少个点是线性卷积的结果。如果改做30点fft，又会出现什么结果？

依据实验结果总结：要用DFT来做线性卷积，DFT的点数必须满足什么要求？ 2、研究高密度频谱与高分辨率频谱 （一）实验例程 设有连续信号

Xa(t)=cos(2π×6.5×10^3t)+cos(2π×7×10^3t)+cos(2π×9×10^3t) 以采样频率fs=32kHz对信号Xa(t)采样,分析下列三种情况的幅频特性。 （1） 采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。

（2） 采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。（3） 采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。 观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。 源程序如下：

fs=32000; %采样频率 N=16; %采集16点 n=0:N-1; %做16点DFT

xa=cos(2*pi*6.5*10^3*n/fs)+cos(2*pi*7*10^3*n/fs)+cos(2*pi*9*10^3*n/fs); Xk1=fft (xa); Am1=abs(Xk1);

n=[0:(length(xa)-1)]; w=(2*pi/length(xa))*n; figure(1)

plot(w,Am1,'b'); grid; title('Magnitude part');

xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

x=[xa(1:16),zeros(1,240)]; %补零到256 Xk2=fft(x); %做256点DFT Am2=abs(Xk2); n=[0:(length(x)-1)]; w=(2*pi/length(x))*n; figure(2)

plot(w,Am2,'b'); grid; title('Magnitude part');

xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

N0=256; %采集256点

18

n0=0:N0-1; %做256点DFT

xa0=cos(2*pi*6.5*10^3*n0/fs)+cos(2*pi*7*10^3*n0/fs)+cos(2*pi*9*10^3*n0/fs);

Xk3=fft(xa0); Am3=abs(Xk3);

n=[0:(length(xa0)-1)]; w=(2*pi/length(xa0))*n; figure(3)

plot(w,Am3,'b'); grid; title('Magnitude part');

xlabel('frequency in radians'); ylabel('|X(exp(jw))|');

运行结果：N=16点DFT幅频特性

图3－2 N=16点的DFT

高密度频谱：

19

图3－3 补零到256点的DFT

高分辨率频谱：

图3－4 256点的DFT

观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。

将16点DFT补零到256点，做256点DFT,由结果可以看出，频率间隔缩小（由2?/16减小为2???256），得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。

第三次做的是256点的DFT，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。

（二）练习

20

练习一：

设 x(n)?cos(0.48?n)?cos(0.52?n)(1)取x(n)(0?n?100)时，求x(n)的FFT变换X(k); (2)将（1）中的 x(n)以补零方式加长到长度为100，求X(k); （3）取x(n)(0?n?100) ，求X(k) 要求：

用三个图形窗口；各图形窗口包括两个子图，分别显示x(n)和X(k)的幅值； 比较三种情况的频谱成分；

说明高密度谱（补0）和高分辨率谱（3）的不同。 尝试N=120,200等，有什么现象？并解释。 练习二： 模拟信号

x(t)?2sin(4?t)?5sin(8?t)?0.8w(t)以t=0.01n（n=0:N-1）进行取样，N=64，

w(n)?randn(1,N)试比较有无噪声时的信号谱。 尝试增加噪声的幅值。 3.重叠保留法做线性卷积

x=0:134; %x为长序列 y=[];

h=[ones(1,8),zeros(1,8)]; %h为单位抽样响应 H=fft(h); for i=0:14

if (i==0) xk=[zeros(1,7),0:8]; %x序列的前面补零M-1点(M-1=7)

else xk=x(9*i-6:9*i+9); %截取xk(n)中的N点(N=16) end

XK=fft(xk);

YK=H.*XK;

yk=ifft(YK); %利用FFT和IFFT做循环卷积 y=[y,yk(8:16)]; %截取yk(n)中的未混叠部分并衔接 end

stem(x,y);

21

图3－5 重叠保留法做线性卷积

利用MATLAB内部函数filter对此结果进行验证, 输入语句

filter(h,[1],x),

则可得到如下曲线. 与重叠保留法所得结果相同.

图3－6 直接卷积结果

三、思考题

1、什么是高密度频谱？什么是高分辨率频谱？实验内容1中为区分三个频率分量至少应采集几个样本？

2、实验内容2中，分别求X1(e)和x2(e)与原X(e)的关系，并说明为什么抽取后能分辨出两个频率分量？有无任何条件的限制？

jw?

j?w

j?w

22

%====================================================== %Hr=Amplitude Response

%w=500 frequencies between [0 pi] over which Hr is computed %c=Type-4 LP filter coefficients %L=Order of Hr

% h=Type-4 LP filter impulse response %

N=length(h);

L=(N)/2 d=2*h(L:-1:1)]; n=[1:L];n=n-0.5;

w=[0:500]'*pi/500;

Hr=cos(w*n)*d';

练习：仿照例程实现以下题目

采用频域抽样方法设计具有下列指标的低通滤波器：

?p?0.2?,Rp?0.25dB,?s?0.3?,As?50dB画出滤波器的频率响应（dB）、理想和实际的冲激响应。 三、思考题

1、你选择“合适窗函数”的依据是什么？ 2、简述利用窗函数法设计滤波器的原理与步骤。

3、利用频率取样法设计FIR滤波器时，为什么要在过渡带增加取样点？

28

实验六 IIR滤波器设计

一、实验目的

1、学习IIR滤波器的设计方法。 2、熟悉MATLAB编程的特点。 二、实验内容

利用Butterworth原型设计具有下列指标的低通数字滤波器 ??0.2?,p

Rp?1dB,?s?0.3?,As?15dB画出滤波器的频率响应（|H|，dB）、相频）。 实验参考例程： %IIR by impulse wp=0.2*pi; ws=0.3*pi; Rp=1; As=15; T=1;

omegaP=wp/T; omegaS=ws/T;

% ep=sqrt(10^(Rp/10)-1);% % Ripple=sqrt(1/(1+ep*ep)); % Attn=1/(10^(As/20)); %-----analog butter--

[cs,ds]=afd_butt(omegaP,omegaS,Rp,As); [b,a]=imp_invr(cs,ds,T); % [C,B,A]=dir2par(b,a); %---------- figure(1);

[db,mag,pha,grd,w]=freqz_m(b,a); subplot(2,2,1); plot(w/pi,mag); subplot(2,2,2); plot(w/pi,db); subplot(2,2,3); plot(w/pi,pha/pi);

29

subplot(2,2,4); plot(w/pi,grd); 各函数文件：

function [b,a]=u_buttap(N,omegaC) [z,p,k]=buttap(N); p=p*omegaC; k=k*omegaC^N; B=real(poly(z)); b0=k; b=k*B; a=real(poly(p));

function [b,a]=afd_butt(wp,ws,Rp,As)

N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws))); omegaC=wp/((10^(Rp/10)-1)^(1/(2*N))); [b,a]=u_buttap(N,omegaC);

function [b,a]=imp_invr(c,d,T) [R,p,k]=residue(c,d); p=exp(p*T); [b,a]=residuez(R,p,k); b=real(b'); a=real(a');

function I=cplxcomp(p1,p2); I=[];

for j=1:length(p2); for i=1:length(p1);

if(abs(p1(i)-p2(j))<0.0001) I=[I,i]; end end end I=I';

function [C,B,A]=dir2par(b,a); M=length(b); N=length(a);

30

[r1,p1,C]=residuez(b,a); p=cplxpair(p1,10000000*eps); I=cplxcomp(p1,p); r=r1(I); K=floor(N/2); B=zeros(K,2); A=zeros(K,3); if K*2==N;

for i=1:2:N-2

Brow=r(1:1:i+1,:);

Arow=p(1:1:i+1,:);

[Brow,Arow]=residuez(Brow,Arow,[ ]); B(fix((i+1)/2),:)=real(Brow'); A(fix((i+1)/2),:)=real(Arow'); end

[Brow,Arow]=residuez(r(N-1),p(N-1),[]); B(K,:)=[real(Brow') 0]; A(K,:)=[real(Arow') 0]; else

for i=1:2:N-1 Brow=r(1:1:i+1,:); Arow=p(1:1:i+1,:);

[Brow,Arow]=residuez(Brow,Arow,[]); B(fix((i+1)/2),:)=real(Brow'); A(fix((i+1)/2),:)=real(Arow'); end end

练习一：仿照例程用双线性变换法实现以上题目要求

练习二：利用函数buttord()和butter()编程实现以上题目要求 三、思考题

1、比较双线性变换法和冲激响应不变法的特点。 2、简述设计IIR滤波器的原理与步骤。

31

附 录

一、变量

MATLAB对使用变量名称的规定：

1、变量名称的英文大小写是有区别的（apple, Apple, AppLe，三个变量不同）。 2、变量的长度上限为 19 个字元。

3、变量名的第一个字必须是一英文字，随后可以掺杂英文字、数字或是底线。 以下列出MATLAB所定义的特别变量及其意义

变量名 help who ans eps pi inf 意义 线上说明, 如 help quit 列出所有定义过的变量名称 预设的计算结果的变量名 MATLAB定义的正的极小值=2.2204e-16 内建的π值 ∞值，无限大 （10) NaN 无法定义一个数目 ??0?? ?0?二、数组和矩阵的产生方法

如果是要个别键入元素，须用中括号[ ] 将元素置于其中。数组为一维元素所构成，而矩阵为多维元素所组成，例如 >> x = [1 2 3] % 一维 1x3 数组

>> x = [1 2 3; 4 5 6] % 二维 2x3 矩阵，以“;”区隔各列的元素 >> x = [1 2 3 % 二维 2x3 矩阵，各列的元素分二行键入

4 5 6]

上面提到数组产生的方式须个别键入其元素，这方法只适用于数组元素很少时。如果要建立的数组的元素多达数百个，则须采用以下的数种方式 >> x=(0:0.02:1) % 以:区隔起始值=0、增量值=0.02、终止值=1

>> x=linspace(0,1,51) % 利用linspace，以区隔起始值=0终止值=1之间的元素数目=51 >> x=(0:0.01:1)*pi % 注意数组外也可作运算 >> a=1:5, b=1:2:9 % 这二种方式更直接

32

a = 1 2 3 4 5 b = 1 3 5 7 9

>> c=[b a] % 可利用先前建立的数组 a 及数组 b ，组成新数组 c =

1 3 5 7 9 1 2 3 4 5

>> d=[b(1:2:5) 1 0 1] % 由数组 b 的三个元素再加上三个元素组成 d =

1 5 9 1 0 1

特殊矩阵的定义，如元素皆为0, 1或是单位矩阵，因为在运算时常会用到，所以在此先介绍。 zeros函数是形成元素皆为0 的矩阵；ones函数是形成元素皆为 1 的矩阵； eye则是产生一个单位矩阵，之所以称为eye是取其发音与原来单位矩阵符号I相同，而又避免与定义复数中的虚部所用的符号i雷同，所以改以eye替代。上述三个函数的使用语法都相似，如zeros(m)可以产生一个m×m的方阵，而zeros(m,n)产生的是m×n的矩阵。也可以使用这三个函数将一个m×n矩阵原来元素全部取代成0, 1 或是单位矩阵的值，不过要加上size指令来指出其矩阵大小是m,n，所以语法为zeros(size(A))，其中A是原来矩阵。

三、数组和矩阵的运算 数组运算功能： + 加 - 减 .* 乘 ./ 左除 .^ 次方 .’转置

MATLAB 在许多运算皆是以数组为对象，即是以数组的元素为对象。因此除了+, - 这二个运算外，其余的运算符号（乘、除、次方）皆须加上“.”，来强调数组之间的运算。例如 >> x = 1.5; % x 是纯量

>> y = exp(x^2); % exp(x^2) 是纯量运算 >> y1 = x/y % x/y 是纯量运算 >> x = 1:0.1:2; % x 是阵列

>> y = exp(x.^2); % exp(x.^2) 是阵列运算 >> y1= x./y % x./y 是阵列运算

如果运算符号（乘、除、次方）未加上“.”则表示矩阵之间的运算。 四、波形产生 单位抽样序列

33

?(n)???1n?0在MATLAB中可以利用zeros( )函数实现:

x?zeros(1,N);x(1)?1;?0 n?0

或者

x=[l zeros(1，N—1)]

如果??n?在时间轴上延迟了k个单位，得到??n?k?即：

?(n?k)???1

n?k?0 n?0

延时K个样本且长度为N的单位抽样序列x[n]，可用如下的MATLAB命令产生，其中K

x=[zeros(1，K) l zeros(1，N-K-1)]；

? 单位阶跃序列

?1u(n)???0n?0n?0

在MATLAB中可以利用ones()函数实现。

x?ones(1,N);

延时单位阶跃序列

?1n?n0u(n?n0)??。

?0n?n0产生延时单位阶跃序列的方法，类似于产生延时单位样本序列的方法。

? 矩形序列

?1RN(n)???00?n?N?1其他

式中的N为序列的长度。矩形序列可以方便地用单位阶跃序列之差来生成：

RN(n)??(n)??(n?N)

? 正弦序列

x(n)?Asin(2?fn/Fs??)

在MATLAB中

n?0:N?1;x?A*sin(2*pi*f*n/Fs?fai); x(n)?ej?n? 复正弦序列

在MATLAB中

34

n?0:N?1;x?exp(j*w*n),

x(n)?a

n? 指数序列

在MATLAB中

n?0:N?1;x?a.^n;

? 随机信号

在区间(0，1)中均匀分布的长度为N的随机信号，可通过如下的MATLAB命令产生： x=rand(1，N)；

同理，使用下面的MATLAB命令，可产生长度为N且具有零均值和单位方差的正态分布的随机信号x[n]

x=randn(1，N); ? 周期序列

如果对所有的n，序列x[n]满足

?[n]?x?[n?kN] x 则x称为周期序列。满足上述关系的最小数N称为基本周期。

用MATLAB把一个周期长的序列{x(n)，O≤n≤N-1}，拓展为有K个周期的序列x，有多种方法可以采用，这里介绍两种。

(1)简单复制法：设x是一个已赋值的行向量，其长度为N=length(x)，把它复制K次，得到的xtide长度为K*N。因此有：

xtide(n)=[x，x,?，x]；n=0：K*N-1 这个方法在K太大时容易数错，不是很好。

(2)求余函数mod法：函数n1=mod(n,N)完成运算nl=(n mod N)。这个算式把大于等于N的n值，减去N的整倍数，使余数n1在0与N-1之间。对于小于等于N的n值，则加以N的整倍数,也使n1在0与N-1之间。把这一运算用到位置向量上，就可以方便地实现有限序列的周期延拓。

先设置位置向量，要复制K个，则新向量的长度应为K*N。假如起始位置为0，则可用下列语句：

nxtide=0：KN-1; ％设置延拓序列的位置向量

xtide=x(mod(nxtide，N)+1); ％确定位置向量各点对应的x值 在最后一条语句中，注意mod(nxtide,N)的值永远在0与N-1之间，而按MATLAB的规则，x的下标为nx=l：N。因此必须把mod函数的结果加1才能与下标规则匹配。 ? sawtooth

功能：产生锯齿波或三角波。 格式：

x=sawtooth(t)

35

x=sawtooth(t，width)

说明：

sawtooth(t)函数类似于sin(t)，它产生周期为2π，幅值从一1到+1的锯齿波，在2π变量的整数倍处，其值为-l，并以1/π的斜率线性上升至+1。

sawtooth(t，width)用于产生三角波，其width(O

x=square(t)

x=square(t，duty) 说明：

square(t)类似于sin(t)，产生周期是2π，幅值为±1的方波，square(t，duty)产生指定周期的方波，其duty用于指定正半周期的比例。

五、MATLAB 的绘图功能

plot是用来划函数x对函数y的二维图，例如要画出 y = sin (x), 0≤x≤2π。plot可以在一个图上划数条曲线，且以不同的符号及颜色来标示曲线，其指令见线上说明help plot。如要在x及y轴及全图加注说明，则可利用指令xlabel, ylabel, title，其指令见线上说明help xlabel, help ylabel, help title。三维图的指令为plot3，其指令见线上说明help plot3。此外二维图及三维图皆可使用指令grid 加上栅格线。

除此之外还有二个指令 text, gtext 可以在图中加上文字用以说明图中的曲线或图形代表什么。text是依据所绘图的座标来放置文字说明，其语法为text(x,y, 'string')，x, y是要放置说明的座标值，string是说明的文字。gtext则是依据滑鼠或上下左右游标键来放置文字说明，其语法为gtext('string')。 我们来看几个例子： >> x=linspace(0,2*pi,30); y=sin(x); z=cos(x); >> plot(x,y,x,z) % 划二条曲线 y=sin(x), z=cos(x)

>> text(2.5,0.7,'sin(x)') % (2.5,0.7)是依据绘图大小的座标值 >> gtext('cos(x)') % 将滑鼠移至适当位置再按滑鼠键

一般的x-y 图在横轴及纵轴皆是以线性尺度来绘图，如果要绘图的数据的 x 或 y 值变化范围太大，就需要改用对数 (log) 尺度来绘图才可得到合理的图。MATLAB 提供三种对数尺度的绘图指令：semilogx,semilogy, loglog， 它们的作用分别是x轴以对数尺度绘图，y 轴以对数尺度绘图，x 和 y 轴以对数尺度绘图。我们来看几个例子， 藉以说明在何种场合须要用对数尺度绘图。 >> y=0:0.1:10; x=10.^y

>> plot(x,y) % 会画出看不出所以然的图

36

>> semilogx(x,y) % 改以对数尺度绘图就清楚多了 >> x=[0 2 5 7 10 12 15 17 20 21];

>> y=[0.1 0.2 0.5 0.6 0.9 1 1.2 1.26 1.22 1.2]; >> plot(x,y) % 先以线性尺度绘图，再分别以三种对数尺度绘 >> semilogx(x,y) >> semilogy(x,y) >> loglog(x,y)

要将数个相关的图画在同一页时，可以用subplot这个指令。其语法为 subplot(m,n,p)，其中 m, n代表绘图成 m x n 个子图，m表示在 y方向有 m 个图， n表示在 x 方向有 n 个图，p 是代表第几个子图。

六、常见数学函数

使用函数须注意几点。首先函数一定出现在计算等式的右边，等式左边是代表这个函数的计算值。此外，一个函数可以被当做另一个函数的引数(argument)。例如：log_x=log(abs(x))其中abs和log皆为内建函数，其意思是先计算abs(x)，所得值再代入log函数。 ? 基本内建函数：

round(x) 将x值进位至最接近的整数 fix(x) 将x值进位至最接近0的整数 floor(x) 将x值进位至最接近-∞的整数

ceil(x) 将x值进位至最接近∞的整数

sign(x) 如果x<0传回值为-1，如果x=0传回值为0，如果x >0传回值为 1 rem(x,y) 传回x/y的余数，例如rem(25,4)的值为1 exp(x) 指数函数

log(x) 以ｅ＝2.718282为底的对数函数，及自然对数 log10(x) 为10底的对数函数

? 三角函数和双曲线函数的使用，和一般数学式相似，其语法也很直接易懂。例如三角函

数有：sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x)。常用到的双曲线函数有：sinh(x), cosh(x), tanh(x), asinh(x), acosh(x), atanh(x)。上述函数用法请参考MATLAB的线上说明或使用手册。

? MATLAB 是以i或j字元来代表虚部，其它的复数相关函数有real, imag, conj, abs,

angle等等，详见线上说明lookfor complex。 ? 如果复数表示为 x?a?bi，则：共轭复数x?a-bi, 复数大小r?a2?b2, 复数

向量的夹角 ??tan?1(b/a) 复数实部a = r cosθ, 复数虚部b = r sinθ, 复数指数表示法 x?rej?

上述各函数对应MATLAB的复数指令为 ：

a=real(x), b=imag(x),x=conj(x), r=abs(x), ?=angle(x), x=r*exp(i*angle(x)) ? 多项式函数：

37

令p(x) 代表一个多项式如下 p(x)?x?4x?7x?10

MATLAB 以一最简便方式代表上述的多项式 p=[1 4 -7 -10]，其中的数值是多项式的各阶项（从高到低）的 各个系数，其实p 也是一个阵列不过是用以代表这个多项式。

为了能直接运用多项式，可以用函数 polyval直接做运算，语法为 polyval(p,x)，其中p 即是代表多项式各阶系数的阵列。因此

>> x=linspace(-1,3); >> p=[1 4 7 -10]; >> v=polyval(p,x);

当二个多项式间要做加减乘除时，加减运算可以直接进行。定义函数conv做乘法运算以及函数deconv做除法运算。当二多项式相乘，在数学上等于二个数组做卷积(convolution)运算（因为我们是以数组来代表一个多项式的各阶系数）， 因此可利用conv函数做乘法运算，其语法为conv(a,b)，其中a, b代表二个多项式的数组。而二多项式相除就相当于反卷积(de-convolution) 运算，因此有 deconv 函数，其语法稍有不同。 [q,r]=deconv(a,b)，其中q,r分别代表整除多项式及余数多项式。

以下就介绍相关范例，来说明二个多项式的加减乘除运算： >> a=[1 2 3 4]; b=[1 4 9 16]; >> c=a+b c = 2 6 12 20 >> d=a-b d = 0 -2 -6 -12 >> e=conv(a,b) e =

1 6 20 50 75 84 64 >> g=e+[0 0 0 c] g =

1 6 20 52 81 96 84 >> [f,r]=deconv(e,b) f = 1 2 3 4 r =

0 0 0 0 0 0 0 % 因为是整除所以余数多项式的各系数皆为零 >> [h,r]=deconv(g,a) f =

1 4 9 18 r =

38

32

0 0 0 0 2 6 12 % 余数多项式为 2*x^2 + 6*x + 12

poly, real，这二个函数的用途是要验算将求解的根展开，求得原多项式。 例如有一个二次方程式的根为 2, 1，则以下式计算原多项式

p(x)?(x?2)(x?1)?x?3x?2

2poly函数就是求出多项式的各阶系数，其语法为 poly(r)，其中 r 是代表根的阵列。而 real则是用来去除因计算时产生的假虚部系数，为何会有此种情形请参考以下的例子。

>> r=[-2 -1];

>> pp=poly(r) % pp=(x+2)(x+1)=x^2+3x+2 pp =

1 3 2

再看另一个多项式 》pp=[1 7 12 9]; % >> r=roots(pp) r = -4.9395

-1.0303 + 0.8721i -1.0303 - 0.8721i

>> pp=poly(r) % 注意因计算的误差会有假虚部产生 pp =

1.0000 7.0000 12.0000 9.0000 + 0.0000i

>> pp=real(pp) % 可用real将假虚部去除，将原多项式还原 pp =

1.0000 7.0000 12.0000 9.0000

? 数据分析函数

最大值max， 最小值min，平均值 mean，一组数据的中位数median，总和值sum，连乘值prod，累积总和值cumsum， 累积连乘值cumprod，排序函数sort。它们的使用方式如下：

max(x) 找出x阵列的最大值

max(x,y) 找出x及y阵列的最大值，会有二个极值分属x及y阵列 [y,i]=max(x) 找出x阵列的最大值以y显示，其在x阵列的位置以i显示 min(x) 找出x阵列的最小值

min(x,y) 找出x及y阵列的最小值，会有二个极值分属x及y阵列 [y,i]=min(x) 找出x阵列的最小值以y显示，其在x阵列的位置以i显示

mean(x) 找出x阵列的平均值 median(x) 找出x阵列的中位数 sum(x) 计算x阵列的总和值

39

prod(x) 计算x阵列的连乘值 cumsum(x) 计算x阵列的累积总和值 cumprod(x) 计算x阵列的累积连乘值 七、关系及逻辑运算 符号 关系的意义 < 小于 <= 小于等于 > 大于 >= 大于等于 == 等于 ~= 不等于 & 逻辑 and | 逻辑 or ~ 逻辑 not

上述的各个运算元须用在二个大小相同的阵列或是矩阵的比较 利用关系及逻辑运算产生信号的例子： >> x=linspace(0,10,100); % 产生数据 >> y=sin(x); % 产生 sine 函数

>> z=(y>=0).*y; % 将 sin(x) 的负值设为零 >> z=z + 0.5*(y<0); % 再将上式的值加上0.5 >> z=(x<8).*z; % 将大于 x=8 以后的值设为零 >> hold on >> plot(x,z)

>> xlabel('x'),ylabel('z=f(x)') >> title('A discontinuous signal') >> hold off

除了上述的运算元之外，尚有以下的逻辑关系函数：xor(x,y), any(x), all(x), isnan(x), isinf(x), finite(x), find(x)，其使用方式详见线上说明。

八、Matlab的程序结构

控制流极其重要，因为它使过去的计算影响将来的运算。MATLAB提供三种决策或控制流结构。它们是：For循环，While循环和If-Else-End结构。由于这些结构经常包含大量的MATLAB命令，故经常出现在M文件中，而不是直接加在MATLAB提示符下。

（1）MATLAB IF-ELSE-END 结构 最简单的If-Else-End结构是：

if expression {commands} end

40

假如有两个选择，If-Else-End结构是： if expression

commands evaluated if True else

commands evaluated if False end

当有三个或更多的选择时，If-Else-End结构采用形式 if expression1

commands evaluated if expression1 is True elseif expression2

commands evaluated if expression2 is True elseif expression3

commands evaluated if expression3 is True elseif expression4

commands evaluated if expression4 is True elseif ?? ． ． ． else

commands evaluated if no other expression is True

end

(2) MATLAB While 循环 While循环的一般形式是：

while expression {commands} end

只要在表达式里的所有元素为真，就执行while和end 语句之间的{commands}。通常，表达式的求值给出一个标量值，但数组值也同样有效。

（3）MATLAB For 循环

For循环允许一组命令以固定的和预定的次数重复。For循环的一般形式是：

for x = array

{commands}

end

在for和end语句之间的{commands}按数组中的每一列执行一次。在每一次迭代中，x

被指定为数组的下一列，即在第n次循环中，x=array(:, n)。 For循环的其它重要方面是：

41

1. For循环不能用For循环内重新赋值循环变量n来终止。 2.在For循环内接受任何有效的MATLAB数组。 3. For循环可按需要嵌套。

4. 当有一个等效的数组方法来解给定的问题时，应避免用For循环。两种方法得出同样的结果，而后者执行更快，更直观，要求较少的输入。

5. 为了得到最大的速度，在For循环(While循环)被执行之前，应预先分配数组。 九、滤波器设计 1、滤波器结构 IIR滤波器结构 IIR滤波器可写成

H(z)=B(z)A(z)?b0+b1z+...+bMz1+a1z+...+aNz-1-1-M-N

即可用滤波器系数bi和ak表示。当aN?0时，滤波器的阶为N。IIR滤波器的差分方程表示为

MNy(n)=?bmx(n-m)-?amy(n-m)

m=0m=1实现IIR滤波器可采用三种不同的结构；直接形式、级联形式和并联形式。 1．直接形式

滤波器差分方程直接用延迟线、乘法器和加法器实现，例如M=4，N=4时，得到如图3所示的结构，其差分方程为

y(n)=b0x(n）＋...+b4x(n-4)-a1y(n-1)-...-a4y(n-4)

这种结构称为直接形式1。

将这种结构进行简化，可得到直接形式2，如图4所示。

在MATLAB中，可利用DSP工具箱函数filter实现IIR的直接形式。

图3 IIR滤波器直接形式1结构(N=4) 图4 IIR滤波器直接形式2结构(N=4)

2．级联形式

下面先假设N为偶数，则IIR可简化

42

H(z)=b0+b1z+...+bNz1+a1z+...+aNzK-1-1-1-1-N-N

=b0?k=11+Bk,1z+Bk,2z1+Ak,1z+Ak,2z-2-2

其中K：N／2，每个二阶形式都可用直接形式2得以实现，例如当N=4时，其级联结构如图1．5所示。

图5 IIR滤波器的级联结构(N=4)

在MATLAB中，给定直接形式滤波器的系数{bk，ak }，可计算出b0，Bk,i，Ak,i这可由扩展函数dir2cas实现。然后利用扩展函数casfiltr实现滤波器的级联形式。

当然，在给定级联形式时，也可得到滤波器的直接形式，这由扩展函数cas2dir实现。 3．并联形式

在级联形式的基础上，利用部分分式展开，可进一步得到(设M≥N)

H(z)=b0+b1z+...+bMz1+a1z+...+aNz-1-1-N-N

=?+b?z-1+...+b?z-Nb01N-11+a1z+...+aNzK-1-NM-N??Ck=0kz-k

??1?Ak?1Bk0?Bk,1zk,1?1?2M-Nz?1?Ak,2z??Ck=0kz-k

其中K=N／2，每一部分均为二阶，它们之间是并联关系。同样以M=4，N=4为例，画出其结构如图6所示。

43

图6 IIR滤波器的并联结构

在MATLAB下，我们设计了将滤波器直接形式转化为并联形式的扩展函数dir2par，其中用到了复共轭对比较函数cplxcomp。这样我们就可利用parfiltr函数实现滤波器的并联形式。最后还设计出了将并联形式转换成直接形式的扩展函数par2dir。

H(z)=b0+b1z+...+bM-1z-1-(M-1)

?bnh(n)=??00?n?M-1其他

y(n)=b0x(n)+b1x(n-1)+...+bM-1x(n?M?1)

图7 FIR滤波器的直接结构形式

H(z)=b0+b1z+...+bM-1z-1-(M-1)

-1-(M-1) =b0+b1z+...+bM-1z

其中K=M/2.以M=7为例，其结构如图8所示。

图8 FIR滤波器级联形式(M=7)

在MATLAB中，可采用dir2cas和cas2dir函数，在直接形式和级联形式之间转换，滤波

44

器实现可通过casfiltr函数实现。

3．线性相位形式线性相位条件

?H(j?)?????,??????

其中?=O或?=±?／2，a为常数。对因果FIR滤波器，设冲激响应处于[O，M-1]区间上，则上述条件就表示下列对称性

(对称冲激响应): h(n)?h(M-1-n),??0,0?n?M-1

(反对称冲激响应): h(n)?-h(M-1-n),????/2,0?n?M-1

以M=7或M=6为例，对称型的结构如图1．9所示。在MATLAB中，实现FIR线性相位形式的函数等同于直接形式，即采用filter函数。

图9 FIR滤波器线性相位形式(M=7或M=6)

4．频率取样形式冲激响应h(n)的M点DFT为H(k)(O≤k≤M一1)，则有

H(z)=Z?h(n)?=Z??IDFT?H(k)???

利用内插公式可得:

?1-z-M?H(Z)=???M?M-1?1-Wk=0H(k)-kM

-1z这表明在这种结构中采用了离散傅里叶变换H(k)，而不是冲激响应h(n)。这种结构分成两部分：第一部分由?1?z?M?/M(构成，第二部分由M个H(k)/1-WMz并联而成。

-k-1在MATLAB中，给定h(n)或H(k)，则可借助于扩展函数dir2fs得到频率取样形式， 2、滤波器分析与实现函数 1．abs

功能：求绝对值(幅值)。 格式： y=abs(x) 说明：

y=abs(x)用于计算x的绝对值。当x为复数时，得到的是复数模(幅值)，即

abs(x)?45

?Re(x)?2??Im(x)? 2

当x为字符串时，abs(x)得到字符串的各个字符的ASCII码，例如x=’123 ’，则abs(x)得到49 50 51。 2．angle

功能：求相角。 格式：

p=angle(h) ， 说明：

p=angle(h)用于求取复矢量或复矩阵h的相角(以弧度为单位)，相角介于-π和π之间。例如对复数h可用两种方法表示

h=x+i y=mej p

则幅值m和相角p可由x+i y格式求出

m=abs(h) p=angle(h)

当然，由m和p也可求取x+i y格式 h=m．*exp(i*p) x=real(h) y=imag(h) 3．conv

功能：求卷积。 格式： C=CONV(a，b)

说明：

conv(a，b)用于求取矢量a和b的卷积，即

N-1c(n+1)=?a(k+1)b(n-k)

k=0其中N为矢量a和b的最大长度。例如，当a=[1 2 3]，b：[4 5 6]时，则 c=conv(a，b) c=

4 13 28 27 18 4.filter

功能：利用递归滤波器(IIR)或非递归滤波器(FIR)对数据进行滤波。 格式：

y=filter(b，a，x) [y，zf]=filter(b，a，x) y=filter(b，a，x，zi)

46

说明：

filter采用数字滤波器对数据进行滤波，其实现采用移位直接型Ⅱ结构，因而适用于IIR和FIR滤波器.滤波器的z域表示为

Y(z)=b(1)+b(2)z+...+b(nb+1)z1+a(2)z+...+a(na+1)z-1-na-1-nbX(z)

这里a(1)=1，如果输入的a(1)≠1，则MATLAB将自动进行归一化系数的操作；如果a(1)=O，则给出出错信息。

y=filter(b，a，x)利用给定的矢量a和b，对x中的数据进行滤波，结果放入y矢量，长度取max(nb，na)。

[y，zf]=filter(b，a，x)除得到结果矢量y外，还得到x的最终状态矢量zf。 y=filter(b，a，x，zi)可在zi中指定x的初始状态。 5. filtic

功能：为filter函数选择初始条件。 格式：

z=filtic(b，a，y，x) z=filtic(b，a，y) 说明：

filtic函数在给定过去输出y和输入x的前提下，为移位直接形式Ⅱ滤波器实现中的延迟找到初始条件z。矢量x和y分别表示过去的输入和输出

x={x(-1)，x(-2)，?，x(-nb)，?} y={y(-1)，y(-2)，?，y(-na)，?}

其中nb=length(b)-l(分子阶数)，na=length(a)-1(分母阶数)。矢量x和y的长度应分别等同于nb和na，如果长度不够，则补O；如果长度超出，则截断处理，使x和y的长度正好为na和nb。输出矢量z的长度为max(na，nb)。

z=filtic(b，a，y，x)可得到给定输入x和输出y时的初始状态；

z=filtic(b，a，y)则假定输入x=O。 6.freqs

功能：模拟滤波器的频率响应。 格式:

h=freqs(b，a，w) l [h，w]=freqs(b，a) I [h，w]=freqs(b，a，n) 说明：

freqs用于计算由矢量a和b构成的模拟滤波器 :

H(s)=47

b(1)snanb+b(2)s?nb?1?+...+b(nb+1)s+a(2)sna?1+...+a(na+1)

本文来源：https://www.bwwdw.com/article/rba8.html