实验五 FIR 滤波器的设计

HUNAN  UNIVERSITY

数字信号

课程实验报告

专业班级：           通信工程一班              

完 成 日  期 ：   　       2014\\06\\01                

实验五 FIR 滤波器的设计

1、实验目的

    认真复习 FIR 数字滤波器的基本概念，线性相位FIR 滤波器的条件和特点、幅度函数

特点、零点位置的基本特点与性质；窗函数设计法的基本概念与方法，各种窗函数的性能和

设计步骤，线性相位FIR 低通、高通、带通和带阻滤波器的设计方法，频率采样设计法的基本概念和线性相位的实现方法。

    掌握几种线性相位的特点，熟悉和掌握矩形窗、三角形窗、汉宁窗、海明窗、布莱克曼

窗、凯塞窗设计IIR 数字滤波器的方法，熟悉和掌握频率抽样设计法的线性相位的设计方法，

并对各种线性相位的频率抽样法的设计给出调整和改进。

熟悉利用 MATLAB 进行各类FIR 数字滤波器的设计方法。

2、实验内容

a. 设线性相位FIR 滤波器单位抽样响应分别为

h(n) ={ -4,1, -1, -2,5,6,5, -2, -1,1, -4}

           ↑

h(n)={ -4,1,- 1, -2,5,6,6,5, -2,-1,1, -4}

           ↑

h(n) ={ -4,1, -1, -2,5,0, -5,2,1, -1,4}

           ↑

h(n) { -4,1, -1, -2,5,6, -6, -5,2,1, -1,4}

          ↑

分别求出滤波器的幅度频率响应 H(ω)，系统函数H(z)以及零极点分布，并绘制相应的

波形和分布图。

1、

代码：

function [H,w,a,L]=H_Type1(h);  %求幅度频率响应

M=length(h);

L=floor((M-1)/2);

a=[h(L+1) 2*h(L:-1:1)];

n=[0:1:L];

w=[0:1:500]'*pi/500;

H=cos(w*n)*a';

主函数

clc;

syms z;

%偶对称单位冲激响应h(n)，N为奇数

h1=[-4,1,-1,-2,5,6,5,-2,-1,1,-4];

N1=length(h1);

n1=[0:1:N1-1];

[H1,w1,a1,L1]=H_Type1(h1);

subplot(221);

stem(h1);  title('h1(n)');

subplot(222);

plot(H1);   title('幅度函数H1(w)');

subplot(223);

zplane(h1);  title('零极点');

HZ1=h1./(z.^n1)                  %系统函数H(z)

结果：

HZ1 =

[      -4,     1/z,  -1/z^2,  -2/z^3,   5/z^4,   6/z^5,   5/z^6,  -2/z^7,  -1/z^8,   1/z^9, -4/z^10]

2、

代码：

function [H,w,b,L]=H_Type2(h); %求幅度频率响应

M=length(h);

L=floor(M/2);

b=2*h(L:-1:1);

n=[1:1:L];

w=[0:1:500]'*pi/500;

H=cos(w*(n-1/2))*b';

主函数

%偶对称单位冲激响应h(n)，N为偶数

figure;

h2=[-4,1,-1,-2,5,6,6,5,-2,-1,1,-4];

N2=length(h2);

n2=[0:1:N2-1];

[H2,w2,a2,L2]=H_Type2(h2);

subplot(221);

stem(h2);   title('h2(n)');

subplot(222);

plot(H2);   title('幅度函数H2(w)');

subplot(223);

zplane(h2);   title('零极点');

HZ2=h2./(z.^n2)                 %系统函数H(z)

结果：

HZ2 =

[      -4,     1/z,  -1/z^2,  -2/z^3,   5/z^4,   6/z^5,   6/z^6,   5/z^7,  -2/z^8,  -1/z^9,  1/z^10, -4/z^11]

3、

代码：

function [H,w,c,L]=H_Type3(h); %求幅度频率响应

M=length(h);

L=floor((M-1)/2);

c=2*h(L:-1:1);

n=[1:1:L];

w=[0:1:500]'*pi/500;

H=sin(w*n)*c';

主函数：

%奇对称单位冲激响应h(n)，N为奇数

figure;

h3=[-4,1,-1,-2,5,0,-5,2,1,-1,4];

N3=length(h3);

n3=[0:1:N3-1];

[H3,w3,a3,L3]=H_Type3(h3);

subplot(221); 

stem(h3);    title('h3(n)');

subplot(222);

plot(H3);   title('幅度函数H3(w)');

subplot(223);

zplane(h3);   title('零极点');

HZ3=h3./(z.^n3)                 %系统函数H(z)

结果：

HZ3 =

[     -4,    1/z, -1/z^2, -2/z^3,  5/z^4,      0, -5/z^6,  2/z^7,  1/z^8, -1/z^9, 4/z^10]

4、

代码：

function [H,w,d,L]=H_Type4(h); %求幅度频率响应

M=length(h);

L=floor(M/2);

d=2*h(L:-1:1);

n=[1:1:L];

w=[0:1:500]'*pi/500;

H=sin(w*(n-1/2))*d';

主函数：

%奇对称单位冲激响应h(n)，N为偶数

figure;

h4=[-4,1,-1,-2,5,6,-6,-5,2,1,-1,4];

N4=length(h4);

n4=[0:1:N4-1];

[H4,w4,a4,L4]=H_Type4(h4);

subplot(221); 

stem(h4);    title('h4(n)');

subplot(222);

plot(H4);    title('幅度函数H4(w)');

subplot(223);

zplane(h4);  title('零极点'); 

HZ4=h4./(z.^n4)               %系统函数H(z)

结果：

HZ4 =

[      -4,     1/z,  -1/z^2,  -2/z^3,   5/z^4,   6/z^5,  -6/z^6,  -5/z^7,   2/z^8,   1/z^9, -1/z^10,  4/z^11]

b. 设计FIR 数字低通滤波器，技术指标为：ωp=0.2π，ωst=0.3π，δ1=0.25dB，δ2=50dB。

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

(3) 选择凯塞窗函数设计该滤波器，并绘制相应的波形图。

代码：

function [db,mag,pha,grd,w]=freqz_m(b,a)

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501));

w=(w(1:1:501));

mag=abs(H);

db=20*log10(mag+eps/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

function hd=ideal_lp(wc,N)

alpha=(N-1)/2;

n=[0:1:(N-1)];

n=n-alpha;

fc=wc/pi

hd=fc*sinc(fc*n)

主函数：

%低通滤波器

clc;

wp=0.2*pi;

ws=0.3*pi;

as=50;

tr_width=ws-wp;

N=ceil((as-7.95)/(2.286*tr_width)+1);

n=[0:1:N-1];

beta=0.102*(as-8.7) 

wc=(ws+wp)/2;

hd=ideal_lp(wc,N);            %理想低通滤波器冲激响应

window=(kaiser(N,beta))';     %凯泽窗

h=hd.*window;                 %实际设计低通滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1)))

As=-round(max(db(ws/delta_w+1:1:501)))

subplot(221)

stem(n,hd);title('理想滤波器单位抽样响应');

axis([0,N-1,-0.1,0.3]);

xlabel('n');ylabel('hd(n)');

subplot(222)

stem(n,window);title('kaiser Window')

axis([0 N-1 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223);

stem(n,h);title('所设计的滤波器脉冲响应');

axis([0 N-1 -0.1 0.3]);xlabel('n');ylabel('h(n)')

subplot(224);

plot(w/pi,db);

title('频率响应');grid

axis([0 1 -100 10]);xlabel('w/pi');ylabel('H(w)');

hold on;

plot([0,0.4],[-50,-50],'r');

text(0.3,-45,num2str(as));

结果

c. 设计FIR 数字带通滤波器，技术指标为：

下阻带边缘：ωst1=0.2π，δs1=60dB，下通带边缘：ωp1=0.35π，δp1=1dB；

上通带边缘：ωp2=0.65π，δp1=1dB，上阻带边缘：ωst2=0.8π，δs2=60dB；

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

clc;

%用布莱克曼窗设计带通滤波器

wp1=0.35*pi;

ws1=0.2*pi;

wp2=0.65*pi;

ws2=0.8*pi;

As=60;

tr_width=min((wp1-ws1),(wp2-ws2));

N=ceil(11*pi/tr_width)+1;

N=-N;

n=[0:1:N-1];

wc1=(ws1+wp1)/2;

wc2=(ws2+wp2)/2;

hd=ideal_lp(wc2,N)-ideal_lp(wc1,N);  %理想带通滤波器冲激响应

window=(blackman(N))';              %布莱克曼窗

h=hd.*window;                     %实际设计带通滤波器冲激响应  

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w));

As=-round(max(db(ws2/delta_w+1:1:501)));

subplot(221);

stem(n,hd);title('理想带通滤波器冲激响应');

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)');

subplot(222);

stem(n,window);title('Blackman Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计带通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224);

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('frequence in \\pi unit');ylabel('decibels');

hold on;

plot([0,0.3],[-74,-74],'r');

text(0.2,-70,num2str(As));

结果：

d. 设计FIR 数字带通滤波器，技术指标为：

下阻带边缘：ωst1=0.2π，δs1=60dB，下通带边缘：ωp1=0.4π，δp1=1dB；

上通带边缘：ωp2=0.6π，δp1=1dB，上阻带边缘：ωst2=0.8π，δs2=60dB；

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

clc;

%用布莱克曼窗设计带通滤波器

wp1=0.4*pi;

ws1=0.2*pi;

wp2=0.6*pi;

ws2=0.8*pi;

As=60;

tr_width=min((wp1-ws1),(wp2-ws2));

N=ceil(11*pi/tr_width)+1;

N=-N;

n=[0:1:N-1];

wc1=(ws1+wp1)/2;

wc2=(ws2+wp2)/2;

hd=ideal_lp(wc2,N)-ideal_lp(wc1,N)  %理想带通滤波器冲激响应

window=(blackman(N))';             %布莱克曼窗

h=hd.*window;                      %实际设计带通滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w));

As=-round(max(db(ws2/delta_w+1:1:501)));

subplot(221);

stem(n,hd);title('理想带通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('Blackman Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计带通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('frequence in \\pi unit');ylabel('decibels');

hold on;

plot([0,0.3],[-74,-74],'r');

text(0.2,-70,num2str(As));

结果：

e. 设计FIR 数字带通滤波器，技术指标为：

下阻带边缘：ωst1=0.2π，δs1=20dB，下通带边缘：ωp1=0.4π，δp1=1dB；

上通带边缘：ωp2=0.6π，δp1=1dB，上阻带边缘：ωst2=0.8π，δs2=20dB；

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

clc;

%用矩形窗设计带通滤波器

wp1=0.4*pi;

ws1=0.2*pi;

wp2=0.6*pi;

ws2=0.8*pi;

As=20;

tr_width=min((wp1-ws1),(wp2-ws2));

N=ceil(11*pi/tr_width)+1;

N=-N;

n=[0:1:N-1];

wc1=(ws1+wp1)/2;

wc2=(ws2+wp2)/2;

hd=ideal_lp(wc2,N)-ideal_lp(wc1,N)  %理想带通滤波器冲激响应

window=(boxcar(N))'; %矩形窗

h=hd.*window;        %实际设计带通滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w));

As=-round(max(db(ws2/delta_w+1:1:501)));

subplot(221);

stem(n,hd);title('理想通带滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('Rectangle Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计通带滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

hold on;

plot([0,0.3],[-21,-21],'r');

text(0.2,-21,num2str(As));

结果：

f. 设计FIR 数字高通滤波器，技术指标为：通带截止频率为ωp=15π/27，阻带截止频率

为ωst=11π/27，通带最大衰减为δ1=2.5dB，阻带最小衰减为δ2=55dB。

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

clc;

%用布莱克曼窗设计高通滤波器

wp=15/27*pi;

ws=11/27*pi;

As=55;

tr_width=wp-ws;

N=ceil(11*pi/tr_width)+1;

n=[0:1:N-1];

wc=(ws+wp)/2;

hd=ideal_lp(pi,N)-ideal_lp(wc,N)  %理想高通滤波器冲激响应

window=(blackman(N))';           %布莱克曼窗

h=hd.*window;                    %实际设计高通滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-min(db(1:1:wp/delta_w+1));

As=-round(max(db(ws/delta_w+1:1:501)));

subplot(221);

stem(n,hd);title('理想高通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('blackman Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计高通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

结果：

g. 设计FIR 数字高通滤波器，技术指标为：通带截止频率为ωp=0.6π，阻带截止频率为

ωst=0.4π，通带最大衰减为δ1=0.25dB，阻带最小衰减为δ2=40dB。

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

clc;

%用海明窗设计带通滤波器

wp=0.6*pi;

ws=0.4*pi;

As=40;

tr_width=wp-ws;

N=ceil(11*pi/tr_width)+1;

n=[0:1:N-1];

wc=(ws+wp)/2;

hd=ideal_lp(pi,N)-ideal_lp(wc,N)  %理想高通滤波器冲激响应

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

h=hd.*window;                 %实际设计高通滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

Rp=-min(db(1:1:wp/delta_w+1));

As=-round(max(db(ws/delta_w+1:1:501)));

subplot(221);

stem(n,hd);title('理想高通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('hamming Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计高通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

结果：

h. 滤波器的技术指标为：通带截止频率为ωp=0.6π，阻带截止频率为ωst=0.4π，通带最

大衰减为δ1=0.25dB，阻带最小衰减为δ2=40dB。

(1) 通过技术指标，选择一种窗函数设计一个具有π/2 相移的FIR 高通滤波器；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

function hd=ideal_pi(wc,N)  %有pi/2相移的滤波器冲激响应

alpha=(N-1)/2;

n=[0:1:(N-1)];

n=n-alpha;

fc=wc/pi

hd=fc*sinc(fc*(n-1/2));

function[db,mag,pha,w]=freqz_m2(b,a)

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501));

w=(w(1:1:501));

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

主函数：

%设计高通滤波器

clc;

Wp=0.6*pi;

Ws=0.4*pi;

tr_width=Wp-Ws;

M=ceil(6.2*pi/tr_width)+1;

n=0:1:M-1;

Wc=(Ws+Wp)/2;

hd=ideal_pi(pi,M)-ideal_pi(Wc,M);  %理想高通滤波器冲激响应

w_ham=(hanning(M))';                %汉宁窗

h=hd.*w_ham;                        %实际设计高通滤波器冲激响应 

[db,mag,pha,w]=freqz_m2(h,[1]);

delta_w=2*pi/1000;

Ap=-(min(db(Wp/delta_w+1:1:501)))

As=-round(max(db(1:1:Ws/delta_w+1)))

subplot(221)

stem(n,hd);title('理想高通滤波器冲激响应')

axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222)

stem(n,w_ham);title('Hanning Window')

axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')

subplot(223)

stem(n,h);title('实际设计高通滤波器冲激响应')

axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');grid

axis([0 1 -100 10]);

xlabel('w/pi');ylabel('幅度');

结果：

i. 设计FIR 数字带阻滤波器，其技术指标为：

低端阻带边缘：ωst1=0.4π，δs1=40dB，低端通带边缘：ωp1=0.2π，δp1=1dB；

高端通带边缘：ωp2=0.8π，δp1=1dB，高端阻带边缘：ωst2=0.6π，δs2=40dB；

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

代码：

function hd=ideal_bs(Wcl,Wch,m);  %带阻滤波器冲激响应

alpha=(m-1)/2;

n=[0:1:(m-1)];

m=n-alpha+eps;

hd=[sin(m*pi)+sin(Wcl*m)-sin(Wch*m)]./(pi*m)

function[db,mag,pha,w]=freqz_m2(b,a)

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501));

w=(w(1:1:501));

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

主函数：

%设计带阻滤波器

clc;

wpl=0.2*pi;

wph=0.8*pi;

wsl=0.4*pi;

wsh=0.6*pi;

tr_width=min((wsl-wpl),(wph-wsh));

M=ceil(6.2*pi/tr_width);        

n=0:1:M-1;

wcl=(wsl+wpl)/2;

wch=(wsh+wph)/2;

hd=ideal_bs(wcl,wch,M);   %理想带阻滤波器冲激响应

w_ham=(hanning(M))';      %汉宁窗

h=hd.*w_ham;               %实际设计带阻滤波器冲激响应 

[db,mag,pha,w]=freqz_m2(h,[1]);

delta_w=2*pi/1000;

Ap=-(min(db(1:1:wpl/delta_w+1)))

As=-round(max(db(wsl/delta_w+1:1:wsh/delta_w+1)))

subplot(221)

stem(n,hd);title('理想带阻滤波器冲激响应')

axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('hd(n)')

subplot(222)

stem(n,w_ham);title('Hanning Window')

axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')

subplot(223)

stem(n,h);title('实际设计带阻滤波器冲激响应')

axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');grid

axis([0 1 -100 10]);

xlabel('w/pi');ylabel('幅度');

结果：

j. 设计FIR 数字带阻滤波器，其频率响应函数为：

其中阻带衰减为 50dB。

(1) 通过技术指标，选择一种窗函数进行设计；

(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。

(3) 用凯塞窗设计设计该滤波器，并绘制相应的波形图。

代码：

%设计带阻滤波器

clc;

M=45; As=50; 

n=[0:1:M-1];

beta=0.1102*(As-8.7);

w_kai=(kaiser(M,beta))';         %凯泽窗

wc1=pi/3; wc2=2*pi/3;

hd=ideal_lp(wc1,M)+ideal_lp(pi,M)-ideal_lp(wc2,M); %理想带阻滤波器冲激响应

h=hd.*w_kai;                             %实际设计带阻滤波器冲激响应 

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(2,2,1); 

stem(n,hd); title('理想带阻滤波器冲激响应')

axis([-1 M -0.2 0.8]);

xlabel('n'); ylabel('hd(n)')

subplot(2,2,2); 

stem(n,w_kai);title('Kaiser Window')

axis([-1 M 0 1.1]); xlabel('n'); ylabel('w(n)')

subplot(2,2,3);

stem(n,h);title('实际设计带阻滤波器冲激响应 ')

axis([-1 M -0.2 0.8]); xlabel('n'); ylabel('h(n)')

subplot(2,2,4);

plot(w/pi,db); axis([0 1 -80 10]);

title('滤波器频谱响应');grid;

xlabel('w/pi');ylabel('幅度');

结果：

k. 设计FIR 数字低通滤波器，其滤波器的截止频率为100Hz，滤波器的阶次为80。分

别用矩形窗、三角形窗、汉宁窗、海明窗、布莱克曼窗设计数字滤波器，并绘制相应的窗函

数波形和幅度频率响应波形。

代码：

%数字低通滤波器

clc;

%用矩形窗设计带通滤波器

wc=100;

N=80;

n=[0:1:N-1];

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

window=(boxcar(N))'; %矩形窗

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('Rectangle Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

%用巴特列特窗设计带通滤波器

figure;

wc=100;

N=80;

n=[0:1:N-1];

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

window=(bartlett(N))'; %巴特列特窗

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('bartlett Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

%用汉宁窗设计带通滤波器

figure;

wc=100;

N=80;

n=[0:1:N-1];

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

window=(hanning(N))'; %汉宁窗

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('hanning Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

%用海明窗设计带通滤波器

figure;

wc=100;

N=80;

n=[0:1:N-1];

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

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

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('hammimg Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

%用布拉克曼窗设计带通滤波器

figure;

wc=100;

N=80;

n=[0:1:N-1];

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

window=(blackman(N))'; %布拉克曼窗

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('blackman Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

%用凯泽窗设计带通滤波器

figure;

wc=100;

N=80;

n=[0:1:N-1];

beta=9;

hd=ideal_lp(wc,N) %理想低通滤波器冲激响应

window=(kaiser(N,beta))'; %凯泽窗

h=hd.*window;    %实际设计通带滤波器冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(221);

stem(n,hd);title('理想低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,window);title('kaiser Window');

axis([-10 N+9 0 1.1]);xlabel('n');ylabel('w(n)');

subplot(223)

stem(n,h);title('实际设计低通滤波器冲激响应')

axis([0 N-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')

subplot(224)

plot(w/pi,db);title('滤波器频谱响应');

grid

axis([0 1 -150 10]);

xlabel('w/pi');ylabel('幅度');

结果：

分析：从上面几个图可以看出，不同的窗函数得出的滤波器的阻带衰减的程度不一样，矩形窗效果差，三角形窗、汉宁窗、海明窗、布莱克曼窗越来越好，而凯泽窗效果最好。

l. .FIR 滤波器的单位抽样响应为

h(n)=1/9* {1,2,3,2,1}

                 ↑

编制MATLAB 程序求系统的频率采样型结构的系数，并画出频率抽样型结构。

代码：

clc;

h=[1,2,3,2,1];

M=5;

n=(M-3)/2;

a=2*h(3-n)

k=-500:500;

K=500;

w=k*pi/K;

p=a*cos(w*n);

plot(w,p);

ylabel('p(w)'); xlabel('w');

 结果：

m. 一个理想差分器的频率响应为：

用长度为 21 的汉宁窗设计一个数字FIR 差分器，并绘制其时域和频率的响应波形。

代码：

function [Hr,w,c,L] = Hr_Type3(h);

M=length(h);

L=(M-1)/2;

c=[2*h(L+1:-1:1)];

n=[0:1:L];

w=[0:1:500]'*pi/500;

Hr=sin(w*n)*c';

主函数：

%差分器

clc;

M=21; 

alpha=(M-1)/2; n=0:M-1;

hd=(cos(pi*(n-alpha)))./(n-alpha); %理想差分器冲激响应

hd(alpha+1)=0;

w_ham=(hamming(M))'; %海明窗

h=hd.*w_ham;         %实际差分器冲激响应

[Hr,w,P,L]=Hr_Type3(h);

subplot(221); 

stem(n,hd); title('理想差分器冲激响应')

axis([-1 M -1.2 1.2]);

xlabel('n'); ylabel('hd(n)')

subplot(222);

stem(n,w_ham);

title('Hamming Window')

axis([-1 M 0 1.2]); xlabel('n'); ylabel('w(n)')

subplot(223); 

stem(n,h);title('实际差分器冲激响应')

axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('h(n)')

subplot(224);

plot(w/pi,Hr/pi);

title('差分器频谱响应');grid;

xlabel('w/pi'); ylabel('幅度');

axis([0 1 0 1]);

结果：

n. 利用汉宁窗设计一个长度为25 的数字希尔伯特变换器，并绘制它的时域和频域的响

应波形。

代码：

%希尔伯特变换器

clc;

M=25;

alpha=(M-1)/2;

n=0:M-1;

hd=(2/pi)*((sin((pi/2)*(n-alpha)).^2)./(n-alpha));

hd(alpha+1)=0;

w_han=(hann(M))';

h=hd.*w_han; 

[Hr,w,P,L]=Hr_Type3(h);

subplot(2,2,1); 

stem(n,hd); title('Ideal Impulse Response')

axis([-1 M -1.2 1.2]);

xlabel('n'); ylabel('hd(n)');

subplot(2,2,2); 

stem(n,w_han);title('Hann Window')

axis([-1 M 0 1.2]); xlabel('n'); ylabel('w(n)')

subplot(2,2,3); 

stem(n,h);title('Actual Impulse Response')

axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('h(n)')

w=w';

Hr=Hr';

w=[-fliplr(w),w(2:501)]; 

Hr=[-fliplr(Hr), Hr(2:501)];

subplot(2,2,4);plot(w/pi,Hr); title('Amplitude Response');grid;

xlabel('frequency in pi units'); ylabel('Hr');

axis([-1 1 -1.1 1.1]);

结果：

p. FIR 数字低通滤波器的技术指标为：ωp=0.2π，ωst=0.3π，δ1=0.25dB，δ2=50dB。

利用频率采样方法设计 FIR 数字滤波器，并绘制滤波器的单位冲激响应、幅度频率响应

的波形。

代码：

function [Hr,w,b,L] = Hr_Type2(h);

M=length(h); L=M/2;

b=2*[h(L:-1:1)]; 

n=[1:1:L]; 

n=n-0.5;

w=[0:1:500]'*pi/500; Hr=cos(w*n)*b';

function [db,mag,pha,grd,w]=freqz_m(b,a)

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501));

w=(w(1:1:501));

mag=abs(H);

db=20*log10(mag+eps/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

主函数：

clc;

M=20; 

alpha=(M-1)/2;

l=0:M-1;

wl=(2*pi/M)*l;

Hrs=[1,1,1,zeros(1,15),1,1];  %Ideal Amp Res sampled

Hdr=[1,1,0,0];

wdl=[0,0.25,0.25,1]; %Ideal Amp Res for plotting

k1=0:floor((M-1)/2); 

k2=floor((M-1)/2)+1:M-1;

angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)];

H=Hrs.*exp(j*angH); 

h=real(ifft(H,M));

[db,mag,pha,grd,w]=freqz_m(h,1); 

[Hr,ww,a,L]=Hr_Type2(h);

subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);

axis([0,1,-0.1,1.1]);

title('Frequency Samples: M=20')

xlabel('frequency in pi units'); ylabel('Hr(k)');

subplot(2,2,2); 

stem(l,h); 

axis([-1,M,-0.1,0.3])

title('Impulse Response');

xlabel('n'); ylabel('h(n)');

subplot(2,2,3);

plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');

axis([0,1,-0.2,1.2]); title('Amplitude Response')

xlabel('frequency in pi units'); ylabel('Hr(w)')

subplot(2,2,4);

plot(w/pi,db); 

axis([0,1,-60,10]); grid

title('Magnitude Response'); xlabel('frequency in pi units');

ylabel('Decibels');

结果：

q. 用窗函数法设计一个线性相位的FIR 数字低通滤波器，其技术指标为：ωp=0.2π，

ωst=0.4π，δ1=0.25dB，δ2=50dB。

代码：

function[db,mag,pha,H,w]=freqz_m3(b,a)

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501))'; w=(w(1:1:501))';

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

主函数：

clc;

wp=0.2*pi;

ws=0.4*pi;

tr_width=ws-wp;

M=ceil(6.6*pi/tr_width)+1;

n=[0:1:M-1];

wc=(ws+wp)/2;

hd=ideal_lp(wc,M);

w_ham=(hamming(M))';   %海明窗

h=hd.*w_ham;

[db,mag,pha,H,w]=freqz_m3(h,[1]);

delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1)));

As=-round(max(db(ws/delta_w+1:1:501)));

figure(1);

subplot(221);

stem(n,hd);title('Ideal Impulse Rresponse')

axis([0 M-1 -0.1 0.3]);

xlabel('n');ylabel('hd(n)')

subplot(222);

stem(n,w_ham);title('Hamming Window');

axis([0 M-1 0 1.1]);

xlabel('n');ylabel('w(n)');

subplot(223);

stem(n,h);title('Actual Impuse Response');

axis([0 M-1 -0.1 0.3]);

xlabel('n');ylabel('h(n)');

subplot(224);

plot(w/pi,db);

title('Magnitude Response in db');grid

axis([0 1 -100 10]);

xlabel('frequence in \\pi unit');ylabel('decibels');

%solution of problem (b)

Rn=ones(1,10);

n=0:1:9;

dw=w(2)-w(1);

w=[-fliplr(w),w(2:501)];

H=[fliplr(H),H(2:501)];

xw=Rn*(exp(-j)).^(n'*w);

answ=H.*xw;

Rnn=answ*exp(-j).^(w'*n)*dw/(2*pi);

figure(2);

subplot(221);

stem(n,Rn);

xlabel('n');ylabel('Rn');

title('input rectangle signal in T-domain');

subplot(222)

stem(w/pi,xw)

xlabel('w/pi');ylabel('Rw');

title('input rectangle signal in F-domain');

subplot(223);

plot(w/pi,abs(answ))

xlabel('w/\\pi');ylabel('Magnitude of H(w)');

title('output signal in F-domain after filteration');

subplot(224);

stem(n,abs(Rnn));

xlabel('n');ylabel('Rnn');

title('output signal in T-domain after filteration');

结果：

实验总结：

    这个实验都是滤波器的知识，通过这个实验让我明白了滤波器是如何设计的。这个实验的内容有点多，在做的过程中花了较多的时间。这个实验也让我知道了用不同的窗设计滤波器的效果是不一样的；知道了如何去使用这些窗函数；针对不同的阻带衰减要用不同的窗函数。凯泽窗的效果是最好的，它可以使阻带衰减从-30到-100，当beta在10左右时，阻带衰减可以到达100dB。数字差分器和希尔伯特变换器也都可以有窗函数来设计。由此可知窗函数的应用是比较广泛的。