数字信号处理 实验报告

1.DFT在信号频谱分析中的应用

1.1设计目的

（1） 熟悉DFT的性质。

（2） 加深理解信号频谱的概念及性质。 （3）了解高密度谱与高分辨率频谱的区别。 1.2设计任务与要求

（1）学习用DFT和补零DFT的方法来计算信号的频谱。

（2）用MATLAB语言编程来实现，在做课程设计前，必须充分预习课本DTFT、DFT及补零DFT的有关概念，熟悉MATLAB语言，编写程序。

1.3设计原理

所谓信号的频谱分析就是计算信号的傅里叶变换。连续信号与系统的傅里叶分析显然不便于直接用计算机进行计算，使其应用受到，而DFT是一种时域和频域均离散化的变换，适合数值运算，成为分析离散信号和系统的有力工具。

工程实际中，经常遇到的连续信号Xa(t),其频谱函数Xa(jW)也是连续函数。数字计算机难于处理，因而我们采用DFT来对连续时间信号的傅里叶变换进行逼近，进而分析连续时间信号的频谱。

1.4设计内容

1.4.1用MATLAB实现DFT与IDTF (1)点序列x(n) 的DFT为：

N1Xkxnwn0knN,wNej2N

点序列x(n) 的IDFT为：

N1xnn0X(k)wknN,wNej2N

(2)N点DFT的矩阵为：

(3)根据DFT公式与矩阵展开，通过MATLAB实现DFT与IDFT，程序如下：

1

Matlab中的内部函数文件fft.m文件:

unction [varargout] = fft(varargin) if nargout == 0

builtin('fft', varargin{:});

else

[varargout{1:nargout}] = builtin('fft', varargin{:}); end

运算量估计：对于N=2M

点序列进行时间抽选奇偶分解FFT计算，需分M级，每级计

N2算N/2个蝶。每一级需N/2次复乘、N次复加，因此总共需要进行：复乘：复加：NMNlog2MN2log2N ；

N直接计算N点的DFT，需要N次复乘、N(N-1)次复加。N值越大，

2时间抽选奇偶分解FFT算法越优越。例如当N=2048点时，时间抽选奇偶分解FFT算法比直接计算DFT速度快300多倍。

1.4.2. 对离散确定信号 x(n)cos(0.48n)cos(0.52n) 作如下谱分析：

（1）截取x(n)使x(n)成为有限长序列N(0nN-1)，(长度N自己选)写程序计算出

x(n)的N点DFT X(k),并画出相应的幅频图X(k)~k。

N=32;

n=0:1:N-1;

xn=cos(0.48*pi*n)+cos(0.52*pi*n); subplot(3,1,1)

stem(n,xn,'.k');

title('时域序列图xn'); xlabel('n');

axis([0,10,-2.5,2.5]); w=2*pi*(0:1:2047)/2048; Xw=xn*exp(-j*n'*w); subplot(3,1,2);

plot(w/pi,abs(Xw));

title('幅频特性曲线X(ejw)'); xlabel('w');

axis([0,1,0,20]);

2

subplot(3,1,3)

Xk=dft(xn,N);

k1=0:1:31;w1=2*pi/32*k1; stem(w1/pi,abs(XK),'.k'); title('频域序列图XK'); xlabel('频率（单位：pi）'); axis([0,1,0,20]);

图1.4.1

由图1可见，由于截断函数的频谱混叠作用，X(k)不能正确分辨w1=0.48π、w2=0.52π这两个频率分量。

（2）将 (1)中x(n)补零加长至M点（长度M自己选），编写程序计算x(n)的M点DFT

X1(k) ,并画出相应的图X1(k)~k，并利用补零DFT计算 (1)中N点有限长序列x(n)频

谱X(ej)并画出相应的幅频图X(ej)~。

N=32; n=0:N-1;

xn=cos(0.48*pi*n)+cos(0.52*pi*n); N1=128;

n1=0:N1-1;

x1=[xn(1:32) zeros(1,96)]; subplot(3,1,1) stem(n1,x1,'.k'); title('时域序列图x1'); xlabel('n');

axis([0,100,-2.5,2.5]); w=2*pi*(0:2047)/2048; X1=x1*exp(-j*n1'*w); subplot(3,1,2);

plot(w/pi,abs(X1));

3

title('幅频特性曲线X(ejw)'); xlabel('w');

axis([0,1,0,20]); subplot(3,1,3); k=0:1:N1-1;

WN=exp(-j*2*pi/N1); nk=n'*k; WNnk=WN.^nk; XK=xn*WNnk; k1=0:1:127; w1=2*pi/128*k1;

stem(w1/pi,abs(XK(1:1:128)),'.k'); title('频域序列图XK');

xlabel('频率（单位：pi）'); axis([0,1,0,20]);

图1.4.2

x(n)补零至点对应的x(n)、X(ejw)、X(k)如图2所示。由图可见，x(n)补零至点，只是改变X(k)的密度，截断函数的频谱混叠作用没有改变，这时的物理分辨率使X(k)仍不能正确分辨w1=0.48π、w2=0.52π这两个频率分量。这说明，补零仅仅是提高了计算分辨率，得到的是高密度频谱，而得不到高分辨率谱。

1.5研究高密度谱与高分辨率频谱。

333对连续确定信号xa(t)cos(26.510t)cos(2710t)cos(2910t)

以采样频率fs=32kHz对信号xa(t)采样得离散信号x(n)，分析下列三种情况的幅频特性。 (1)采集数据x(n)长度取N=16点，编写程序计算出x(n)的16点DFTX(k),并画出相应的幅频图。

4

N=16;%length of sequence

n=0:1:N-1;%time sample T=1/32000; t=n*T;

xn=cos(2*pi*6500*t)+cos(2*pi*7000*t)+cos(2*pi*9000*t);% subplot(3,1,1) stem(n,xn,'.k'); title('时域序列图x1'); xlabel('n');

axis([0,100,-2.5,2.5]); w=2*pi*(0:2047)/2048; X1=xn*exp(-j*n'*w); subplot(3,1,2);

plot(w/pi,abs(X1));

title('幅频特性曲线X(ejw)'); xlabel('w'); axis([0,1,0,20]); subplot(3,1,3); k=0:1:N-1;

WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; XK=xn*WNnk; k1=0:1:15;

w1=2*pi/16*k1;

stem(w1/pi,abs(XK(1:1:16)),'.k'); title('频域序列图XK'); xlabel('频率（单位：pi）'); axis([0,1,0,20]);

图1.5.1

5

(2)采集数据x(n)长度N=16点，补零加长至M点(长度M自己选)，利用补零DFT计算 x(n)的频谱X(ej)并画出相应的幅频图X1(ej)~。选取M=128 N=16;%length of sequence n=0:1:N-1;%time sample T=1/32000; t=n*T;

xn=cos(2*pi*6500*t)+cos(2*pi*7000*t)+cos(2*pi*9000*t);% N1=128; n1=0:N1-1;

x1=[xn(1:16) zeros(1,112)]; subplot(3,1,1) stem(n1,x1,'.k'); title('时域序列图x1'); xlabel('n');

axis([0,100,-2.5,2.5]); w=2*pi*(0:2047)/2048; X1=x1*exp(-j*n1'*w); subplot(3,1,2); plot(w/pi,abs(X1));

title('幅频特性曲线X(ejw)'); xlabel('w'); axis([0,1,0,20]); subplot(3,1,3); k=0:1:N1-1;

WN=exp(-j*2*pi/N1); nk=n'*k;

WNnk=WN.^nk; XK=xn*WNnk;

k1=0:1:127;

w1=2*pi/128*k1;

stem(w1/pi,abs(XK(1:1:128)),'.k'); title('频域序列图XK'); xlabel('频率（单位：pi）'); axis([0,1,0,20]);

6

图1.5.2

(3) 采集数据x(n)长度取为M点（注意不是补零至M），编写程序计算出M点采集数据

x(n)的的频谱X(ej)并画出相应的幅频图X2(ej)~。

N=128;%length of sequence n=0:1:N-1;%time sample T=1/32000;

t=n*T;

xn=cos(2*pi*6500*t)+cos(2*pi*7000*t)+cos(2*pi*9000*t);% subplot(3,1,1) stem(n,xn,'.k');

title('时域序列图x1'); xlabel('n');

axis([0,100,-2.5,2.5]); w=2*pi*(0:2047)/2048; X1=xn*exp(-j*n'*w);

subplot(3,1,2);

plot(w/pi,abs(X1));

title('幅频特性曲线X(ejw)'); xlabel('w'); axis([0,1,0,70]); subplot(3,1,3); k=0:1:N-1;

WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; XK=xn*WNnk;

7

k1=0:1:127;

w1=2*pi/128*k1;

stem(w1/pi,abs(XK(1:1:128)),'.k'); title('频域序列图XK');

xlabel('频率（单位：pi）'); axis([0,1,0,70]);

图1.5.3

1.6.思考题及解答

(1)对比设计内容2中图1.5.1、图1.5.2、图1.5.5，说明补零DFT的作用。

由图1-图1.5.1、图1.5.2、图1.5.52可知DFT是有限长序列的频谱等间隔采样所得到的样本值，这就相当于透过一个栅栏去观察原来信号的频谱，因此必然有一些地方被栅栏所遮挡，这些被遮挡的部分就是未被采样到的部分，这种现象称为栅栏效应。如下图

图1.6.1

由于栅栏效应总是存在的，因而可能会使信号频率中某些较大的频率分量由于被“遮挡”二无法得到反映。此时，通常在有限长序列的尾部增补若干个零值，介意改变原序列的长度。这样加长的序列作DFT时，由于点数增加就相当于调整了原来栅栏的间隙即间隔频率，可以使得原来的不到反映的那些较大的频率分量落在采样点上而得到反映。但要注意，由于栅栏效应，使得被分析的频谱变得较为稀疏，为此，在采样样本序列x(n)后补零，在数据长度Tp不变的情况下，可以改变频谱的频率取样密度，得到高密度频谱。

(2)解释设计内容3中

X1(ej)~图和

X2(ej)~图有什么区别？补零DFT能否提高

信号的频谱分辨率，说明提高频谱密度、频谱分辨率的措施各是什么？

8

图1.5.2在图1.5.1的基础上提高了频谱密度和计算分辨率，但是频谱的包络没有发生变化所以

X1(ej)~图和

X2(ej)~图没有区别。补零作用不能提高信号的频谱分辨率，

因在x(n)后面补零并没有增加新的信息量，改善的仅是栅栏效应，所以补零是不能提高频率分辨率的，即得不到高分辨率谱。这说明，补零仅仅是提高了计算分辨率，得到的是高密度频谱，而要得到高分辨率谱，则要通过增加数据的记录长度Tp来提高物理分辨率。

2. 有噪声情况下信号幅度谱的研究

2.1设计目的

（1）了解并掌握白噪声的产生方法。

（2）了解并掌握正弦信号及白噪声信号的相关函数求法。 2.2设计任务与要求

（1）仿真在正弦信号加白噪声情况下，求其幅度谱及相关函数。

（2）用MATLAB语言编程来实现，在做课程设计前，应查阅信号去噪的相关理论知识，熟悉MATLAB语言，编写程序。 2.3.设计内容

3.1.编写产生均匀分布白噪声序列的M函数文件drand.m 。 A=6; x0=1; M=255; f=2;

N=100; %初始化； x0=1; M=255;

for k=1: N %乘同余法递推100次；

x2=A*x0; %分别用x2和x0表示xi+1和xi-1；

x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（xi）中；

v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中； v(:,k)=(v1-0.5 )*f; %将v1中的数（）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化； x0=x1; % xi-1= xi；

v0=v1;

end %递推100次结束； v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中； k1=k;

%grapher %以下是绘图程序； k=1:k1;

plot(k,v,k,v,'b'); xlabel('k');

9

ylabel('v');

title(' (-1,+1)均匀分布的白噪声')

图2.3.1

2.3.2.编写计算序列x(n) 正弦信号加白噪声的自相关序列的M函数文件dcor.m。 （1）设定正选信号的频率为10HZ，抽样频率为100HZ； （2）设定N(0,0.25)高斯白噪声，及噪声功率为0.25W； （3）最后将噪声叠加到正弦信号上，观察其三者时域波形。 fs=100; fc=10;

x=(0:1/fs:2);

n=201;

y1=sin(2*pi*fc*x); %原正弦信号，频率为10 a=0;b=0.5; %均值为a，方差为b^2 subplot(3,2,1); plot(x,y1,'b');

title('y=sin(20pi*x)'); ylabel('y');

xlabel('x/20pi');

y2=a+b*randn(1,n); %高斯白噪声 subplot(3,2,2); plot(x,y2,'b');

title('N(0，0.25)的高斯白噪声'); ylabel('y');

xlabel('x/20pi');

y=y1+y2; %加入噪声之后的信号 subplot(3,2,3);

plot(x,y,'b');

title('叠加了高斯白噪声的sinx'); ylabel('y');

10

xlabel('x/20pi');

2.3.3.编写m程序文件，分析含噪信号的相关函数及功率密度谱。

首先对原正弦信号直接进行FFT，得出其频谱；其次求白噪声的自相关函数，随机序列自相关函数的无偏估计公式为：

^Nm1（1）rxx(m)1Nmx(n)x(nm) 0mN1

n0^^（2）rxx(m)rxx(m)

0mN1

FY=fft(y); %傅里叶变换得出频谱函数 FY1=fftshift(FY); %频谱校正 f=(0:200)*fs/n-fs/2; subplot(3,2,4);

plot(f,abs(FY1),'b'); title('函数频谱图'); ylabel('F(jw)'); xlabel('w');

%求高斯白噪声的自相关函数 m=50;

i=-0.49:1/fs:0.49; for j=1:m

R(j)=sum(y2(1:n-j-1).*y2(j:199),2)/(n-j);%无偏自相关函数的估计 Rx(49+j)=R(j); Rx(51-j)=R(j); end

subplot(3,2,5); plot(i,Rx,'b');

title('白噪声自相关函数图'); ylabel('Rx'); xlabel('x');

Fy2=fft(Rx); %傅里叶变换得出白噪声功率谱函数 Fy21=fftshift(Fy2); %功率谱校正 f=(0:98)*fs/99-fs/2; subplot(3,2,6);

plot(f,abs(Fy21),'b'); axis([-50 50 -0.5 1]); title('白噪声功率谱函数图'); ylabel('F(Rx)'); xlabel('w');

11

图2.3.2.1

2.4.简要回答如下思考题： （1）分析白噪声的特点，白噪声有哪些主要参数?如何调整随机序列的输出平均率和平均值? 白噪声是一种功率频谱密度为常数的随机信号或随机过程。换句话说，此信号在各个频段上的功率是一样的。相对的，其他不具有这一性质的噪声信号被称为有色噪声。 󰀀󰀀理想的白噪声具有无限带宽，因而其能量是无限大，这在现实世界是不可能存在的。实际上，我们常常将有限带宽的平整讯号视为白噪音，因为这让我们在数学分析上更加方便。然而，白噪声在数学处理上比较方便，因此它是系统分析的有力工具。一般，只要一个噪声过程所具有的频谱宽度远远大于它所作用系统的带宽，并且在该带宽中其频谱密度基本上可以作为常数来考虑，就可以把它作为白噪声来处理。例如，热噪声和散弹噪声在很宽的频率范围内具有均匀的功率谱密度，通常可以认为它们是白噪声。白噪声的主要参数为均值和方差，其中方差代表其平均功率。

(2) 计算正弦信号的平均功率、功率密度谱和自相关函数?当截取的点数N不为正弦信号周期的整数倍时，会有什么结果?

ˆf一个长L的信号xLn的PSD的周期图估计是PxxXLffsL2，这里XLf运用的是

L1matlab里面的fft的定义不带归一化系数，所以要除以L，其中XLfn0xLne2jfn/fs实际对XLf的计算可以只在有限的频率点上执行并且使用FFT。大多数周期图法的应用

ˆf都计算N点PSD估计PxxkL1XLfk2fsL，fkkfsN,k0,1,,N1，其中

XLfkn0xLne2jkn/N。选择N是大于L的下一个2的幂次是明智的，要计算XLfk我

12

们直接对xLn补零到长度为N。假如L>N，在计算XLfk前，我们必须绕回xLn模N。 （3）设计中应当如何正确选择点数?

采样点数可以采用N= length(x)；来取，x是采样数据；󰀀采样频率fs = 1/Ts 即采样时间的倒数，也就是采样信号中两个数据点的时间间隔的倒数；󰀀采样频率一定时，采样点数越多越好，即采样时间越长越好，这样fs/N就越小，频域的频率分辨率越大，FFT结果就越准确，采样点数最好是2的整数次幂，可以加快FFT运算。在实际应用时，由于受内存计算等的要求，采样点数满足FFT计算的一定精度要求就行了，不必太多。补零可以使N满足2的整数次幂，可以增加频谱谱线的光滑性，但不能太多。

3. 离散时间系统频域分析

3.1设计目的

(1)学习离散时间系统频率特性的计算方法。

(2)深刻理解离散时间系统频率特性与滤波特性的关系。

(3)掌握离散时间系统的系统参数、系统零极点及系统频率特性间的关系。 3.2设计任务与要求 (1)用MATLAB语言编程分析数字滤波器的各种滤波特性及其与滤波特性相关的参数。 (2)在做课程设计前，必须充分预习相关理论知识，熟悉MATLAB语言，编写程序。 3.3设计原理

一个时域离散系统或网络的表示方法有三种：

MNii1． 差分方

y(n)bx(ni)ai0My(ni)

i12． 系统函数H(z)Y(z)X(z)i0bizNi

i11ai1zi3． 单位脉冲响应 h(n)ZT3.4.设计内容

[H(z)]

(1) 用MATLAB语言编写计算N阶差分方程所描述系统频响函数H(efr.m。

function [H]=fr(b,a,w);

%计算N阶差分方程所描述系统频响函数 %w 采样频率

%b 系统函数H(z)的分子项(对FIR，b=h)

%a 系统函数H(z)的分母项(对FIR，a=1)

m=0:length(b)-1; %length是返回矩阵最长维的的长度 l=0:length(a)-1;

num=b*exp(-j*m'*w);%exp是以e为底的指数函数,m'是m的转置

j)的m函数文件

13

den=a*exp(-j*l'*w);

H=num./den; %./是向量右除

(2) 根据频响特性与系统零极点的关系，自己构造一个N阶差分方程，使该差分方程为 数字低通滤波器。利用MATLAB程序画出相应的幅频图H(ej)~。 n=0:0.01:2;

N=5

[z,p,k]=buttap(N); [b,a]=zp2tf(z,p,k); [H,w]=freqs(b,a,n); magH=(abs(H)).^2; plot(w,magH); axis([0 2 0 1]); xlabel('w/wc');

ylabel('|H(jw)|^2');

title('Butterworth analog filter prototype');

图3.4.1 低通滤波器

(3)改变2.中差分方程的系数，使该差分方程分别为数字高通及全通滤波器。利用MATLAB程序画出相应的幅频图H(ej)~。 N=5; Ws=50; Fs=100;

[b,a]=butter(N,Ws/Fs,'high'); [z,p,k]=butter(N,Ws/Fs,'high'); [H,W]=freqz(b,a);

plot(W*Fs/(2*pi),abs(H));grid; xlabel('Hz');

14

ylabel('Db');

N=5; Ws=50; Fs=100;

[b,a]=butter(N,Ws/Fs,'high'); [z,p,k]=butter(N,Ws/Fs,'high'); [H,W]=freqz(b,a);

plot(W*Fs/(2*pi),abs(H)); xlabel('Hz'); ylabel('Db');

图3.4.2 高通滤波器

Wp1=500; Wp2=800; Ws1=600; Ws2=700;

Wp=[Wp1 Wp2]; Ws=[Ws1 Ws2]; Rp=0.1; Rs=50; Fs=2000;

[N,Wn]=ellipord(2*Wp/Fs,2*Ws/Fs,Rp,Rs); [num,den]=ellip(N,Rp,Rs,Wn,'stop'); [H,W]=freqz(num,den); plot(W*Fs/(2*pi),abs(H)); xlabel('HZ'); ylabel('DB');

15

图3.4.3 阻带滤波器

Wp=[60 200]/500;

Ws=[50 250]/500; Rp=3;Rs=40;

[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs) [b,a]=cheby1(n,Rp,Wn); [H,W]=freqz(b,a); plot(W/(2*pi),abs(H)); grid;

xlabel('Hz'); ylabel('DB');

title('Chebysheve I Bandpass Filter')

图3.4.4 带通滤波器

16

（5）FIR滤波器

function[hd]=ideal(wc,N) %求理想低通滤波器单位抽样响应 %wc理想低通滤波器的截止频率 %N滤波器长度 q=(N-1)/2; n=0:N-1;

m=n-q+eps; %eps是浮点相对经度＝2^-52 hd=sin(wc*m)./(pi*m);

构造一个N阶差分方程，使该差分方程为数字低通滤波器。利用MATLAB程序画出相应的幅频图H(ej)~。 p=0.2*pi; ws=0.3*pi;

width=ws-wp %确定过度带宽

N=ceil(6.6*pi/width)+1 % 确定滤波器阶数 ，ceil是向上取整函数 n=0:N-1;

wc=(ws+wp)/2 %理想低通的截止频率 hd=ideal(wc,N);

wn=(hamming(N))'; %海明窗

h=hd.*wn; %截取得到实际的单位抽样响应

k=0:500;w=(pi/500)*k;

[H]=fr(h,1,w); %计算实际滤波器的幅度响应 mag=abs(H);

db=20*log10((mag+eps)/max(mag));%eps是浮点相对经度＝2^-52 wth=pi/500;

Ap=-(min(db(1:1:wp/wth+1))) %实际通带纹波 As=-round(max(db(ws/wth+1:1:500))) %实际阻带纹波 %画图

subplot(2,2,1); stem(n,hd,'fill'); title('理想抽样响应'); axis([0 N-1 -0.1 0.3]); ylabel('hd(n)'); subplot(2,2,2); stem(n,wn); title('海明窗'); axis([0 N-1 0 1.1]); ylabel('wn'); subplot(2,2,3); stem(n,h,'fill'); title('实际抽样响应'); axis([0 N-1 -0.2 0.3]); xlabel('n');

17

ylabel('h(n)'); subplot(2,2,4); plot(w/pi,db);

title('幅度响应(dB)'); axis([0 1 -100 10]); grid;

xlabel('以pi为单位的频率') ylabel('分贝数');

图3.4.5

p=0.2*pi; ws=0.3*pi;

width=ws-wp %确定过度带宽

N0=ceil(6.6*pi/width) ; % 确定滤波器阶数 ，ceil是向上取整函数N=N0+mod(N+1,2); n=0:N-1;

wc=(ws+wp)/2 %理想低通的截止频率 hd=ideal(wc,N);

wn=fir1(N-1,wc,'high',hanning(N)); %海明窗 h=hd.*wn; %截取得到实际的单位抽样响应 k=0:500;w=(pi/500)*k;

[H]=fr(h,1,w); %计算实际滤波器的幅度响应 mag=abs(H);

db=20*log10((mag+eps)/max(mag));%eps是浮点相对经度＝2^-52 wth=pi/500;

Ap=-(min(db(1:1:wp/wth+1))) %实际通带纹波 As=-round(max(db(ws/wth+1:1:500))) %实际阻带纹波 %画图

subplot(2,2,1);

stem(n,hd,'fill');grid title('理想抽样响应'); axis([0 N-1 -0.1 0.3]); ylabel('hd(n)');

18

subplot(2,2,2); stem(n,wn);

title('海明窗');grid axis([0 N-1 0 1.1]); ylabel('wn'); subplot(2,2,3); stem(n,h,'fill');

title('实际抽样响应');grid axis([0 N-1 -0.2 0.3]); xlabel('n'); ylabel('h(n)'); subplot(2,2,4); plot(w/pi,db);

title('幅度响应(dB)'); axis([0 1 -100 10]); grid;

xlabel('以pi为单位的频率') ylabel('分贝数');

图3.4.6

3.5 简要回答如下思考题

（1）你所构造的数字滤波器是IIR还是FIR?试画出该滤波器的运算结构图。

IIR数字滤波器设计。IIR网络的特点是信号流图中含有反馈支路，即含有环路，其单位脉冲响应是无限长的。基本网络结构有三种，即直接型、级联型和并联型。 a. 直接型

将N阶差分方程重写如下：

MNiiy(n)bx(ni)ai0i1y(ni)

设M=N=2,其系统函数如下：

19

2H(z)i0bizi112H1(z)H2(z)iai1zi

图3.5.1

按照差分方程可以直接画出网络结构如图3.5.1(a)所示。图中第一部分系统函数用H1(z)表示，第二部分用H2(z)表示，那么b.级联型

系统函数H(z)中，分子、分母均为多项式，且多项式的系数一般为实数。现将分子、

分母多项式分别进行因式分解，得到：

MH(z)H1(z)H2(z)。

H(z)A(1Cr1Nrz)1

r(1dr1z)1式中A是常数，Cr和dr分别表示零点和极点。由于多项式的系数是实数，Cr和dr是实数或者是共轭成对的复数，将共轭成对的零点（极点）放在一起，形成一个二阶多项式，其系数仍为实数；再将分子、分母均为实系数的二阶多项式放在一起，形成一个二阶网络Hj(z)。Hj(z)如下式：

20

Hj(z)0j1jz11jz112jz222jz

式中，0j、1j、2j、1j和2j均为实数。这样H(z)就分解成一些一阶或二阶数字网络的级联形式，如下式：

H(z)H1(z)H2(z)...Hk(z)

式中Hi(z)表示一个一阶或二阶的数字网络的系统函数，每个Hi(z)的网络结构均采用前面介绍的直接型网络结构，如图3.5.2所示。

图3.5.2

利用双线性变换设计IIR滤波器（上述程序利用巴特沃斯数字低通滤波器的设计），首先要设计出满足指标要求的模拟滤波器的传递函数所要设计的IIR滤波器的系统函数H(z)。 c.并联型

如果将级联形式的H(z)展成部分分式形式，则得到IIR并联型结构。

H(z)H1(z)H2(z)...Hk(z)

Ha(s)，然后由

Ha(s)通过双线性变换可得

式中，Hi(z)通常为一阶网络或二阶网络，网络系统均为实数。二阶网络的系统函数一般为 Hi(z)0i1iz11iz1122iz

式中，0i、1i、1i和2i都是实数。如果2i0，则构成一阶网络。由(6.2.4)式，其输出Y(z)表示为:Y(z)H1(z)X(z)H2(z)X(z)...Hk(z)X(z)

上式表明将x(n)送入每个二阶（包括一阶）网络后，将所有输出加起来得到输出y(n)。 如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率

wp和ws的转换，对

p和s指标不作变化。边界频率的转换关系为

21

2Ttan(12w)。接着，按照模拟低通滤波器的技术指标根据相应设计公式求出滤波器

c的阶数N和3dB截止频率；根据阶数N查巴特沃斯归一化低通滤波器参数表，得到归

psc一化传输函数

Ha(p)；最后，将代入

Ha(p)去归一，得到实际的模拟滤波器传输

11函数

Ha(s)s21zT1z。之后，通过双线性变换法转换公式，得到所要设计的IIR滤波器

的系统函数H(z)。 类型 低通 映射 1设计参数 z1za1asin[(cc)/2]sin[(cc)/2]''1az c' ：目的滤波器截止频率 cos[(cc)/2]cos[(cc)/2]'' 高通 z1z1a1a1az c' ：目的滤波器截止频率 cos[(cc)/2]cos[(cc)/2]'' 带通 z1z2[2ab/(b1)]z21a[(b1)/(b1)]1 [(b1)/(b1)]z[2ab/(b1)]z1 bcot[(c1c2)/2]tan(c/2)c1c2：带通DF通带下限截止频率 ：带通DF通带上限截止频率 cos[(cc)/2]cos[(cc)/2]'' 带阻 z1[2a/(b1)]z21z2a[(1b)/(1b)]1 [(1b)/(1b)]z[2a/(b1)]z1 btan[(c1c2)/2]tan(c/2)c1c2：带通DF通带下限截止频率 ：带通DF通带上限截止频率 （2）解释系统参数、系统零极点与系统频响特性的关系。

22

在z平面上直接设计IIR数字滤波器，以滤波器响应作为依据，直接在z平面上，通过多次选定极点和零点位置逼近该响应。 即在单位园内设置一对共轭极点 ，频响在w0处就有一峰值。r越近于1，极点位置越接近单位园，则峰值就越尖锐。同理，若在单位园内设置一对零点 ,频响就会在w1处出现零值，即可实现陷波。如特性还达不到要求，可再移动零、极点，这样作二、三次调整后，就可以获得一些简单的DF。这种方法，可以设计一些简单阶数很低（1~2阶)的DF。

（3）说明在数字域上如何定义高通、低通、带通、全通滤波器。

低通滤波器容许低频信号通过, 但减弱(或减少)频率高於截止频率的信号的通过。高通滤波器容许高频信号通过, 但减弱(或减少)频率低于於截止频率的信号的通过。带通滤波器容许一定频率范围信号通过, 但减弱(或减少)频率低于於下限截止频率和高于上限截止频率的信号的通过。带阻滤波器减弱(或减少)一定频率范围信号, 但容许频率低于於下限截止频率和高于上限截止频率的信号的通过。

4.设计体会

通过该课程设计，我们受益匪浅。Matlab是一个基于矩阵运算的软件，它的运算功能非常强大，编程效率高，强大而智能化的作业图功能，可扩展性强。在这段时间里，我们主要学习MATLAB的工具的使用，熟悉其最基础的功能，锻炼了我的实际动手能力。通过用MATLAB实现DFT与FFT，我熟练了MATLAB的基本函数与用法，会通过矩阵、数组进行的复杂运算，并且更加熟练了用图形展示实验结果。Help是MATLAB中最有效的命令。遇到问题，通常都可以借助help解决问题。学好MATLAB是不容易的，这是一件需要持之以恒的事，必须要坚持不懈的学习，需要我们勤于思考，勤于记忆，勤于动手。程序设计是实践性和操作性很强的事情，需要我们亲自动手。因此，我们应该经常自己动手实际操作设计程序，熟悉MATLAB的操作，这对提高我们的操作能力非常有效。

离散傅里叶变换建立了有限长序列与其近似频谱之间的联系，在理论上具有重要意义。离散傅里叶变换DFT在数字通信、语音处理、图像处理、谱估计、仿真、系统分析、雷达、光学影像、地震等各个领域得到广泛应用，但是这都是以卷积和相关运算，对连续信号和序列进行谱分析为基础的。通过本次课程设计，我们对DFT在进行频谱的分析上有了根深刻的理解和掌握。DFT实现了频域采样，同时DFT存在快速算法FFT，所以在实际应用中，可以利用计算机，用DFT来逼近连续时间信号的傅里叶变换，进而分析连续时间信号频谱。同时知道了补零点的作用，其仅仅是提高了计算分辨率，得到的是高密度频谱，并不能得到高分辨率谱，要提高频率分辨率，则要通过增加数据记录长度来提高物理分辨率。

在编程实现中，遇到了一些问题，为此我翻阅一些了参考书。期间我们不仅学到了许多课本上的知识，还有课本以外的内容，学到了许多课本上所没提到的东西，这些东西都让我们耳目一新，开阔了视野，拓宽了知识面。从以前仅仅掌握离散傅里叶变换的概念，到现在渐渐领悟到离散傅里叶变换的一些实际应用，更明白它在实际设计中的作用，从理论到实践的逐步过渡，增了动手能力。 通过课程设计，我更了解了自己学习过程中的不足，这主要体现在做课程设计的过程中，我们深深感觉到自身所学知识的有限，书本上没有提及的环节，我们基本都没有去研究过，做的过程有时突然间觉得有点茫然，虽然通过查阅可以解决问题，但还是浪费了许多时间，这一点是我们在以后的学习中必须加以改进的地方，同时在以后的学习过程中也要督促自己不断地完善自我，超越自我。

23