布莱克曼BlackmanFIR滤波器课程设计要点.docx

1、布莱克曼BlackmanFIR滤波器课程设计要点成绩 布莱克曼（Blackman）FIR滤波器课程设计目录一、设计任务与要求 2二、总体方案设计 2三、单元电路设计与参数计算 33.1 FIR滤波器 33.2 窗口设计法 33.3 布莱克曼窗 4四、性能测试与分析 54.1语音信号录制 54.2 滤波器设计 74.3信号滤波处理 84.4结果分析 9五、结论与心得 10六、参考文献 10七、附录 117.1 布莱克曼窗设计FIR滤波器的MATLAB主程序 117.2 调用函数ideal_lp.m的MATLAB源代码 137.3 调用函数freqz_m.m的MATLAB源代码 13布莱克曼窗FI

2、R滤波器的设计一、设计任务与要求1滤波器指标必须符合工程实际。； 2设计完后应检查其频率响应曲线是否满足指标。；3. 独立完成课程设计并按要求编写课程设计报告书。二、总体方案设计本课程设计是对录制的语音信号进行加噪处理并分析加噪前后语音信号的时域图与频域图，再用布莱克曼窗设计一个FIR滤波器，而滤波器的设计必须符合其相应的指标，否则将不能滤掉加入的噪声。最后将滤波前后的波形图进行比较看是否相同。课程设计流程图如图1所示。图1：课程设计流程图三、单元电路设计与参数计算3.1 FIR滤波器FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中

3、最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。FIR数字滤波器以其良好的线性特性被广泛应用于现代电子通信系统中，是数字信号处理的重要内容之一。在实际信号处理中，往往要求系统兼具实时性和灵活性，而已有的一些软件或硬件实现方案(如DSP)则难以同时达到这两方面的要求。使用具有并行处理特性的FPGA来实现FIR滤波器，既有很强的实时性，又兼顾了灵活性，为数字信号处理提供了一种很好的解决方案。设FIR滤波器的单位冲激响应h(n)为一个N点序列，0nN1，则滤波器的系统函数为： (2-1)就是说，它有阶极点在处，有个零点位于有限

4、z平面的任何位置因此是永远稳定的。稳定和相位特性是FIR滤波器突出的优点。FIR滤波器有以下几种基本结构：横截型（卷积型、直接型）、级联型、频率抽样型、快速卷积结构。 3.2 窗口设计法 窗函数设计法的基本思路是用FIRDF逼近希望的滤波特性。设希望逼近的滤波器的频率响应函数为，其单位脉冲响应用表示。为了设计简单方便，通常选择为具有片段常数特性的理想滤波器。因此是无限长非因果序列，不能直接作为FIRDF的单位脉冲响应。窗函数设计法就是截取为有限长的一段因果序列，并用合适的窗函数进行加权做为FIRDF的单位脉冲响应。下面介绍窗函数设计法的基本设计过程。窗口设计法的主要工作是计算和，但当较为复杂时

5、，就不容易由反傅里叶变换求得。这时一般可用离散傅里叶变换代替连续傅里叶变换，求得近似值。窗口法的设计步骤如下：（1）通过傅里叶变换换的理想滤波器的单位脉冲响应。 （2）根据指标选择窗口形状、大小和位置。确定窗口类型的主要依据是过渡带宽和阻带最小衰耗的指标。 （3）给定理想频响由和，加窗得。 （4）检验滤波器性能。由求是否在误差容限之内。窗口函数对理想特性的影响：改变了理想频响的边沿特性，形成过渡带，宽为,等于的主瓣宽度；过渡带两旁产生肩峰和余振（带内、带外起伏），取决的旁瓣，旁瓣多，余振多；旁瓣相对值大，肩峰强，与N无关；N增加,过渡带宽减小,肩峰值不变。因主瓣附近(2-2)其中,所以N的改变

6、不能改变主瓣与旁瓣的比例关系，只能改变的绝对值大小和起伏的密度，当N增加时，幅值变大，频率轴变密，而最大肩峰永远为8.95%，这种现象称为吉布斯（Gibbs）效应。3.3 布莱克曼窗布莱克曼窗的时域形式可表示为： (2-3)它的频域特性为： (2-4)其中为矩形窗函数的幅度频率特性。增加一个二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，为。加N可减少过渡带。布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。四、性能测试与分析4.1语音信号录制 通过手机录制一段长度24秒的语音文件，如图2所示。图2：语音信号录制录制好语音信号后，打开MATLAB软件平台，利用函数aud

7、ioread对语音信号进行采样，再调用函数sound此时可以听见录制的语音。采样完后再语音信号中加入一个单频噪声，单频的噪声的频率可以自己设置。按照加入噪声后的采样频率调用sound函数，这时可以明显的听见播放的语音信号中有尖锐的单频啸叫声。下面是调用该语言信号以及加入噪声的程序：x,fs=audioread(test_yfw.mp3);% 使用audioread读取音频文件的采样率fssound(x,fs); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围

8、，样本数除以采样频率x=x(:,1);y=x+0.1*sin(fn*2*pi*t);sound (y,fs); %明显听出有尖锐的单频啸叫声现在是对加入噪声前后的语音信号进行频谱分析，先对原始和加噪后的语音信号进行傅里叶变换，再计算频谱的频率范围和谱线间隔。最后就可以画出未加入噪声和加入噪声后的时域图和频域图。所有未加和加入噪声的时域图和频域图如图3。下面是对未加和加入噪声的频谱分析的程序：X=abs(fft(x);Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换X=X(1:N/2);Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=

9、0:deltaf:fs/2-deltaf; % 计算频谱频率范围 图3：语音信号的时域图和频率图4.2 滤波器设计本课程设计就是要设计一个滤波器虑掉加入的噪声，使其恢复原始的语音信号。而设计滤波器的方法有很多，例如：窗函数法、频率采样法、脉冲响应不变法和双线性变换法。而本课程设计采用的是窗函数法设计FIR滤波器。而FIR滤波器的设计也有很多方法。在Matlab中，可以利用矩形窗、三角窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗等设计FIR滤波器。而本次采用的是布莱克曼窗来设计滤波器。在用布莱克曼窗设计滤波器的时候，首先要确定滤波器的性能指标。从六种窗函数的基本参数中我们可以得到旁瓣峰值n=-57，过

10、度带宽,最小阻带衰减，这就表明在设置这些值时其参数必须不大于这些值。而其它带阻滤波器的设计指标则要根据加入噪声的频率来确定。若不能按照这些来设计滤波器则不可能虑掉噪声。当所有的指标都设置完后，可以用这些数字来计算上下边带的中心频率和频率间隔，并计算布莱克曼窗设计该滤波器所需要的阶数和产生几阶的布莱克曼窗。当所有的准备工作完成后就可以调用自编的函数计算理想带阻滤波器的脉冲响应和用窗函数法计算实际的滤波器的脉冲响应。最后调用freqz函数得到滤波器的频率特性。从画出的图中可以清楚的看见滤波器的幅频和相频特性。下面是用布莱克曼窗设计滤波器的整个程序：fpd=1800;fsd=2050;fsu=195

11、0;fpu=2000;Rp=1;As=70; % 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi; M=ceil(10*pi/dw)+1; % 计算布莱克曼窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_black=blackman(M); % 产

12、生M阶的布莱克曼窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_black.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性图4：滤波器的频率特性4.3信号滤波处理对语音信号信号进行滤波处理主要是滤掉加入的噪声。不同的滤波器利用不同的函数对语音信号进行滤波，FIR滤波器利用函数fftfilt对信号进行滤波。对信号进行滤波处理要先利用函数filter对y进行滤波，然后对y进

13、行傅里叶变换。而画频谱时只取前面一半。最后在同一个图中画出原始信号的、加入噪声的语音信号以及滤波后语音信号的频域图和时域图。这样便于将所有的图进行对比和分析，而且还可以直观的观察该课程设计是否成功。当将设计好的滤波器滤掉噪声后我们也可以再一次调用函数sound，听此时的声音是否与原始语音信号基本一样，若没有单频啸叫声则说明此次设计是成功的，否则应重新设置指标。下面是对语音信号进行滤波的程序：y_fil=filter(h_bs,1,y); %用设计好的滤波器对y进行滤波Y_fil=fft(y_fil); %对y进行傅里叶变化Y_fil=Y_fil(1:length(Y_fil)/2); % 计算

14、频谱取前一半图5：滤波前后语音信号对比4.4结果分析要确定本课程设计是否成功就得看原始信号的频域图和时域图与经过滤波器后的语音信号的频域图和时域图是否一样，若一样则表示该设计是成功的，否则是不成功的。在第一个图中：第一幅图和第二幅图是原始语音信号的时域图和频域图,第三幅图和第四幅图是加入频率为2100的噪声。从图中可以看出，第一图和第三图相比因为加入噪声的缘故所以第三图y轴的幅度要比第一图要大，但其形状还是基本没有改变。而第二图与第四图相比较在频率f=2100时多了一个尖锐脉冲。说明原始语音信号加入噪声是成功的。在滤波器频率特性的图中可以看到：第一个图是以db为单位的幅频特性，第二图是幅频特性

15、，第三个图是滤波器的相频特性，最后一个图是滤波器的脉冲响应。从图中可以清楚的了解滤波器的幅频和相频特性。在滤波前后信号比较的图中我们可以得到：原始的语音信号与滤波后的信号的图基本一样，只是滤波后的图在原始信号的基础上有所延迟。所以用布莱克曼窗设计的滤波器是符合要求的，也就是说该课程设计是成功的。五、结论与心得为期两周的数字信号处理课程设计已经结束了，但在这次设计中我学到了许多的东西。通过这次的设计，不仅加深了我对课本基础理论知识的理解，而且增强了我的实践能力，同时更加认识到理论知识和实践结合的重要性。首先，更加深入理解了滤波器设计的各个关键环节，包括在什么情况下使用哪种方法设计FIR滤波器最好

16、以及在选择特定的窗函数进行滤波器的设计时我们应该怎样确定其性能指标；其次，更加深刻的认识了语音原始信号与加噪后语音信号的波形及频谱；再次，较大地提高了综合运用专业基础知识及软件设计能力，在一定程度上对自己的动手能力有很大的帮助。虽然这次课程设计已经完成了，但是遇到的困难也是很多的。其中最主要的问题要属怎样设置滤波器的指标问题，如果指标的设置有问题那么后续的工作就不可能得到原始的语音信号。在设置过程中有很多次因为设置的参数不合适而导致设计的滤波器不能虑出单频噪声信号。所以在设计指标问题时一定要结合布莱克曼本身的特点还要考虑加入噪声的频率。其次就是一些函数的细节问题。虽然在这次课程设计中遇到很多的

17、困难，但通过自己查找有关资料以及老师和同学的帮助下都一一解决了，而且在与同学交流的过程中使同学之间的感情更进一步。在此向帮助我的老师及热心同学表示忠心的感谢！希望今后还能参加更多的课程设计，以锻炼自己在各个方面的能力，尤其是综合运用专业基础知识和实践结合的能力。设计的过程中，我通过查阅大量有关资料，与同学交流经验和自学，并向老师请教等方式，使我学到了不少的东西，虽然有许多的辛酸，但是看到自己课程设计完成后心中的那份激动是无法用言语来形容的。六、参考文献1 康华光模拟电子技术基础（第五版）M.北京：高等教育出版社，20062 童诗白，华成英，模拟电子技术基础（第三版）M.高等教育出版社，2006

18、3 何方白，张德民，阳莉，李强.数字信号处理M.高等教育出版社，20114 王靖,李永全.椭圆滤波器Matlb设计与实现J.现代电子技术,2007七、附录7.1 布莱克曼窗设计FIR滤波器的MATLAB主程序x,fs=audioread(test_yfw.mp3); % 使用audioread读取音频文件的采样率fssound(x,fs); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1);y=x+0.1*sin(fn*2*

19、pi*t);sound (y,fs); %明显听出有尖锐的单频啸叫声X=abs(fft(x);Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换X=X(1:N/2);Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围figuresubplot(2,2,1);plot(t,x);xlabel(时间（单位：s）);ylabel(幅度);title(原始语音信号);subplot(2,2,2);plot(f,X);xlabel(频率（单位：Hz）);ylabel(幅度谱);tit

20、le(语音信号幅度谱);subplot(2,2,3);plot(t,y);xlabel(时间（单位：s）);ylabel(幅度);title( 加入单频干扰后的语音信号);subplot(2,2,4);plot(f,Y);xlabel(频率（单位：Hz）);ylabel(幅度谱);title(加入单频干扰后的语音信号幅度谱);fpd=1800;fsd=2050;fsu=1950;fpu=2000;Rp=1;As=70; % 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率，和频率间隔

21、wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(10*pi/dw)+1; % 计算布莱克曼窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_black=blackman(M); % 产生M阶的布莱克曼窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_black.*hd_bs; % 用窗口法

22、计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性figuresubplot(2,2,1);plot(w/pi,db);title(滤波器的db);xlabel(w/pi);ylabel(db);subplot(2,2,2);plot(w/pi,mag);title(滤波器的幅频特性);xlabel(w/pi);ylabel(mag);subplot(2,2,3);plot(w/pi,pha);title(滤波器的相频特性);xlabel(w/pi);ylabel(pha);subplot(2,2,4);plot(h_

23、bs);title(滤波器的脉冲响应 );xlabel(w/pi);ylabel(h_bs);y_fil=filter(h_bs,1,y); %用设计好的滤波器对y进行滤波Y_fil=fft(y_fil); %对y进行傅里叶变化Y_fil=Y_fil(1:length(Y_fil)/2); % 计算频谱取前一半figuresubplot(3,2,1);plot(t,x);title(未加噪声的时域图);xlabel(t);ylabel(x);subplot(3,2,2);plot(f,X);title(未加噪声的频域图);xlabel(f);ylabel(X);subplot(3,2,3);p

24、lot(t,y);title(加噪声后的时域图);xlabel(t);ylabel(y);subplot(3,2,4);plot(f,Y);title(加噪声后的频域图);xlabel(f);ylabel(Y);subplot(3,2,5);plot(t,y_fil);title(加噪声后的频域图);xlabel(f);ylabel(Y);title(滤去噪声后的时域图);xlabel(t);ylabel(y_fil);subplot(3,2,6);plot(f,Y_fil);title(滤去噪声后的频域图);xlabel(f);ylabel(Y_fil);7.2 调用函数ideal_lp.m

25、的MATLAB源代码function hd=ideal_lp(wc,M);% Ideal LowPass filter computation% -% hd = ideal_lp(wc,M);% hd = ideal impulse response between 0 to M-1% wc = cutoff frequency in radians% M = length of the ideal filter%alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha +eps; % add smallest number to avoi divided by

26、zerohd = sin(wc*m)./(pi*m);7.3 调用函数freqz_m.m的MATLAB源代码function db,mag,pha,grd,w = freqz_m(b,a)% Modified version of freqz subroutine% -% db,mag,pha,grd,w = freqz_m(b,a);% db = Relative magnitude in dB computed over 0 to pi radians% mag = absolute magnitude computed over 0 to pi radians % pha = Phase

27、 response in radians over 0 to pi radians% grd = Group delay over 0 to pi radians% w = 501 frequency samples between 0 to pi radians% b = numerator polynomial of H(z) (for FIR: b=h)% a = denominator polynomial of H(z) (for FIR: a=1)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);% pha = unwrap(angle(H); grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0 grd(1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;