基于窗函数法的FIR数字高通滤波器设计doc.docx
《基于窗函数法的FIR数字高通滤波器设计doc.docx》由会员分享,可在线阅读,更多相关《基于窗函数法的FIR数字高通滤波器设计doc.docx(19页珍藏版)》请在冰豆网上搜索。
![基于窗函数法的FIR数字高通滤波器设计doc.docx](https://file1.bdocx.com/fileroot1/2023-2/4/aa0885cc-3313-4212-a2a5-95ce9c9fe196/aa0885cc-3313-4212-a2a5-95ce9c9fe1961.gif)
基于窗函数法的FIR数字高通滤波器设计doc
课程设计说明书
题目:
基于窗函数法的FIR数字高通滤波器设计
姓名:
院(系):
专业班级:
学号:
指导教师:
成绩:
时间:
课程设计任务书
题目基于窗函数法的FIR数字高通滤波器设计
主要内容、基本要求、主要参考资料等:
主要内容:
利用MATLAB软件读取一段含有噪声的.wav格式的语音信号,然后基于FFT对该信号进行频谱分析;基于含噪语音信号的频谱确定滤波器的参数,利用窗函数法设计一个FIR数字高通滤波器,并利用所设计的滤波器对信号进行滤波处理。
比较滤波前后语音信号的时域波形及频谱,分析滤波前后的语音变化。
基本要求:
1、基于含噪语音信号的频谱确定滤波器的参数;
2、分别采用矩形窗、汉明窗和布莱克曼窗设计FIR数字高通滤波器;
3、掌握利用wavread函数读取、播放.wav格式语音信号的方法;
4、对语音信号进行滤波,绘制滤波前后信号的时域波形及频谱;
5、回放语音信号,分析滤波前后的语音变化。
主要参考资料:
1、从玉良.数字信号处理原理及其MATLAB实现[M].北京:
电子工业出版社.2009.7
2、胡广书.数字信号处理理论、算法与实现[M].北京:
清华大学出版社.2003,8
完成期限:
一
指导教师签名:
课程负责人签名:
摘要1
1数字滤波器简介2
2FIR数字滤波器设计2
2.1FIR数字滤波器的原理2
2.2设计工具2
2.3用窗函数法设计FIR数字滤波器3
2.3.1常用窗函数3
2.4FIR数字滤波器的一般设计步骤:
5
3MATLAB仿真设计6
4主要程序7
5仿真结果图11
6总结13
参考资料14
基于窗函数法的FIR数字高通滤波器设计
摘要:
数字滤波器在图像处理、数据传输等场合具有广泛应用,其设计是信号处理的核心问题之一,软件实现数字滤波优势体现在滤波器参数的改变伴随滤波器性能的改变.
阐述了数字滤波器的设计方法,讨论了线性相位的条件和幅度特性,并以窗函数在MATLAB软件中实现了滤波器的仿真设计。
仿真表明,借助计算机软件设计数字滤波器的过程简便易行,滤波性能指标均能达到指定的要求!
对数字滤波器的推广应用有一定的促进作用。
关键词:
数字滤波器;MATLAB;窗函数;高通滤波
1数字滤波器简介
数字信号处理(DigitalSignalProcessing,简称DSP)是一门涉及许多学科而又广泛应用于许多领域的新兴学科。
20世纪60年代以来,随着计算机和信息技术的飞速发展,数字信号处理技术应运而生并得到迅速的发展。
在过去的二十多年时间里,数字信号处理已经在通信等领域得到极为广泛的应用。
数字滤波器的种类很多,从功能上我们可以将其分为低通、高通、带通和带阻等滤波器(输入输出都为数字的滤波器称之为数字滤波器)数字滤波器是通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。
(数字滤波器从单位脉冲响应分类,可以分为无限脉冲响应(HR)滤波器和有限脉冲响应(FIR)滤波器(由数字信号处理的一般理论可知,IIR滤波器的特征是具有无限持续时间的冲激响应,而FIR滤波器的冲激响应只能持续一定的时间。
2FIR数字滤波器设计
2.1FIR数字滤波器的原理
FIR滤波器的传递函数为:
H(Z)=KQ=£;i(")zF(2T)
X(Z)企
可得FIR滤波器的系统系数差分方程,当FIR滤波器的系数满足下列中心对称条件时:
b(n)=b(N-1-n)或b(n)=-b(N一1一n)(2一2)
线性相位滤波器的相位滞后和延迟在整个频带上是相等且不变的。
对于一个N阶的线性相位FIR滤波器,群延迟为常数,即滤波后的信号简单地延迟常数个时间步长。
这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。
2.2设计工具
在MATLAB中有三种方法可以实现设计FIR数字滤波器:
程序设计法、FDATool设计法和SPTool设计法,程序设计法具有设计简便、可移植性强等特点,故在本设计中采用的是程序设计法。
而在程序设计方法中,又有很多方法可以实现FIR滤波器的设计,如窗函数设计法、频率采样设计法和最优化设计法等,这些方法都比较快捷方便,其中最常用的是窗函数设计法。
2.3用窗函数法设计FIR数字滤波器
设所希望设计的滤波器传输函数为Hid"),h(n)是与其对应的单位脉冲响应,因此
.7T
hd(n)=v-1HX^e^dw
由已知的求出h(n),经过Z变换可得到滤波器的系统函数。
但一般情况下,
H(elw)是逐段稳定的,在边界频率处有不连续点,因而h(n)是无限时宽的,且是非因
果序列(但是从实现的角度来说,我们希望得到一个长度为Z的线性相位滤波器,因此只能通过对h(n)进行加窗得出。
窗函数设计滤波器的基本思想,就是根据给定的滤波器技术指标,选择滤波器的阶数N和合适的窗函数w(n)。
即用一个有限长度的窗口函数序列w(n)来截取一个无限长的序列h(n)获得一个有限长序列h(n),即H(n)=h(n)*w(n),这样我们用一个有限长的序列H(n)去代替h(n),肯定会引起误差,表现在频域就是通常所说的吉布斯效应该效应引起通带内和阻带内的波动性,尤其使阻带的衰减小,影响滤波器的性能。
为了减小吉布斯效应,从原理上来说可以通过加大N,但实验表明,加大N虽然可以使H(w)过渡带变窄,但对带内波动以及阻带衰减并没有多大的影响。
因此,我们希望能找到不同的窗函数,使其主瓣包含更多的能量,同时旁瓣幅度减小,从而使通带%阻带波动减小,加大阻带衰减。
工程实际中常用的窗函数有6种,即矩形窗、三角形窗、汉宁窗、哈明窗、布莱克曼窗和凯泽窗,不同窗函数的过渡带和阻带衰减特性有所不同。
2.3.1常用窗函数
1汉宁(Hanning)窗
汉宁窗函数的时域形式表示为:
(3T)w(n)=0.5-0.5cos(^^),/?
=0,1,2,•••,2V-1N-1
频域形式为
(3-2)
W(J®)«^0.5Ws((y)+0.25Ws(®-—)+T¥s(®+—)je2
汉宁窗函数的最大旁瓣值比主瓣值低31dB,但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍,为8ji/No
hanning函数:
生成汉宁窗
调用方式:
(1)w=hanning(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
注意:
此函数不返回是零点的窗函数的首尾两个元素。
(2)w=hanning(n,'symmetric'):
与上面相类似。
(3)w=hanning(n,'periodic'):
此函数返回包括为零点的窗函数的首尾两个元素。
图2.3.1汉宁窗及其频谱特性
2汉明(Hamming)窗:
函数的时域形式可以表示为
w伙)=0.54-0.46cos^27t^k=\,2,…,N(3-3)
频域形式为:
=0.54IVr(o)+0.23WR(a>--]+%仏+二叮(3-4)
其中,叫⑷为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB,但它和汉宁窗函数的主瓣宽度是一样大的。
Hamming函数:
生成海明窗
调用方式
(1)w=hamming(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2)w=hamming(n,sflag):
参数sflag用来控制窗函数首尾的两个元素值;其取值
为symmetric或periodic;默认值为symmetrico
3布莱克曼窗函数的时域形式可以表示为
w(灯=0.42-0.5cos]2兀一)+0.08cos[4兀一jk=\,2,…,N(3-5)
它的频域特性为
=0.42WR(a))+0.25WR
27T
CD
N—1
+啲+語
+0.04叭
4兀
G)
N
(3-6)
其中,肥⑷为矩形窗函数的幅度频率特性函数。
布莱克曼窗函数的最大旁瓣值比主瓣值低57dB,但是主瓣宽度是矩形窗函数的主瓣宽度的3倍,为12jr/No
Blackman函数:
生成海明窗
调用方式
(1)w=blackman(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2)w=blackman(n,sflag):
参数sflag用来控制窗函数首尾的两个元素值;其取值为symmetric或periodic;默认值为symmetric□
图2.3.2布莱克曼窗及其频谱特性
2.4FIR数字滤波器的一般设计步骤:
(1)按照任务要求确定滤波器的指标;
(2)用因果稳定的离散线性时不变系统的系统函数逼近这一性能指标要求;
(3)利用有限精度算法实现该系统函数,包括选择运算结构、适当的字长和有效
位处理的方法;
(4)进行实际技术实现,常采用软件、硬件或者软硬件结合的方法实现。
3MATLAB仿真设计
本文的设计利用MATLAB的FDAtool采用窗函数法来设计语音高通数字滤波器。
在MATLABcommand窗口下,执行fdatool命令,将会出现如下窗口
MfilterDesign&AnalysisTool-[untitled.fda]
FileEditAnalysisTargetsViewWindowHelp
□0$Q3Qi炉0QXOE4-Ej-^^IOSS粉
—FrterSpecifications.
ResponseType
_FilterOrder
Q|Lowpass
CSpecifyorder:
pO~
r|Highpass
CBandpass
CBandstop
(•hArwnumorder
lOptions
广|DifferertiatorZalgnMother!
DensityFactor:
[20
|Butterworth
QFIR|Equiripple
FrequencySpecifications
_MagnitudeSpecifications
5it$:
|H2
Units:
|dB
Fs:
(48000
Apass|1
Fpass
|9600
Astop
Fstop
[12000-
步骤如下:
选择一段时间较短的单声道音频,用计算机windows附件中的录音机将其他格式的语音转换成。
Wav格式并命名为gaot保存到MATLAB的work文件夹中,也可以自行创建work文件夹。
用MATLAB中的wavread()函数将.wav语音转换成数组格式数据,用wavwrite()函数将数组格式数据转换成.wav格式语音,利用MATLAB的FDAtool工具设计所需FIR滤波器系数h(n)。
由于本文实验所用音频录制时的抽样频率为44.lkH,所以在设计h(n)系数E1寸应采用44.lkH,其他参数设置为,Response选择Highpass,Design选择FIR的window通过频率Fp=44.5kH,截止频率Fs=20kH.将设计好的高通滤波器系数导出到工作空间,命名为higho
4主要程序
[xl,fsl,bitsl]=wavread(,F:
\MATLAB课程设计\work\SHE.wav,);
%soundsc(xl,fsl);
[x2,fs2,bits2]=wavread(,F:
\MATLAB课程设计\work\Highpass.wav');
%soundsc(x2,fs2);
figure
(1)
subplot(2,1,1);
plot(xl)%做原始语音信号以44.Ik采样后的时域图形
titleC原始语音采样后时域信号');
xlabel('时间轴n,);
ylabel(,幅值A');
subplot(2,1,2);
plot(x2)%做加噪声语音信号以44.Ik采样后的时域图形
title('加噪声语音采样后时域信号');
xlabel('时间轴n,);
ylabel(,幅值A');
pause;
Xl=fft(xl,4096);
figure
(2)
subplot(2,1,1);
f=fsU(0:
2047)/4096;
plot(f,abs(XI(1:
2048)));
titlef原始信号频谱');
X2=fft(x2,4096);
subplot(2,1,2);
plot(f,abs(X2(1:
2048)));
titleC加噪声信号频谱');pause;
%fs3=44100;
%fpl=250;fp2=7500;
%fsl=180;fs2=7570;
%As=60;
%Wsl=(fpl+fsl)/fs3;Ws2=(fp2+fs2)/fs3;
%w=(fpl-fsl)/fs3;
%M=ceil((As-7.95)/(14.36*w));
%window=boxcar(M+l);
%%window=hamming(M+l);
%%window=blackman(M+l);
%[b,a]=firl(M,[Wsl,Ws2],window):
%figure(3)
%plot(b);
%设定采样频率
%第一截止频率
%第二截止频率
%最小阻带衰减
%截止频率归一化处理
%求归一化过渡带
%计算所需滤波器的阶数
%生成长度为M+1的矩形窗
%生成长度为M+1的汉宁窗
%生成长度为M+1的布莱克曼窗
%生成设计的fir滤波器
%pause:
%figure(4)
%freqz(b);%绘制幅频和相频响应曲线
%wpl=0.007*pi;wsl=0.017*pi;
%tr_width=wsl-wpl;
%M=ceil(6.6*pi/trwidth)+1;
%wcl=(wsl+wpl)/2;%求截止频率
%求过渡带宽度
%求得所需窗函数的长度
%hd=ideal_hpl(wcl,M);
%w_ham=(hamming(M))';
%h=hd.*wham:
%求得理想带阻的冲击响应
%得到长度为M的汉宁窗
%利用窗函数截短
%figure(3)%freqz(hd);
%figure(4)
%freqz(h);
%理想高通滤波器的幅频响应频谱
%加窗之后的幅频响应频谱
%pause;
wpl=0.007*pi;wsl=0.017*pi;
tr_width二wsl—wpl;
M=ceil(6.6*pi/tr_width)+1;wcl=(wsl+wpl)/2;hd=ideal_hpl(wcl,M);w_ham=(blackman(M));
h=hd.*w_ham;
%求过渡带宽度
%求得所需窗函数的长度
%求截止频率
%求得理想带阻的冲击响应
%得到长度为M的布莱克曼窗
%利用窗函数截短
figure(3)freqz(hd);
figure(4)freqz(h);
%理想高通滤波器的幅频响应频谱
%加窗之后的幅频响应频谱
pause;
f2=filter(h,1,x2);
FO=fft(f2,4096);
figure(5)
subplot(3,1,1);
plot(f,abs(XI(1:
2048)));
titlef原始信号频谱');subplot(3,1,2);
plot(f,abs(X2(1:
2048)));title('滤波前的频谱')xlabel('Hz,);
ylabel('fuzhi');
subplot(3,1,3)
plot(f,abs(F0(1:
2048)));titlef滤波后的频谱')xlabel('Hz,);
ylabel('fuzhi');
pause;
figure(6)subplot(2,1,1);
plot(xl)%做原始语音信号以44.Ik采样后的时域图形
titleC原始语音采样后时域信号');
xlabel('时间轴n,);
ylabel(,幅值A');
subplot(2,1,2);
plot(f2)%做滤波之后语音信号以44.Ik采样后的时域
图形
titleC滤波之后语音采样后时域信号');
xlabel('时间轴n,);
ylabel(,幅值A');
pause;
sound(xl,44100);
sound(f2,44100);%播放滤波后的语音信号
functionhd二ideal_hpl(Wc,N)
%computetheidealhighpassfiterunitpulserespondencehd(n)
%wc:
cutofffrequency
%N:
windowlength
%hd:
unitpulserespondence
alpha=(N-l)/2;
n=0:
1:
N-1;
m=n-alpha+eps;
hd=[sin(pi*m)-sin(Wc*ni)]・/(pi*ni);
5仿真结果图
原始信号频谱
0.511.5
2
8060
40
20
加噪声信号频谱
2.5
4
x10
60
40
20
0.51
1.5
2.5
4
x10
20
0.20.30.40.50.60.70.80.91
NormalizedFrequency(xjirad/sample)
oo
-2bp)apn-Ebew
40
(sasBap)as_2d
o-6
0.20.30.40.50.60.70.80.91
NormalizedFrequency(xtirad/sample)
原始信号频谱
50
o
0.5
1.5
滤波前的频谱
2.5
1034
00o
5C
0.5
11.5
Hz
滤波后的频谱
2.5
4
x10
50
一llza
0.5
1
1.5
2
2.5
x104
o
Hz
原始语音采样后时域信号
505
o.0<1
0.51
1.52
时间轴n
2.5
50.5
o.0<1
0.511.522.533.5
时间轴n%
滤波之后语音采样后时域信号
6总结
一周的课程设计,感觉收获很多,刚开始接到设计任务书,不懂得从何处下手,同时发现自己以前学的知识基本忘了差不多,比如FIR滤波器的具体概念,还有高通的意田竺竺
JHl'-vj*O
面对这么多问题,我首先认真回忆了这些知识在哪学过,然后找到相应的书籍认真的复习的一遍,当然还不够理解这些概念,然后在网上查一些资料,借鉴了FIR滤波器设计的基本过程,加深了自己的理解。
通过这个过程,我深刻认识了读书心存志远,实践悟出真知。
只有自己去思考了,动手了才能对这个知识的深入理解。
当然,在此次课程设计还有很多知识是我新接触的,这个过程我成长了很多。
这些知识不仅在课堂上有效,对以后的学习也有很大的指导意义,在日常生活中更是有着现实意义;也对自己的动手能力是个很大的锻炼。
实践出真知,纵观古今,所有发明创造无一不是在实践中得到检验的。
没有足够的动手能力,就奢谈在未来的科研尤其是实验研究中有所成就。
总的来说,对这次的课程设计很满意,因为我真正的学到了很多,对自己也认识了其中的不足,知道接下去该如何补充自己了,永远不要停止学习的脚步。
参考资料
1、从玉良.数字信号处理原理及其MATLAB实现[M].北京:
电子工业出版社.2009.7
2、胡广书.数字信号处理理论、算法与实现[M].北京:
清华大学出版社.2003,8
33.5
x105