序列逆Z变换的Matlab实现.ppt
一、序列逆Z变换的Matlab实现,函数residuez:适合计算离散系统有理函数的留数和极点，可以用于求解序列的逆Z变换。,函数residuez基本调用方式:r,p,c=residuez(b,a);输入参数:b=b0，b1，bM为分子多项式的系数,a=a0，a1，aN为分母多项式的系数，这些多项式都按z的降幂排列 输出参数:r是极点的留数，p是极点，c是无穷项多项式的系数项，仅当MN时存在。,例：计算逆Z变换,例：计算 的逆Z变换。,解:有理分式X(z)分子和分母多项式都按z的降幂排列。,b=0,1;a=2,-3,1;%多项式的系数r,p,c=residuez(b,a);%求留数、极点和系数项disp(留数:);disp(r);%显示输出参数disp(极点:);disp(p);disp(系数项:);disp(c);,程序运行结果为留数:1-1极点:1.0000 0.5000系数项:,X(z)的部分分式形式为,逆Z变换为,二、DFT物理意义的Matlab实现,序列的N点DFT的物理意义:对X(ej)在0，2上进行N点的等间隔取样。函数fft用于快速计算离散傅里叶变换，调用方式为y=fft(x);y=fft(x,N);y=fft(x)利用FFT算法计算序列x的离散傅里叶变换。当x为矩阵时，y为矩阵x每一列的FFT。当x长度为2的整数次幂时，函数fft采用基-2的FFT算法，否则采用混合基算法。y=fft(x,N)采用N点FFT。当序列x长度小于N时，函数fft自动对序列尾部补零，构成N点数据；当x长度大于N时，函数fft自动截取序列前面N点数据进行FFT。,函数ifft用于快速计算向量或矩阵的离散傅里叶逆变换，与函数fft的调用规则基本相同。调用方式为y=ifft(x);y=ifft(x,N);,例：利用FFT实现线性卷积,例：利用FFT实现线性卷积。已知序列x(n)=R4(n)，求：（1）用conv函数求x(n)与x(n)的线性卷积y(n)，并绘出图形；（2）用FFT求x(n)与x(n)的4点循环卷积y1(n)，并绘出图形；（3）用FFT求x(n)与x(n)的8点循环卷积y2(n)，并将结果与（1）比较，说明线性卷积与循环卷积之间的关系。,解 程序如下：,N1=4;N2=8;n1=0:1:N1-1;n2=0:1:N2-1;x=1,1,1,1;%构造序列x(n)x1=1,1,1,1,0,0,0,0;%在序列x(n)后补4个零figure(1)subplot(2,2,1)stem(n1,x),grid on;title(序列x(n)y1=conv(x,x);%y1为x(n)与x(n)的线性卷积subplot(2,2,2)stem(0:1:length(y1)-1,y1),grid on;title(x(n)与x(n)线性卷积),X2=fft(x);%计算x(n)与x(n)的4点循环卷积Y2=X2.*X2;y2=ifft(Y2);subplot(2,2,3)stem(n1,y2),grid on;title(x(n)与x(n)的4点循环卷积)X3=fft(x1);%计算x(n)与x(n)的8点循环卷积Y3=X3.*X3;y3=ifft(Y3)subplot(2,2,4)stem(n2,y3),grid on;title(x(n)与x(n)的8点循环卷积),程序运行结果图,三、频域取样定理的Matlab实现,例：设x(n)=(0.7)nu(n)，在单位圆上以M=5和M=20，对其Z变换取样，研究时域信号受M变化的影响。（1）对x(n)进行Z变换；（2）对X(z)进行等角取样，取样点数为M，求X(k)；（3）对X(k)进行IDFT变化，得到M点序列，请比较几个序列，并作分析。,解 x(n)=(0.7)nu(n)的Z变换为,程序清单,n=0:19;x=0.7.n;na=0:4;za=exp(j*2*pi*na/5);%在z平面的单位圆上对其进行5点的等角距取样Xa=za./(za-0.7);xa=abs(ifft(Xa);nb=0:19;zb=exp(j*2*nb*pi/20);%在z平面的单位圆上对其进行20点的等角距取样Xb=zb./(zb-0.7);xb=abs(ifft(Xb);,figure(1)subplot(2,2,1);%画出原始时域信号stem(n,x)title(时域信号x(n)subplot(2,2,2);xa=xa,xa,xa,xa;stem(n,xa)title(5点取样恢复的序列)subplot(2,2,3);stem(n,xb)title(20点取样恢复的序列),程序运行结果,四、用FFT进行谱分析的Matlab实现,设模拟信号，以 t=0.01n(n=0:N-1)进行取样，试用fft函数对其做频谱分析。N分别为：(1)N=45；(2)N=50；(3)N=55；(2)N=60。,程序清单如下,%计算N=45的FFT并绘出其幅频曲线N=45;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y)title(FFT N=45),程序清单,%计算N=50的FFT并绘出其幅频曲线N=50;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y)title(FFT N=50),%计算N=55的FFT并绘出其幅频曲线N=55;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,3)plot(q,abs(y)title(FFT N=55),%计算N=60的FFT并绘出其幅频曲线N=60;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,4)plot(q,abs(y)title(FFT N=60),程序运行结果,从图中可以看出，这几种情况下均有较好的精度。,程序运行结果分析,分析：由t=0.01n进行取样可得，采样频率fs=100Hz。而连续信号的最高模拟角频率为8，由2 f可得，最高频率为8/2=4Hz。因此，满足采样定理的要求。采样序列为,即,为周期序列，周期N=50。,将程序中plot改为stem函数，则可以更清楚地看出频谱。,修改程序运行结果,6.5 IIR数字滤波器的Matlab仿真实现,IIR数字滤波器设计 模拟滤波器到数字滤波器的转换,五、IIR数字滤波器设计,设数字滤波器系统函数为模拟滤波器的系统函数为函数butter和cheby1可以确定Butterworth和Chebyshev I型滤波器的系统函数。,函数butter的调用格式,函数butter的调用格式为b,a=butter(n,Wc,)%设计数字Butterworth滤波器b,a=butter(n,Wc,ftype)%设计模拟Butterworth滤波器其中，n为滤波器阶数，Wc为截止频率。,函数cheby1的调用格式,函数cheby1的调用格式为b,a=cheby1(n,Rp,Wc)%设计数字Chebyshev滤波器b,a=cheby1(n,Rp,Wc,ftype)%设计模拟Chebyshev滤波器其中，n为滤波器阶数，Rp为通带内的纹波系数，Wc为截止频率。,例：设计butterworth低通滤波器,设计一模拟butterworth低通滤波器，通带截止频率300Hz，通带最大衰减2dB，阻带截止频率800Hz，阻带最小衰减30dB。,解 滤波器的阶数和截止频率可由式和确定，程序段为Wp=2*pi*300;Ws=2*pi*800;Rp=2;Rs=30;N=ceil(log10(10(0.1*Rs)-1)/(10(0.1*Rp)-1)/(2*log10(Ws/Wp);Wc=Wp/(10(Rp/10)-1)(1/(2*N);b,a=butter(N,Wc,s);freqs(b,a),程序运行结果,运行程序，得到N=4，Wc=2.0157e+003。幅频特性和相频特性如图7.16所示。,模拟滤波器到数字滤波器的转换,设模拟滤波器系统函数为数字滤波器的系统函数为从模拟滤波器到数字滤波器的转换有两种方法，即脉冲响应不变法和双线性变换法。(b,a分别为模拟滤波器的系统函数的系数，bz,az分别为数字滤波器的系统函数的系数),脉冲响应不变法,脉冲响应不变法:用代换Ha(s)中的(s-sk)即可得到H(z)，从而将模拟滤波器转换为数字滤波器格式。可用函数impinvar实现，调用格式为bz,az=impinvar(b,a,fs)其中，fs为取样频率。,双线性变换法,双线性变换法:用代换Ha(s)中的s即可得到H(z)，从而将模拟滤波器转换为数字滤波器格式。可用函数bilinear实现，调用格式为zd,pd,kd=bilinear(z,p,k,fs)其中，z，p，k和zd，pd，kd分别为s域和z域系统函数的零点、极点和增益。,例：模拟滤波器转换数字滤波器,利用impinvar将一模拟低通滤波器变换成数字滤波器(取样频率为10Hz)，,程序段为b,a=butter(4,.3,s);bz,az=impinvar(b,a,10);程序运行结果为bz=1.0e-006*-0.0000 0.1324 0.5192 0.1273 0az=1.0000-3.9216 5.7679-3.7709 0.9246,六、FIR数字滤波器的Matlab仿真实现,窗函数法设计FIR滤波器FIR滤波器的优化设计,窗函数法设计FIR滤波器,窗函数法是通过对理想滤波器的单位取样响应加窗来逼近理想滤波器的。函数fir1用于设计标准的低通、带通、高通和带阻滤波器。,函数fir1的调用格式,函数fir1的调用格式为 b=fir1(n,Wc,ftype,Windows)其中，n为滤波器阶数，Wc为截止频率ftype决定滤波器类型，ftype=high，设计高通FIR滤波器，ftype=stop，设计带阻FIR滤波器。Windows指定窗函数类型，默认为Hamming窗；可选Hanning、Hamming、Blackman、triangle、bartlett和boxcar窗，每种窗都可以由Matlab的相应函数生成。,例：设计butterworth低通滤波器,设计一个15阶FIR低通滤波器，截止频率为0.2。,b=fir1(15,0.2);freqz(b,1,512);函数freqz(b,a,N)用于计算由a和b构成的数字滤波器的频率响应，并用图形方式分别表示其幅度响应和相位响应。,程序运行结果,程序运行结果如图所示。,FIR滤波器的优化设计,前面介绍了FIR滤波器的优化设计方法，通过迭代的方法求解FIR滤波器，过程十分复杂。在Matlab中，可以调用函数remez实现滤波器的设计。b=remez(n,f,m)函数remez采用Parks-McClellan算法设计线性相位FIR滤波器，n为滤波器阶数，其幅频特性由f和m指定。,例：设计带通滤波器,采用Parks-McClellan算法设计一个17阶带通滤波器，并画出期望的幅频特性曲线和实际的幅频特性曲线，,程序段为f=0 0.3 0.4 0.6 0.7 1;m=0 0 1 1 0 0;b=remez(17,f,m);h,w=freqz(b,1,512);plot(f,m,w/pi,abs(h);,程序运行结果,程序运行结果如图所示。,
