数字信号处理实验报告

实验报告

课程名称： 数字信号处理 实验题目： 1、用FFT作谱分析

2、用窗函数法设计FIR数字滤波器

院 系： 电子与信息工程学院 班 级： 姓 名： 学 号： 指导教师： 实验时间： 2013年11月

哈尔滨工业大学

实验一、用FFT作谱分析

1、实验目的

(1)、进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。 (2)、熟悉FFT算法原理和FFT子程序的应用。 (3)、学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。

2、实验原理

1、离散傅立叶变换(DFT)

离散傅立叶变换是周期序列，有N个的数值，所以它的许多特性可以通过有限长序列延拓来得到。对于一个长度为N的有限长序列x(n)，n只在0~(N-1)各点上有非零值，即

x(n),0nN1 x(n)0,其他x(n)，则有 把序列x(n)以N为周期进行周期延拓得到周期序列~x(n),0nN1 ~x(n)0,其他所以，有限长序列x(n)的离散傅立叶变换(DFT)为

knX(k)DFT[x(n)]x(n)WN,0kN1

n0N12、用FFT对信号作频谱分析

对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是N/2，因此要求2Π/N<=D。可以根据此式选择FFT的变换区间N。

误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。 周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

如果给出的是模拟信号xa(t)，则首先要根据其最高频率确定抽样频率fs以及由频率分辨率选择抽样点数N，然后对其进行软件抽样(即计算 x(n)=xa(nT)，0≤n≤N-1)，产生对应序列 x(n)。再利用MATLAB所提供的库函数fft(n,x)进行FFT计算。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

DFT的运算量是与N2成正比的，所以N越小对计算越有利，因而小点数序列的DFT比大点数序列的DFT运算量要小。FFT正是基于这样的基本思路而发展起来的，她的算法基本上可分成两大类，即按时间抽取法和按频率抽取法。 我们最常用的是N=2M的情况，该情况下的变换成为基2快速傅里叶变换。 完成一次完整的FFT计算总共需要(N/2)log2N次复数乘法运算和Nlog2N次复数加法运算。很明显，N越大，FFT的优点就越突出。

3、实验步骤与实验结果

3.1、实验步骤：

(1)、复习DFT的定义、 性质和用DFT作谱分析的有关内容。

(2)、复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图。 (3)、编写信号产生程序， 产生以下典型信号供谱分析用：

x1(n)R4(n)

n1,0n3

x(n) 2  8  n 4  n7  其它n

0

4n0n3

x3(n)n34n7

 0 其它n

 x4(n)cosn，0n19 4 x5(n)sinn，0n19 8x6(t)cos8tcos16tcos20t

(4)、按实验内容要求，上机实验，并写出实验报告。

(5)、上机实验内容：

①、画出(3)中所给出的信号，并逐个进行谱分析（即画出幅频特性）。下面给出针对各信号的FFT变换区间N以及对连续信号x6(t)的采样频率fs，供实验时参考。 x1(n), x2(n), x3(n), x4(n), x5(n): N=8, 16 x6(t): fs=(Hz), N=16, 32, (n=0:1:69)

②、令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换，X(k)=DFT[x(n)]。 ③、令x(n)=x4(n)+jx5(n)，重复②。

3.2、实验结果

(1)、对x1(n)分别做8点和16点的DFT（使用FFT实现），x1(n)的波形图和所得频谱图如下图所示：

图1、x1(n)的8点和16点的频谱图

(2)、对x2(n)分别做8点和16点的DFT（使用FFT实现），x2(n)的波形图和所得频谱图如下图所示：

图2、x2(n)的8点和16点的频谱图

(3)、对x3(n)分别做8点和16点的DFT（使用FFT实现），x3(n)的波形图和所得频谱图如下图所示：

图3、x3(n)的8点和16点的频谱图

(4)、对x4(n)分别做8点和16点的DFT（使用FFT实现），x4(n)的波形图和所得频谱图如下图所示：

图4、x4(n)的8点和16点的频谱图

(5)、对x5(n)分别做8点和16点的DFT（使用FFT实现），x5(n)的波形图和所得频谱图如下图

所示：

图5、x5(n)的8点和16点的频谱图

(6)、对x6(t)抽样后的序列x6(n)分别做16点、32点和点的DFT（使用FFT实现），x6(t)和抽样后的x6(n)的波形图和所得频谱图如下图所示：

图6、x6(t)抽样后的16点、32点和点的频谱图

(7)、对x(n)=x4(n)+x5(n)分别做8点和16点的DFT（使用FFT实现），x(n)的波形图和所得频谱图如下图所示：

图7、x(n)=x4(n)+x5(n)的8点和16点的频谱分析图

(8)、对xj(n)=x4(n)+j*x5(n)分别做8点和16点的DFT（使用FFT实现），x(n)的波形图和所得频谱图如下图所示：

图7、xj(n)=x4(n)+j*x5(n)的8点和16点的频谱分析图

4、实验主要结论、误差分析和参数选取方法

1、本实验主要是求x1(n)，x2(n)，x3(n)，x4(n)，x5(n)，x6(t)，x(n)=x4(n)+x5(n)，xj(n)=x4(n)+j*x5(n)的DFT。实验中的到图像和理论值基本相符。N取值的不同对应的频谱图也不同，N较小时，有些点的值比较小，区间选取小不能完全的反应原函数的频谱特征；增大N值可以得到更能完整反映原信号的频谱。

离散傅立叶变换可以看作是X(z)在z=e(j2Π)/N时的Z变换，即表明x(n)的N点DFT是x(n)的Z变换在单位圆上的N点等间隔采样。DFT也可以看作X(k)在ω=(2Π)/N时的傅立叶变换，即表明X(k)可以看作x(n)的傅立叶变换X(ejω)在区间[0,2Π]上的N点等间隔采样。

2、 ①、序列x1(n)理论FFT频谱为一个抽样信号，由图1看出，随着采样率的提高，得到的FFT频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。

②、序列x2(n)理论FFT频谱为一个抽样函数平方的信号，由图2看出，随着采样率的提高，得到的FFT频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。

③、序列x3(n)由图3看出，随着采样率的提高，得到的FFT频谱分辨率就越高，但对于这个函数来说，N仍然过小，因此幅频特性曲线与理论结果有较大差距。

④、序列x4(n)的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。由图4可以看出，当N=8,16时，满足抽样定理。

⑤、序列x5(n)的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。由图5可以看出，当N=8时，不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。

⑥、x6(t)抽样后的序列x6(n)的理论FFT抽样频谱为有三个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。由图6可以看出，当N=16时，不满足抽样定理，因此产生了混叠失真，当N=32,时，满足抽样定理，因此可得到与理论曲线一致的图形。

⑦、序列x(n)=x4(n)+x5(n)的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。由图7可以看出，当N=8时，不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。

⑧、序列xj(n)=x4(n)+j*x5(n)的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。由图8可以看出，当N=8时，

不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。结合x4(n)和x5(n)的频谱分析，从以上的图可以看出，将原信号的DFT变换分为共轭对称部分合共个反对称部分，则可以得出原信号的实部对应离散傅里叶变换的共轭对称部分，原信号的虚部对应离散信号的共轭反对称部分。

误差原因：因为实际应用与软件相比，软件抽样点与实际相比，抽样点是有限的。 3、N的选取：如果信号长度为M，N要比信号的长度M大，这样才能由频率抽样序列X(k)不失真的回复X(jω)，增大N可以更有效地反映原信号的信息。

5、思考题

(1) 、在N=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?

答：相同。在N=8的时候，x3(n)只是x2(n)以它的长度8为周期，将其延拓成周期 序列然后加以移位，最后截取主值区间(n=0到8)的序列值，因此x3(n)是x2(n) 进行循环位移后的结果。由于序列的循环位移性质

km DFT[x(nm)]WNDFT[x(n)]

因此，循环位移只影响到DFT后的相频特性，而不影响幅频特性，因此x2(n) 和x3(n)的幅频

特性会相同。

当N=16时，对x2(n)和x3(n)做DFT就是为它们分别在末尾补零之后再进行的DFT，则此时两个原序列经过补零以后的新序列就不满足循环位移的性质，因此x2(n)和x3(n)的幅频特性就会发生变化。

(2)、如果周期信号的周期预先不知道，如何用FFT进行谱分析?

答：先对信号进行N1点的FFT，观察并记录其包络；取N2>N1，再观察包络。比较二者结果，若二者的差别满足分析误差的要求，就可以近似得到该信号的频谱；假若不满足误差要求则继续加长截取的长度进行FFT，直到满足条件，得到较完整幅频特性。

6、程序代码：

%序列x1(n)的频谱分析 n=0:1:19;

y=[ones(1,4),zeros(1,16)]; subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x1(n)');

title('序列x1(n)的波形'); grid on; N=8;

y1=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,2);

stem(n,abs(y1),'.'); xlabel('n'); ylabel('幅值');

title('序列x1(n)的N=8的FFT频谱');

grid on;

N=16;

y2=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,3); stem(n,abs(y2),'.'); xlabel('n'); ylabel('幅值');

title('序列x1(n)的N=16的FFT频谱'); grid on;

%序列x2(n)的频谱分析 n=0:1:19;

y=(n+1).*(n>=0 & n<=3)+(8-n).*(n>=4 & n<=7); subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x2(n)');

title('序列x2(n)的波形'); grid on; N=8;

y1=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,2);

stem(n,abs(y1),'.'); xlabel('n'); ylabel('幅值');

title('序列x2(n)的N=8的FFT频谱'); grid on;

N=16;

y2=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,3); stem(n,abs(y2),'.'); xlabel('n'); ylabel('幅值');

title('序列x2(n)的N=16的FFT频谱'); grid on;

%序列x3(n)的频谱分析 n=0:1:19;

y=(4-n).*(n>=0 & n<=3)+(n-3).*(n>=4 & n<=7); subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x3(n)');

title('序列x3(n)的波形'); grid on; N=8;

y1=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,2); stem(n,abs(y1),'.'); xlabel('频率(Hz)'); ylabel('幅值');

title('序列x3(n)的N=8的FFT频谱'); grid on;

N=16;

y2=fftshift(fft(y,N)); n=0:N-1;

subplot(3,1,3); stem(n,abs(y2),'.'); xlabel('频率(Hz)'); ylabel('幅值');

title('序列x3(n)的N=16的FFT频谱'); grid on;

%序列x4(n)的频谱分析 n=0:1:19;

y=cos(pi*0.25*n); subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x4(n)');

title('序列x4(n)的波形图'); grid on;

N=8; n=0:N-1; n1=0:7;

x1=cos(pi*(n1)/4); n2=8:15;

x2=cos(pi*(n2)/4);

n3=16:19;

x3=cos(pi*(n3)/4); y11=fft(x1,N); y12=fft(x2,N); y13=fft(x3,N);

y1=fftshift(abs(y11+y12+y13)); subplot(3,1,2); stem(n,y1,'.'); xlabel('n'); ylabel('幅值');

title('序列x4(n)的N=8的FFT频谱'); grid on;

N=16; n=0:N-1; n1=0:15;

x1=cos(pi*n1/4); n2=16:19;

x2=cos(pi*n2/4); y21=fft(x1,N); y22=fft(x2,N);

y2=fftshift(abs(y21+y22)); subplot(3,1,3); stem(n,y2,'.'); xlabel('n'); ylabel('幅值');

title('序列x4(n)的N=16的FFT频谱'); grid on;

%序列x5(n)的频谱分析 n=0:1:19;

y1=sin(pi*0.125*n); subplot(3,1,1); stem(n,y1,'.'); xlabel('n'); ylabel('x5(n)');

title('序列x5(n)的波形'); grid on;

N=8; n=0:N-1; n1=0:7;

x1=sin(pi*n1/8); n2=8:15;

x2=sin(pi*n2/8); n3=16:19;

x3=sin(pi*n3/8); y11=fft(x1,N); y12=fft(x2,N); y13=fft(x3,N);

y1=fftshift(abs(y11+y12+y13)); subplot(3,1,2); stem(n,y1,'.'); xlabel('n'); ylabel('幅值');

title('序列x5(n)的N=8的FFT频谱'); grid on;

N=16; n=0:N-1; n1=0:15;

x1=sin(pi*n1/8); n2=16:19;

x2=sin(pi*n2/8); y21=fft(x1,N); y22=fft(x2,N);

y2=fftshift(abs(y21+y22)); subplot(3,1,3); stem(n,y2,'.'); xlabel('n'); ylabel('幅值');

title('序列x5(n)的N=16的FFT频谱'); grid on;

%连续函数x6(t)的抽样频率为Hz的抽样信号x6(n)的频谱分析 x=0:1/:69/;

y=cos(pi*8*x)+cos(pi*16*x)+cos(pi*20*x); subplot(2,3,1); plot(x,y); xlabel('x'); ylabel('y');

title('x6(t)的波形');

m=0:1:69; n=m/;

yn=cos(pi*8*n)+cos(pi*16*n)+cos(pi*20*n); subplot(2,3,2); stem(x,yn,'.');

xlabel('n'); ylabel('yn');

title('x6(t)的抽样频率为Hz的抽样后的波形'); grid on;

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

y1=fftshift(abs(fft(yn,16)+fft(yn(17:32),16)+fft(yn(33:48),16)+... fft(yn(49:),16)+fft(yn(65:69),16))); subplot(2,3,4); stem(n,y1,'.'); xlabel('n'); ylabel('幅值');

title('序列x6(n)的N=16的FFT频谱'); grid on;

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

y2=fftshift(abs(fft(yn,32)+fft(yn(33:),32)+fft(yn(65:69),32))); subplot(2,3,5); stem(n,y2,'.'); xlabel('n'); ylabel('幅值');

title('序列x6(n)的N=32的FFT频谱'); grid on;

N=; n=0:N-1;

y3=fftshift(abs(fft(yn,)+fft(yn(:69),))); subplot(2,3,6); stem(n,y3,'.'); xlabel('n'); ylabel('幅值');

title('序列x6(n)的N=的FFT频谱'); grid on;

%x(n)的频谱分析(x(n)=x4(n)+x5(n)) n=0:1:19;

y=cos(pi*0.25*n)+sin(pi*0.125*n); subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x(n)');

title('序列x(n)的波形图');

N=8; n=0:N-1; n1=0:7;

x1=cos(pi*n1/4)+sin(pi*n1/8); y1=fft(x1,N); n2=8:15;

x2=cos(pi*n2/4)+sin(pi*n2/8); y2=fft(x2,N); n3=16:19;

x3=cos(pi*n3/4)+sin(pi*n3/8); y3=fft(x3,N);

y1=fftshift(abs(y1+y2+y3)); subplot(3,1,2); stem(n,y1,'.'); xlabel('n'); ylabel('幅值');

title('序列x(n)的N=8的FFT频谱'); grid on;

N=16; n=0:N-1; n1=0:15;

x1=cos(pi*n1/4)+sin(pi*n1/8); y1=fft(x1,N); n2=16:19;

x2=cos(pi*n2/4)+sin(pi*n2/8); y2=fft(x2,N);

y2=fftshift(abs(y1+y2)); subplot(3,1,3); stem(n,y2,'.'); xlabel('n'); ylabel('幅值');

title('序列x(n)的N=16的FFT频谱'); grid on;

%xj(n)的频谱分析(xj(n)=x4(n)+j*x5(n)) n=0:1:19;

y=cos(pi*0.25*n)+j*sin(pi*0.125*n); subplot(3,1,1); stem(n,y,'.'); xlabel('n'); ylabel('x(n)');

title('序列xj(n)的波形图');

N=8; n=0:N-1; n1=0:7;

x1=cos(pi*n1/4)+j*sin(pi*n1/8); y1=fft(x1,N); n2=8:15;

x2=cos(pi*n2/4)+j*sin(pi*n2/8); y2=fft(x2,N); n3=16:19;

x3=cos(pi*n3/4)+j*sin(pi*n3/8); y3=fft(x3,N);

y1=fftshift(abs(y1+y2+y3)); subplot(3,1,2); stem(n,y1,'.'); xlabel('n'); ylabel('幅值');

title('序列xj(n)的N=8的FFT频谱'); grid on;

N=16; n=0:N-1; n1=0:15;

x1=cos(pi*n1/4)+j*sin(pi*n1/8); y1=fft(x1,N); n2=16:19;

x2=cos(pi*n2/4)+j*sin(pi*n2/8); y2=fft(x2,N);

y2=fftshift(abs(y1+y2)); subplot(3,1,3); stem(n,y2,'.'); xlabel('n'); ylabel('幅值');

title('序列xj(n)的N=16的FFT频谱'); grid on;

实验二： 用窗函数法设计FIR数字滤波器

1、实验目的

(1)、熟悉矩形窗、汉宁窗、海明窗和布莱克曼窗。

(2)、掌握用上述窗函数法设计FIR数字滤波器的原理和方法。 (3)、熟悉线性相位FIR数字滤波器特性。 (4)、了解各种窗函数对滤波特性的影响。

2、实验原理

滤波器的理想频率响应函数为Hd(ejω)，则其对应的单位脉冲响应为

1hd(n) =

2Hd(ej)ejnd

窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数。w(n)将hd(n)截断，并进行加权处理：h(n)=hd(n)*ω(n) h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数H(ejω)为

H(ejω) =

h(n)en0N1jn

如果要求线性相位特性， 则h(n)还必须满足：h(n)=h(N-1-n)或h(n)=-h(N-1-n)

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。

3、实验步骤与实验结果

3.1、实验步骤

(1)、复习用窗函数法设计FIR数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。 (2)、编写程序，其中幅度特性要求用dB表示。备注：不许用 freqz() 这个函数 设

H(k)DFT[h(n)]

H(k)Hk(k)jHI(k)

2H(k)HR(k)HI2(k)

2画图时，用 H ( k ) 打印幅度特性。第k点对应的频率为 w k  k 为使曲线包络更接20logNjw近 H (e ) 的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h(n)末尾补零的方法，使长度变为，再进行点DFT，则可得到更精确的幅度衰减特性曲线。 (3)、上机实验内容：

设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即

jwaN1e,wwc jwHde其中a 20,wcw 11wcjwajwnjwjwnhnHeedweedwdd 22wc sinwcna nan=15，n=33，ωc=/Π4，用四种窗函数设计线形相位低通滤波器。要求在两种窗口长度下，绘制相应的幅频和相频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。

3.2、实验结果

(1)、加Boxcar窗后的滤波器的幅频特性和相频特性(N=15和33)

图9

(2)、加Hanning窗后的滤波器的幅频特性和相频特性(N=15和33)

图10

(3)、加Hamming窗后的滤波器的幅频特性和相频特性(N=15和33)

图11

(4)、加Blackman窗后的滤波器的幅频特性和相频特性(N=15和33)

图12

4、实验结果分析

1、由图像观察得过渡带宽大小依次为：布莱克曼窗>海明窗>汉宁窗>矩形窗。

比较不同窗函数的幅频特性图像发现滤波特性有如下关系：布拉克曼窗>汉宁窗>三角窗>矩形窗。

实验结果和理论基本符合。其中，改变窗口长度N，可以改变窗口频谱的主瓣宽度，从而改变系统函数的过渡带的宽度，增大N，系统函数的过渡带变窄。但是，N值得改变不能改变其主瓣与旁瓣电平的比值。要想获得比较低的旁瓣和主瓣的比值，从而获得较高的阻带衰减，就应该选择旁瓣衰减大的窗函数。比如，Hanning窗的阻带衰减为-44dB，而Hamming窗的阻带衰减为-53dB，比Hanning窗就要好。 2、窗口长度N和窗函数类型对滤波特性的影响：

增加窗的长度，可以减少窗的主瓣宽度，从而减少H(ω)过渡带的带宽，但增加N并不能减少带内波动以及加大阻带，这两个指标只能从窗函数的形状上找解决方法。通过对窗函数选择，可以获得符合要求的阻带衰减指标。 3、窗函数法特点：

窗口法设计的主要优点是简单，使用方便，且相位为线性变化。窗口函数大多有封闭的公式可循，性能、参数都已有表格、资料可供参考，计算程序简便，所以很实用。 缺点是通带和阻带的截止频率不易控制

5、思考题

(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器?写出设计步骤。

答：①、根据给定的阻带衰减选择合适的窗口ω(n)，在保证阻带衰减的前提下，尽量选择过渡带较窄的窗函数。

②、根据所选的窗函数的D值和给定的过度带宽求出窗口的长度N;

③、由给定的指标给定希望逼近的频率响应函数Hd(jω)，包括幅频特性与相频特性，同时相频特性应该满足线性相位的要求；

④、根据给出的Hd(jω)求出相应的Hd(n)。可以利用如下公式计算

1hd(n) =

2Hd(ej)ejnd

⑤、计算所设计的滤波器的单位冲激响应为h(n)=hd(n)*ω(n)；

⑥、验证计算结果，计算出H(jω)，核对其是否与要求的指标一致，如果不符合则重复上面步骤直至符合为止；

⑦、根据单位冲激响应计算其系统函数：

H(z)N10h(n)znn

⑧、根据系统函数实现FIR滤波器系统。

(2)、如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd(n) .

答：只需将一个高通滤波器器和一个低通滤波器相加，即可得到带通滤波器，具体设计过程与低通滤波器相同。

理想带通滤波器的频率响应是

ejHd(e)0j

012其他n

其中求此滤波器的单位冲激响应hd(n)，即

hd(n)12N12

Hd(ej)ejd

nn1sinn2sinn1(n)1(21)

根据对通带、阻带衰减的要求以及对过渡带宽的要求，可选择窗函数及窗的点数N。由此可求得所需的线性相位带通滤波器的单位冲激相应h(n)=hd(n)*ω(n)。

由h(n)求Hd(ejω)，检验各项指标是否满足要求。如不满足要求，则要改变N，或改变窗形状，然后重新计算。

6、程序代码

%boxcar n=0:1:14;

wr1=ones(1,15);

hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps)); h1=hd.*wr1; N=;

H1=fft(h1,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,1);

plot(w,(20*log((abs(H1))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('boxcar窗，幅频特性（n=15）'); subplot(2,2,3);

plot(w,unwrap(angle(H1))); xlabel('w/rad');

ylabel('相位');

title('boxcar窗，相频特性（n=15）'); grid on; n=0:1:32;

wr2=ones(1,33);

hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps)); h2=hd.*wr2; N=;

H2=fft(h2,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,2);

plot(w,(20*log((abs(H2))))); grid on;

xlabel('w/rad')

ylabel('20lg|H(jw)|/dB')

title('boxcar窗，幅频特性（n=33）'); subplot(2,2,4);

plot(w,unwrap(angle(H2))); xlabel('w/rad');

ylabel('相位');

title('boxcar窗，相频特性（n=33）'); grid on;

%hanning窗 n=0:1:14;

wh1=0.5*(1-cos(2*pi*n/14));

hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps)); h1=hd.*wh1; N=;

H1=fft(h1,N); n=0:N-1;

w=2*pi*n/; subplot(2,2,1);

plot(w,(20*log((abs(H1))))); grid on;

xlabel('w/rad')

ylabel('20lg|H(jw)|/dB');

title('hanning窗，幅频特性（n=15）'); subplot(2,2,3);

plot(w,unwrap(phase(H1))); xlabel('w/rad');

ylabel('相位');

title('hanning窗，相频特性（n=15）'); grid on;

n=0:1:32;

wh2=0.5*(1-cos(2*pi*n/32));

hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps)); h2=hd.*wh2; N=;

H2=fft(h2,N); n=0:N-1;

w=2*pi*n/; subplot(2,2,2);

plot(w,(20*log((abs(H2))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('hanning窗，幅频特性（n=33）'); subplot(2,2,4);

plot(w,unwrap(angle(H2))); xlabel('w/rad'); ylabel('相位');

title('hanning窗，相频特性（n=33）'); grid on;

%hamming窗 n=0:1:14;

wH=0.54-0.46*cos(2*pi*n/(14+eps));

hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps)); h1=hd.*wH; N=;

H1=fft(h1,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,1);

plot(w,(20*log((abs(H1))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('hamming窗，幅频特性（n=15）'); subplot(2,2,3);

plot(w,unwrap(angle(H1))); xlabel('w/rad');

ylabel('相位');

title('hamming窗，相频特性（n=15）'); grid on;

n=0:1:32;

wH=0.54-0.46*cos(2*pi*n/(32+eps));

hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps)); h1=hd.*wH; N=;

H1=fft(h1,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,2);

plot(w,(20*log((abs(H1))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('hamming窗，幅频特性（n=33）'); subplot(2,2,4);

plot(w,unwrap(angle(H1))); xlabel('w/rad');

ylabel('相位');

title('hamming窗，相频特性（n=33）');

grid on;

%blackman窗 n=0:1:14;

wb1=0.42-0.5*cos(2*pi/(14+eps)*n)+0.08*cos(4*pi/(14+eps)*n); hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps)); h1=hd.*wb1; N=;

H1=fft(h1,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,1);

plot(w,(20*log((abs(H1))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('blackman窗，幅频特性（n=15）'); subplot(2,2,3);

plot(w,unwrap(angle(H1))); xlabel('w/rad');

ylabel('相位');

title('blackman窗，相频特性（n=15）'); grid on;

n=0:1:32;

wb2=0.42-0.5*cos(2*pi/(32+eps)*n)+0.08*cos(4*pi/(32+eps)*n); hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps)); h2=hd.*wb2; N=;

H2=fft(h2,N); n=0:N-1;

w=2*pi/*n; subplot(2,2,2);

plot(w,(20*log((abs(H2))))); grid on;

xlabel('w/rad');

ylabel('20lg|H(jw)|/dB');

title('blackman窗，幅频特性（n=33）'); subplot(2,2,4);

plot(w,unwrap(angle(H2))); xlabel('w/rad');

ylabel('相位');

title('blackman窗，相频特性（n=33）'); grid on;