1、2018/10/18,1,MATLAB在通信中的应用,一、 MATLAB基础知识 二、 MATLAB与信号处理 三、 MATLAB与通信仿真 四、 Simulink应用,2018/10/18,2,第二部分 信号处理MATLAB仿真,信号及其运算表示,线性时不变系统的时域分析,信号与系统的频谱分析,数字滤波器设计,2018/10/18,3,连续信号 严格来说，MATLAB并不能处理连续信号，而是用等时间间隔点的样值来近似表示连续信号。当取样时间间隔足够小时，这些离散的样值就能较好地近似连续信号。在MATLAB可视化绘图中 对于以t为自变量的连续信号，在绘图时统一用plot函数； 而对n为自变量的
2、离散序列，在绘图时统一用stem函数。,一、信号及其运算表示,2018/10/18,4,一、 几种常用信号的表示,1单位冲击函数、序列,2018/10/18,5,单位冲激函数、单位冲激序列,t = -5:0.01:5; y = (t=0); subplot(1,2,1); plot(t, y, r); n = -5:5; x = (n=0); subplot(1,2,2); stem(n, x);,2018/10/18,6,1、单位冲激信号x=zeros(1,N); x(1)=1; %注：Matlab下标从1开始。,2018/10/18,7,2单位阶跃函数、序列,(一）几种常用信号的表示,20
3、18/10/18,8,单位阶跃函数、单位阶跃序列,t =-5:0.01:5; y=(t=0); subplot(1,2,1); plot(t, y, r)n = -5:5; x = (n=0); subplot(1,2,2); stem(n, x);,2018/10/18,9,2、单位阶跃信号x=ones(1,N); %注：Matlab下标从1开始。,2018/10/18,10,3实指数序列,直接实现：n=ns:nf; x=a.n;,4复指数序列,直接实现：n=ns:nf; x=exp(sigema+jw)*n);,5正(余)弦序列,直接实现：n=ns:nf; x=cos(w*n+sita);
4、,2018/10/18,11,生成上述三种信号,t =-5:0.01:5; subplot(2,2,1); a=2 y1=a.t plot(t, y1, r) subplot(2,2,2); a=2; theat=pi/3; y2=sin(2*pi*t+theat) plot(t, y2) subplot(2,1,2); w=4; y3=exp(a+j*w)*t); plot(t, y3, y),2018/10/18,12,工具箱中的信号产生函数,2018/10/18,13,产生锯齿波或三角波 调用格式：sawtooth(x,y)例： 产生f=50Hz的锯齿波、三角波 Fs=10000; %采
5、样频率 t=0:1/Fs:0.1; %采样间隔1/Fs f=50; %50Hz x1=sawtooth(2*pi*50*t,0); x2=sawtooth(2*pi*50*t,1); x3=sawtooth(2*pi*50*t,0.5);subplot(311); plot(t,x1); subplot(312); plot(t,x2); subplot(313); plot(t,x3);,1、sawtooth函数,2018/10/18,14,(周期)矩形脉冲矩形波调用格式： square(x,y)例：产生50Hz占空比分别为20 和50的矩形波。 Fs=10000; %采样频率 t=0:1/
6、Fs:0.1; %采样间隔1/Fs f=50; %50Hz x1=square(2*pi*50*t,20); x2=square(2*pi*50*t,50);subplot(211); plot(t,x1); subplot(212); plot(t,x2);,2、square函数,2018/10/18,15,产生矩形波 调用格式为：x = rectpuls (t, width)例：生成幅度为2，宽度T = 4、 中心在t = 0的矩形波x(t)以及x(t-T/2),t = -4 : 0.0001 : 4; T = 4; x1 = 2*rectpuls(t, T); subplot(121);
7、plot(t, x1); title(x(t); axis(-4 6 0 2.2); x2 = 2*rectpuls(t-T/2,T); subplot(122);plot(t, x2); title(x(t-T/2); axis(-4 6 0 2.2);,3、 rectpuls函数,2018/10/18,16,产生sinc波形或sin(t)/(t)波形例： t=linspace(-10,+10,200); x=sinc(t);plot(t,x);,4、sinc函数,2018/10/18,17,x=rand(1,N); %产生0,1上均匀分布的随机信号。,5、均匀分布的随机信号,2018/10
8、/18,18,x=randn(1,N); %产生均值为0，方差为1的高斯分布随机信号（即白噪声信号）。,6、高斯分布的随机信号,2018/10/18,19,（二） 信号的基本运算,1 信号的相加与相乘,2 信号的移位与圆周移位,3 信号的翻转与累加,4 信号的卷积与相关,5 信号的能量与功率,2018/10/18,20,（二） 信号的基本运算,1信号的相加与相乘 y(n)=x1(n)+x2(n) MATLAB实现：y=x1+x2y(n)=x1(n)x2(n)MATLAB实现：y=x1.*x2,2018/10/18,21,t=-5:0.01:5; y1=jieyue(t,-2); y2=jiey
9、ue(t,2); f=(1+t./2).*(y1-y2) plot(t,f),function y=jieyue(t,t0) %产生一个阶跃序列 %t0代表起始时刻 y=(t-t0)=0),2018/10/18,22,2序列移位与周期延拓移位运算,序列移位：y(n)=x(n-m)。 MATLAB实现：y=x; ny=nx+m序列圆周延拓：y(n)=x(n)M， MATLAB实现：ny=nxs:nxfy=x(mod(ny,M)+1) 序列圆周延拓移位：y(n)=x(n-m)M， MATLAB实现：ny=nxs:nxfy=x(mod(ny-m,M)+1),2018/10/18,23,clc;clf
10、; N=24;M=7;m=3; n1=0:8; x1=(0.8).n1; %生成指数序列 x2=x1; n2=n1+m; subplot(4,1,1);stem(n1,x1);ylabel(x(n);axis(0 N-1 0 1); subplot(4,1,2);stem(n2,x2);ylabel(x(n-3);axis(0 N-1 0 1); n=0:N-1; if Mlength(x1)for i=1:length(x1)-Mx1(i)=x1(i)+x1(length(x1)-i)end end xc=x1(mod(n,M)+1); %产生x(n)周期延拓x(n)8 xcm=x1(mod
11、(n-m,M)+1); %产生移位序列x(n-3)周期延拓x(n-3)8 subplot(4,1,3);stem(n,xc);ylabel(x(n)_8);axis(0 N-1 0 max(xc); subplot(4,1,4);stem(n,xcm);ylabel(x(n-3)_8);axis(0 N-1 0 max(xcm);,2018/10/18,24,3 信号翻褶与累加运算,序列翻褶：y(n)=x(-n)。MATLAB可实现: y=fliplr(x),序列累加的数学描述为：,MATLAB实现：y=cumsum(x),2018/10/18,25,4 序列的卷积与相关运算,两序列卷积运算：
12、,MATLAB实现： y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。,两序列相关运算：,MATLAB实现： y=xcorr(x1,x2)。,2018/10/18,26,5 信号的能量和功率,1.信号能量,数字定义：,MATLAB实现: E=sum(x.*conj(x);或 E=sum(abs(x).2);,数字定义：,2. 信号功率,MATLAB实现: P=sum(x.*conj(x)/N; 或 E=sum(abs(x).2)/N;,2018/10/18,27,已知某线性移不变系统单位抽样相应为,系统输入信号为 求系统输出,n=0:5 x1=(n-2)=0); x2=(n
13、=0)-(n=3); y=conv(x1,x2) n1=0:10 subplot(311);stem(n,x1); subplot(312);stem(n,x2); subplot(313);stem(n1,y),举例,2018/10/18,28,f1(t)=e-t(t)- (t-2)、 f2(t) = (t)- (t-3)，编制一个m文件，绘出f1(t)、f2(t)和f(t)=f1(t)* f2(t)的波形。,p=0.001; k1=0:p:2 f1=exp(-k1) k2=0:p:3; f2=rectpuls(k2-1.5,3); f=conv(f1,f2); %计算序列f1,f2的卷积和
14、f f=f*p; k0=k1(1)+k2(1); %计算序列f非零样值的起点位置 k3=length(f1)+length(f2)-2; %计算卷积和f的非零样值的宽度 k=k0:p:k3*p; %确定卷积和f非零样值的时间向量 subplot(2,2,1); plot(k1,f1); title(f1(t); xlabel(t); ylabel(f1(t); axis(0 5 0 1.2); subplot(2,2,2); plot(k2,f2); title(f2(t); xlabel(t); ylabel(f2(t); axis(0 5 0 1.2); subplot(2,1,2); p
15、lot(k,f); %画出卷积f(t)的波形 title(f(t)=f1(t)*f2(t); xlabel(t); ylabel(f(t); axis(0 5 0 1.2);,2018/10/18,29,二、线性时不变系统时域分析,系统描述系统表示方法转换时域响应 连续系统 离散系统,2018/10/18,30,1 系统的描述,1常系数线性微分/差分方程,2系统传递函数：tf,3零极点增益模型zpk,连续系统：,连续系统：,离散系统：,离散系统：,2018/10/18,31,4极点留数模型,离散系统：,连续系统：,5二次分式模型sos,连续系统：,离散系统：,6状态空间模型ss,连续系统：,离
16、散系统：,1 系统的描述,2018/10/18,32,在Matlab中，传递函数用分子、分母两个多项式的系数表示，系数为降幂排列。分子（Numerator） ：B=b(1) b(2) b(m+1) 分母（Denominator）：A=a(1) a(2) b(n+1)例：num=1 0.2 1; den=1 0.5 1;,(1)传递函数表示法,2018/10/18,33,在Matlab中，增益系数、零点向量、极点向量用三个列向量表示。增益系数（Gain） ：k 零点向量（Zero）：z=z1 z2 zn 极点向量（Pole）：p=p1 p2 pnsys=zpk(z,p,k) %获得零-极点模型表
17、达式,(2)零极点模型表示法zpk,zplane：离散系统零极点分布图,pzmap：连续系统零极点分布图,2018/10/18,34,连续系统状态空间方程：离散系统状态空间方程：状态向量：x 输出向量：y 激励向量（输入向量）：u在Matlab中，用矩阵A、B、C、D表示系统的状态空间模型。,(3)状态空间模型表示法：ss,2018/10/18,35,2 系统模型的转换函数,在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零极点增益模型。其中sos表示二次分式，g为比例系数，sos为L6的矩阵，即,（415）,1ss2tf函数 格式：num, de
18、n=ss2tf(A,B,C,D,iu) 功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函数模型num,den。,2zp2tf函数 格式：num,den=zp2tf(z,p,k) 功能：将给定系统的零极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。,2018/10/18,36,2 系统模型的转换函数,2018/10/18,37,3、线性时不变系统时域响应,连续系统 离散系统 冲激响应、阶跃响应 零状态响应、零输入响应、全响应,2018/10/18,38,连续系统的时域响应,1、impulse函数 求连续系统的单位冲击响应。,调用格式：impulse
19、(sys,t)：绘出系统在 0t 时间范围内冲激响应的时域波形,2018/10/18,39,连续系统的时域响应,2、step函数 求连续系统的单位阶跃响应。,调用格式： step (sys,t)： 绘出系统在 0t 时间范围内冲激响应的时域波形,2018/10/18,40,连续系统的时域响应,3、 initial函数 求连续系统的零输入响应。,调用格式：initial(sys,t,x0)：绘出连续时间LTI系统由于初始状态x0所引起的零输入响应,2018/10/18,41,格式：y,x=lsim(a,b,c,d,u,t) 功能：返回连续LTI系统,对任意输入时系统的输出响应y和状态记录x，其中
20、u给出每个输入的时序列，一般情况下u为一个矩阵；t用于指定仿真的时间轴，它应为等间隔。,4.对任意输入的连续LTI系统响应函数lsim( ),2018/10/18,42,2、当输入 时，系统的输出。,已知描述某连续系统的微分方程为： 1、试用Matlab绘出该系统冲激响应和阶跃响应。,3、y(0+)=-1，y(0+)=0时，再求解上问。,连续时间系统时域响应,举例,2018/10/18,43,*1、impz函数 求离散系统（数字滤波器）的单位冲击响应。 （注：Matlab 7.0不再支持dimpulse函数）,离散时间系统响应,2018/10/18,44,离散系统的时域响应,2、stepz函数
21、 求离散系统的单位阶跃响应。, b=1 2; a=1 2 1; stepz(b,a) %阶跃响应,2018/10/18,45,离散系统的时域响应,3、 dinitial函数 求离散系统的零输入响应。,调用格式：dinitial(sys,t,x0)绘出连续时间LTI系统由于初始状态x0所引起的零输入响应,2018/10/18,46,4、对任意输入的离散LTI系统响应函数dlsim( ) 格式：y,x=dlsim(a,b,c,d,u) 功能：返回离散LTI系统,对输入序列u的响应y和状态记录x。,离散系统的时域响应,2018/10/18,47,4滤波函数filter 格式：y=filter(B,A
22、,x),例 设系统差分方程为,MATLAB源程序为： B=1; A=1,-0.8; N=0:31; x=0.8.n; y=filter(B,A,x); subplot(2,1,1);stem(x) subplot(2,1,2);stem(y),，求该系统对信号,的响应。,离散系统的时域响应,2018/10/18,48,试画出系统的零极点分布图、 求系统的单位脉冲与阶跃响应 判断系统是否稳定,举例：,离散系统的时域响应,2018/10/18,49,三、频域分析,傅里叶分析 连续时间傅里叶变换 离散傅里叶变换 FFT与卷积 系统频域响应,2018/10/18,50,傅里叶(Fourier)变换,1
23、 连续时间、连续频率傅里叶变换,2 连续时间、离散频率傅里叶级数,正变换：,逆变换：,正变换：,逆变换：,2018/10/18,51,3 时间离散、连续频率序列傅里叶变换,4 离散时间、离散频率离散傅里叶级数,5离散时间、离散频率离散傅里叶变换(DFT),正变换：,逆变换：,正变换：,逆变换：,正变换：,逆变换：,傅里叶(Fourier)变换,2018/10/18,52,信号的傅里叶分析,Matlab的符号运算工具箱（Symbolic Math Toolbox）提供了能直接求解傅里叶变换和逆变换的符号运算函数fourier()和ifourier()。函数的调用格式如下F=fourier(f)f
24、=ifourier(F),2018/10/18,53,求单边指数函数 的傅里叶变换，画出其幅频特性和相频特性图,syms t w f f=exp(-2*t)*sym(Heaviside(t); F=fourier(f) subplot(3,1,1);ezplot(f,0:2,0:1.2); subplot(3,1,2);ezplot(abs(F),-10:10); subplot(3,1,3);ezplot(angle(F),-10:10),信号的傅里叶分析-举例,2018/10/18,54,计算的时域波形x(t) =e-tu(t) 及其频谱,T=0.01; t=-1:T:10; x=exp(
25、-t).*(t0); subplot(3,1,1); plot(t,x);%时域波形 axis(-1 10 -0.1 1.1); xlabel(timesec); w=-40:0.01:40; X=x*exp(-j*t*w)*T; subplot(3,1,2); plot(w,abs(X);%频域幅度谱 subplot(3,1,3); plot(w,angle(X);%频域相位谱,连续傅里叶变换-举例,2018/10/18,55,傅里叶级数,MATLAB实现：ak = x1*exp(-j*k*w0*t)*dt/T;,2018/10/18,56,傅里叶级数,clc; dt=0.01; T=2;t
26、=-4:dt:4;w0=2*pi/T; x1=rectpuls(t-0.5,1); x=0; for m = -1:1 % Periodically extend x1(t) to form a periodic signalx = x+rectpuls(t-0.5-m*T),1); end subplot(2,1,1) plot(t,x); %line(0 0,-0.1 1.1); axis(-4 4 0 1.1); title(x(t); xlabel(Time index t); N=10; for k=-N: N; ak(N+1+k) = x1*exp(-j*k*w0*t)*dt/T;
27、 end k=-N: N; subplot(2,1,2) stem(k,abs(ak),k.); title(The Fourier series coefficients); xlabel(Frequency index k);,2018/10/18,57,DFS、DFT,function n,X=dft(x,N)n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;x=x zeros(1,N-length(x);X=x*WNnk;,2018/10/18,58,clear;close all; %计算8点DFT n=0:7; x=cos(n*
28、pi/6); n1,X1=dft(x,8); subplot(2,1,1); stem(n1,abs(X1),.); ylabel(|X(k)|); title(8点DFT); %计算16点DFT n2,X2=dft(x,16); subplot(2,1,2); stem(n2,abs(X2),.); ylabel(|X(k)|); title(16点DFT);,2018/10/18,59,1一维快速正傅里叶变换函数fft 格式：X=fft(x, N) 功能：采用FFT算法计算序列向量x的N点DFT变换， 当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计
29、算，否则用混合算法。2一维快速逆傅里叶变换函数ifft 格式：x=ifft(X, N) 功能：采用FFT算法计算序列向量X的N点IDFT变换。,FFT,2018/10/18,60,计算的时域波形x(t) =e-tu(t) 及其频谱,T=0.01; t=-1:T:10; x=exp(-t).*(t0); subplot(3,1,1); plot(t,x);%时域波形 axis(-1 10 -0.1 1.1); xlabel(timesec); w=-40:0.01:40; X=x*exp(-j*t*w)*T; subplot(3,1,2); plot(w,abs(X);%频域幅度谱 subplo
30、t(3,1,3); plot(w,angle(X);%频域相位谱,2018/10/18,61,k=0:15;f=cos(pi*k./2);F_16=fft(f);F_512=fft(f,512);L=0:511;stem(L/512,abs(F_512);hold on;stem(k/16,abs(F_16),r:);grid on;xlabel(Normalized frequency);ylabel(Magnirude);hold off,FFT-举例,2018/10/18,62,MATLAB程序(部分)： %线性卷积 xn= sin(0.4*1:15); %对序列x(n)赋值, M=15
31、 hn= 0.9.(1:20); %对序列h(n)赋值, N=20 yn=conv(xn,hn); % 直接调用函数conv计算卷积%园周卷积 L=pow2(nextpow2(M+N-1); Xk=fft(xn,L); Hk=fft(hn,L); Yk=Xk.*Hk; yn=ifft(Yk,L);,例 用快速傅里叶变换FFT计算下面两个序列的卷积。,FFT-举例,2018/10/18,63,连续傅里叶变换,上式用MATLAB表示为：X=x*exp(-j*t*w)*T,2018/10/18,64,x=1 2 0 1; h=2 2 1 1;,x=1 2 0 1;h=2 2 1 1;L=length
32、(x)+length(h)-1;XE=fft(x,L);HE=fft(h,L);y1=ifft(XE.*HE);k=0:L-1;subplot(2,1,1);stem(k,real(y1);axis(0 6 0 7);title(Result of Linear Convolution);xlabel(Time index k);ylabel(Amplitude);y2=conv(x,h);error=y1-y2;subplot(2,1,2);stem(k,abs(error);xlabel(Time index k);ylabel(Amplitude);title(Error Magnitu
33、de);,2018/10/18,65,S变换与Z变换,S变换与系统频响Z变换与系统频响,2018/10/18,66,1求连续系统的H(S),MATLAB的符号数学工具箱(Symbolic Math Tools)提供了计算拉普拉斯变换的函数：正S变换： laplace(s)函数 逆S变换：ilaplace(W)函数,2018/10/18,67,求以下函数的拉普拉斯变换。,syms t f1=sym(exp(-2*t)*Heaviside(t); F1=laplace(f1) %求f1(t)的拉普拉斯变换 f2=sym(t*exp(-t)*Heaviside(t); F2=laplace(f2),
34、2018/10/18,68,2求连续系统Ha(s)的频率响应函数freqs( ),格式：Hfreqs(B,A,W) 功能：计算由向量W(rad/s)指定的频率点上系统函数Ha(s)的频率响应Ha(j)，结果存于H向量中。,例 已知某线性时不变的系统函数,求该系统的频率响应。,MATLAB源程序如下。 B=1;A=1 2.6131 3.4142 2.6131 1; W=0:0.1:2*pi*5; freqs(B,A,W),2018/10/18,69,3求离散系统的H(z),MATLAB的符号数学工具箱(Symbolic Math Tools)提供了计算函数：z正变换：ztrans(x)函数 逆z
35、变换：iztrans(Z)函数,2018/10/18,70,例 已知某离散系统函数为,求该系统的频率响应,4求离散系统H(z)的频率响应函数freqz( ) 格式：H=freqz(B,A,W)功能：计算由向量W(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw )。,MATLAB源程序为： B=1 0 0 0 0 0 0 0 1; A=1; freqz(B,A),2018/10/18,71,四、滤波器设计,IIR滤波器设计 FIR滤波器设计,2018/10/18,72,滤波器结构,在MATLAB中，提供了以下几种扩展函数来实现这些结构形式的转换及实现。 dir2cas变直接型为
36、级联型、 cas2dir 变级联型为直接型、 dir2par 变直接型为并联型、 par2dir 变并联型为直接型、 dir2fsm 变直接型为频率采样型、 filter滤波器的直接实现 casfiltr滤波器的级联实现、parfiltr 滤波器的并联实现。,2018/10/18,73,IIR数字滤波器的设计方法,1. 数字滤波器的频率响应函数,幅度响应：,相位响应：,图 理想低通、高通、带通、带阻数字滤波器幅度特性,2018/10/18,74,2. 滤波器的技术指标幅度响应指标、相位响应指标,图4.38 数字低通滤波器的幅度特性,通带要求：,阻带要求：,通带最大衰减：,阻带最小衰减：,201
37、8/10/18,75,MATLAB的直接设计,cheby1函数：Chebyshev 型滤波器设计 （通带等纹波）cheby2函数：Chebyshev 型滤波器设计 （阻带等纹波）,n,Wn = buttord(Wp,Ws,Rp,Rs); b,a=butter(n,Wn,bandpass);,调用格式,2018/10/18,76,例：从多种频率成分叠加的信号中，提取出某一频带的信号。例如：从500Hz、1kHz、2kHz三种频率成分叠加的信号中，提取出频带为800Hz1.5kHz的信号。要求：阻带衰减不低于10dB，通带衰减不高于3dB。t=0:1/(40000):1/200; x1=1.8*s
38、in(2*pi*500*t+pi/6); x2=2.0*sin(2*pi*1000*t-pi/10); x3=2.2*sin(2*pi*2000*t+pi/12); y=x1+x2+x3; %先选择滤波器阶数 Fs=5000; %Fs2*2kHz Wp=800 1500/(Fs/2); Ws=500 2000/(Fs/2); Rp=3; Rs=10; n,Wn = buttord(Wp,Ws,Rp,Rs); %再进行滤波器设计 b,a=butter(n,Wn,bandpass); z=filter(b,a,y);,subplot(511);plot(x1); subplot(512);plot
39、(x2); subplot(513);plot(x3); subplot(514);plot(y); subplot(515);plot(z);,butter函数：Butterworth滤波器设计,2018/10/18,77,例 试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1650Hz、wp2850Hz；阻带的截止频率：ws1700Hz、ws2800Hz；通带内的最大衰减为rp0.1dB；阻带内的最小衰减为rs50dB；采样频率为Fs2000Hz。 MATLAB源程序设计如下：wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs
40、=2000;wp=wp1,wp2/(Fs/2);ws=ws1,ws2/(Fs/2); %利用Nyquist频率频率归一化N,wc=ellipord(wp,ws,rp,rs,z); %求滤波器阶数num,den=ellip(N,rp,rs,wc,stop); %求滤波器传递函数H,W=freqz(num,den); %绘出频率响应曲线plot(W*Fs/(2*pi),abs(H);grid;xlabel(频率/Hz);ylabel(幅值),ellip函数：椭圆滤波器设计,2018/10/18,78,冲激响应不变法,2.MATLAB信号处理工箱中的专用函数impinvar( )： 格式：BZ,AZ
41、 =impinvar(B,A,Fs) 功能：把具有B,A模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的 传递函数模型BZ,AZ。采样频率Fs的默认值为Fs=1。,1. 冲激响应不变法设计IIR数字滤波器的基本原理：,例4 MATLAB源程序如下： num=1; %模拟滤波器系统函数的分子 den=1,sqrt(5),2,sqrt(2),1; %模拟滤波器系统函数的分母 num1,den1=impinvar(num,den) %求数字低通滤波器的系统函数 程序的执行结果如下： num1 = -0.0000 0.0942 0.2158 0.0311 den1 = 1.0000 -
42、2.0032 1.9982 -0.7612 0.1069,2018/10/18,79,MATLAB信号处理工具箱中的专用双线性变换函数bilinear( ) 格式：numd,dendbilinear(num,den,Fs) 功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。,双线性变换法,双线性变换利用频率变换关系：,例 MATLAB源程序如下：num=1; %模拟滤波器系统函数的分子den=1,sqrt(3),sqrt(2),1; %模拟滤波器系统函数的分母num1,den1=bilinear(num,den,1) %求数字滤波器的传递函数 运算的结果如下： num1 = 0.
43、0533 0.1599 0.1599 0.0533 den1 = 1.0000 -1.3382 0.9193 -0.1546,2018/10/18,80,IIR数字滤波器的频率变换设计法,1. IIR数字滤波器的频率变换设计法的基本原理根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转换为相应的模拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法，将模拟滤波器数字化成相应的数字滤波器。,2018/10/18,81,1MATLAB的典型设计,(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标； (2)根据转换后的技术指标使用滤波器阶数函数，确定滤波
44、器的最小阶数N和截止频率Wc； (3)利用最小阶数N产生模拟低通滤波原型； (4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器； (5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。,2018/10/18,82,MATLAB源程序设计如下：%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;N,Wc=buttord(wp,ws,rp,rs,s); %选择滤波器的最小阶数Z,P,K=buttap(N); %创建Butterworth低通滤波器原型A,B,C,D=zp2s
45、s(Z,P,K); %零极点增益模型转换为状态空间模型AT,BT,CT,DT=lp2hp(A,B,C,D,Wc); %实现低通向高通的转变num1,den1=ss2tf(AT,BT,CT,DT); %状态空间模型转换为传递函数模型%运用双线性变换法把模拟滤波器转换成数字滤波器num2,den2=bilinear(num1,den1,100); H,W=freqz(num2,den2); %求频率响应plot(W*Fs/(2*pi),abs(H);grid; %绘出频率响应曲线xlabel(频率/Hz); ylabel(幅值),例 设计一个数字信号处理系统，它的采样率为Fs100Hz，希望在该系
46、统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。,2018/10/18,83,FIR数字滤波器设计,格式：w = boxcar(M) 功能：返回M点矩形窗序列。 MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1( )。 格式：Bfir1(N,wc) 功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B 为滤波器的系数(单位冲激响应序列)，其长度为N+1。,1 窗函数设计法 窗函数设计的基本原理：h(n)=w(n)hd(
47、n)w(n)为窗函数， hd(n)理想数字滤波器的单位冲激响应。在MATLAB信号处理工具箱中为用户提供了Boxcar (矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这些窗函数的调用格式相同。,FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称h(n)=h(N-n-1) 或 h(n)=-h(N-n-1) FIR数字滤波器具有线性相位：,或,2018/10/18,84,窗函数,1.矩形窗： w = boxcar ( N ) 2.巴特利特窗：w = bartlett ( N ) 3.汉宁窗： w = hanning ( N ) 4.汉明窗： w = hamming (
48、N ) 5.布莱克曼窗： w = blackman ( N ) 6.凯塞（Kaiser）窗： w = kaiser (N, beta),2018/10/18,85,例 用矩形窗设计线性相位FIR低通滤波器。该滤波器的通带截止频率wc=pi/4，单位脉冲响h(n)的长度M=21。并绘出h(n)及其幅度响应特性曲线。M=21; wc=pi/4; n=0:M-1; r=(M-1)/2; nr=n-r+eps*(n-r)=0); hdn=sin(wc*nr)./(pi*nr); if rem(M,2)=0, hdn(r+1)=wc/pi; end wn1=boxcar(M); hn1=hdn.*wn1; subplot(2,1,1);stem(n,hn1,.); line(0,20,0,0); xlabel(n),ylabel(h(n),title(矩形窗设计的h(n); hw1=fft(hn1,512); w1=2*0:511/512; subplot(2,1,2), plot(w1,20*log10(abs(hw1) xlabel(w/pi), ylabel(幅度(dB); title(幅度特性(dB);,
