采用频率采样法的FIR滤波器

吹管乐滤波去噪——基于频率采样法的FIR滤波器学生姓名：焦阳指导老师：胡双红摘要本课程设计主要内容是设计利用频率采样法设计一个FIR滤波器，对一段吹管乐进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。

本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。

首先在网上找到一段笛子独奏，加入一单频噪声，对信号进行频谱分析以确定所加噪声频率，设计滤波器进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。

由分析结果可知，滤波器后的音频信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。

关键词滤波去噪；FIR滤波器；频率采样法；MATLAB1 引言滤波去噪[1]是信号处理中一种非常基本但十分重要的技术。

利用滤波可以从复杂的信号中提取所需的信号，一直不需要的信号。

滤波器就是这样一种可以在时域和频域对信号进行滤波处理的系统。

通常情况下，有用信号和干扰信号是在不同频段上的，于是通过对滤波器的频率特性精心设计就能达到滤波的目的。

本课程设计是采用频率采样法设计频率抽样型滤波器，从而对吹管乐信号滤波去噪。

通过对比滤波前后的波形图及回放滤波前后的吹管乐信号，来判断滤波器对噪声信号确实有滤除作用。

1.1 课程设计目的（1）熟悉使用MATLAB；（2）了解FIR滤波器原理及结构；（3）利用所学数字信号处理想干知识用MATLAB设计一个FIR滤波器；（4）提高自己动手能力；（5）对加噪声的语音信号进行滤波去噪处理，比较滤波前后的时域波形和频谱并进行分析；1.2 课程设计要求（1）滤波器指标必须符合工程设计；（2）设计完后应检查其频率响应曲线是否满足指标；（3）处理结果和分析结论应该一致，而且应符合理论；（4）独立完成课程设计并按要求编写课程设计报告；1.3 设计平台本课程设计仿真平台为MATLAB7.0。

MATLAB的名称源自Matrix Laboratory，1984年由美工Mathworks公司推向市场。

它是一种科学计算软件，专门以矩阵的形式处理数据。

MATLAB将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信号处理等领域的分许、仿真和设计工作。

1993年MathWorks公司从加拿大滑铁卢大学购得MAPLE软件的使用权，从而以MAPLE为“引擎”开发了符号数学工具箱（Symbolic Math Toolbox）[2]。

2 设计原理用网上找一段吹管乐，绘制波形并且观察其频谱，给定相应技术指标，用频率采样法设计的一个满足指标的频率采样型FIR滤波器，对该信号进行滤波去噪处理，比较滤波前后的波形和频谱进行分析。

2.1 FIR滤波器的设计FIR（Finite Implse Response）[3]滤波器：有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，他可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。

因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

有限长单位冲激响应（FIR）滤波器有以下特点：（1）系统的单位冲激响应h（n）在有限个n值处不为0；（2）系统函数H（z）在|z|>0处收敛，极点全部在z=0处（因果系统）；（3）结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。

2.2 频率采样型结构把一个有限长序列（长度为N点）的z变换H（z）在单位圆上作N等分抽样，就得到H （k ），其主值序列就等于h(n)的离散傅里叶变换H （k ）。

那里也说到用H （k ）表示的H （z ）的内插公式为----==--∑1101()()(1)1N Nk k N H k H z z N W z （2-1） 这个公式就为FIR 滤波器提供了另外一种结构，这种结构由两部分级联组成。

-==∑11()()()N c k k H z H z H z N (2-2) 其中级联的第一部分为梳状滤波器，其结构如图所示：-=-()(1)N c H z z （2-3）-Z -N 图2-1 梳状滤波器结构图第二部分由N 各谐振器组成的谐振柜。

它是由N 个一阶网络并联而成，而这每一个一阶网络都是一个谐振器--=-1()()1K K N H k H z W z （2-4）其结构如下图所示：H （k ） H k (z)W -A图2-2 一阶谐振器频率抽样型结构特点：（1）它的系数H （k ）直接就是在滤波器在π=2k w k N处得频率响应。

因此，控制 得频率响应是很直接得。

（2）结构有两个主要缺点：a ．所有的相乘系数及H （k ）都是复数，应将它们先化成二阶实数，这样乘起来较复杂，增加乘法次数，存储量。

b ．所有谐振器的极点都是在单位圆上，由-k N w 决定考虑到系数量化的影响，当系数量化时，极点会移动，有些极点就不能被梳状滤波器的零点所抵消。

（零点由延时单元决定，不受量化的影响）系统就不稳定了。

（3）将一阶网络合并为二阶网络：a ．第k 和第N-k 个谐振器合并为一个实系数的二阶网络，因为h(n)是实数，他的DFT 也是圆周共轭对称的。

=-*()()H k H N k =⋯-1,2,3,1k N （2-5） 因此，可以将第k 和第N-k 个谐振器合并为一个二阶网络。

----------≈+=+----**1()111()()()()1111k k N k k K N N N N H k H N k H k H k H rW z rW z rW z rW z ---------β+β==π-++-+**101122122()21[]12cos()k k K k k k N N N N z H k k z W W r r W W z z r r z N（2-6） 其中：β=2β01k k k N Ｒe[H(k)],=-2r Ｒe[H(k)Ｗ]b ．第k 和第N-k 个谐振器合并为一个二阶网络的极点在单位圆内，而不是在单位圆上，因而从频率响应的几何解释可知，它相当于一个有限Q 的谐振器。

其谐振频率为π=2k w k N。

图2-3 二阶网络结构图除了共轭复根外，还有实根。

当N=偶数时，有一对实根，它们分别是k=0，k=N/2两个点。

-=-01(0)()1H H z rz 和-=+12()2()1N N H H z rz （2-7） 当N=奇数时，只有一个实根z=r(k=0)，即只有H0(z)。

c ．修正频率抽样结构流图（N=偶数）图2-4 修正频率抽样结构流图（N=偶数） -------=β+β=-⨯++π-+-+∑11201111221()1(0)2()(1)[]21112cos()NN N k k k N H z H H z r z N rz rz z r k r z N （2-8） 修正频率抽样结构流图（N=奇数）图2-5 修正频率抽样结构流图（N=奇数） -------=β+β=-++π-+-+∑11201111221()1(0)2()(1)[]21112cos()N N N k k k N H z H H z r z N rz rz z r k r z N（2-9）2.3 频率采样法2.3.1 设计思路：这种设计方法是从频域进行设计的一种方法，首先给定一个希望逼近的频率响应。

πω==|2()()jwd k N H k He =⋯-1,2,3,1k N （2-10）知道H （k ）后，由IDFT 定义，可以用这N 个采样值H （k ）来唯一确定有限长序列h(n)，即 --==∑101()()N nk N k h n H k W N =⋯-0,1,2,1n N （2-11） --==∑10()()N n n H z h n z （2-12）-ωω==∑10()()N j j n k H e h n e （2-13） 内插公式： -ω-ωΦ(ω=ω(1)/2sin(/2))sin(/2)j N N e N （2-15） 四种线性相位的FIR 滤波器如下表2-1所示：表2-1 四种线性相位的FIR 滤波器2.3.2 逼近误差及其改进措施：（1）采样点上滤波器的实际响应是严格地和理想频率响应数值相等的；（2）在采样点之间的频率响应则是由各采样点的加权内插函数的延伸叠加而成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状。

（3）理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小；（4）如果采样点之间的理想频率特性变化越陡，则内插值与理想值的误差就越大，因而在理想频率特性的不连续点附近，就会产生尖峰和起伏。

2.3.3 滤波器性能的改善：（1）增加过滤带采样点，它可以大大减少震荡，阻带衰减也可以得到进一步改善。

一般一点到二点的过滤带采样即可得到满意的结果。

（2）增加采样点密度，过渡带的宽度与采样点数N 成反比。

但N 值意味着或长度的增加，滤波器运算量必然增大[4]。

3 设计步骤3.1 设计流程图图3-1 流程图3.2 音频信号获取在网上下载一段吹管乐，并截取适当片段。

把音频放到电脑录音机里面，把音频属性设置为8000Hz。

如图3-2所示：图3-2音频信号设置然后在MATLAB软件平台下，利用函数wavread对音频信号进行采样，源程序为：[x,fs,bits]=wavread(‘梅花三弄.wav’)，记住采样频率和采样点数，MATLAB实现得：fs=8000;bits=8。

3.3 音频信号的频谱分析[x,fs,bits]=wavread('e:\梅花三弄.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。

sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放n=length(x); % 计算信号x的长度fn=3200; % 单频噪声频率，此参数可改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可以看出，在频域为3200Hz处加入一个单频噪声，而加入噪声之后，时域的波形图出现了明显失真，通过听取原音频信号x和加噪声音频信号y，可以明显听到y 音频信号中有一明显尖锐噪声。