数字信号处理滤波器设计

滤波器设计课程设计题目一、FIR（选取合适的窗函数及窗长度）通带0~6000Hz±0.01通带波动阻带8500Hz~Fs/2-70dB44100Hz阻带波动抽样频率二、IIR（可选巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型）1、设计数字低通滤波器，满足以下指标：通带边带频率通带波动2000Hz0.01阻带边界频率2700Hz0.0511025Hz阻带波动抽样频率要求：1、求出h(n)和H(ejω)。(a)用H(ejω)直接画出幅度响应曲线（dB）及相位响应曲线;(b)做h(n)的FFT,有此画出幅度响应曲线（dB）及相位响应曲线；(c)并两种方法的结果进行比较。(d)验证你的设计结果满足设计指标。FIR和IIR滤波器设计过程FIR滤波器的设计：1.设计步骤：（1）写出理想低通滤波器的传输函数和单位脉冲响应。11sin⁡(𝑤𝑐𝑛)𝑗𝑤𝑗𝑤𝑗𝑤ℎ𝑑n=H𝑒𝑒𝑑𝑤=𝑒𝑑𝑤=2𝜋2𝜋𝜋𝑛−𝜋−𝑤𝑐𝜋𝑤𝑐（2）根据题目已有的模拟滤波器频率进行归一化求得数字滤波器的频率，并对题目所给的通带波动和阻带波动进行转换，可求得通带截止频率𝑤𝑝=0.8549𝐻𝑧，阻带截止频率𝑤𝑠=1.2110𝐻𝑧，通带波动𝑎𝑝=0.174𝑑𝐵，阻带衰减𝑎𝑠=70𝑑𝐵，（3）选择用布莱克曼窗函数来设计FIR滤波器，通过查阅资料可知道布莱克曼窗的过渡带宽为11πM，最小衰减为74𝑑𝐵，则可以通过matlab中ceil函数向上取整的方式以及阻带频率和通带频率之差所求的过渡带频率来确定滤波器的长度,在经过调试可以取得N=99。当然，在设计中可以之间用blackman()这个函数，但是在这里我们会利用布莱克曼窗函数的公式构造布莱克曼窗：wn=[0.42−0.5cos2𝜋𝑛4𝜋𝑛+0.08cos⁡()]𝑅(𝑛)𝑁−1𝑁−1𝑁（4）通过求得的矩形窗函数和单位脉冲响应ℎ𝑑n进行时域上的乘积即可求得FIR滤波器的时域响应，再利用DTFT直接求出H(𝑒𝑗𝑤)，并且可以画出幅度响应曲线（dB）及相位响应曲线。∞X𝑒𝑗𝑤=DTFTxn=𝑥(𝑛)𝑒−𝑗𝑤𝑛𝑛=−∞（5）根据时域所求的函数做FFT并由此画出幅度响应曲线（dB）及相位响应曲线，对这两种方法的结果进行比较并验证设计结果是否满足设计指标。3.设计过程：（1）求出h(n)，并做h(n)的DTFT从而直接求出H(𝑒𝑗𝑤)的幅度曲线响应（dB）和相位响应曲线。由设计步骤中可以知道𝑑n=sin⁡(𝑤𝑐𝑛)𝜋𝑛以及自己通过公式所构造的布莱克曼窗函数即可求得h(n)。然后将h(n)经过DTFT即可求得H(𝑒𝑗𝑤)，进而求得幅度曲线响应和相位响应曲线。波形如下：图一布莱克曼窗函数的波形图二理想滤波器冲击响应hd[n]的波形图三h(n)的波形𝐡(𝐧)的前100个值Columns1through13-0.0000-0.0000-0.0000-0.00000.00000.00010.0000-0.0001-0.0002-0.00010.00020.00040.0003Columns14through26-0.0003-0.0008-0.00050.00050.00140.0010-0.0008-0.0023-0.00160.00120.00360.0026-0.0016Columns27through39-0.0053-0.00410.00210.00770.0062-0.0026-0.0109-0.00930.00310.01540.0139-0.0035-0.0220Columns40through52-0.02110.00390.03270.0340-0.0043-0.0549-0.06480.00450.13910.27290.32880.27290.1391Columns53through650.0045-0.0648-0.0549-0.00430.03400.03270.0039-0.0211-0.0220-0.00350.01390.01540.0031Columns66through78-0.0093-0.0109-0.00260.00620.00770.0021-0.0041-0.0053-0.00160.00260.00360.0012-0.0016Columns79through91-0.0023-0.00080.00100.00140.0005-0.0005-0.0008-0.00030.00030.00040.0002-0.0001-0.0002Columns92through99-0.00010.00000.00010.0000-0.0000-0.0000-0.0000-0.0000图四H(𝑒𝑗𝑤)的幅度响应曲线图四H(𝑒𝑗𝑤)的相位响应曲线（2）根据时域所求的函数h(n)做FFT并由此画出幅度响应曲线（dB）及相位响应曲线。图五H(𝑒𝑗𝑤)的幅度响应曲线图六H(𝑒𝑗𝑤)的相位响应曲线由题目可知滤波器的指标：通带波动为0.174dB，阻带衰减为-70dB，通带截止频率为0.8549Hz，阻带截止频率为1.2110Hz。通过DTFT直接求H(𝑒𝑗𝑤)所得到的FIR滤波器的指标是：通带波动为0.001dB，阻带衰减为-73.14dB。通过FFT求H(𝑒𝑗𝑤)所得到的FIR滤波器的指标是：通带波动为0.0014dB，阻带衰减为-75dB。图七H(𝑒𝑗𝑤)的幅度响应曲线(|H𝑒𝑗𝑤|)由图六可以知道这个滤波器两个正负尖峰频率差值为：1.2612-0.80381=0.45739Hz设计的的FIR滤波器的过渡带宽为：1.2110-0.8549=0.3561Hz图八DTFT和FFT所求的幅度响应(dB)曲线由该图可以看出DTFT所求幅度响应（dB）的曲线和FFT所求的幅度响应（dB）的曲线是基本一致的。图八DTFT和FFT所求的相位响应曲线综上所述：1.所设计的布莱克曼窗FIR滤波器在阻带衰减上接近了-70(dB)的指标，但是在通带波动上0.02(dB)略小于所要求的0.174(dB)的指标。2.FIR滤波器的过渡带宽小于两个正负尖峰频率差值，满足要求。3.DTFT和FFT两种方法所设计的滤波器结果及指标基本没有区别，但是DTFT的计算量会比FFT相对大，而FFT所求的幅度和相位响应的波形是对称的，需要经过处理才能达到要求。4.DTFT所算出来的指标会相对来说会更接近于题目中要求的指标。IIR滤波器的设计：1.设计步骤：（1）选择采用冲激响应不变法设计巴特沃斯低通滤波器，已知通带边带频率为2000Hz，阻带边带频率为2700Hz，通带波动为0.01，阻带波动为0.05，抽样频率为11025Hz，根据对频率的归一化以及对阻带波动通带波动转换成（dB）的形式可得，Wp=1.1398Hz，Ws=1.5387Hz，Rp=0.087（dB）,As=26(dB)。（2）计算样本模拟滤波器的阶次100.1𝐴𝑠−1𝛺𝑠N≥（lg0.1𝑅𝑝)/[2lg⁡()]𝛺𝑝10−1可求得N=17。（3）发现Rp≠3(dB)，𝛺𝑝≠𝛺𝑠，由于巴特沃斯滤波器归一化低通原型的通带截止频率𝛺𝑐=1，去归一化成任意截止频率时必须是3（dB）衰减处的𝛺𝑐，利用通带或者阻带满足衰减指标的公式去求得3(dB)衰减的截止频率。𝛺𝑐≥𝛺𝑝/100.1𝑅𝑝−1𝛺𝑐≤𝛺𝑠/100.1𝑅𝑝−1在这里选择的是阻带满足衰减指标的公式，可求得𝛺𝑐=3.8003e+03。（4）根据𝛺𝑐和阶数N求模拟低通滤波器的系统函数分子分母系数。（5）冲激响应不变法利用模拟低通滤波器的系统函数和采样频率转换成数字低通滤波器的系统函数。（6）利用数字低通滤波器的系统函数的分子分母系数求得冲激响应即数字低通滤波器的时域响应。（7）根据所求的时域响应作FFT以求得在频率上的响应。（8）检测步骤（5）中的求得的模拟滤波器系统函数分子分母系数以及数字滤波器系统函数分子分母系数来检验衰减指标。3.设计过程：（1）将要求的各种指标代入巴特沃斯函数中以求得阶数和3（dB）衰减的截止2𝑁2𝑁频率，并由此确定模拟滤波器系统函数的分子分母系数，再利用impinvar将模拟滤波器系统函数的分子分母系数转换成数字滤波器系统函数的分子分母系数。然后利用freqz（）函数求出数字滤波器频域上的幅度响应曲线（dB）和相位响应曲线。图一H(𝑒𝑗𝑤)的幅度响应曲线图二H(𝑒𝑗𝑤)的相位响应曲线（2）将（1）中所求的数字滤波器系统函数的分子分母的系数利用imz函数即可求得系统的冲激响应即时域函数h(n)。图三h(n)的波形𝐡(𝐧)的前100个值Columns1through13-0.0000-0.00000.00000.00000.00040.00510.03090.10830.24340.36100.33520.1315-0.1050Columns14through26-0.1764-0.05330.09240.0972-0.0112-0.0796-0.03560.04150.0494-0.0062-0.0414-0.01650.0239Columns27through390.0250-0.0060-0.0226-0.00670.01440.0124-0.0050-0.0123-0.00220.00860.0060-0.0037-0.0066Columns40through52-0.00040.00500.0028-0.0026-0.00350.00030.00290.0012-0.0017-0.00180.00040.00160.0005Columns53through65-0.0010-0.00090.00040.00090.0002-0.0006-0.00040.00030.00050.0000-0.0004-0.00020.0002Columns66through780.0003-0.0000-0.0002-0.00010.00010.0001-0.0000-0.0001-0.00000.00010.0001-0.0000-0.0001Columns79through91-0.00000.00000.0000-0.0000-0.000000000000Columns92through1040000000000000（3）根据时域所求的函数h(n)做FFT并由此画出幅度响应曲线（dB）及相位响应曲线。图四H(𝑒𝑗𝑤)的幅度响应曲线图五H(𝑒𝑗𝑤)的相位响应曲线由题目可知滤波器的指标：通带波动为0.087dB，阻带衰减为26dB，通带截止频率为1.1398Hz，阻带截止频率为1.5387Hz。由于求利用impz函数求h(n)只是数字滤波器的单位冲激响应逼近模拟滤波器的单位冲激响应的过程，再进行FFT变换所得到的数字滤波器的指标误差相对较大，所以在这里只检验直接求频域响应所得到的数字滤波器的指标。利用书上检验衰减指标的代码并将参数代入可求得：通带波动为0.0634dB，阻带衰减为26dB，非常接近题目中要求的指标。图六直接求和FFT所求的幅度响应（dB）曲线图七直接求和FFT所求的相位响应曲线由以上两幅图可以知道直接求所得到的幅度响应曲线及相位响应曲线在w/pi=0~0.7这个范围内的波形基本保持一致，但是在w/pi＞0.7之后FFT所得到的幅度响应的曲线的谐波较多，而且呈现的也不是线性相位。综上所述，直接求H(𝑒𝑗𝑤)所设计的低通滤波器基本能够满足题目要求的指标，在比较了直接求H(𝑒𝑗𝑤)和FFT设计的数字低通滤波器后，认为在设计低通巴特沃斯滤波器上还是用直接求H(𝑒𝑗𝑤)的方式较为精确。设计过程中遇到的问题和解决方法：1.FIR滤波器中在求窗长度的时候11pi/(ws-wp)得不到整数。解决方法：用ceil函