使用布莱克曼窗设计的FIR滤波器正文之欧阳学创编.docx
《使用布莱克曼窗设计的FIR滤波器正文之欧阳学创编.docx》由会员分享，可在线阅读，更多相关《使用布莱克曼窗设计的FIR滤波器正文之欧阳学创编.docx（12页珍藏版）》请在冰豆网上搜索。
使用布莱克曼窗设计的FIR滤波器正文之欧阳学创编
音乐信号滤波去噪
时间：
2021.03.03
创作：
欧阳学
——使用布莱克曼窗设计的FIR滤波器
摘要本课程设计主要是用麦克风采集一段语音信号，绘制其波形并观察其频谱。
然后在该语言信号中加一个噪音，利用布莱克曼窗设计一个FIR滤波器，对该语音信号进行虑噪处理，然后比较滤波前后的波形与频谱。
在本课程设计中，是用MATLAB的集成环境完成一系列的设计。
首先对加噪的语音信号进行虑波去噪处理，再比较滤波前后的频率响应曲线，若一样则满足所设计指标，否则不满足。
也可以调用函数sound听滤波前后其语音信号是否带有噪声。
若无噪声也说明该滤波器的设置也是成功的。
关键词音乐信号；MATLAB；FIR滤波器；滤波去噪
1引言
人们在语音通信的过程中将不可避免的会受到来自周围环境的干扰，例如传输媒介引入的噪声，通信设备内部的电噪声，乃至其他讲话者的话音等。
正因为有这些干扰噪声的存在，接受者接受到的语音已不是原始的纯净语音信号，而是受噪声干扰污染的带噪声语音信号。
而本课程设计就是利用MATLAB集成环境用布莱克曼窗的方法设计一个FIR滤波器，对语音信号进行滤波去噪处理，并将虑噪前后的频谱图进行对比。
1.1课程设计目的
数字信号处理课程设计是数字信号处理课程的重要实践性环节，是学生在校期间一次较全面的工程师能力训练，在实现学生总体培养目标中占有重要地位。
综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为编程工具进行计算机实现，从而复习巩固了课堂所学的理论知识，提高了对所学知识的综合应用能力，并从实践上初步实现了对数字信号的处理。
本课程设计能使学生对通信工程领域各种技术的DSP实现的设计有较熟练的掌握。
且通过自身的实践，对DSP的设计程序、内容和方法有更深入的掌握，提高实际运用的能力。
并可综合运用这些知识解决一定的实际问题，使学生在所学知识的综合运用能力上以及分析问题、解决问题能力上得到一定的提高。
1.2课程设计的要求
（1）滤波器指标必须符合工程实际。
（2）设计完后应检查其频率响应曲线是否满足指标。
（3）处理结果和分析结论应该一致，而且应符合理论。
（4）独立完成课程设计并按要求编写课程设计报告书。
1.3设计平台
MATLAB名称是有两个英文单词Matrix和Laboratory的前三个字母组成。
MATLAB7.0是美国MathWorks公司开发的优秀计算软件MATLAB的最新版本。
MATLAB自20世纪80年代面世以来，以其强大的数值计算能力、优秀的绘图功能以及与其他软件良好的交互功能在众多的数学计算软件中独领风骚，特别是它源代码的开放性使用户可以二次开发，受到了广大使用者的格外赞赏。
MATLAB是一个为科学和工程计算机专门设计的交互式大型软件，是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。
它集图和精确计算与一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到了广泛应用。
它不仅是一个在各类工程设计中便于使用的计算工具，在世界各地的高等院校中十分流行，在各类工业应用中更有不俗的表现。
MATLAB可以几乎所有的PC机和大型计算机上运行，适用于Window、UNIX等多种系统平台。
本课程设计我们就可以直接诶使用MATLAB提供的模块，实现模拟通信系统的仿真。
MATLAB软件有很强的开放性和适应性。
在保持内核不变的情况下，MATLAB可以针对不同的应用学科推出相应的工具箱，目前已经推出了图像处理工具箱、信号处理工具箱、小波工具箱、神经网络工具箱以及通信工具箱等多个学科的专用工具箱，极大的方便了不同学科的研究工作。
国内已有越来越多的科研和技术人员认识到MABLAB的强大作用，并在不同领域内使用MATLAB来快速实现科研构想和提高工作效率。
2设计原理
2.1FIR滤波器
FIR（FiniteImpulseResponse）滤波器：
有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。
因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。
FIR数字滤波器以其良好的线性特性被广泛应用于现代电子通信系统中，是数字信号处理的重要内容之一。
在实际信号处理中，往往要求系统兼具实时性和灵活性，而已有的一些软件或硬件实现方案（如DSP）则难以同时达到这两方面的要求。
使用具有并行处理特性的FPGA来实现FIR滤波器，既有很强的实时性，又兼顾了灵活性，为数字信号处理提供了一种很好的解决方案。
FIR滤波器系数计算较为繁琐，在设计时借助Matlab工具箱，选择合适的窗函数，可以方便地计算滤波器系数，并分析其幅频、相频特性。
有限长单位冲激响应（FIR）滤波器有以下特点：
（1）、系统的单位冲激响应h（n）在有限个n值处不为零；
（2）、系统函数H（z）在|z|>0处收敛，极点全部在z=0处（因果系统）；
（3）、结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。
设FIR滤波器的单位冲激响应h（n）为一个N点序列，0≤n≤N—1，则滤波器的系统函数为：
H（z）=
-n（2-1）
就是说，它有（N—1）阶极点在z=0处，有（N—1）个零点位于有限z平面的任何位置因此H（z）是永远稳定的。
稳定和相位特性是FIR滤波器突出的优点。
FIR滤波器有以下几种基本结构：
横截型（卷积型、直接型）、级联型、频率抽样型、快速卷积结构。
FIRDF的设计方法主要分为两类：
第一类是基于逼近理想滤波器特性的方法，包括窗函数法、频率采样法和等波纹最佳逼近法；第二类是最优设计法。
2.2窗口设计法
数字信号处理的主要数学工具是博里叶变换．而傅里叶变换是研究整个时间域和频率域的关系。
不过，当运用计算机实现工程测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。
做法是从信号中截取一个时间片段，然后用观察的信号时间片段进行周期延拓处理，得到虚拟的无限长的信号，然后就可以对信号进行傅里叶变换、相关分析等数学处理。
无线长的信号被截断以后，其频谱发生了畸变，原来集中在f（0）处的能量被分散到两个较宽的频带中去了（这种现象称之为频谱能量泄漏）。
为了减少频谱能量泄漏，可采用不同的截取函数对信号进行截断，截断函数称为窗函数，简称为窗。
窗函数设计法的基本思路是用FIRDF逼近希望的滤波特性。
设希望逼近的滤波器的频率响应函数为hd（ejω），其单位脉冲响应用hd（n）表示。
为了设计简单方便，通常选择hd（ejω）为具有片段常数特性的理想滤波器。
因此hd（n）是无限长非因果序列，不能直接作为FIRDF的单位脉冲响应。
窗函数设计法就是截取hd（n）为有限长的一段因果序列，并用合适的窗函数进行加权做为FIRDF的单位脉冲响应h（n）。
下面介绍窗函数设计法的基本设计过程。
窗口设计法的主要工作是计算hd（n）和w（n），但当Hd（ejω）较为复杂时，hd（n）就不容易由反付里叶变换求得。
这时一般可用离散付里叶变换代替连续付里叶变换，求得近似值。
窗口法的设计步骤如下：
（1）、通过傅里叶变换忽的理想滤波器的单位脉冲响应hd（n）。
（2）、根据指标选择窗口形状、大小和位置。
确定窗口类型的主要依据是过渡带宽和阻带最小衰耗的指标。
（3）、给定理想频响由hd（ejω）和hd（n），加窗得h（n）=w（n）hd（n）。
（4）、检验滤波器的性能。
由h（n）求H（ejω）是否在误差容限之内。
如果不满足，则返回第
（2）步。
以上步骤中hd（n）、H（ejω）的计算可采用傅氏变换的现成公式和程序，窗函数w（n）也是现成的。
但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算出。
整个设计可用计算机编程来做。
窗口法的优点是简单，有闭合的公式可用，性能及参数都有表格资料可查，计算程序简单，较为实用。
缺点是当Hd（ejω）较为复杂时，hd（n）就不容易由反付里叶变换求得。
边界频率因为加窗的影响而不易控制。
窗口函数对理想特性的影响：
改变了理想频响的边沿特性，形成过渡带，宽为
等于WR（ω）的主瓣宽度；过渡带两旁产生肩峰和余振（带内、带外起伏），取决于WR（ω）的旁瓣，旁瓣多，余振多；旁瓣相对值大，肩峰强，与N无关；N增加,过渡带宽减小,肩峰值不变。
因主瓣附近
（2-2）
其中x=Nω/2,所以N的改变不能改变主瓣与旁瓣的比例关系，只能改变WR（ω）的绝对值大小和起伏的密度，当N增加时，幅值变大，频率轴变密，而最大肩峰永远为8.95%，这种现象称为吉布斯（Gibbs）效应。
肩峰值的大小决定了滤波器通带内的平稳程度和阻带内的衰减，所以对滤波器的性能有很大的影响。
改变窗函数的形状，可改善滤波器的特性，窗函数有许多种，但要满足以下两点要求:
窗谱主瓣宽度要窄，以获得较陡的过渡带；相对于主瓣幅度，旁瓣要尽可能小，使能量尽量集中在主瓣中，这样就可以减小肩峰和余振，以提高阻带衰减和通带平稳性。
但实际上对同样长度的窗这两点不能兼得，一般总是通过增加主瓣宽度来换取对旁瓣的抑制。
2.3布莱克曼窗
布莱克曼窗的时域形式可表示为：
N（n）（2-3）
它的频域特性为：
W
R
R
R
R
R
（2-4）
其中
为矩形窗函数的幅度频率特性。
增加一个二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，为
。
加N可减少过渡带。
布莱克曼窗函数的最大旁瓣之比主瓣值低57db，但是主瓣宽度是矩形窗函数的主瓣宽度的三倍。
布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。
3设计步骤
3.1设计流程图
本课程设计是对录制的语音信号进行加噪处理并分析加噪前后语音信号的时域图与频域图，再用布莱克曼窗设计一个FIR滤波器，而滤波器的设计必须符合其相应的指标，否则将不能滤掉加入的噪声。
最后将滤波前后的波形图进行比较看是否相同。
下面是整个课程设计的流程图如图3.1所示：
图3.1语音信号的整个流程图
3.2录制语音信号
点击桌面上的“开始”菜单，再选择“程序”中的“附件”，在“附件“的菜单栏中选择“娱乐”，最后点击“录音机”。
就可以得到如图3.2的图。
然后点击开始录制语音信号，时间大约在2~3秒之间。
图3.2录制语音的录音机
录制好语音信号后，打开MATLAB软件平台，利用函数waveread对语音信号进行采样，记住采样频率的和采样点数。
再调用函数sound此时可以听见录制的语音。
采样完后再语音信号中加入一个单频噪声，单频的噪声的频率可以自己设置。
按照加入噪声后的采样频率调用sound函数，这时可以明显的听见播放的语音信号中有尖锐的单频啸叫声。
下面是调用该语言信号以及加入噪声的程序：
[x,fs,bits]=wavread（'e:
yuyin.wav'）;%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数
sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放
N=length（x）;%计算信号x的长度
fn=2100;%单频噪声频率，此参数可改
t=0:
1/fs:
（N-1）/fs;%计算时间范围，样本数除以采样频率
x=x（:
1）';
y=x+0.1*sin（fn*2*pi*t）;
sound（y,fs,bits）;%明显听出有尖锐的单频啸叫声
现在是对加入噪声前后的语音信号进行频谱分析，先对原始和加噪后的语音信号进行傅里叶变换，再计算频谱的频率范围和谱线间隔。
最后就可以画出未加入噪声和加入噪声后的时域图和频域图。
将所有未加和加入噪声的时域图和频域图画在同一个图中，便于比较和分析。
下面是对未加和加入噪声的频谱分析的程序：
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;%计算频谱频率范围
图3.3语音信号的时域图和频率图
3.3滤波器设计
本课程设计就是要设计一个滤波器虑掉加入的噪声，使其恢复原始的语音信号。
而设计滤波器的方法有很多，例如：
窗函数法、频率采样法、脉冲响应不变法和双线性变换法。
而本课程设计采用的是窗函数法设计FIR滤波器。
而FIR滤波器的设计也有很多方法。
在Matlab中，可以利用矩形窗、三角窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗等设计FIR滤波器。
而本次采用的是布莱克曼窗来设计滤波器。
在用布莱克曼窗设计滤波器的时候，首先要确定滤波器的性能指标。
从六种窗函数的基本参数中我们可以得到旁瓣峰值
n=-57，过度带宽
B=11
最小阻带衰减
s=74db，这就表明在设置这些值时其参数必须不大于这些值。
而其它带阻滤波器的设计指标则要根据加入噪声的频率来确定。
若不能按照这些来设计滤波器则不可能虑掉噪声。
当所有的指标都设置完后，可以用这些数字来计算上下边带的中心频率和频率间隔，并计算布莱克曼窗设计该滤波器所需要的阶数和产生几阶的布莱克曼窗。
当所有的准备工作完成后就可以调用自编的函数计算理想带阻滤波器的脉冲响应和用窗函数法计算实际的滤波器的脉冲响应。
最后调用freqz函数得到滤波器的频率特性。
从画出的图中可以清楚的看见滤波器的幅频和相频特性。
下面是用布莱克曼窗设计滤波器的整个程序：
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））;%计算上下边带中心频率，和频率间隔
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;%用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=freqz_m（h_bs,1）;%调用自编函数计算滤波器的频率特性
图3.4滤波器的频率特性
3.4信号滤波处理
对语音信号信号进行滤波处理主要是滤掉加入的噪声。
不同的滤波器利用不同的函数对语音信号进行滤波，FIR滤波器利用函数fftfilt对信号进行滤波。
对信号进行滤波处理要先利用函数filter对y进行滤波，然后对y进行傅里叶变换。
而画频谱时只取前面一半。
最后在同一个图中画出原始信号的、加入噪声的语音信号以及滤波后语音信号的频域图和时域图。
这样便于将所有的图进行对比和分析，而且还可以直观的观察该课程设计是否成功。
当将设计好的滤波器滤掉噪声后我们也可以再一次调用函数sound，听此时的声音是否与原始语音信号基本一样，若没有单频啸叫声则说明此次设计是成功的，否则应重新设置指标。
下面是对语音信号进行滤波的程序：
y_fil=filter（h_bs,1,y）;%用设计好的滤波器对y进行滤波
Y_fil=fft（y_fil）;%对y进行傅里叶变化
Y_fil=Y_fil（1:
length（Y_fil）/2）;%计算频谱取前一半
图3.5滤波前后语音信号对比
sound（x，fs，bits）;可以感觉滤波前后的声音有变化。
3.5结果分析
要确定本课程设计是否成功就得看原始信号的频域图和时域图与经过滤波器后的语音
信号的频域图和时域图是否一样，若一样则表示该设计是成功的，否则是不成功的。
在第一个图中：
第一幅图和第二幅图是原始语音信号的时域图和频域图,第三幅图和第四幅图是加入频率为2100的噪声。
从图中可以看出，第一图和第三图相比因为加入噪声的缘故所以第三图y轴的幅度要比第一图要大，但其形状还是基本没有改变。
而第二图与第四图相比较在频率f=2100时多了一个尖锐脉冲。
说明原始语音信号加入噪声是成功的。
在滤波器频率特性的图中可以看到：
第一个图是以db为单位的幅频特性，第二图是幅频特性，第三个图是滤波器的相频特性，最后一个图是滤波器的脉冲响应。
从图中可以清楚的了解滤波器的幅频和相频特性。
在滤波前后信号比较的图中我们可以得到：
原始的语音信号与滤波后的信号的图基本一样，只是滤波后的图在原始信号的基础上有所延迟。
所以用布莱克曼窗设计的滤波器是符合要求的，也就是说该课程设计是成功的。
4出现的问题及解决方法
虽然课程设计已经完成了，但是在设计的过程中还是遇到了许多的问题。
总结起来大概有以下几个方面：
首先，最主要的是要把设计滤波器的参数设置正确才能滤除语音信号中的噪声。
有几次因为前面的噪声频率设置为f=2100，而后面的带阻滤波器设计指标没有与前面的频率相匹配导致画出来的图怎么也不能滤掉噪音。
后来慢慢的改变指标使其在噪声频率的左右，这时才能滤掉噪声得到原始的语音信号。
还有用布莱克曼窗设计的滤波器时期过渡带宽度和阻带最小衰减必须符合布莱克曼窗的基本参数否则也不可能得到所期望的结果。
其次，在对加入噪声后的信号进行滤波时没有使用正确的形式也不能滤去噪声。
后来在同学的帮助下解决了这些问题。
再次，在利用自编函数计算理想带阻滤波器的脉冲响应，如果在MATLAB软件中的Work下没有定义，那么不能调用自编的函数，否则将会报错。
最后，在MATLAB软件下编程时最好新建一个File文档。
因为在编程的过程过有可能出现错误，如果建一个文档有助于程序出现错误时可以在文档中直接修改，这样可以省很多的时间，又这个课程设计的程序多而繁杂，一不小心就有可能写错，如果在工作环境下修改这样利于将所有的程序复制在课程设计中，而且还要对复制后的程序进行删除。
5结束语
为期两周的数字信号处理课程设计已经结束了，但在这次设计中我学到了许多的东西。
通过这次的设计，不仅加深了我对课本基础理论知识的理解，而且增强了我的实践能力，同时更加认识到理论知识和实践结合的重要性。
首先，更加深入理解了滤波器设计的各个关键环节，包括在什么情况下使用哪种方法设计FIR滤波器最好以及在选择特定的窗函数进行滤波器的设计时我们应该怎样确定其性能指标；其次，更加深刻的认识了语音原始信号与加噪后语音信号的波形及频谱；再次，较大地提高了综合运用专业基础知识及软件设计能力，在一定程度上对自己的动手能力有很大的帮助。
虽然这次课程设计已经完成了，但是遇到的困难也是很多的。
其中最主要的问题要属怎样设置滤波器的指标问题，如果指标的设置有问题那么后续的工作就不可能得到原始的语音信号。
在设置过程中有很多次因为设置的参数不合适而导致设计的滤波器不能虑出单频噪声信号。
所以在设计指标问题时一定要结合布莱克曼本身的特点还要考虑加入噪声的频率。
其次就是一些函数的细节问题。
虽然在这次课程设计中遇到很多的困难，但通过自己查找有关资料以及老师和同学的帮助下都一一解决了，而且在与同学交流的过程中使同学之间的感情更进一步。
这次设计不仅让我学会如何独立完成一项工作，而且提高了独立解决问题的能力，为以后的课程设计打下良好的基础。
在此向帮助我的老师及热心同学表示忠心的感谢！
希望今后还能参加更多的课程设计，以锻炼自己在各个方面的能力，尤其是综合运用专业基础知识和实践结合的能力。
设计的过程中，我通过查阅大量有关资料，与同学交流经验和自学，并向老师请教等方式，使我学到了不少的东西，虽然有许多的辛酸，但是看到自己课程设计完成后心中的那份激动是无法用言语来形容的。
参考文献
[1]丁玉美,高西全，阔永红.数字信号处理.第一版.西安：
西安电子科技出版社,2001年
[2]陈后金.数字信号处理.第三版.北京：
高等教育业出版社，2004年
[3]程佩青.数字信号处理教程.第四版.北京：
清华大学出版社，2002年
[4]刘敏，魏玲.Matlab通信仿真与应用.第二版.北京：
国防工业出版社，2001年
[5]张圣勤.MATLAB7.0实用教程.第三版.北京：
机械工业出版社，2006
附录：
用布莱克曼窗设计FIR滤波器的整个源程序
[x,fs,bits]=wavread（'e:
yuyin.wav'）;%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数
sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放
N=length（x）;%计算信号x的长度
fn=2100;%单频噪声频率，此参数可改
t=0:
1/fs:
（N-1）/fs;%计算时间范围，样本数除以采样频率
x=x（:
1）';
y=x+0.1*sin（fn*2*pi*t）;
sound（y,fs,bits）;%明显听出有尖锐的单频啸叫声
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;%计算频谱频率范围
subplot（2,2,1）;plot（t,x）;xlabel（'时间（单位：
s）'）;ylabel（'幅度'）;title（'原始语音信号'）;
subplot（2,2,2）;plot（f,X）;xlabel（'频率（单位：
Hz）'）;ylabel（'幅度谱'）;title（'语音信号幅度谱'）;
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））;%计算上下边带中心频率，和频率间隔
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'
升级会员 