FIR音频滤波.docx

FIR音频滤波.docx

《FIR音频滤波.docx》由会员分享，可在线阅读，更多相关《FIR音频滤波.docx（22页珍藏版）》请在冰豆网上搜索。

FIR音频滤波

音频滤波：

去噪

FIR系统的系统函数

系数b0,b1,…,bM即为系统的单位抽烟响应h（0），h

（1），…，h（M），且当n>M时，h（n）=0。

FIR滤波器的设计主要是建立在对理想滤波器频率特性作某种近似的基础上，这些近似方法有窗函数法、频率抽样法以及等波纹切比雪夫逼近法。

窗函数法设计的基本思想是把给定的频率响应通过，求得脉冲响应，然后利用加窗函数对它进行截断和平滑，以实现一个物理可实现且具有线性相位的FIR数字滤波器的设计目的。

其核心是从给定的频率特性，通过加窗确定有限长单位取样响应h（n）。

频率采样法设计的基本思想是，把给出的理想频率响应进行取样，通过IDFT从频谱样点直接求得有限脉冲响应。

等波纹切比雪夫逼近法，利用MATLAB提供的remez函数实现Remez算法，设计滤波器逼近理想频率响应。

一、窗函数法

1.1设计原理

设计FIR数字滤波器的最简单的方法是窗函数法，通常也称之为傅立叶级数法。

窗函数法设计的基本思想是，加窗截取理想滤波器的单位脉冲响应的一段作为所要设计的滤波器的脉冲响应。

其核心是从给定的频率特性，通过加窗确定有限长单位取样响应h（n）。

设计过程如下：

加窗的作用是通过把理想滤波器的无限长脉冲响应hd（n）乘以窗函数w（n）来产生一个被截断的脉冲响应，即h（n）=hd（n）w（n），并且对频率响应进行平滑。

用窗函数法设计FIRDF的具体设计步骤如下：

（1）选择窗函数类型和长度，写出才窗函数的表达式。

根据阻带最小衰减选择窗函数的的类型，再根据过渡带宽度确定所选窗函数的长度。

用窗函数法设计的FIRDF通带波纹幅度近似等于阻带波纹幅度。

一般阻带最小衰减达到40dB以上，则通带最大衰减就小于0.1dB。

所以用窗函数法设计FIRDF时，通常只考虑阻带最小衰减就可以了。

（2）构造希望逼近的频率响应函数。

根据设计需要，一般选择线性相位理想滤波器（理想低通、理想高通、理想带通、理想带阻）。

理想滤波器的截止频率近似为最终设计的FIRDF的过渡带中心频率，幅度函数衰减一半（约-6dB）。

所以一般取

，

和

分别为通带边界频率和阻带边界频率。

（3）计算

。

（4）加窗得到设计结果：

。

常见的窗函数有矩形窗、bartlett窗、三角窗、Hamming窗、Hanning窗、Blackman窗、Chebyshey窗、Kaiser窗等等。

矩形窗具有最窄的主瓣，但也有最大的边瓣峰值和最慢的衰减速度。

Hamming窗和Hanning窗主瓣稍宽，但有较小的边瓣和较大的衰减速度，是较为常用的窗函数。

1.2实验流程

用麦克风采集一段wav音乐或下载wav音乐

加入单频噪声

对语音信号进行频谱分析，画出时域和频域波形图

设计FIR滤波器

画出其频率响应

用FIR滤波器对语音信号进行滤波

画出语音信号滤波前后波形并且进行比较分析

开始

结束

1.2.1获取wav文件并加噪

wav音乐文件的获取，可以自行采集，也可以从网上下载。

然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使用，让我们很快理解了采样频率、采样位数等概念。

采集完成后在信号中加入一个单频噪声。

（对应程序见附录5.1）

接下来，输出显示原始信号和加噪信号的时域图和频谱图。

（程序见附录5.2）

图1-1原始信号和加噪信号的时域图和频谱图

1.2.2带阻滤波器概念

根据前面所加的噪声知道，这里需要使用带阻滤波器来滤除加上的噪声。

数字带阻滤波器也具有频率响应的周期性，频率变量以数字频率ω来表示（ω=ΩT=Ω/fs，Ω为模拟角频率，T为抽样时间间隔，fs为抽样频率）。

所以，数字滤波器设计中，必须给出抽样频率。

下图显示数字带阻滤波器的理想幅频响应，这样的理想频率响应是不太可能实现的，原因是频带之间幅度响应是突变的，因而其单位抽样响应是非因果的。

一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。

以低通滤波器为例，如下图所示，频率响应有通带、过渡带和带阻三个范围（非理想）。

在通带内幅度响应以误差α1逼近于1，即

在带阻中，幅度响应以误差α2逼近于0，即

其中ωc，ωst分别为通带截止频率和阻带截止频率。

为了逼近理想低通滤波器的特性，还必须有一个非零宽度（ωc-ωst）的滤波带，在这个过渡带内的频率响应平滑地从通带下降到阻带。

1.2.3滤波器设计

在该滤波器的设计中，我们给出该滤波器的性能指标如下：

fpd=1800;fsd=1850;fsu=1950;fpu=2000;Rp=1;As=40;

截止频率也可以任意自选，在单频噪声干扰附近即可。

在这里，很重要的是通带截止db值的设置。

这个值一定要根据我们使用的设计滤波器的方法来设定。

因为我使用的是NUTTALLWIN窗法，NUTTALLWIN窗函数中，滤波器的过渡带宽为15.4，最小阻带衰减为108db。

所以，一定要将通带截止db值设置的小于108，所以，我将其设置为100db。

在这里我是使用窗函数法设计上面要求的FIR滤波器。

在Matlab中，利用NUTTALLWIN窗设计FIR滤波器，利用Matlab中的函数freqz画出该滤波器的频率响应。

首先，我们利用数字信号处理里面学过的知识，根据自己选定的参数，用指定的方法设计FIR滤波器，得到FIR滤波器的阶数M。

随后调用nuttallwin（M）函数产生M阶的NUTTALLWIN窗。

然后，调用自编ideal_lp函数计算理想带阻滤波器的脉冲响应。

最后，再调用自编freqz_m（h_bs,1）函数即可计算得到该滤波器的频率特性。

（程序见附录5.5）

在将滤波器设计好，频率特性求出来之后，编写程序，来画出该带阻滤波器的幅频和相频响应曲线以及滤波器单位脉冲响应曲线。

如下图所示：

图1-2频率响应曲线

图1-3单位脉冲响应曲线

从图中可以看出，幅频响应曲线满足带阻滤波器特性，我们计算的该滤波器的阶数是3397，非常的大。

与下图300阶的滤波器幅频响应曲线相比，要好很多，前者从通道到阻带和从阻带到通带过渡带要债一些，衰减也大一些。

另外，从相频响应曲线中可以看出，该滤波器满足线性相位。

图1-4300阶滤波器的频率响应曲线

1.2.4滤波以及结果分析

在将滤波器设计好后，我们用自己设计的带阻滤波器对采集的语音信号进行滤波。

在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。

在将加噪信号滤波之后，我们将滤波前后语音信号的波形及频谱图相互比较。

在同一张大图里分别绘制原始信号x，加噪信号y，滤波去噪信号y_fil的时域波形和频谱，以便比较和分析。

（程序见附录5.4）

三个信号的时域波形和频谱图如下图所示：

图1-5三个信号的时域波形和频谱图

从图中我们可以看出，原信号与滤波去噪信号的时域图基本相似，只有边缘部分有点差异；原信号与滤波去噪信号的频谱图波形也大致相似。

通过观察可以看到，加噪信号的时域图中部分被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在1900Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。

二、频率抽样法

2.1设计思路

频率抽样法的设计思路，是对理想频率响应等间隔抽样，作为实际FIR数字滤波器的频率特性的抽样值。

k=0,1,…,N-1

知道H（k）后，由IDFT定义，可以用这N个采样值H（k）来唯一的确定有限长序列h（n），即

n=0,1,2,…,N-1

h（n）为待设计的滤波器的单位脉冲响应。

其系统函数H（z）为

以上就是频率抽样法设计滤波器的基本原理，即通过有限的频率特性取样值去逼近理想滤波特性，然后由有限的频率特性取样值（如系统冲击响应的DFT）取得系统函数。

窗函数设计法是从时域出发，把理想的单位脉冲响应hd（n）用一定形状的窗函数截取成有限长的h（n），以h（n）来近似hd（n），从而使得频率响应函数

近似理想频率响应

。

频率抽样法则是从频域出发，对理想的频率响应

进行等间隔取样，以有限的频率响应

去逼近理想频率响应。

2.2滤波器设计

在Matlab中，fir2函数可以实现任意FIR数字滤波器的频率取样法设计，可以设计任意形状频率响应的滤波器。

（Matlab中fir2函数的功能说明是“Frequencysampling-basedfiniteimpulseresponsefilterdesign”，即基于频率抽样的有限响应FIR滤波器设计）格式如下：

b=fir2（N,F,A）或b=fir2（N,F,A,window）

输出参数：

b为FIR数字滤波器的N+1个系数构成的矩阵。

输入参数：

N为滤波器的阶数。

F指定归一化的各频带边界频率，从0到1递增，1对应于fsam/2，即数字频率π。

A指定各频带边界频率处的幅度响应，因此F和A的长度相等。

window指定窗函数，若不指定，则默认为哈明窗。

在将滤波器设计好，频率特性求出来之后，编写程序，来画出该带阻滤波器的幅频和相频响应曲线以及滤波器单位脉冲响应曲线。

如下图所示：

图2-1滤波器频率响应曲线

图2-2滤波器单位脉冲响应曲线

同样是3397阶的滤波器，上图中的幅频响应曲线，其衰减要小于窗函数法中设计的滤波器，即频率抽样法设计的滤波器性能不及窗函数法。

该滤波器同样满足线性相位。

2.3滤波结果分析

原始信号和加噪信号和之前几种方法中所使用的相同。

在将加噪信号滤波之后，我们将滤波前后语音信号的波形及频谱图相互比较。

在同一张大图里分别绘制原始信号x，加噪信号y，滤波去噪信号y_fil的时域波形和频谱，以便比较和分析。

（程序见附录5.4）

三个信号的时域波形和频谱图如下图所示：

图2-3三个信号的时域波形和频谱图

从图中我们可以看出，原信号与滤波去噪信号的时域图基本相似，只有边缘部分有点差异；原信号与滤波去噪信号的频谱图波形也大致相似。

通过观察可以看到，加噪信号的时域图中部分被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在1900Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。

2.4改善滤波器性能的措施

如果给出的理想低通滤波器在通带的频谱

等于1而阻带为0，则不论样点N取的如何密，在临界频率处总有两个幅度突变的样点，它们之间的落差为1。

于是阻带边缘产生反冲和阻尼振荡，其最大幅度取决于sinc函数，是个固定的值。

这样设计出来的滤波器的阻带最小衰耗固定为-20dB，与矩形窗一样。

增加采样点数N并不能改善阻带最小衰耗，改善阻带衰耗的唯一办法是加宽过渡带。

具体方法是，在通、阻带交界处，人为地安排一到几个过渡点，其值介于0和1之间，这样可以减小样点间的落差，使过渡平滑，反冲减小，阻带最小衰耗增大。

经验表明，每增大一个过渡点，过渡带增宽2π/N，最优情况下阻带衰耗可增大20~30dB。

兼顾过渡带宽和阻带最小衰减，需要增加采样点，同时在通带、阻带交界处安排过渡点。

频率采样法特别适用于设计窄带选频滤波器。

因为只有少数几个非零值的H（k），计算量大为降低。

但是，由于频率抽样点的分布必须符合一定的规律，在规定通阻带截止频率方面不够灵活。

比如，当截止频率ωc不是2π/2N整数倍时，会产生较大的逼近误差，因而该方法的应用不及窗口法普遍。

另外，在FDATool中，带阻FIR滤波器的设计方法只有窗函数法、等波纹逼近法这一类、最小二乘法这一类。

没有单独列出频率抽样法，是因为频率抽样法和窗函数法是很相似了，前一个是在频域加窗截取，后一个是在时域加窗截取。

三、等波纹切比雪夫逼近法

3.1滤波器设计

窗口法设计和频率采样法设计都存在某些缺陷。

首先，在设计中不能将边缘频率ωp和ωs精确地给定。

其次，不能够同时标定波纹因子δ1和δ2。

最后，近似误差（即理想响应和实际响应之差）在频率区间上不是均匀分布的。

而等波纹切比雪夫逼近法，即最优等波纹设计法能解决上面三个问题。

对于线性相位FIR滤波器来说，有可能导得一组条件，对这组条件能够证明，在最大近似误差最小化的意义下这个设计解是最优的。

具有这种性质的滤波器就是通过最优等波纹设计法设计得到的等波纹滤波器。

切比雪夫近似问题现在能定义为：

确定这组系数

或

或

或

[或等效为a（n）或b（n）或c（n）或d（n）]以使在通带和阻带内E（w）的最大绝对值最小，即

（2..2）

交错点定理：

设S是闭区间[0,π]内任意闭合子集，为使P（ω）是在S上对

的唯一最大值最小近似，其必要与充分条件是E（ω）在S内至少呈现出（L+2）个“交错点”或极值频率；这就是说，在S内一定存在（L+2）个频率

使之有

E（

）=-E（

）=

（2.3）

将这个定理与前面的结论结合在一起，表明最优等波纹滤波器在S内它的误差函数不是有（L+2）个就是有（L+3）个交错点。

Parks-McClellan算法：

假定滤波器长度M（或L）和比值

已知，选取加权函数，正确的选定阶M，当这个解得到时就有

。

M和

是互为关联的；M愈大，

就愈小。

近似的M可由下面这个公式得到。

（2.4）

Parks-McClellan算法从估算（L+2）个极值频率{

}开始并估计出在这些频率上的最大误差

。

然后通过（2.3）式给出的点拟合一个L阶多项式。

在一个很细的密度上确定局部最大误差，并在这些新的极值上调整极值频率{

}。

通过这些新的极值频率又拟合出一个新的L阶多项式，这个过程一直重复下去。

这一迭代过程一直持续到最优一组频率{

}和全局最大误差

被找到为止。

最后脉冲响应被计算出来。

由于是对M的近似，最大误差

可能不等于

。

若是这样，那么必须增大M（若

≥

），或者减小M（若

≤

），再求出一个新的

。

重复这个过程直到

≤

为止。

这样，最优等波纹FIR滤波器就确定了。

Parks-McClellan算法在MATLAB中作为一个称为firpm的函数是可以得到的，这个函数一种的句法是

[h]=firpm（N,f,m,weights）；

设计一N阶（注意，滤波器长度是M=N+1）FIR数字滤波器，其频率响应由数组f和m给定。

数组f包含以π为单位的频带边缘频率，也即0.0

数组m包含有在f所给定频率上的期望幅度响应。

数组weights给出了再每个频带内的加权函数。

为了估计出滤波器的阶N，SP工具箱提供了函数firpmord。

这个函数也能估计在firpm函数中用到的其他参数。

这个基本语法是

[N,f0,m0,weights]=firpmord（f,m,delta）；

这个函数计算窗的阶N，在f0中的归一化频带，m0中的振幅响应，以及在weights中的频带加权值。

向量f是归一化频带边缘频率向量，m是由f定义的频带上给定期望振幅值的向量。

F的长度小于两倍m的长度，也即f不包括0或1。

向量delta给定在每个频带内的容度。

估计出的参数现在能用于firpm函数。

（程序见附录5.7）

在将滤波器设计好，频率特性求出来之后，编写程序，来画出该带阻滤波器的幅频和相频响应曲线以及滤波器单位脉冲响应曲线。

如下图所示：

图3-1滤波器频率响应曲线

图3-2滤波器单位脉冲响应

该滤波器阶数为1024。

与前面两种方法设计的滤波器相比，本方法设计的滤波器衰减特性更好，衰减更快，过渡带更窄。

该滤波器同样满足线性相位。

3.2滤波结果分析

原始信号和加噪信号和之前几种方法中所使用的相同。

在将加噪信号滤波之后，我们将滤波前后语音信号的波形及频谱图相互比较。

在同一张大图里分别绘制原始信号x，加噪信号y，滤波去噪信号y_fil的时域波形和频谱，以便比较和分析。

（程序见附录5.4）

三个信号的时域波形和频谱图如下图所示：

图3-3三个信号的时域波形和频谱图

从图中我们可以看出，原信号与滤波去噪信号的时域图基本相似；原信号与滤波去噪信号的频谱图波形也大致相似。

通过观察可以看到，加噪信号的时域图中部分被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在1900Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。

与之前不同的是，滤波信号低频部分的幅度比原始信号高了一些，也就是说加强了低频分量。

四、结论

比较窗函数法、频率取样法、等纹切比雪夫逼近法，可以看出：

（1）窗函数法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，基本可以达到阻带衰耗要求；

（2）频率抽样法特别适用于设计窄带选频滤波器，因为这时只要少数几个非0值的H（k），计算量大为降低。

但是，频率抽样法偏离设计指标最明显，阻带衰减最小，而且设计比采用窗函数法复杂，应用不及窗口法普遍。

只有适当选取过渡带样点值，才会取得较好的衰耗特性；

（3）利用等波纹切比雪夫逼近法则的设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。

五、附录

5.1音乐采集和加噪

[x,fs,bits]=wavread（'test01.wav'）;%输入参数为文件的路径和文件名，输出参数分布是样本值、采样率、编码位数

sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放

N=length（x）%计算信号x的长度

fn=1900;

t=0:

1/fs:

（N-1）/fs;%计算时间范围，样本数除以采样频率

x=x（:

1）;

x=x';

y=x+0.1*sin（fn*2*pi*t）;

sound（y,fs,bits）;%应该可以听出有尖锐的单频啸叫声

5.2原始信号和加噪信号时域和频谱图显示

X=abs（fft（x））;Y=abs（fft（y））;%对原始信号和加噪信号进行fft变换

X=X（1:

length（X）/2）;Y=Y（1:

length（Y）/2）;%截取前半部分

deltaf=fs/2/length（X）;%计算频谱的谱线间隔

f=0:

deltaf:

fs/2-deltaf;%计算频谱频率范围

figure;

subplot（2,2,1）;plot（t,x）;gridon;

xlabel（'时间'）;ylabel（'原始信号'）;

subplot（2,2,2）;plot（f,X）;gridon;

xlabel（'频率'）;ylabel（'原始信号频谱'）;

subplot（2,2,3）;plot（t,y）;gridon;

xlabel（'时间'）;ylabel（'加入的噪声信号'）;

subplot（2,2,4）;plot（f,Y）;gridon;

xlabel（'频率'）;ylabel（'加入的噪声信号频谱'）;

5.3输出滤波器的频率特性

%h_bs为滤波器的单位脉冲响应

figure;

freqz（h_bs）;

figure;

plot（h_bs）;

5.4加噪信号滤波并显示时域波形和频谱

y_fil=fftfilt（h_bs,y）;%用设计好的滤波器对y进行滤波

Y_fil=abs（fft（y_fil））;Y_fil=Y_fil（1:

length（Y_fil）/2）;%计算频谱取前一般

figure;

subplot（3,2,1）;plot（t,x）;gridon;

xlabel（'时间'）;ylabel（'时域原始信号'）;title（'时域原始信号'）;

subplot（3,2,2）;plot（f,X）;gridon;

xlabel（'频率'）;ylabel（'原始信号频谱'）;title（'原始信号频谱'）;

subplot（3,2,3）;plot（t,y）;gridon;

xlabel（'时间'）;ylabel（'加噪信号'）;title（'加噪信号'）;

subplot（3,2,4）;plot（f,Y）;gridon;

xlabel（'频率'）;ylabel（'加噪信号频谱'）;title（'加噪信号频谱'）;

subplot（3,2,5）;plot（t,y_fil）;gridon;

xlabel（'时间'）;ylabel（'滤波信号'）;title（'滤波信号'）;

subplot（3,2,6）;plot（f,Y_fil）;gridon;

xlabel（'频率'）;ylabel（'滤波信号频谱'）;title（'滤波信号频谱'）;

sound（y_fil,fs,bits）;%应该可以听到与原语音信号基本相似的语音

5.5nuttallwin窗设计滤波器

fpd=1800;fsd=1850;fsu=1950;fpu=2000;Rp=1;As=100;%带阻滤波器设计指标

fcd=（fpd+fsd）/2;fcu=（fpu+fsu）/2;

df=min（（fsd-fpd）,（fpu-fsu））;%计算上下带中心频率和频率间隔

%将Hz为单位的模拟频率换算为rad为单位的数字频率

wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;

wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;

M=ceil（15.4*pi/dw）+1;%计算nuttallwin窗设计该滤波器时需要的阶数

n=0:

M-1;%定义时间范围

w_ham=nuttallwin（M）;%产生M阶的nuttallwin窗

%调用自编函数计算理想带阻滤波器的脉冲响应

hd_bs=ideal_lp（wcd,M）+ideal_lp（pi,M）-ideal_lp（wcu,M）;

%用窗口法计算实际滤波器脉冲响应

h_bs=w_ham'.*hd_bs;

5.6频率取样法设计滤波器

fpd=1800;fsd=1850;fsu=1950;fpu=2000;Rp=1;As=100;%带阻滤波器设计指标

fcd=（fpd+fsd）/2;fcu=（fpu+fsu）/2;

df=min（（fsd-fpd）,（fpu-fsu））;%计算上下带中心频率和频率间隔

%构造带阻滤波器

M=3396;L=（M+mod（M,2））/2;

wsd=fsd/fs;wsu=fsu/fs;

F=[0:

1/L:

1];%F长度是length（L）+1

A=[ones（1,fix（wsd*M））,zeros（1,fix（wsu*M）-fix（wsd*M））,ones（1,L+1-fix（wsu*M））];

%fir2函数——频率抽样法构造带阻滤波器

B=fir2（M,F,A）;

5.7等波纹切比雪夫逼近法设计滤波器

%fs=22050;%设定采用频率，此处fs需通过采样wav音乐得到

rp1=0.5;%通带波纹

rp2=1;

rs=60;%阻带波纹

fv=[1800185019502000];%截止频率

a=[101];%期望幅度

dev=[（10^（rp1/20）-1）/（10^（rp1/20）+1）10^（-rs/20）（10^（rp2/20）-1）/（10^（rp2/20）+1）];

%FIR滤波器有B个频带时，f,a,dev是分布有2B-2,