数字信号处理实验指导书通信

1、数字信号处理实验徐俊2016年8月目 录实验一离散时间系统及系统响应2实验二离散傅立叶变换及其特性验证8实验三时域采样与频域采样17实验四冲激响应不变法IIR数字滤波器设计24 实验一 离散时间系统及系统响应一、实验目的1、掌握求解离散时间系统脉冲响应和阶跃响应的方法；2、掌握用线性卷积求解离散时间系统响应的基本方法。二、实验原理与设计方法1、用impz和dstep函数求解离散系统的单位脉冲响应和阶跃响应【例1-1】已知某因果系统的差分方程为yn+0.5yn-1=xn+2x(n-2)系统为零状态，求系统的脉冲响应和阶跃响应。解：该系统是一个2阶系统，列出bm和ak系数为a0=1，a1=0.5，

2、a2=0，b0=1，b1=0，b2=2MALAB程序如下（取16点作图）：a=1,0.5,0;b=1,0,2;n=16;hn=impz(b,a,n); %脉冲响应gn=dstep(b,a,n); %阶跃响应subplot(1,2,1),stem(hn,'k');title('系统的单位脉冲响应');ylabel('h(n)');xlabel('n');axis(0,n,1.1*min(hn),1.1*max(hn);subplot(1,2,2),stem(gn,'k');title('系统的单位阶跃响应&#

3、39;);ylabel('g(n)');xlabel('n');axis(0,n,1.1*min(gn),1.1*max(gn);结果如下图所示：2、用conv函数进行卷积计算求系统响应【例1-2】某离散时间系统的脉冲响应为hb(n)=(n)+2.5(n-1)+2.5(n-2)+(n-3)激励信号为xt=Ae-nTsin0nT 0n<50设A=444.128，=502，0=502。试求该系统在输入信号激励下的响应。解：MATLAB程序如下：n=1:50; %定义序列的长度是50hb=zeros(1,50); %注意：MATLAB中数组下标从1开始hb(1)

4、=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;close all;subplot(3,1,1);stem(hb);title('系统hn');m=1:50;T=0.001; %定义序列的长度和采样率A=444.128;a=50*sqrt(2.0)*pi; %设置信号有关的参数w0=50*sqrt(2.0)*pi;x=A*exp(-a*m*T).*sin(w0*m*T);subplot(3,1,2);stem(x);title('输入信号xn');y=conv(x,hb);subplot(3,1,3);stem(y);title('输出信号

5、yn');结果如下图所示：3、用filter函数求系统响应线性常系数差分方程是描述离散时间LTI系统的另一个时域模型，即系统的输入信号xn输出信号yn关系可以用下面的差分方程来表达 为了求得系统响应信号的显式表达式（Explicit expression），必须对差分方程求解。在MATLAB中，是用差分方程的系数来描述一个离散时间LTI系统的。例如，一个LTI离散时间系统的差分方程为 yn+yn-1-6yn-2=xnMATLAB则用两个系数向量num = 1和den = 1 1 -6来描述该系统，其中num和den分别表示系统差分方程右边和左边的系数，按照差分运算的递减排列。在用num

6、和den定义了一个离散时间LTI系统之后，可以利用MATLAB来求解系统响应。求解离散时间系统的响应的一个非常有用的MATLAB函数就是filter()。它可以用来求解系统的在输入信号xn作用下的零状态响应，也可以求解系统在这个输入信号作用下的完全响应。其用法描述如下：y = filter(num, den, x)，求系统在输入x作用下的零状态响应y，x与y的长度相同。y = filter(num, den, x, ic)，求系统在输入x作用下的完全响应y，x与y的长度相同。其中ic为系统的初始条件（Initial Condition），即ic = y-1, y-2, y-3, ., y-N。

7、【例1-3】给定一个离散时间LTI系统，它的差分方程为yn+0.1yn-1-0.06yn-2=xn若输入信号为xn = 0.5nun，系统的初始条件为ic = 0, 1，编写程序，绘制输入信号xn、系统的零状态响应yzsn和系统的完全响应信号yn的波形图。MATLAB程序如下：clear,close all,num = 1; den = 1 .1 -0.06;n = 0:20; x = 0.5.n.*ones(1,21); ic = 0 .9;yzs = filter(num,den,x); y = filter(num,den,x,ic); subplot(2,2,1)stem(n,x,&#

8、39;.'), title('The input sequence xn'), axis(0,20,-0.5,1.5)subplot(2,2,2)stem(n,yzs,'.'), title('The zero-state response yzsn'), axis(0,20,-0.5,1.5), xlabel('Time index n')subplot(2,2,3)stem(n,y,'.'), axis(0,20,-0.5,1.5), title('The total solution yn&#

9、39;)xlabel('Time index n')结果如下图所示：三、实验内容1、分别用impz和dstep函数求解下面离散时间系统的脉冲响应和阶跃响应。（1）系统的差分方程为y(n)=0.8y(n-1)-0.64y(n-2)+0.866x(n)（2）系统的系统函数为Hz=1-0.5z-11-z-1+z-22、利用第1题求得的系统的脉冲响应求解系统在激励x(n)=u(n-3)下的响应。3、利用filter函数求解第1题系统在激励x(n)=u(n-3)下的响应。四、实验参考MATLAB基础操作1、 矩阵输入Matlab有很强的数值矩阵处理能力。实际上，矩阵和矢量这两个词经常换用

10、。矩阵是以实数或者复数为元素的长方形矢量。在输入矩阵时，应遵循下列规则：l 整个矩阵用中括号括起来；l 同一行的数据用空格或逗号隔开；l 不同行用分号隔开。在Matlab中，列矢量可被当作只有一列的矩阵；行矢量也可被当作是只有一行的矩阵；标量可被当作只有一列的矩阵。在MATLAB中，用一个列向量来表示一个有限长序列，由于一个列向量并不包含位置信息，因此需要用表示位置的n和表示量值的x两个向量来表示任意一个序列，如：【例】n=-3,-2,-1,0,1,2,3,4;x=2,1,-1,0,1,4,3,7;plot(n,x);运行结果：如果不对向量的位置进行定义，则MATLAB默认该序列的起始位置为n

11、=0。由于内存有限，MATLAB不能表示一个无限序列。2、 利用函数输入矩阵例如：zeros：生成一个元素全部为0的矩阵；ones：生成一个元素全部为1的矩阵；eye：生成一个单位矩阵。A=zeros(2,3) ：生成一个2行3列的全0矩阵；A= 0 0 0 0 0 03、 利用plot，stem进行波形绘制l 绘图命令plot绘制x-y坐标图;l 在绘图过程中，经常要把几个图形在同一个图形窗口中表现出来，而不是简单地叠加这就用到函数subplot其调用格式如下：subplot(m,n,p)。subplot函数把一个图形窗口分割成m×n个子区域，用户可以通过参数p调用个各子绘图区域进

12、行操作子绘图区域的编号为按行从左至右编号l stem绘制离散序列图，常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。实验二 离散傅立叶变换及其特性验证一、实验目的1、掌握离散时间傅立叶变换（DTFT）的计算方法和编程技术。2、掌握离散傅立叶变换（DFT）的计算方法和编程技术。3、掌握用FFT对信号作频谱分析。二、实验原理与设计方法1、离散时间傅立叶变换如果序列x(n)满足绝对可和的条件，即，则其离散时间傅立叶变换定义为： （1）如果x(n)是无限长的，则不能直接用MATLAB由x(n)计算X(ejw)，但可以用它来估计X(

13、ejw)表达式在0,频率区间的值并绘制它的幅频和相频（或实部和虚部）曲线。如果x(n)是有限长的，则可以用MATLAB对任意频率w处的X(ejw)进行数值计算。如果要在0,间按等间隔频点估计X(ejw)，则（1）式可以用矩阵向量相乘的运算来实现。假设序列x(n)在(即不一定在0, N-1)有N个样本，要估计下列各点上的X(ejw)：它们是0,之间的(M+1)个等间隔频点，则（1）式可写成： （2）将x(nl)和X(ejwk)分别排列成向量x和X，则有： X=Wx （3）其中W是一个(M+1)×N维矩阵：将k和n排成列向量，则在MATLAB中，把序列和下标排成行向量，对（3）式取转置得

14、：其中nTk是一个N×(M+1)维矩阵。用MATLAB实现如下：k=0:M; n=n1:n2;X=x*(exp(-j*pi/M).(n*k);2、离散傅立叶变换一个有限长序列的离散傅立叶变换对定义为： （4） （5）以列向量x和X形式排列x(n)和X(k)，则式（4）、（5）可写成：X=WNx其中矩阵WN由下式给出：可由下面的MATLAB函数dft和idft实现离散傅立叶变换运算。function Xk = dft(xn,N)% Computes Discrete Fourier Transform% -% Xk = dft(xn,N)% Xk = DFT coeff. array

15、over 0 <= k <= N-1% xn = N-point finite-duration sequence% N = Length of DFT%n = 0:1:N-1; % row vector for nk = 0:1:N-1; % row vecor for kWN = exp(-j*2*pi/N); % Wn factornk = n'*k; % creates a N by N matrix of nk valuesWNnk = WN . nk; % DFT matrixXk = xn * WNnk; % row vector for DFT coeffi

16、cientsfunction xn = idft(Xk,N)% Computes Inverse Discrete Transform% -% xn = idft(Xk,N)% xn = N-point sequence over 0 <= n <= N-1% Xk = DFT coeff. array over 0 <= k <= N-1% N = length of DFT%n = 0:1:N-1; % row vector for nk = 0:1:N-1; % row vecor for kWN = exp(-j*2*pi/N); % Wn factornk =

17、 n'*k; % creates a N by N matrix of nk valuesWNnk = WN . (-nk); % IDFT matrixxn = (Xk * WNnk)/N; % row vector for IDFT values【例2-1】对以下基本序列进行16点DFT。要求画出原序列、DFT变换后的频谱。x1n=e-j8n,x2n=cos(8n)解：N=16;%产生序列x1(n)，计算DFTx1(n)n=0:N-1;x1n=exp(j*pi*n/8);%产生x1(n)X1k=dft(x1n,N);%计算N点DFTx1(n)%产生序列x2(n)，计算DFTx2(n

18、)x2n=cos(pi*n/8);%产生x2(n)X2k=dft(x2n,N);%计算N点DFTx2(n)%绘图subplot(2,2,1);stem(n,x1n,'.');axis(0,20,-1,1);title('序列x1(n)');xlabel('n');ylabel('x1(n)');subplot(2,2,2);stem(n,x2n,'.');axis(0,20,-1,1);title('序列x2(n)');xlabel('n');ylabel('x2(n)

19、9;);subplot(2,2,3);stem(n,abs(X1k),'.');axis(0,20,0,20);title('16点DFTx1(n)');xlabel('k');ylabel('|X1(K)|');subplot(2,2,4);stem(n,abs(X2k),'.');axis(0,20,0,8);title('16点DFTx2(n)');xlabel('k');ylabel('|X2(K)|');结果如下图所示：3、按时间抽取基-2FFT算法对于有限

20、长离散数字信号xn，0 £ n £ N-1，其离散谱xk可以由离散付氏变换（DFT）求得。DFT的定义为由DFT的定义可以看出，在xn为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。MATLAB中的函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R<N，则x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算

21、R点的谱序列的R点IDFT值；而ifft(X,N)同fft(x,N)的情况。【例2-2】已知x(n)主值序列x=1 2 1 0 0 0 0 1 2 1，周期N=10，用MATLAB计算周期序列x(n)的DFS并用IDFS恢复x(n)的主值序列。解：MATLAB程序如下：x=1 2 1 0 0 0 0 1 2 1;N=10;n=0:9;k=0:9;X=fft(x,N);x1=ifft(X,N);subplot(3,2,1);stem(n,x);title('x');%xlabel('n');ylabel('x');subplot(3,2,2);st

22、em(k,abs(X);title('abs(X)');%xlabel('k');ylabel('abs(X)');subplot(3,2,3);stem(k,real(X);title('real(X)');%xlabel('k');ylabel('real(X)');subplot(3,2,4);stem(k,imag(X);title('imag(X)');%xlabel('k');ylabel('imag(X)');subplot(3,2,5)

23、;stem(n,real(x1);title('real(x1)');%xlabel('n');ylabel('real(x1)');subplot(3,2,6);stem(n,imag(x1);title('imag(x1)');%xlabel('n');ylabel('imag(x1)');结果如下所示：【例2-3】用MATLAB求序列x(n)的离散傅里叶变换（N=11），xn=0.85n,0n100, 其他并计算通过补零后（N=41）及增加x(n)的有效长度（N=41）相应的离散傅里叶变换，并

24、比较三种频谱的分辨率。增加x(n)的有效长度（N=41）如下所示：xn=0.85n,0n400, 其他解：MATLAB程序如下：N1=11;n1=0:N1-1;k1=0:N1-1;x1=(0.85).n1;X1=fft(x1,N1);dw1=(2*pi)/N1;N2=41;n2=0:N2-1;k2=0:N2-1;x2=(0.85).n1,zeros(1,N2-N1);X2=fft(x2,N2);dw=(2*pi)/N2;x3=(0.85).n2;X3=fft(x3,N2);figure(1);subplot(2,1,1);stem(n1,abs(X1);title('N=11,abs(

25、X1)');xlabel('k');ylabel('abs(X1)');subplot(2,1,2);stem(n1,angle(X1);title('N=11,angle(X1)');xlabel('k');ylabel('angle(X1)');figure(2);subplot(2,1,1);stem(n2,abs(X2);title('N=41,abs(X2),补30个零');xlabel('k');ylabel('abs(X2)');subplot(

26、2,1,2);stem(n2,angle(X2);title('N=41,abs(X2),补30个零');xlabel('k');ylabel('angle(X2)');figure(3);subplot(2,1,1);stem(n2,abs(X3);title('N=41,abs(X2),增加有效长度');xlabel('k');ylabel('abs(X1)');subplot(2,1,2);stem(n2,angle(X3);title('N=41,angle(X2),增加有效长度&#

27、39;);xlabel('k');ylabel('angle(X3)');结果如下所示：三、实验内容：1、已知有限长序列x(n)=7,6,5,4,3,2，求x(n)的DFT和IDFT。要求：（1）画出序列DFT对应的|X(k)|和argX(k)的图形。（2） 画出原信号与IDFTX(k)的图形，并进行比较。2、分别计算两个实数序列的128点FFT。3、利用DFT的方式计算下面两序列的线性卷积：gn=3, 4, -2, 0, 1, -4,hn=1, -3, 0, 4, -2, 34、用MATLAB求序列x(n)的离散傅里叶变换（N=12），xn=0.6n+0.7n

28、,0n110, 其他并计算通过补零后（N=32）及增加x(n)的有效长度（N=32）相应的离散傅里叶变换，并比较三种频谱的分辨率。实验三 时域采样与频域采样一、 实验目的1、 掌握模拟信号采样前后频谱的变化，以及如何选择采样频率才能使采样后的信号不丢失信息；2、 掌握频率域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。二、 实验原理与设计方法1、 时域采样在模拟信号数字化的过程中，若采样速率不满足采样定理的要求，采样之后信号的频谱将发生混叠现象，则原信号将无法恢复。因此信号在进行A/D转换时，采样速率的确定十分重要。要使得有限带宽信号xa(t)被采样后能不失真

29、的恢复出原始信号，采样信号的采样频率fs必须符合奈圭斯特定理。假设xa(t)的最高频率为fm，那么采样频率应满足fs2fm。满足奈圭斯特采样定理的采样信号xa(t)，只要经过一个理想的低通滤波器，将原始信号有限带宽以外的部分滤去，保留基带信号，就可重建原模拟信号。信号重建一般采用内插公式完成，即用时域信号与理想低通滤波器的单位脉冲响应进行卷积积分来求解。理想低通滤波器的频率特性为一矩形，其单位脉冲响应为ht=12-H(j)ejtd=sin(t/T)t/T采样信号xa(t)通过滤波器输出，其结果为xa(t)和h(t)的卷积积分，即yat=xat=xat*ht=-xaht-d=n=-xa(nT)s

30、in(t-nT)/T(t-nT)/T上面的公式成为内插公式。由上式可知，模拟信号可以由其采样值及内插函数重构。MATLAB提供了sinc函数，可以很方便的实现内插重构。【例3-1】已知某模拟信号ft=sin2f0t+13sin6f0t，其中f0=1Hz，取最高有限带宽频率fm=5f0。（1） 分别显示原连续信号波形和fs>2fm、fs=2fm、fs<2fm三种情况下采样信号的波形以及采样信号的频谱。（2） 用内插公式实现原始信号的恢复。解：MATLAB程序如下：dt=0.01;f0=1;T0=1/f0;%输入基波频率、周期fm=5*f0;Tm=1/fm;%最高频率为基波频率的5倍t

31、=-2:dt:2;f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);%建立原连续信号subplot(4,1,1),plot(t,f,'k');axis(min(t) max(t) 1.1*min(f) 1.1*max(f);title('原连续信号和采样信号');for i=1:3 fs=i*fm;Ts=1/fs;%确定采样频率和周期 n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n);%生成采样信号 subplot(4,1,i+1);stem(n,f,'filled','k

32、'); axis(min(t) max(t) 1.1*min(f) 1.1*max(f);endfigureN=length(t);%时间轴上采样点数f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);wm=2*pi*fm;k=0:N-1;w1=k*wm/N;%频率轴上采样点F1=f*exp(-j*t'*w1)*dt;%对原始信号进行傅里叶变换subplot(4,1,1),plot(w1/(2*pi),abs(F1),'k');axis(0 max(4*fm) 1.1*min(abs(F1) 1.1*max(abs(F1);title(

33、9;原连续信号和采样信号的频谱');for i=1:3 if i<=2 c=0,else c=1,end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); N=length(n); wm=2*pi*fs; k=0:N-1; w=k*wm/N; F=f*exp(-j*n'*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F),'k'); axis(0 max(4*fm) 1.1*min(abs(F) 1.1*max(abs(F);end

34、t=0:dt:3*T0;x= sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);figure;subplot(4,1,1),plot(t,x,k);axis(min(t) max(t) 1.1*min(x) 1.1*max(x);title('用内插公式重建采样信号');for i=1:3fs=i*fm;Ts=1/fs;n=0:(3*T0)/Ts;t1=0:Ts:3*T0;x1= sin(2*pi*n*f0/fs)+1/3*sin(6*pi*n*f0/fs);%生成采样信号T_N=ones(length(n),1)*t1-n*Ts*ones(1.length(

35、t1);xa=x1*sinc(fs*pi*T_N);%内插subplot(4,1,i+1),plot(t1,xa,k);axis(min(t1) max(t1) 1.1*min(xa) 1.1*max(xa);end运行结果：2、频域采样在单位圆上对任意序列x(n)的z变换X(z)进行N点的等间隔采样可得Xk=X(z)z=ej2Nnk k=0,1,N-1上式实现了序列在频域的采样。那么，由频域采样获得的频谱序列能否和连续信号时域采样一样，在一定的条件下能不失真的恢复出原离散序列呢？频域采样定理给出了答案。频域采样定理由下列公式表述xn=r=-x(n+rN)上式表明，对一个频谱采样后经IDTF生

36、成的周期序列xn是原非周期序列x(n)的额周期延拓序列，其时域周期等于频域采样点数N。假设有限长序列x(n)的长度是M，其频域采样点数为N，则原时域离散信号能够不失真的有频域采样恢复的条件如下。如果x(n)为无限长序列，则必然造成混叠现象，无法无失真恢复。如果x(n)为有限长序列，且频域采样点数N小于序列长度M，则x(n)以N为周期进行延拓也将产生混叠，从xn中不能无失真的恢复出x(n)。如果x(n)为有限长序列，且频域采样点数N大于或等于序列长度M，则x(n)以N为周期进行延拓也将产生混叠，从xn中能够无失真的恢复出x(n)，即xNn=xnRNn=r=-xn+rNRNn=x(n)【例3-2】

37、已知序列x(n)=1,1,1;n=0,1,2，对其频谱进行采样，分别取N=2、3、5，观察频域采样造成的混叠现象。解：MATLAB程序如下：Ts=1;N1=2;D1=2*pi/(Ts*N1);kn1=floor(-(N1-1)/2:-1/2);kp1=floor(0:(N1-1)/2);w1=kp1,kn1*D1;X1=1+1*exp(-j*w1)+exp(-j*2*w1);%2点频谱采样结果n=0:N1-1;x1=abs(ifft(X1,N1);%由2点频谱采样结果恢复的时域信号figuresubplot(1,3,1);stem(n*Ts,x1,'filled','k

38、');N2=3;D2=2*pi/(Ts*N2);kn2=floor(-(N2-1)/2:-1/2);kp2=floor(0:(N2-1)/2);w2=kp2,kn2*D2;X2=1+1*exp(-j*w2)+exp(-j*2*w2);%3点频谱采样结果n=0:N2-1;x2=abs(ifft(X2,N2);%由3点频谱采样结果恢复的时域信号subplot(1,3,2);stem(n*Ts,x2,'filled','k');N3=5;D3=2*pi/(Ts*N3);kn3=floor(-(N3-1)/2:-1/2);kp3=floor(0:(N3-1)/2

39、);w3=kp3,kn3*D3;X3=1+1*exp(-j*w3)+exp(-j*2*w3);%5点频谱采样结果n=0:N3-1;x3=abs(ifft(X3,N3);%由5点频谱采样结果恢复的时域信号subplot(1,3,3);stem(n*Ts,x3,'filled','k');运行结果：三、 实验内容1、已知连续信号xt=sinc(t)，取最高有限带宽频率fm=1Hz。（1）分别显示原信号波形和fs=fm、fs=2fm、fs=3fm三种情况下采样信号的波形。（2）求解原连续信号和采样信号所对应的频谱。（3）用内插公式重建信号。2、已知一个时间序列x(n)

40、=2,4,6,4,2;n=0,1,2,3,4，先求其频谱，再分别取频域采样点数N=3、5、10，用IFFT计算并求出其时间序列x(n)，用图形显示各时间序列。观察时域序列是否存在混叠，有何规律？实验四 冲激响应不变法IIR数字滤波器设计一、实验目的1、掌握构成一个频率响应与给定的滤波特性相接近的模拟滤波器的设计原理。2、掌握用冲激响应不变法设计IIR数字滤波器的基本原理和算法。3、了解数字滤波器和模拟滤波器的频率响应特性，掌握相应的计算方法，分析用冲激响应不变法获得的数字滤波器频率响应特性中出现的混叠现象。二、实验原理与设计方法、冲激响应不变法设计IIR数字滤波器的基本原理和算法采用冲激响应不

41、变法设计数字滤波器，就是使其单位样值响应与相应的模拟滤波器的冲激响应在抽样点处的量值相等，即 (1)其中T为抽样周期。因此用冲激响应不变法设计IIR数字滤波器的基本步骤，就是首先根据设计要求确定相应的模拟滤波器的传递函数，经Laplace反变换求出冲激响应后，对它进行抽样得到的等于数字滤波器的单位样值响应，再经z变换所得就是数字滤波器的传递函数。如果模拟滤波器的传递函数的N个极点都是单极点，则可以将写成部分分式展开的形式 (2)那么，经Laplace反变换求出的模拟滤波器的冲激响应为相对应的数字滤波器的单位抽样响应为对上式作z变换，得 (3)由上面的推导可见，只要模拟滤波器的传递函数的N个极点

42、都是单极点，当已经求出各个极点值和部分分式的系数后，则可以从模拟滤波器的传递函数的表达式（2）直接得到数字滤波器的传递函数的表达式（3）。冲激响应不变法设计数字滤波器的步骤：（1） 设数字滤波器的技术指标为wp、ws、Rp、Rs，选择采用周期T，将wp、ws、Rp、Rs转换为模拟滤波器参数p、s、Rp、Rs。其中p=wp/T，s= ws/T。（2） 选择一种模拟低通滤波器原型（巴特沃斯、切比雪夫）设计模拟低通滤波器。（3） 用冲激响应不变法转换为数字滤波器。、Butterworth模拟滤波器的设计方法（1）Butterworh原型MATLAB提供了函数z,p,k=buttap(N)用来设计N阶

43、归一化（即c=1）的Butterworth模拟低通滤波器；其中N为巴特沃斯低通滤波器的阶数，输出值z、p、k分别对应零点数组、极点数组及增益数组。如果c1，则需去归一化处理，极点数组变为cp，零点数组变为cz，分子也要乘以cN（k=cNk）（2）按给定技术指标设计Butterworth模拟低通滤波器N,Wn=buttord(Wp,Ws,Rp,Rs,s)函数实现指定参数的Butterworth模拟低通滤波器的设计。其中输入值Wp、Ws为通带及阻带截止频率，单位为rad/s，Rp,Rs为通带最大衰减及阻带最小衰减，单位为dB，s说明设计的是模拟滤波器。输出值N为滤波器阶数，Wn为滤波器的截止频率。

44、【例4-1】设计巴特沃斯模拟低通滤波器：通带截止频率fp=10Hz，通带最大衰减Rp=2dB，阻带截止频率fs=20Hz，阻带最小衰减Rs=12dB。要求：（1） 求滤波器的阶次N及截止频率Wn；（2） 画出系统的幅频特性。解：fp=10;fs=20;Rp=2;Rs=12;Wp=2*pi*fp;Ws=2*pi*fs;N,Wn=buttord(Wp,Ws,Rp,Rs,'s');z0,p0,k0=buttap(N); %计算归一化参数p=Wn*p0; %计算极点z=Wn*z0; %计算零点k=WnN*k0;%计算增益b=real(poly(z);b=k*b;%计算零极点所对应的系统

45、函数的分子的系数a=real(poly(p);%计算零极点所对应的系统函数的分子母的系数H,w=freqs(b,a);plot(w/(2*pi),20*log10(abs(H)/max(abs(H); %画系统的幅频特性图title('20*log10(abs(H)/max(abs(H)');xlabel('w/Hz');ylabel('20*log10(abs(H)');grid on;axis(0 60 -30 0);运行结果：3、切比雪夫型（Chebyshev）模拟滤波器的设计方法（1）切比雪夫型原型MATLAB提供了函数z,p,k= Ch

46、eb1ap (N,Rp)用来设计N阶归一化（即c=1）的切比雪夫型模拟低通滤波器；其中N为切比雪夫型低通滤波器的阶数，Rp为通带波动值，输出值z、p、k分别对应零点数组、极点数组及增益数组。如果c1，则需去归一化处理，极点数组变为cp，零点数组变为cz，分子也要乘以cN（k=cN*k）（2）按给定技术指标设计切比雪夫型模拟低通滤波器N,Wn= Cheb1 ord(Wp,Ws,Rp,Rs,s)函数实现指定参数的切比雪夫型模拟低通滤波器的设计。其中输入值Wp、Ws为通带及阻带截止频率，单位为rad/s，Rp,Rs为通带最大衰减及阻带最小衰减，单位为dB，s说明设计的是模拟滤波器。输出值N为滤波器阶

47、数，Wn为滤波器的截止频率。【例4-2】设计一个设计切比雪夫型模拟低通滤波器：通带截止频率fp=10Hz，通带波纹Rp=2dB，阻带截止频率fs=20Hz，阻带最小衰减Rs=12dB。要求：（1） 求滤波器的阶次N及截止频率Wn；（2） 画出系统的幅频特性。解：fp=10;fs=20;Rp=2;Rs=12;Wp=2*pi*fp;Ws=2*pi*fs;N,Wn=cheb1ord (Wp,Ws,Rp,Rs,'s');z0,p0,k0=cheb1ap(N,Rp); %计算归一化参数p=Wn*p0; %计算极点z=Wn*z0; %计算零点k=WnN*k0;%计算增益b=real(poly(z);b=k*b;%计算零极点所对应的系统函数的分子的系数a=real(poly(p);%计算零极点所对应的系统函数的分子母的系数H,w=freqs(b,a);plot(w/(2*pi),20*log10(abs(H)/max(abs(H); %画系统的幅频特性图title('20*log10(abs(H)/max(abs(H)');xlabel('w/Hz');ylabel('20*log10(abs(H)');g