窗函数设计低通滤波器 电信课设.docx

窗函数设计低通滤波器 电信课设.docx

《窗函数设计低通滤波器 电信课设.docx》由会员分享，可在线阅读，更多相关《窗函数设计低通滤波器 电信课设.docx（23页珍藏版）》请在冰豆网上搜索。

窗函数设计低通滤波器电信课设

XXXX大学

课程设计报告

学生：

xxx学号：

xxx

专业班级：

电子信息工程

课程名称：

数字信号处理课程设计

学年学期20XX——20XX学年第X学期

指导教师：

xxx

2014年6月

课程设计成绩评定表

学生

XXX

学号

XXXXXX

成绩

专业班级

XXXXX

起止时间

20XX-X-X至20XX-X-XX

设计题目

1.窗函数设计低通滤波器

2.用哈明窗设计FIR带通数字滤波器

指

导

教

师

评

语

指导教师：

年月日

1.窗函数设计低通滤波器

1.1设计目的……………………………………………………………………1

1.2设计原理推导与计算………………………………………………………1

1.3设计容与要求……………………………………………………………2

1.4设计源程序与运行结果……………………………………………………3

1.5思考题……………………………………………………………………10

2.用哈明窗设计FIR带通数字滤波器

2.1设计要求……………………………………………………………………14

2.2设计原理和分析……………………………………………………………14

2.3详细设计……………………………………………………………………15

2.4调试分析及运行结果……………………………………………………15

2.5心得体会…………………………………………………………………17

参考文献………………………………………………………………………17

1.窗函数设计低通滤波器

1.1设计目的

1.熟悉设计线性相位数字滤波器的一般步骤。

2.掌握用窗函数法设计FIR数字滤波器的原理和方法。

3.熟悉各种窗函数的作用以及各种窗函数对滤波器特性的影响。

4.学会根据指标要求选择合适的窗函数。

1.2设计原理推导与计算

如果所希望的滤波器的理想的频率响应函数为

，则其对应的单位脉冲响应为

（4.1）

窗函数设计法的基本原理是设计设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数

，即

其中

用有限长单位脉冲响应序列

逼近

。

由于

往往是无限长序列，而且是非因果的，所以用窗函数

将

截断，并进行加权处理，得到：

（4.2）

就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数

为

（4.3）

式中，N为所选窗函数

的长度。

用窗函数法设计的滤波器性能取决于窗函数

的类型及窗口长度N的取值。

设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。

各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表

（一）。

窗函数

旁瓣峰值幅度/dB

过渡带宽

阻带最小衰减/dB

过渡带带宽

矩形窗

-13

4π/N

-12

0.9

三角形窗

-25

8π/N

-25

2.1

汉宁窗

-31

8π/N

-44

3.1

哈明窗

-41

8π/N

-53

3.3

不莱克曼窗

-57

12π/N

-74

5.5

凯塞窗（

=7.865）

-57

10π/N

-80

5

表

（一）各种窗函数的基本参数

这样选定窗函数类型和长度N之后，求出单位脉冲响应

，并按照式（4.3）求出

。

是否满足要求，如果

不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。

1.3设计容与要求

（一）设计要求：

1.学会计算滤波器各项性能指标及如何来满足给定的指标要求。

2.用MATLAB语言编程实现给定指标要求的滤波器的设计。

3.熟悉MATLAB语言，独立编写程序。

4.设计低通FIR滤波器的指标：

通带最大波动

，

阻带最小衰减

，

（二）、设计容：

1.熟悉各种窗函数,在MATLAB命令窗下浏览各种窗函数，绘出（或打印）所看到的窗函数图。

2.编写计算理想低通滤波器单位抽样响应hd（n）的m函数文件ideal.m。

3.编写计算N阶差分方程所描述系统频响函数

的m函数文件fr.m。

4.根据指标要求选择窗函数的形状与长度N。

（至少选择两种符合要求的窗函数及其对应的长度）。

5.编写.m程序文件，通过调用ideal.m和fr．m文件，计算你设计的实际低通FIR滤波器的单位抽样响应h（n）和频率响应

，打印在频率区间[O，π]上的幅频响应特性曲线

，幅度用分贝表示。

6.验证所设计的滤波器是否满足指标要求。

1.4设计的源程序及运行结果：

1、利用MATLAB窗口观察各种窗函数：

%巴特利特窗

w=bartlett（20）;

subplot（3,2,1）;

plot（w）;

stem（w,'y'）;%'y'表示黄色

%stem表示以离散图输出

title（'巴特利特床窗'）;

xlabel（'n'）;%横坐标为n

ylabel（'w（n）'）;%纵坐标为w（n）

%布莱克曼窗

w=blackman（20）;

subplot（3,2,2）;

plot（w）;

stem（w,'b'）;%'b'表示蓝色

title（'布莱克曼窗'）;

xlabel（'n'）;

ylabel（'w（n）'）;

%矩形窗

w=boxcar（20）;

subplot（3,2,3）;

plot（w）;

stem（w,'r'）;

title（'矩形窗'）;

xlabel（'n'）;

ylabel（'w（n）'）;

%海明窗

w=hamming（20）;

plot（w）;

stem（w,'m'）;%'m'表示紫色

title（'海明窗'）;

xlabel（'n'）;

ylabel（'w（n）'）;

%汉宁窗

w=hanning（20）;

subplot（3,2,5）;

plot（w）;

stem（w,'g'）;%'g'表示绿色

title（'汉宁窗'）;

xlabel（'n'）;

ylabel（'w（n）'）;

%凯泽窗

beta=5.6533;

w=kaiser（20,beta）;

subplot（3,2,6）;

plot（w）;

stem（w,'k'）;%'k'表示黑色

title（'凯泽窗,beta=5.6533'）;

xlabel（'n'）;

ylabel（'w（n）'）;

常用窗函数的图形

2、理想低通滤波器单位抽样响应hd（n）的m函数文件ideal.m。

functionhd=ideal（wc,M）

%理想低通滤波器计算

%hd为0到M-1之间的理想脉冲响应

%wc为截止频率

%M为理想滤波器的长度

alpha=（M-1）/2;

n=0:

M-1;

m=n-alpha+eps;

hd=sin（wc*m）./（pi*m）;

3、N阶差分方程所描述的系统频响函数的m函数文件fr.m。

function[db,mag,pha,gfd,w]=fr（b,a）

%求解系统响应

%db为相位振幅（db）

%mag为绝对振幅

%pha为相位响应

%grd为群延时

%w为频率样本点矢量

%b为Ha（z）分析多项式系数（对FIR而言，b=h）

%a为Hz（z）分母多项式系数（对FIR而言，a=1）

[H,w]=freqz（b,a,1000,'whole'）;

H=（H（1:

501））';

w=（w（1:

501））';

mag=abs（H）;

db=20*log10（（mag+eps）/max（mag））;

pha=angle（H）;

gfd=grpdelay（b,a,w）;

4、实际低通滤波器FIR：

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

wp=0.2*pi;

ws=0.3*pi;

tr_width=ws-wp;

disp（['海明窗设计低通滤波器参数:

']）;

M=ceil（6.6*pi/tr_width）+1;

disp（['滤波器的长度为',num2str（M）]）;

n=0:

M-1;

wc=（ws+wp）/2;%理想LPF的截止频率

hd=ideal（wc,M）;

w_ham=（hamming（M））';

h=hd.*w_ham;

[db,mag,pha,gfd,w]=fr（h,[1]）;

delta_w=2*pi/1000;

Rp=-（min（db（1:

1:

wp/delta_w+1）））;%求出实际通带波动

disp（['实际带通波动为',num2str（Rp）]）;

As=-round（max（db（ws/delta_w+1:

1:

501）））;%求出最小阻带衰减

disp（['最小阻带衰减为-',num2str（As）],’db’）;

%绘图

subplot（1,1,1）

subplot（2,6,1）

stem（n,hd）;

title（'理想冲击响应'）;

axis（[0M-1-0.10.3]）;

ylabel（'hd（n）'）;

subplot（2,6,2）

stem（n,w_ham）;

title（'海明窗'）;

axis（[0M-101.1]）;

ylabel（'w（n）'）;

subplot（2,6,7）

stem（n,h）;

title（'实际冲激响应'）;

axis（[0M-1-0.10.3]）;

xlabel（'n'）;

ylabel（'h（n）'）;

subplot（2,6,8）

plot（w/pi,db）;

title（'幅度响应（db）'）;

axis（[01-10010]）;

grid;

xlabel（'以pi为单位的频率'）;

ylabel（'分贝数'）;

图

（1）海明窗设计的FIR

海明窗设计低通滤波器参数:

滤波器的长度为67

实际带通波动为0.03936

最小阻带衰减为-52db

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

wp=0.2*pi;

ws=0.3*pi;

tr_width=ws-wp;

disp（['布莱克曼窗设计低通滤波器的参数:

']）;

M=ceil（11.0*pi/tr_width）+1;

disp（['滤波器的长度为',num2str（M）]）;

n=0:

M-1;

%理想LPF的截止频率

wc=（ws+wp）/2;

hd=ideal（wc,M）;

w_bla=（blackman（M））';

h=hd.*w_bla;

[db,mag,pha,gfd,w]=fr（h,[1]）;

delta_w=2*pi/1000;

Rp=-（min（db（1:

1:

wp/delta_w+1）））;%求出实际通带波动

disp（['实际带通波动为',num2str（Rp）]）;

As=-round（max（db（ws/delta_w+1:

1:

501）））;%求出最小阻带衰减

disp（['最小阻带衰减-',num2str（As）],’db’）;

%绘图

subplot（2,6,3）

stem（n,hd）;

title（'理想冲击响应'）;

axis（[0M-1-0.10.3]）;

ylabel（'hd（n）'）;

subplot（2,6,4）

stem（n,w_bla）;

title（'布莱克曼窗'）;

axis（[0M-101.1]）;

ylabel（'w（n）'）;

subplot（2,6,9）

stem（n,h）;

title（'实际冲激响应'）;

axis（[0M-1-0.10.3]）;

xlabel（'n'）;

ylabel（'h（n）'）;

subplot（2,6,10）

plot（w/pi,db）;

title（'幅度响应（db）'）;

axis（[01-10010]）;

grid;

xlabel（'以pi为单位的频率'）;

ylabel（'分贝数'）;

图

（2）布莱克曼窗设计的FIR

布莱克曼窗设计低通滤波器的参数:

滤波器的长度为111

实际带通波动为0.0033304

最小阻带衰减为-73db

5、技术指标比较：

（1）海明窗设计低通滤波器参数:

滤波器的长度为67

实际带通波动为0.03936

最小阻带衰减为-52db

（2）布莱克曼窗设计低通滤波器的参数:

滤波器的长度为111

实际带通波动为0.0033304

最小阻带衰减为-73db

在相同的技术指标下用布莱克曼窗设计的低通滤波器实际带通波动实际带通波动最小，最小阻带衰减，滤波器的长度最大；海明窗和凯泽窗最小阻带衰减差不多，滤波器的长度页差不多，但是海明窗实际波动小于凯泽窗；所以用布莱克曼窗用设计的FIR最逼近理想单位冲击响应。

这三个窗设计的低通滤波器都符合要求。

1.5思考题：

1.设计线性相位数字滤波器的一般步骤：

（1）、给定所要求的频率响应函数Hd（

）以及技术指标阻带衰减

;

（2）、求hd=IDTFT[Hd（

）];

由过渡带带宽及带阻最小衰减的要求，利用六种常见的窗函数基本参数的比较表或凯泽窗设计FIR滤波器的经验公式，选择窗函数的形状及N的大小（一般N要通过几次试探而最后确定）；

（3）求得所设计的FIR滤波器的单位抽样响应：

h（n）=hd（n）w（n）,N=0,1,2,3…,N-1

（4）、求Hd（

）=DIFT[h（n）],校验是否满足设计要求，如果不满足，则重新设计。

2.窗函数有哪些指标要求?

对给定指标要求的低通滤波器，理论计算所需窗函数的长度N。

答:

窗函数的指标要求：

主瓣宽度，旁瓣峰值。

海明窗设计低通滤波器的长度为67，凯泽窗设计低通滤波器的长度为60，布莱克曼窗设计低通滤波器的长度为111。

3.用窗函数法设计FIR滤波器，滤波器的过渡带宽度和阻带衰减与哪些因素有关？

答：

过渡带宽度与窗函数的形状和窗的宽度有关；阻带衰减只有窗函数的形状决定，不受N的影响。

4、计算理想带通滤波器单位抽样响应hd（n）时取N为奇数和N为偶数有没有区别？

你计算时所用的方法是仅适合于N为奇数或偶数的一种还是两种都可以用？

答：

以海明窗为例说明：

滤波器的长度为67

实际带通波动为0.03936

最小阻带衰减为52

滤波器的长度为68

实际带通波动为0.036424

最小阻带衰减为53

滤波器的长度为64

实际带通波动为0.068677

最小阻带衰减为45

N取奇偶数时，低通滤波器的幅度函数是不同的，如上图所示（海明窗设计低通滤波器N取奇偶数时的图），通过比较，当滤波器的长度大于技术指标要求的长度时，选择偶数也是满足要求的，所以海明窗奇偶都满足。

同理，N取奇偶数，选择其他的窗函数也满足设计低通滤波器的指标要求。

比较所选窗长N相同但窗形状不同对滤波器设计结果的影响以及选同一种窗函数但窗长N不同时对滤波器设计结果的影响，将结论写在报告中。

海明窗

布莱克曼窗

凯泽窗

以海明窗设计的低通滤波器的长度N=67，用布莱克曼窗和凯泽窗设计低通滤波器

滤波器的长度为67

实际带通波动0.03936

最小阻带衰减为52

滤波器的长度为67

实际带通波动0.29523

最小阻带衰减为30

滤波器的长度为67

实带通波动0.047919

最小阻带衰减为51

线性斜率beta=4.5513

以凯泽窗设计的低通滤波器的长度N=60，用布莱克曼窗和海明窗窗设计低通滤波器

滤波器的长度为60

实际带通波动为0.12166

最小阻带衰减为39

滤波器的长度为60

实际带通波动为0.47079

最小阻带衰减为26

滤波器的长度为60

实际带通波动0.053711

最小阻带衰减为51

线性斜率beta=4.5513

以布莱克曼窗设计的低通滤波器的长度N=111，用凯泽窗和海明窗窗设计低通滤波器

滤波器的长度为111

实际带通波动为0.02894

最小阻带衰减为59

滤波器的长度为111

实际带通波动0.0033304

最小阻带衰减为73

滤波器的长度为111

实际带通波动0.034426

最小阻带衰减为58

线性斜率beta=4.5513

结论：

1、当以海明窗设计的低通滤波器的长度N，用布莱克曼窗和凯泽窗设计低通滤波器，即N=67时，三个窗函数设计的低通滤波器中布莱克曼窗设计的不满足要求，凯泽窗设计的带通波动大于海明窗设计的低通滤波器带通波动，最小阻带衰减小于海明窗设计的低通滤波器。

凯泽窗设计的结果没有用海明窗设计结果好。

当N取67、60、111时，用海明窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<67时不满足设计技术指标。

2、以凯泽窗设计的低通滤波器的长度N，用布莱克曼窗和海明窗窗设计低通滤波器，即N=60时，海明窗设计的低通滤波器的最小阻带衰减不满足指标要求，布莱克曼窗设计的低通滤波器的最小阻带衰减和带通波动都不满足要。

当N取67、60、111时，用布莱克曼窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<111时不满足设计技术指标。

3、以布莱克曼窗设计的低通滤波器的长度N，用凯泽窗和海明窗窗设计低通滤波器，即N=111时，三个都满足要求。

但是用布莱克曼窗函数设计的低通滤波器的性能更好。

当N取67、60、111时，用凯泽窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N>60时都满足设计技术指标。

2、用哈明窗设计FIR带通数字滤波器

2.1设计要求

针对一个含有5Hz、15Hz和30Hz的混合正弦波信号

，设计滤波系统，滤除5Hz和30Hz的正弦分量，阻带的最小衰减不小于50dB。

①确定采样频率

、采样长度N，对

进行采样；

②选择合适的窗函数，设计FIR带通数字滤波器；

③通过计算机仿真对滤波器的性能进行分析。

2.2设计原理分析和设计

此题要求也是滤波，但是要求滤掉5Hz和30Hz的信号，而保留15Hz的信号，因此必须采用带通滤波器。

要求阻带的最小衰减不低于50dB，查表可知可选哈明窗。

分析题目可知：

f1=5Hz，T1=1/5s；

f2=15Hz，T2=1/15s；

f3=30Hz，T3=1/30s；

则采样时长tT至少应为0.2s，取tT=1s；

采样频率fs≥2fc；可取fs=150；则采样间隔T=1/fs；

所以采样长度N=tT/T=150；

数字技术指标可取

wp1=2*pi*14/fs;

wp2=2*pi*23/fs;

ws1=2*pi*7/fs;

ws2=2*pi*27/fs;

滤波器设计好后，利用卷积和可得输出信号y（n）=x（n）*h（n），卷积和长度N=N1+N2-1;

2.3详细设计

源程序代码：

tT=1;%采样时长

fs=150;%采样频率

wp1=2*pi*14/fs;

wp2=2*pi*20/fs;

ws1=2*pi*7/fs;

ws2=2*pi*27/fs;

trwidth=wp1-ws1;%过渡带带宽

N1=ceil（8*pi/trwidth）%计算滤波器阶次

n=0:

N1-1;

wc1=（ws1+wp1）/2;wc2=（ws2+wp2）/2;%计算3dB截止频率

alpha=（N1-1）/2;%单位响应的对称中心

m=n-alpha;

hd=sin（wc2*m）./（pi*m）-sin（wc1*m）./（pi*m）;%理想带通滤波器的单位响应

wham=（hamming（N1））';

h=hd.*wham;%实际带通滤波器单位响应

w=0:

0.01:

pi;

H=freqz（h,1,w）;%单位滤波器幅频响应

T=1/fs;%对输入信号进行采样

N2=fs*tT;

t=（0:

N2-1）*T;

x=sin（2*pi*5*t）+sin（2*pi*15*t）+sin（2*pi*30*t）;

subplot（411）%绘制数字滤波器幅频响应

plot（w/pi,abs（H））;

gridon;title（'数字滤波器幅频响应'）;

subplot（412）%绘制数字滤波器单位响应

stem（n,h,'.'）;

gridon;title（'数字滤波器单位响应'）

subplot（413）%绘制输入信号

plot（t,x）;

axis（[00.4-44]）

gridon;title（'输入信号'）

subplot（414）%绘制输出信号

y=conv（x,h）;%用卷积求输出信号

N=N1+N2-1;%计算卷积和长度

tt=（0:

N-1）*T;

plot（tt,y）;

axis（[00.8-11]）

gridon;title（'输出信号'）

2.4调试分析及运行结果

代码完成后，滤波器幅频响应、单位响应和输入信号能够正常输出，但是绘制输出信号时用的是plot（y），并未指明自变量，结果绘制的图（如下）是以n为自变量，看起来非常不明了也不清楚。

后将绘图语句直接改为plot（t,y），结果出现错误

提示t和y矩阵长度不匹配，翻阅《信号与系统》上册课本重新复习了离散卷积和部分，得知两信号卷积后长度变化规律，于是将绘制输出信号部分代码改为：

y=conv（x,h）;

N=N1+N2-1;

tt=（0:

N-1）*T;

plot（tt,y）;

然后得到以时间为横坐标的输出信号。

运行结果如下：

2.5心得体会：

通过运用MATLAB设计低通滤波器，学习了MATLAB语言的编程以及熟悉MATLAB软件环境，会利用MATLAB编写基本程序，熟悉绘制图形命令的正确使用，以及掌握MATLAB处理数字信号的有关命令。

在命令窗口运行编写好的m文件，学习M文件调用其他的M文件，文件之间的调用必须满足这两个文件在一个文件夹里。

学习生成常用窗函数的MATLAB语言函数，并通过调用生成各种窗函数MATLAB部函数，以及编写窗函数M文件，在figure窗口观察常用窗函数的图形。

掌握FIR带通滤波器设计原理，并自己查阅有关带通滤波器的设计的书籍，根据设计指标要求利用MATLAB软件设计出符合要求的带通滤波器。

另外，掌握设计FIR滤波器的一般步骤，以及掌握窗函数的选择。

在MATLAB窗口观察用不同窗函数设计的带通滤波器的图形，比较在相同的指标条件下不同的窗函数设计的带通滤波器自之间的差别，以及N取不同值时，同一个窗函数设计的FIR的带通滤波器的技术指标不同。

参考文献

1、王华、有军、建存，《MATLAB电子仿真与应用教程》，国防科技大学，2007

2、闻新、周露、鸿《MATLAB科学图形构建基础与应用》，科技，2002

3、程佩青，数字信号处理及应用，清华大学，2009

4、王宏，MATLAB及其在信号处理中的应用，清华大学，2004