第章 离散时间信号与系统 例1-1 用MATLAB计算序列-2 0 1 1 3和序列1 2 0 -1的离散卷积。 解 MATLAB程序如下： a=-2 0 1 -1 3; b=1 2 0 -1; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel(n); ylabel(幅度); 图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：-2 -4 1 3 1 5 1 -3。 例1-2 用MATLAB计算差分方程 当输入序列为 时的输出结果 。 解 MATLAB程序如下： N=41;a=0.8 -0.44 0.36 0.22;b=1 0.7 -0.45 -0.6;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y)xlabel(n);ylabel(幅度) 图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 例1-3 用MATLAB计算例1-2差分方程 所对应的系统函数的DTFT。 解 例1-2差分方程所对应的系统函数为： 其DTFT为 用MATLAB计算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle(实部)xlabel(omega/pi);ylabel(幅度)subplot(2,2,2);plot(w/pi,imag(h);gridtitle(虚部)xlabel(omega/pi);ylabel(Amplitude)subplot(2,2,3);plot(w/pi,abs(h);gridtitle(幅度谱)xlabel(omega/pi);ylabel(幅值)subplot(2,2,4);plot(w/pi,angle(h);gridtitle(相位谱)xlabel(omega/pi);ylabel(弧度) 第2章 离散傅里叶变换及其快速算法 例2-1 对连续的单一频率周期信号 按采样频率 采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。 解 此时离散序列 ，即k=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下： k=8; n1=0:1:19; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel(t/T);ylabel(x(n); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel(k);ylabel(X(k); n2=0:1:15; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel(t/T);ylabel(x(n); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel(k);ylabel(X(k); 计算结果示于图2.1，(a)和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。 例2-2 用FFT计算两个序列 的互相关函数 。 解 用MATLAB计算程序如下： x=1 3 -1 1 2 3 3 1; y=2 1 -1 1 2 0 -1 3; k=length(x); xk=fft(x,2*k);yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k); m=(-k+1):(k-1); stem(m,rm) xlabel(m); ylabel(幅度); 其计算结果如图2.2所示。 例2-3计算两个序列的的互相关函数，其中 x(n)=2 3 5 2 1 1 0 0 12 3 5 3 0 1 2 0 1 2 y(n)=x(n-4)+e(n), e(n)为一随机噪声，在MATLAB中可以用随机函数rand产生 解 用MATLAB计算程序如下： x=2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2; y=0 0 0 0 2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2; k=length(y); e=rand(1,k)-0.5; y=y+e; xk=fft(x,2*k); yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k); m=(-k+1):(k-1); stem(m,rm)xlabel(m); ylabel(幅度); 计算结果如图2.3(a)，我们看到最大值出现在m=4处，正好是y(n)对于x(n)的延迟。2. 3(b)是x(n)的自相关函数，他和y(n)的区别除时间位置外，形状也略不同，这是由于y(n)受到噪声的干扰。第3章 无限长单位脉冲响应(IIR)滤波器的设计方法 例3-1 设采样周期T=250s（采样频率fs =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc =1kHz。 B,A=butter(3,2*pi*1000,s); num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,s); num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),-.,f,abs(h2),-); grid; xlabel(频率/Hz ) ylabel(幅值/dB) 程序中第一个butter的边界频率21000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图3.1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。同时，也看到双线性变换法，在z=-1即=或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在=处的三阶传输零点通过映射形成的。 例3-2 设计一数字高通滤波器，它的通带为400500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。 wc=2*1000*tan(2*pi*400/(2*1000); wt=2*1000*tan(2*pi*317/(2*1000); N,wn=cheb1ord(wc,wt,0.5,19,s); B,A=cheby1(N,0.5,wn,high,s); num,den=bilinear(B,A,1000); h,w=freqz(num,den); f=w/pi*500; plot(f,20*log10(abs(h); axis(0,500,-80,10); grid; xlabel() ylabel(幅度/dB) 图3.2给出了MATLAB计算的结果，可以看到模拟滤波器在=处的三阶零点通过高通变换后出现在=0（z=1）处，这正是高通滤波器所希望得到的。 例3-3 设计一巴特沃兹带通滤波器，其dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3 = 120kHz处的最小衰减大于dB，采样频率fs=400kHz。 w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0 wr,3,10,s); B,A=butter(N,wn,s); num,den=bilinear(B,A,400); h,w=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel(频率/kHz) ylabel(幅度/dB) 图3.3给出了MATLAB计算的结果，可以看出数字滤波器将无穷远点的二阶零点映射为z=1的二阶零点，数字带通滤波器的极点数是模拟低通滤波器的极点数的两倍。 例3-4 一数字滤波器采样频率fs = 1kHz，要求滤除100Hz的干扰，其dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为 w1=95/500; w2=105/500; B,A=butter(1,w1, w2,stop); h,w=freqz(B,A); f=w/pi*500; plot(f,20*log10(abs(h); axis(50,150,-30,10); grid; xlabel(频率/Hz) ylabel(幅度/dB) 图3.4为MATLAB的计算结果 第4章 有限长单位脉冲响应(FIR)滤波器的设计方法 例2