FIR数字滤波器的(海明)窗函数法设计

课程设计说明书NO.1FIR数字滤波器的（海明）窗函数法设计1．课程设计目的（1）熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数。（2）加深对FIR数字滤波器设计的理解，并用窗函数法进行FIR数字滤波器的设计。（3）将设计出来的FIR数字滤波器利用MATLAB进行仿真。（4）对一段音频文件进行加入噪声处理，对带有噪声的文件进行滤波处理。2.设计方案论证2.1Matlab语言概述MATLAB是一种以矩阵运算为基础的交互式程序语言，专门针对科学、工程计算及绘图的需求。随着版本的不断升级，内容不断扩充，功能更加强大，从而被广泛应用于仿真技术、自动控制和数字信号处理领域。此高级语言可用于技术计算此开发环境可对代码、文件和数据进行管理交互式工具可以按迭代的方式探查、设计及求解问题数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数积分等二维和三维图形函数可用于可视化数据各种工具可用于构建自定义的图形用户界面各种函数可将基于MATLAB的算法与外部应用程序和语言（如C、C++、Fortran、Java、COM以及MicrosoftExcel）集成不支持大写输入，内核仅仅支持小写2.2声音处理语音是人类获取信息的重要来源和利用信息的重要手段。语音信号处理是一门发展十分迅速，应用非常广泛的前沿交叉学科，同时又是一门跨学科的综合性应用研究领域和新兴技术。声音是一种模拟信号，而计算机只能处理数字信息0和1。因此，首先要把模拟的声音信号变成计算机能够识别和处理的数字信号，这个过程称为数字化，也叫“模数转换”。在计算机对数字化后的声音信号处理完课程设计说明书NO.2后，得到的依然是数字信号。必须把数字声音信号转变成模拟声音信号，然后再输出到扬声器，这个过程称为“数模转换”。2.3数字滤波器的介绍数字滤波器(digitalfilter)是由数字乘法器、加法器和延时单元组成的一种装置。其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。数字滤波器是一个离散时间系统（按预定的算法，将输入离散时间信号转换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1／2抽样频率点呈镜像对称。为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。2.4FIR滤波器基本结构FIR滤波器的数学表达式为：式中：N为FIR滤波器的抽头数；x(n)为第n时刻的输入样本；h(i)为FIR滤波器第i级抽头系数。其相应的z变换为：式中：z-i为N-1阶多项式。在自适应处理、数据通信等领域中往往要求信号在传输过程中不能有明显的相位失真，FIR滤波器可以做到线性相位满足此要求。F1R滤波器实质上是一个分节的延迟线，把每一节的输出加权累加，得到滤波器的输出。对于FIR滤波器的单位脉冲响应h(i)只要满足以下2个条件之一，则为线性相位滤波器。线性相位的FIR滤波器具有中心对称的特性，其对称中心在N／2处。课程设计说明书NO.3(2)由性能指标确定窗函数w(n)和窗口长度N(3)求得实际滤波器的单位脉冲响应h(n)(4)检验滤波器性能。设计常用的窗函数有矩形窗、汉宁窗、海明窗、凯撒窗等。其中：海明窗的旁瓣峰值小于主瓣峰值的1％，99．963％的能量集中在主瓣内．且通过海明窗设计的FIR滤波器在较少的阶数下可以得到较小通带纹波，非常适合工程设计。2.5利用Windows进行语音信号的采集利用windows下的录音机，按照（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8000KHz，8位，单声道）的顺序操作，录制一段自己的语音，录制时间为5秒，如图1，图2所示，将自己录好的语音文件保存为“zf.wav”。图1选择windows下的录音机课程设计说明书NO.4图2用8000Hz采样录音2.6语音信号的分析(1)将上一步骤中保存下来的语音信号文件“zf.wav”复制到计算机装有Matlab软件的磁盘中相应Matlab目录中的“work”文件夹中:(C:\ProgramFiles\MATLAB71\work)。(2)双击桌面上Matlab软件的快捷图标，打开Matlab软件。(3)在Matlab菜单栏中选择“FilenewM-File”或是点击快捷按钮，打开m文件编辑器。(4)在m文件编辑器中输入相应的指令将自己的语音信号导入Matlab工作台。(5)画出原始语音信号s的波形，由于原始语音信号开始一段会是无用的语音信号，因此要截取掉，截取的一段语音信号为1至1+fs-1，即从1到7999，画出截取原始语音信号s1的波形，代码如下，波形如图3，所示：closeallclearallclc[s,fs,bits]=wavread('C:\ProgramFiles\MATLAB71\work\zf.wav');s1=s(1:8000);sound(s1,fs,bits);figure(1);课程设计说明书NO.5subplot(211)plot(s)title('原始语音信号')subplot(212)plot(s1)title('截短语音信号');会得到下图：图3原始语音信号和截短语音信号(6)对语音信号进行频谱分析，在Matlab中，利用函数FFT()对信号进行快速傅里叶变换，得到信号的频谱特性，如图4所示wavwrite(s1,fs,'s1.wav');%将被处理信号s1输出为语音文件“s1.wav”S1=fft(s1);figure(2)subplot(311);plot(s1);课程设计说明书NO.6title('截短预处理语音信号')subplot(312)plot(abs(S1))title('预处理语音信号频谱');subplot(313);k=0:4000;plot(k(1:4000)*1,abs(S1(1:4000)));title('预处理语音信号单边带频谱')图4截短预处理语音信号如图4所示，从右向左看，第一个较大的波峰所在的频率即为3db截止频率，第二个波峰所对应的频率为通带截止频率pf，在图中可以读出pf=610Hz；一般在3db截止频率右侧的波谷位置选择阻带截止频率stf=750Hz。通带截止频率课程设计说明书NO.7pf、阻带截止频率stf数值的确定，就可以确定滤波器的基本指标。图4里第二个图是信号)(nx的FFT结果，即)]([)(nxDFTkX是信号)(nx的实际频谱)]([)(nxDFTejXw采样，本设计中信号)(nx的长度取的是L=8000点，所以，图中，)(kX的每两个相邻点之间的频率间隔大小，即频率分辨率：180008000信号的长度采样频率f（Hz），所以，根据它的放大图，即图4中的第三个图的放大图，可以确定HzffHzffHL760760,240240点点。2.7滤波器的设计2.7.1滤波器的参数设定本次课设我设计的是一个线性FIR低通滤波器，利用的窗函数是hamming窗，如上（图4）所述：所以：HzfHzfHzfsstp8000,750,610通带截止频率为4791.0800061022sppffw(rad/sample)阻带截止频率为589.0800075022sststffw(rad/sample)1099.04791.0589.0pst(rad/sample)由于海明窗过渡带满足：Nw23.3求得滤波器阶数stpstpcff22212117.0)(sstpsccffffw189N9421N课程设计说明书NO.8（1）给定所要求的频率响应函数)(jdeH.||0||)(ccjjdeeH（2）求单位采样响应)(nhd)]([)(jwddeHIDFTnhdweeHjwnjw)(21)]([nwSawcc（3）海明窗)()]12cos(46.054.0[)(nRNnnwN（4）滤波器的单位采样响应：)()()(nwnhnhd)()]12cos(46.054.0)][([)(nRNnnwSawnhNcc2.7.2滤波器的MATLAB仿真在M文件中继续编写代码，把计算出来的参数带入代码中。代码如下：%%%%%%%%%加噪声完成信号截取[s,fs,bits]=wavread('C:\ProgramFiles\MATLAB71\work\zf.wav');s1=s(1:8000);sound(s1,fs,bits);figure(1);subplot(211)plot(s)title('原始语音信号')subplot(212)plot(s1)title('截短语音信号');wavwrite(s1,fs,'s1.wav');S1=fft(s1);figure(2)subplot(311);)()]94cos(46.054.0)][94([17.0189nRnnwSac课程设计说明书NO.9plot(s1);title('截短预处理语音信号')subplot(312)plot(abs(S1))title('预处理语音信号频谱');subplot(313);k=0:4000;plot(k(1:4000)*1,abs(S1(1:4000)));title('预处理语音信号单边带频谱')s2=awgn(s1,15);%%%%%%%%%%%%%%%%%%%%%%%%%%完成加噪！15dbwavwrite(s2,fs,'s2.wav');figure(3);subplot(211);plot(s2);title('加噪后语音信号');subplot(212);S2=fft(s2);plot(abs(S1));title('加噪后信号频谱');figure(4)subplot(211);plot(s1);title('语音信号');subplot(212);plot(s2);title('加噪后语音信号');%%%%%%%%%滤波器完成相关参数配置wp=610*2*pi/8000;wst=750*2*pi/8000;wc=(wp+wst)/2;N=ceil(3.3*2*pi/(wst-wp))+1;r=(N-1)/2;hn1=fir1(N-1,wc/pi,'low',hamming(N));%s3=conv(s2,hn1);wavwrite(s3,fs,'s3.wav');S3=fft(s3);figure(5)freqz(hn1);title('滤波器幅频特性与相频特性')figure(6)subplot('111')plot(hn1);title('滤波器系统函数');课程设计说明书NO.10figure(7)subplot(211)plot(s3)title('滤波器处理之后信号图')subplot(212);plot(abs(S3));title('滤波器处理之后频谱');figure(8)subplot(211)plot(s2);title('加噪后语音信号');subplot(212);plot(s3);title('滤波器处理之后信号图');%%%%%%%%%%%%%%求信噪比snr在Workspace中体现！s4=conv(s1,hn1);%p1=sum(s1.^2);%p2=sum(s2.^2)-sum(s1.^2);%SNR1=10*log10(p1/p2);p1=sum(s4.^2)/8000;p2=sum(s3.^2)/8000-sum(s4.^2)/8000;SNR2=10*log10(p1/p2);运行代码效果如图5所示。的噪声生成叠加指令为：awgn，所加的噪声为15