数字信号处理书上实验1.2.3.4

实验一熟悉Matlab环境
一、实验目的
1.熟悉MATLAB的主要操作命令。

2.学会简单的矩阵输入和数据读写。

3.掌握简单的绘图命令。

4.用MATLAB编程并学会创建函数。

5.观察离散系统的频率响应。

二、实验内容
认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。

在熟悉了MATLAB基本命令的基础上，完成以下实验。

上机实验内容：
（1）数组的加、减、乘、除和乘方运算。

输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。

clear all;
a=[1 2 3 4];
b=[3 4 5 6];
c=a+b;
d=a-b;
e=a.*b;
f=a./b;
g=a.^b;
n=1:4;
subplot(4,2,1);stem(n,a);
xlabel('n');xlim([0 5]);ylabel('A');
subplot(4,2,2);stem(n,b);
xlabel('n');xlim([0 5]);ylabel('B');
subplot(4,2,3);stem(n,c);
xlabel('n');xlim([0 5]);ylabel('C');
subplot(4,2,4);stem(n,d);
xlabel('n');xlim([0 5]);ylabel('D');
subplot(4,2,5);stem(n,e);
xlabel('n');xlim([0 5]);ylabel('E');
subplot(4,2,6);stem(n,f);
xlabel('n');xlim([0 5]);ylabel('F');
subplot(4,2,7);stem(n,g);
xlabel('n');xlim([0 5]);ylabel('G');
（2）用MATLAB实现下列序列：
a) x(n)=0.8n0≤n≤15
b) x(n)=e(0.2+3j)n0≤n≤15
c) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15
d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。

e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。

clear all;
N=0:15;
% a) x(n)=0.8n 0≤n≤15
xa=0.8.^N;
figure;subplot(2,1,1);stem(N,xa); xlabel('n');xlim([0 16]);ylabel('xa');
% b) x(n)=e(0.2+3j)n 0≤n≤15
xb=exp((0.2+3*j)*N);
subplot(2,1,2);stem(N,xb);
xlabel('n');xlim([0 16]);ylabel('xb');figure;
% c) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15
xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);
subplot(3,1,1);stem(N,xc);xlabel('n');xlim([0 16]);ylabel('xc');
% d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。

k=0:3;m=0;
for i=1:4
for j=1:16
m=m+1;
n(m)=N(j)+16*k(i);
x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi);
end
end
subplot(3,1,2);stem(n,x16);xlabel('n');ylabel('x16');
% e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。

for j=1:10
x10(j)=x16(j);
end
for i=1:3
for m=1:10
x10(i*10+m)=x10(m);
end
end
n=1:40;
subplot(3,1,3);stem(n,x10); xlabel('n');ylabel('x10');
（3）x(n)=[1,-1,3,5],产生并绘出下列序列的样本: a) x 1(n)=2x(n+2)-x(n-1)-2x(n)
b) ∑=-=5
1k 2)k n (nx (n) x
clear all n=1:4; T=4;
x=[1 -1 3 5]; x(5:8)=x(1:4);
subplot(2,1,1);stem(1:8,x);grid;
for i=1:4 if i-1<0
x1(i)=2*x(i+2)-x(i-1)-2*x(i); else
x1(i)=2*x(i+2)-x(i-1+T)-2*x(i); end end
x1(5:8)=x1(1:4);
subplot(2,1,2);stem(1:8,x1);grid;
（4）绘出下列时间函数的图形，对x 轴、y 轴以及图形上方均须加上适当的标注： a) x(t)=sin(2πt) 0≤t ≤10s
b) x(t)=cos(100πt)sin(πt) 0≤t ≤4s
ta=0:0.05:10; xa=sin(2*pi*ta);
subplot(2,1,1);plot(ta,xa); xlabel('t');ylabel('幅度'); tb=0:0.01:4;
xb=cos(100*pi*tb).*sin(pi*tb); subplot(2,1,2);plot(tb,xb); xlabel('t');ylabel('幅度');
（5）编写函数stepshift(n0,n1,n2)实现u(n-n0),n1<n0<n2,绘出该函数的图形，起点为n1，
终点为n2。

n0=5;ns=1;nf=10;%ns 为起点;nf 为终点；在=n=n0处生成单位阶跃序列 n=[ns:nf];
x=[(n-n0)>=0]; stem(n,x);
（6）给一定因果系统)0.9z 0.67z -1)/(1z 2(1H(z)-2-1-1+++=求出并绘制H(z)的幅频响应
与相频响应。

clear all;
b=[1,sqrt(2),1]; a=[1,-0.67,0.9]; [h,w]=freqz(b,a);
am=20*log10(abs(h)); subplot(2,1,1);plot(w,am);
ph=angle(h);
subplot(2,1,2);plot(w,ph);
（7）计算序列{8 -2 -1 2 3}和序列{2 3 -1 -3}的离散卷积，并作图表示卷积结果。

clear all;
a=[8 -2 -1 2 3];
b=[2 3 -1 -3];
c=conv(a,b); %计算卷积
M=length(c)-1;
n=0:1:M;
stem(n,c);
xlabel('n');ylabel('幅度');
（8）求以下差分方程所描述系统的单位脉冲响应h(n),0≤n≤50
y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1)
clear all;
N=50;
a=[1 -2];
b=[1 0.1 -0.06];
x=[1 zeros(1,N-1)];
k=0:1:N-1;
y=filter(a,b,x);
stem(k,y);
xlabel('n');ylabel('幅度');
实验二快速Fourier变换(FFT)及其应用
一、实验目的
1．在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉FFT子程序。

2．熟悉应用FFT对典型信号进行频谱分析的方法。

3. 了解应用FFT进行信号频谱分析过程中可能出现的问题以便在实际中正确应用FFT。

4．熟悉应用FFT实现两个序列的线性卷积的方法。

5．初步了解用周期图法作随机信号谱分析的方法。

二、实验原理与方法
在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。

这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：
反变换为：
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。

FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。

它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。

常用的FFT 是以2为基数的，其长度。

它的效率高，程序简单，使用非常方便，当要变换的序
列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。

(一)、在运用DFT进行频谱分析的过程中可能产生三种误差：
(1) 混叠
序列的频谱时被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。

避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。

(2) 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。

泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。

为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。

(3) 栅栏效应
DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。

减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖
桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。

(二)、用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。

在一定的条件下，可以使圆周卷积等于线性卷积。

一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT 的长度
N≥N1＋N2
对于长度不足N的两个序列，分别将他们补零延长到N。

当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。

有两种方法：
重叠相加法。

将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。

重叠保留法。

这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。

(三)、用周期图法(平滑周期图的平均法)对随机信号作谱分析
实际中许多信号往往既不具有有限能量，由非周期性的。

无限能量信号的基本概念是随机过程，也就是说无限能量信号是一随机信号。

周期图法是随机信号作谱分析的一种方法，它特别适用于用FFT直接计算功率谱的估值。

将长度为N的实平稳随机序列的样本x(n)再次分割成K段，每段长度为L,即L=N/K。

每段序列仍可表示为：
xi(n)=x(n+(i-1)L)，0≤n≤L-1，1≤i≤K
但是这里在计算周期图之前，先用窗函数w(n)给每段序列xi(n)加权，K个修正的周期图定义为
其中U表示窗口序列的能量，它等于
在此情况下，功率谱估计量可表示为
三、实验内容及步骤
实验中用到的信号序列：
a) Gaussian序列
b) 衰减正弦序列
c) 三角波序列
d) 反三角波序列
上机实验内容：
(1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。

% %定义高斯序列
% function [Xa,Fa] =gauss(p,q)
% n=[0:15];
% Xa(n+1)=exp(-(n+1-p).^2./q);
% F=fft(Xa);
% Fa=abs(F);
clear all;
%%%%%%% p=8,q=2 %%%%%%%%%%%%
[Xa1,Fa1]= gauss(8,2);
k=0:15;
subplot(5,2,1);plot(k,Xa1);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=2');
subplot(5,2,2);plot(k,Fa1);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=2');
%%%%%%% p=8,q=4 %%%%%%%%%%%%
[Xa2,Fa2]= gauss(8,4);
subplot(5,2,3);plot(k,Xa2);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=4');
subplot(5,2,4);plot(k,Fa2);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=4');
%%%%%%% p=8,q=8 %%%%%%%%%%%%
[Xa3,Fa3]= gauss(8,8);
subplot(5,2,5);plot(k,Xa3);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=8');
subplot(5,2,6);plot(k,Fa3);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=8');
%%%%%%% p=13,q=8 %%%%%%%%%%
[Xa4,Fa4]= gauss(13,8);
subplot(5,2,7);plot(k,Xa4);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=13,q=8');
subplot(5,2,8);plot(k,Fa4);
xlabel('n');ylabel('幅频特性');text(8,3,'p=13,q=8');
%%%%%%% p=14,q=8 %%%%%%%%%%
[Xa5,Fa5]= gauss(14,8);
subplot(5,2,9);plot(k,Xa5);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=14,q=8');
subplot(5,2,10);plot(k,Fa5);
xlabel('n');ylabel('幅频特性');text(8,3,'p=14,q=8');
(2)、观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。

% 定义衰减正弦序列
% function [Xb,Fb] = downsin(a,f)
% n=[0:15];
% Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536
% F = fft(Xb);
% Fb=abs(F);
clear all;
k=0:15;
%%%%%%%%% a=0.1,f=0.0.0625 %%%%%%%%%%%
[Xb,Fb]=downsin(0.1,0.0625);
subplot(3,2,1); plot(k,Xb);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.0625');
subplot(3,2,2); plot(k,Fb);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.0625');
%%%%%%%%% a=0.1,f=0.4375 %%%%%%%%%%%
[Xb1,Fb1]=downsin(0.1,0.4375);
subplot(3,2,3); plot(k,Xb1);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.4375');
subplot(3,2,4); plot(k,Fb1);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.4375');
%%%%%%%%% a=0.1,f= 0.5625 %%%%%%%%%%
[Xb2,Fb2]=downsin(0.1,0.5625);
subplot(3,2,5); plot(k,Xb2);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.5625');
subplot(3,2,6); plot(k,Fb2);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.5625');
(3)、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。

在xc(n)和xd(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？
clear all;
n=[0:3];k=[1:8];
%定义三角波序列
Xc(n+1) = n;Xc(n+5) =4-n;
%定义反三角波序列
Xd(n+1) = 4-n;Xd(n+5) =n;
%%%%%%%%%%% 三角波特性%%%%%%%%%%%%%
subplot(2,2,1);plot(k-1,Xc);
xlabel('n');ylabel('时域特性');text(1,3,'三角波');
subplot(2,2,2);plot(k-1,abs(fft(Xc)));
xlabel('k');ylabel('幅频特性');text(4,10,'三角波'); %%%%%%%%%%% 反三角波特性%%%%%%%%%% subplot(2,2,3);plot(k-1,Xd);
xlabel('n');ylabel('时域特性');text(3,3,'反三角波');
subplot(2,2,4);plot(k-1,abs(fft(Xd)));
xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');
%末尾补0，计算32点FFT
Xc(9:32)=0;Xd(9:32)=0;k=1:32;figure; %%%%%%%%%%% 三角波特性%%%%%%%%%%%%% subplot(2,2,1);plot(k-1,Xc);
xlabel('n');ylabel('时域特性');text(1,3,'三角波');
subplot(2,2,2);plot(k-1,abs(fft(Xc)));
xlabel('k');ylabel('幅频特性');text(4,10,'三角波'); %%%%%%%%%%% fan三角波特性%%%%%%%%%% subplot(2,2,3);plot(k-1,Xd);
xlabel('n');ylabel('时域特性');text(3,3,'反三角波');
subplot(2,2,4);plot(k-1,abs(fft(Xd)));
xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');
(4)、一个连续信号含两个频率分量，经采样得
x(n)=sin2π*0.125n+cos2π*(0.125+Δf)n n=0,1……,N-1
已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？
clear all;
%%%%%%% N = 16 %%%%%%%%%%%%
N=16;detf=1/16;n=[0:N-1];
x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
%%%%%%% N = 16,detf = 1/16 %%%%%%%%%%%%
subplot(2,2,1);stem(n,x1);hold;plot(n,x1);
xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/16');
subplot(2,2,2);stem(n,abs(fft(x1)));
xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16');
%%%%%%% N = 16,detf = 1/64 %%%%%%%%%%%%
subplot(2,2,3);stem(n,x2);
xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/64');
subplot(2,2,4);stem(n,abs(fft(x2)));
xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/64');
%%%%%%% N = 128 %%%%%%%%%%%%
N=128;detf=1/16;n=[0:N-1];
x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
detf = 1/64;
x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
%%%%%%% N = 128,detf = 1/16 %%%%%%%%%%%%
figure;subplot(2,2,1);stem(n,x3);
xlabel('n');ylabel('时域特性');axis([0 128 -2 2]);text(6,1.5,'N=128,detf=1/16');
subplot(2,2,2);stem(n,abs(fft(x3)));
xlabel('n');ylabel('幅值特性');axis([0 128 -10 70]);text(40,60,'N=128,detf=1/16'); %%%%%%% N = 128,detf = 1/64 %%%%%%%%%%%%
subplot(2,2,3);stem(n,x3);
xlabel('n');ylabel('时域特性');axis([0 128 -2 2]);text(6,1.5,'N=128,detf=1/16'); subplot(2,2,4);stem(n,abs(fft(x4)));
xlabel('n');ylabel('幅值特性');axis([0 128 -10 70]);text(40,60,'N=128,detf=1/16');
(5)、用FFT分别实现xa(n)（p＝8，q＝2）和xb(n)（a＝0.1，f＝0.0625）的16点圆周卷积和线性卷积。

clear all;
N=16;
n=0:N-1;
p=8;q=2;
Xa(n+1)=exp(-(n-p).^2./q);
a=0.1;f=0.0625;
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);
%16点循环卷积
Fa=fft(Xa); Fb=fft(Xb);
Fx=Fa.*Fb;
X51=ifft(Fx);
stem(n,X51);
%16点线性卷积
Xa(N+1:2*N-1)=0;
Xb(N+1:2*N-1)=0;
Fa=fft(Xa); Fb=fft(Xb);
Fc=Fa.*Fb;
X52=ifft(Fc);
figure;stem(1:2*N-1,X52);
（7）用FFT分别计算x a(n)(p=8,q=2)和x b(n)(a=0.1,f=0.0625)的16点循环相关和线性相关，问一共有多少种结果，他们之间有何异同点。

clear all;
N=16;
n=0:N-1;
p=8;q=2;
Xa(n+1)=exp(-(n-p).^2./q);
a=0.1;f=0.0625;
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);
N=length(Xa);
%16点循环相关
Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);
Fx=conj(Fa).*Fb;
X71=real(ifft(Fx));
X71=[X71(N+2:2*N) X71(1:N)];
n=(-N+1):(N-1); stem(n,X71);
%16点线性相关
Xa(N+1:2*N-1)=0;
Xb(N+1:2*N-1)=0;
Fa=fft(Xa); Fb=fft(Xb);
Fc=conj(Fa).*Fb;
X72=real(ifft(Fc));
figure;stem(1:2*N-1,X72);
（8）用FFT分别计算x a(n)(p=8,q=2)和x b(n)(a=0.1,f=0.0625)的自相关函数。

clear all;
N=16;
n=0:N-1;
p=8;q=2;
Xa(n+1)=exp(-(n-p).^2./q);
a=0.1;f=0.0625;
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536
N=length(Xa);
% Xa(n) 16点自相关
Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);
F1=conj(Fa).*Fa;
X81=real(ifft(F1));
X81=[X81(N+2:2*N) X81(1:N)];
n=(-N+1):(N-1);
subplot(2,1,1);stem(n,X81);
xlabel('n'); ylabel('幅度');
% Xb(n) 16点自相关
Fb=fft(Xb,2*N);
F2=conj(Fb).*Fb;
X82=real(ifft(F2));
X82=[X82(N+2:2*N) X82(1:N)];
% n=(-N+1):(N-1);
subplot(2,1,2);stem(n,X82);
xlabel('n'); ylabel('幅度');
实验三 IIR 数字滤波器的设计
(1)kHz f c 3.0=，dB 8.0=δ，kHz f r 2.0=，dB At 20=,ms T 1=;设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。

分析：f=200Hz 时阻带衰减大于30dB ，通过修改axis([0,fs/2,-80,10])为axis([200,fs/2,-1,1]) 发现通带波动rs 满足<0.8。

bz =[0.0262 -0.1047 0.1570 -0.1047 0.0262] az =[1.0000 1.5289 1.6537 0.9452 0.2796] 系统函数为：
43214
3212796.09452.06537.15289.110262.01047.01570.01047.0-0262.0)(H --------+++++-+=
z z z z z z z z z
（2）kHz f c 2.0=，dB 1=δ，kHz f r 3.0=，dB At 25=,ms T 1=;分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和
衰减量，检查是否满足要求。

比较这两种方法的优缺点。

bz1 =[0.0000 0.0002 0.0153 0.0995 0.1444 0.0611 0.0075 0.0002 0.0000 0] az1 =[1.0000 -1.9199 2.5324 -2.2053 1.3868 -0.6309 0.2045 -0.0450 0.0060 -0.0004] 因此脉冲响应不变法的系统函数为:
1234567
123456789
0.00020.01530.09950.14440.06110.00750.0002()1 1.9199 2.5324 2.2053 1.38680.63090.20450.04500.00600.0004imp z z z z z z z H z z z z z z z z z z ----------------+++---=
-+-+-+-+-
bz2 =[0.0179 0.1072 0.2681 0.3575 0.2681 0.1072 0.0179] az2 =[1.0000 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.0025] 因此双线性变换法的系统函数为:
123456
123456
0.01790.10720.26810.35750.26810.10720.0179()10.60190.91300.29890.15010.02080.0025bil z z z z z z H z z z z z z z
------------++++++=-+-+-+ 分析：
脉冲响应不变法的N=9，双线性变换法的N=6，由图知它们都满足要求，但脉冲响应的衰减较快，双线性变换的过渡带窄一些，且阶数比脉冲小，容易实现。

实验四 FIR 数字滤波器的设计
一、实验目的
1. 掌握用窗函数法，频率采样法及优化设计法设计FIR 滤波器的原理及方法，熟悉响
应的计算机编程；
2. 熟悉线性相位FIR 滤波器的幅频特性和相频特性；
3. 了解各种不同窗函数对滤波器性能的影响。

二、实验原理与方法
线性相位实系数FIR 滤波器按其N 值奇偶和h(n)的奇偶对称性分为四种： 1、h(n)为偶对称，N 为奇数
H(ejω)的幅值关于ω=0，π，2π成偶对称。

2、h(n)为偶对称，N为偶数
H(ejω)的幅值关于ω=π成奇对称，不适合作高通。

3、h(n)为奇对称，N为奇数
H(ejω)的幅值关于ω=0，π，2π成奇对称，不适合作高通和低通。

4、h(n)为奇对称，N为偶数
H(ejω) ω=0、2π＝0，不适合作低通。

(一) 窗口法
窗函数法设计线性相位FIR滤波器步骤
确定数字滤波器的性能要求：临界频率{ωk}，滤波器单位脉冲响应长度N；
根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ejω)的幅频特性和相频特性；
求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ejω)按M(M远大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)；
选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应；求H(ejω)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。

窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。

W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有：
矩形窗w(n)=RN(n)；
Hanning窗；
Hamming窗；
Blackmen窗；
Kaiser窗。

式中Io(x)为零阶贝塞尔函数。

（二）频率采样法
频率采样法是从频域出发，将给定的理想频率响应Hd(ejω)加以等间隔采样
然后以此Hd(k)作为实际FIR数字滤波器的频率特性的采样值H(k),即令
由H(k)通过IDFT可得有限长序列h(n)
将上式代入到Z变换中去可得
其中Φ(ω)是内插函数
（三）FIR滤波器的优化设计
FIR滤波器的优化设计是按照最大误差最小化准则，使所设计的频响与理想频响之间的最大误差，在通带和阻带范围均为最小，而且是等波动逼近的。

为了简化起见，在优化设计中一般将线性相位FIR滤波器的单位脉冲响应h(n)的对称中心置于n=0处，此时，线性相位因子α=0。

当N为奇数，且N=2M+1，则
如希望逼近一个低通滤波器，这里M，ωp和ωs固定为某个值。

在这种情况下有
定义一逼近误差函数：
E(ω)为在希望的滤波器通带和阻带内算出的误差值，W(ω)为加权函数。

k应当等于比值δ1/δ2，δ1为通带波动，δ2为阻带波动。

在这种情况下，设计过程要求
|E(ω)|在区间的最大值为最小，它等效于求最小δ2。

根据数学上多项式逼近连续函数的理论，用三角多项式逼近连续函数，在一定条件下存在最佳逼近的三角多项式，而且可以证明这个多项式是唯一的。

这一最佳逼近定理通常称作交替定理。

在逼近过程中，可以固定k，M，ωp，ωs而允许改变δ2，按照交替定理，首先估计出(M+2)个误差函数的极值频率点{ωi}，i=0,1,...,M+1，共计可以写出(M+2)个方程
式中ρ表示峰值误差。

一般仅需求解出ρ，接着便可用三角多项式找到一组新的极值频率点，并求出新的峰值误差ρ。

依此反复进行，直到前、后两次ρ值不变化为止，最小的ρ即为所求的δ2。

这一算法通常称作雷米兹(Remez)交替算法。

三、实验内容及步骤
(1)N=15, 。

用Hanning窗设计一线性相位带通滤波器，观察它
的实际3dB和20dB带宽。

N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N变化的影响；
(2)分别改用矩形窗和Blackman窗，设计(1)中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点；。