基于matlab的布莱克曼窗函数法设计的低通滤波器

基于matlab的布莱克曼窗函数法设计的低通滤波器
第一章设计要求 1.1基本要求 (1)理解FIR数字低通滤波器的作用及应用领域,掌握布菜克曼窗函数法设计 FIR数字低通滤波器的原理及其在 matlab仿真中的实现方法 (2)掌握 matlab的编程方法; (3)通过脚本编程或 SIMULINK实现FIR数字LPF;使用布莱克曼窗函数;M-1, n=[0:1:M-1;Wc=0.2米pi; (4)完成课程设计报告 1.2提高要求 (1)实现Wc和M均可交的布莱克曼窗 (2)用设计岀的滤波器对声音信号加噪声后进行滤波,对滤波前后信号进行对 比,并在时域和频域进行分别进行分析。 第二章设计原理 2.1FIR滤波器 FIR( Finite Impulse response)滤波器:有限长单位冲激响应滤波器,又称 为非递归型滤波器,是数字信号处理系统中最基本的元件,它可以在保证幅频特 性的冋时具有严格的线性相频特性,同时其单位抽样响应是有限长的,因而滤波 器是稳定的系统。因此,FIR滤波器在通信、图像处理、模式识别等领域都有着 广泛的应用。 有限长单位冲激响应(FIR)滤波器有以下特点 (1)、系统的单位冲激响应h(n)在有限个n值处不为零: (2)、系统函数H(z)在z>0处收敛,极点全部在z=0处(因果系统) (3)、结构上主要是非递归结构,没有输出到输入的反馈,但有些结构中(例 如频率抽样结构)也包含有反馈的递归部分。 FIR低通滤波器的设计方法主要分为两类:第一类是基于逼近理想滤波器特 性的方法,包括窗函数法、频夲采样法和等波纹最佳逼近法:第二类是最优设计 法。实验采用布莱克曼窗函数法设计。 2.2窗函数设计法 窗函数设计法的基木思路是用FIR低通滤波器逼近希望的滤波特性。设希 望逼近的滤波器的频率响应函数为Hd(e"),其单位脉冲响应用hd(n)表示。其 中 Hd(e/w)=y hd(n e/ wn hd()=1∫ Hd(elwewndw WC hd(n)是无限时宽的,是非因果序列。窗函数设计法就是截取hdan)为有限 长的一段因果序列,并用合适的窗函数进行加权做为FIR低通滤波器的单位脉 冲响应h(n)。 窗口法的设计步骤如下: (1)、通过傅里叶变换后的理想滤波器的单位脉冲响应hd(n)。 (2)、根据指标选择窗口截止频率Wc和阶数M。 (3)、给定理想频响由Hd(e")和hd(n),加窗得h(n)=w(n)hd(n) 窗函数又包含有多种,包括矩形窗,三角窗,汉明窗,布莱克曼窗等,实验 采用布莱克曼窗 2.3布莱克曼窗函数 布莱克曼窗的时域形式可表示为: wvn()=(042-05c0s2m+0.08c0sm)R2n(n) 它的频域特性为: Wa(ey)=042(y)-02WR(e0x-)+Wk{(e)+ 0.04[WR(ejw-w-1)+WR(e(wt N-1) (2-4) 其幅度函数为 Wn(o)=042Wn()-0.25(0-2x)+W/(+2)1+ 00(-)+Wm(+) (2-5) 增加一个二次谐波余弦分量,可进一步降低旁瓣,但主瓣宽度进一步增加,为 12π∧N。加N可减少过渡带。布莱克曼窗函数的最大旁瓣之比主鐮值低57db, 但是主瓣宽度是矩形窗函数的主瓣宽度的三倍。布莱克曼窗主瓣宽,旁瓣小,频 率识别精度最低,但幅值识别精度最高 第三章设计步骤及结果分析 实验分基础和提髙两部分,基础部分实现窗函数阶数与截止频率均固定的布 莱克曼窗,提高部分实现窗函数阶数与截止频率均可变的布莱克曼窗并测试通过 不同频率杂音信号,实现滤除高频分量,并且显示其时域及频域图 3.1基础部分 3.1.1设计流程图 输入M,Wc 产生单位康 向 生布菜克曼 加窗了 频率响应 绘圈 结束 图3.1基础部分流程图 由流程图,基础部分通过输入固定M,Wc,调用已有函数进行布菜克曼滤波 器设计,然后绘制波形图进行观察。 3.1.2程序设计 (1) i deal_1p函数 该函数用于产生理想低通脉冲响应 function hd=ideal lp(wc, N) alpha(N-1)/2 n=[0:1:N-1] m=n-alpha hd-sin(wc*m). /(pi*m) en (2)窗函数阶数与截止频率均固定的布莱克曼窗程序 clc: clear all 4 M-11; n=[0:l:M-1]; Wc=0.2*p1; hd= ideal_lp(We,M);%产生单位脉冲响应 Wbla=( blackman(M)’;%产生布莱克曼窗 h-hd.*Wbla;%加窗 [H,W]= frege(h,I);%频率响应 If 1figure: set(Hf 1, NumberTitle', off, Name, M-16, Wc-0 2*pi 1i3k 克曼窗’);%设置标题 subplot(2,2,1);stem(n,hd); title('理想脉冲响应’);axis(0,M1,0.1,0 3]);xlabel(M); ylabel(hd (M)); grid subplot(2,2,2);stem(n,wbla); title('布莱克曼窗’);axis([0,M-1,0,1,1]); xlabel('w): ylabel(wBl(M)); grid subplot(2,2,3);stem(m,h); Litle(实际脉冲响应’);axis([0,M-1,-0.1,0.3]); xlabel(M); ylabel(h(M))); grid subplot(2, 2,1); plot(w/pi, 20*log10(abs(h))); axis([o, 1,-200, 10]); title (幅度响应/dB'); xlabe1(o/π’); y labe(20lg|Ig(o)');grid 3.1.3实验结果及其分析 M=11c=02pi菜克窗 File Edit View Insert Tools Desktop Window Help 固出冷气國函鸟日图日回 NaN+0000000coo000000i Nan NeNi 理想脉冲响应 布菜克曼窗 Nan+ NaN Nan+ NaNi NaN十 NaNI Nan+ NaNI m0.5 Nan Nani Nan+ NaNi NaN+ NaN↓ NaNI NaN十 NeNi 实际尿冲响应 幅度响应d NaN十 Neni 33 Nan+ NaNi Nan+ NaN Nan+ NaNi Nan NeNi Nan+ NaN Nan+ NaN NaN+ NaN WTT NaN十 NENI 图3.2M=11时的窗函数设计 由图,布莱克曼窗已经显示出来,但由于阶数太小,幅度响应并不能显示(见 图32左),因为经过 fretz函数后,幅度是非数(见图32右)。经测试,当M>=16 时,幅度响应才能够正确显示。可以看出,效果并不是很好(见图3.3)。 M-16W-02*pi莱克曼窗 File Edit View Insert Tools Desktop Window Hell 理想脉 布菜克曼窗 C.3 .2 .1 0.5 5 15 15 实际脉冲响应 幅度响应dB .1 15 10 15 图3.3M=16时的窗函数设计 3.2提高部分 3.2.1设计流程图 本课程设计是对录制的语音信号进行加噪处理并分析加噪前后语音信号的 时域图与频域图,再用布莱克曼窗设计一个FIR滤波器,而滤波器的设计必须 符合其相应的指标,否则将不能滤掉加入的噪声。最后将滤波前后的波形图进行 比较看是否相同。下面是整个设计的流程图如图3.4所示: 开始 原始信号时域图←读取语音信号}原始信号频谱 加噪信号时城图←指音信号加噪}“加信号谱 否 是否输 入MWc 是 默以。个用读器 匚滤波信号时域图←「滤波器滤除噪声滤波信号烦谱 绘图 结束 图3.4语音信号的整个流程图 由图知,系统首先通过手机录制声音信号,然后用 matlab软件读取,再对 语音信号加噪声,经过设计好的滤波器进行滤波,最后比铰通过滤波器前后的声 音信号的时域及频域的区别。 由于整个系统过于复杂,因此分为以下几个子系统分别进行实现 3.2.2读取语音信号 打开 MATLAB软件平台,将用手机录制好的语音信号(为了区别低频高频, 这里通过讲话声和口哨声区分)利用函数 audiohead对语音信号进行采样,记住 采样频率的和釆样点数。再凋用函数 sound函数听录制的语音。采样完后再语音 信号中加入一个单频噪声,单频的噪声的频率可以自己设置。按照加入噪声后的 采样频率再次调用 sound函数,这时可以明显的听见播放的语音信号中有尖锐的 噪声。 fs=44100 x1- audiohead(C: USers、Al1 bymyself\ Desktop\Fir.aac;%输入参数为文件的 全路径和文件名,输出参数是每个样木的值 t-0:1/44100:(1 ength(x1)-1)/44100;%计算时间范围,样本数除以采样频率 yl=fft(x1,441000); f=fs米(0:220499)/441000 % sound(x1,44100);%原始声音信号 Au=0.01; d=LAU*COsS(2*p*50000kt)”;%单频噪声,频率参数可改 x2=x1+d; % sound(x2,44100);%明显听出有尖锐的单频啸叫声 y2=fft(x2,111000 3.2.3低通滤波器设计 (1) black man函数设计 该函数用于产生布莱克曼窗,需要注意的是这个函数与自带函数b1 ackman产生 列向量不同,它产生的是行向量 function[x]= black man (N) n-0:N-1; x=(0.42-0.5米COS((2来pi*n)/(N-1)+0.08*COs(4**n)/(N-1)) end 7 (2) blackman2 kebian数设计 该函数为频率和阶数均可变的布莱克曼窗函数此函数设计了缺省参数,当不 输入参数,或只输入一个参数时,函数会设置默认参数,并且把Ⅵ和Wc显示在 左上角,方便观察 function [ blackman2kebian(, Wc) %设置缺省参数,当输入一个参数,默认Wc=0.2*pi,当输入零个参数,默认M-=11 Wc=0. 2*pi if nargin--1 Wc=0. 2*pi; end if nargin==0 M=ll; Wc=0. 2*pi; end %设置标题 a b=num2str(m) c=strcat(a, b) d=,wC b=num2str(Wc) d=strcat(a, b) e=strcat(c, d) e-strcat( blackman', e) e= sTreat(e,”)’) Hf 1=figure; set(Hf 1, NumberTitle', off, Name', e n-[0:1:M-1; hd= ideal_lp(We,M;%理想低通滤波器脉冲响应 wbla=( black man(M));%调用自编 black man窗 hn=hd.水wbla; H,W]= fretz(hn,1);%频率参数 绘图 subplot(2,2,1);stem(n,hd); title('理想脉冲响应’);axis(0,M1,0.1,0. 3); xlabel(M); ylabel(hd(M))); grid subplot(2,2,2);stem(m,wbla); lille(布莱克曼窗’);axis([0,M-1,0,1.1]) xlabel(M); ylabel(’wB1(M’);grid subplot(2,2,3);stem(m,hn); title('实际脉冲响应');axis(0,M-1,-0.1,0