序列逆Z变换的Matlab实现

一、序列逆Z变换的Matlab实现•函数residuez:适合计算离散系统有理函数的留数和极点，可以用于求解序列的逆Z变换。101111001()()(2.98)()1MNMNkMkkNkkNkbbzbzRBzXzCzaazazAzpz•函数residuez基本调用方式:[r,p,c]=residuez(b,a);–输入参数:b=[b0，b1，…，bM]为分子多项式的系数,a=[a0，a1，…，aN]为分母多项式的系数，这些多项式都按z的降幂排列–输出参数:r是极点的留数，p是极点，c是无穷项多项式的系数项，仅当M≥N时存在。例：计算逆Z变换例：计算的逆Z变换。解:有理分式X(z)分子和分母多项式都按z的降幂排列。2()231zXzzz12120()23123zzXzzzzz•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.00000.5000•系数项:X(z)的部分分式形式为1111()110.5Xzzz逆Z变换为()()(0.5)()nxnunun二、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),gridon;title('序列x(n)')y1=conv(x,x);%y1为x(n)与x(n)的线性卷积subplot(2,2,2)stem(0:1:length(y1)-1,y1),gridon;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),gridon;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),gridon;title('x(n)与x(n)的8点循环卷积')程序运行结果图00.511.522.5300.20.40.60.81序列x(n)012345601234x(n)与x(n)线性卷积00.511.522.5301234x(n)与x(n)的4点循环卷积0246801234x(n)与x(n)的8点循环卷积三、频域取样定理的Matlab实现例：设x(n)=(0.7)n·u(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)n·u(n)的Z变换为()0.7zXzz程序清单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点取样恢复的序列')程序运行结果0510152000.20.40.60.81时域信号x(n)0510152000.511.55点采样恢复的序列0510152000.511.520点采样恢复的序列四、用FFT进行谱分析的Matlab实现设模拟信号，以t=0.01n(n=0:N-1)进行取样，试用fft函数对其做频谱分析。N分别为：(1)N=45；(2)N=50；(3)N=55；(2)N=60。()2sin(4)5cos(8)xttt程序清单如下%计算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('FFTN=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('FFTN=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('FFTN=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('FFTN=60')程序运行结果02468050100150FFTN=4502468050100150FFTN=5002468050100150FFTN=5502468050100150FFTN=60从图中可以看出，这几种情况下均有较好的精度。程序运行结果分析分析：由t=0.01n进行取样可得，采样频率fs=100Hz。而连续信号的最高模拟角频率为Ω＝8π，由Ω＝2πf可得，最高频率为8π/2π=4Hz。因此，满足采样定理的要求。采样序列为()2cos(4)5cos(8)xnTnTn48()2cos5cos100100xnnn即为周期序列，周期N=50。将程序中plot改为stem函数，则可以更清楚地看出频谱。修改程序运行结果02468050100150FFTN=4502468050100150FFTN=5002468050100150FFTN=5502468050100150FFTN=606.5IIR数字滤波器的Matlab仿真实现•IIR数字滤波器设计•模拟滤波器到数字滤波器的转换五、IIR数字滤波器设计•设数字滤波器系统函数为•模拟滤波器的系统函数为•函数butter和cheby1可以确定Butterworth和ChebyshevI型滤波器的系统函数。11()(1)(2)(1)()()1(2)(1)nnBzbbzbnzHzAzazanz11()(1)(2)(1)()()(2)(1)nnannBsbsbsbnHsAssasan函数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分别为数字滤波器的系统函数的系数)11()(1)(2)(1)()()(1)(2)(1)NNaNNBsbsbsbNHsAsasasaN11(1)(2)(1)()(1)(2)(1)NNNNbzsbzsbzNHzazsazsazN脉冲响应不变法•脉冲响应不变法:用代换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