基于频率抽样法二维非递归滤波器设计及其简单应.docx

基于频率抽样法二维非递归滤波器设计及其简单应.docx

《基于频率抽样法二维非递归滤波器设计及其简单应.docx》由会员分享，可在线阅读，更多相关《基于频率抽样法二维非递归滤波器设计及其简单应.docx（20页珍藏版）》请在冰豆网上搜索。

基于频率抽样法二维非递归滤波器设计及其简单应

一、开展本课题的意义及工作内容：

数字滤波器被广泛应用于语音、图像、无线电等领域，具有广阔的发展空间，特别是在图像处理中的应用非常广泛。

随着数字信号处理系统的大量运用及数字信号处理技术的不断发展，人们对数字滤波器的要求也越来越多样化。

本设计主要根据FIR滤波器性能要求，运用最小二乘法设计一个一维优化的FIR滤波器，然后通过一定的映射关系将一维滤波器转换成二维数字滤波器，实现对数字图像的滤波，并分析滤波器性能。

二、课题工作的总体安排及进度：

Forpersonaluseonlyinstudyandresearch;notforcommercialuse

总体安排：

2007-2008学年第一学期19-21周确定毕业设计课题；2007~2008学年第二学期1-7周在教师指导下进行做毕业论文（设计），第8周，交毕业设计（论文）成果，院组织专门人员根据规范要求对毕业设计（论文）进行形式审查。

第9周，指导教师评定本科毕业生的毕业论文（设计），毕业论文交叉评阅，学生答辩资格审定及组织学生答辩。

进度：

2007-2008学年第一学期19-21周确定毕业设计课题，在寒假期间及2007-2008学年第二学期1-2周准备好毕业设计的相关文献资料。

3-7周进行毕业设计及论文的编写。

其中在第6周基本完成预期工作，论文初稿及仿真程序提交审批。

04月22日对照审批意见对初稿进行修改并在24好最终定稿。

三、课题预期达到的效果：

设计基本达到预期的成果，成功设计了一个一维优化的FIR滤波器，然后通过一定的映射关系将一维滤波器转换成二维数字滤波器，实现对数字图像的滤波。

基于频率抽样法二维非递归滤波器设计

及其简单应用

摘要

随着数字信号处理系统的大量运用及数字信号处理技术的不断发展，人们对数字滤波器的要求也越来越多样化。

数字滤波器在数字信号处理的各种应用中发挥着十分重要的作用，它是通过对采样数据信号进行数学运算处理来达到频域滤波的目的。

数字滤波器是提取有用信息非常重要、非常灵活的方法，是现代信号处理的重要内容，它被广泛应用于语音、图像、无线电、医疗器件等领域，具有广阔的发展空间，特别是在医疗图像处理中的应用非常广泛。

本设计利用最小二乘优化算法设计一个一维FIR数字滤波器，再根据频率抽样法设计二维FIR滤波器，并在Matlab中进行仿真。

利用所设计出的滤波器对一幅X光儿童手腕骨医学图像滤波，为X光图像的后续处理作准备。

1前言

1前言

1.1课题研究的目的和意义

根据结构的不同，数字滤波器分为无限冲击响应（IIR）数字滤波器和有限冲击响应（FIR）数字滤波器。

其中，由于FIR数字滤波器总是稳定的，而且可以做到严格的线性相位，故而人们对它的设计问题进行了大量而深入的研究。

[1]

FIR数字滤波器的设计指标通常有最小二乘指标、Chebyshev指标、Lp指标及混合指标等。

因为最小二乘指标设计简单，Chebyshev指标物理意义明确，所以此两种指标最为常用。

70年代McClellan和Parks设计出基于Remez算法的计算机程序—McClellan-Parks（MP）程序—能高效的解决线性相位FIR数字滤波器数字滤波器的Chebyshev设计问题，这也使得Chebyshev误差指标在FIR滤波器的设计中被广泛的采用。

[2]

数字滤波器被广泛应用于语音、图像、无线电等领域，具有广阔的发展空间，特别是在图像处理中的应用非常广泛。

随着数字信号处理系统的大量运用及数字信号处理技术的不断发展，人们对数字滤波器的要求也越来越多样化。

比如在通讯系统中要求尽量减小阻带能量，在图像处理中要求保持信号的直流分量不变，在许多场合下都要求去掉50Hz的电源干扰信号，而为达到这些目的就需要在滤波器的设计中加入各种约束。

1.2二维数字滤波器的研究现状

滤波器的种类很多，分类方法也不同，如可以从功能上分，也可以从实现方法上分，或从设计方法上来分等等。

经典滤波器是假定输入信号x（n）中的有效信号和噪声（或干扰）信号成分各在不同的频带，当x（n）通过一个线性滤波系统后，可以将噪声信号成分有效地去除。

可是，如果有效信号和噪声信号的频率带相互重叠，那么经典的滤波器将无能为力。

现在的地质雷达信号处理中的滤波器主要采用经典的滤波器进行处理。

因此有时滤波效果较好，有时较差。

现代滤波理论研究的主要内容是从含有噪声的数据记录（又称为时间序列）中估计出信号的某些特征或信号本身。

一旦信号被估计出，那么估计出的信号将比原信号会有高的信噪比。

现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征（如自相关函数、功率谱函数等等）导出一套最佳的估值算法，然后用硬件和软件实现。

目前现代滤波器主要有：

维纳滤波器、卡尔曼滤波器、线性预测器、自适应滤波器等，很多专家将基于特征分解的频率估计及奇异值分解算法都归入现代滤波器的范畴。

[2]

1.3本课题研究的内容

本设计主要根据FIR滤波器性能要求，运用最小二乘法设计一个一维优化的FIR滤波器，然后通过一定的映射关系将一维滤波器转换成二维数字滤波器，实现对数字图像的滤波，并分析滤波器性能。

2FIR滤波器设计技术及滤波特性

2.1传统FIR滤波器设计方法

FIR滤波器设计最常用的方法是窗函数设计法，其它常用的方法有频率抽样法，频率变换法，切比雪夫逼近法等，以下对这些常用方法作简单介绍。

窗函数法又称傅里叶级数法，是最常用的设计方法，它以hd（n）为媒介的时域设计法,而滤波器指标往往是在频域给出的,为此,要由Hd（ejω）算出hd（n）,加窗后又从h（n）算出hd（n）来检验。

窗函数设计方法首先根据要求选择一个适当的理想滤波器,由于理想滤波器的脉冲响应是非因果且无限长,用适当的窗函数来截取它的脉冲响应,从而得到线性相位和因果的FIR滤波器。

频率抽样法先对理想频响Hd（ejω）采样,得到样值H（k）,再利用插值公式直接求出系统转换函数H（z）；或者求出频响H（ejω）,以便与理想频响进行比较。

在[0,2π]区间上对Hd（ejω）进行N点采样,等效于时域以N为周期延拓。

频率抽样法可以看做为插值法,它在采样的ωi上保证

等于

；而在非插值点（采样点）上,H（ejω）是插值函数的线性叠加。

为了提高逼近的质量，减少逼近误差，可以采用人为的扩展过渡带的方法，即在频率相应的过渡带内插入一个或多个比较连续的采样点，使得过渡带比较连续，从而使得通带和阻带之间变法比较缓慢，使得设计得到的滤波器对理想滤波器的逼近误差较小[3]。

切比雪夫方法从数值逼近的理论来看,对某个函数f（x）的逼近一般有三个方法[10,11]：

插植法；最小平方逼近法；最佳一致逼近法。

所谓插值,即寻找一n阶多项式（或三角多项式）p（x）,使它在n+1个点x0,x1…xn处满足p（xk）=f（xk）,k=0,1,…,n而在非插值点上,p（x）是f（xk）的某种组合。

当然,在非插植点上,p（x）和f（xk）存在一定的误差。

频率抽样法可以看作为插值法,它在抽样点ωk上保证了

=

而在非抽样点上,H（ejω）是插值函数S（ω,k）的线性组合,其权重是

。

这种设计方法的缺点是通带和阻带的边缘不容易精确的确定。

这种设计方法是着眼于使整个区间[a,b]的总误差为最小,但它并不一定能保证在每个局部位置误差都最小。

在某些位置上,有可能存在着较大的误差。

实际上,傅立叶级数法就是一种最小平方逼近法。

该方法在间断点处出现了较大的过冲（Gibbs现象）。

为了减少这种过冲和欠冲,采用了加窗口的方法,当然,加窗口以后的设计方法,已不在是最小平方逼近。

2.2改进的频率抽样法

频率抽样法的特点是:

（1）在采样频率上的逼近误差为零，也就是理想和实际响应的差为零；

（2）其余频率上的逼近误差取决于理想响应的形状；理想响应的轮廓越陡，则逼近误差越大；

（3）靠近带的边缘的误差大，在带内的误差小。

频率抽样法目前有两种设计方法，第一种直接用上面的基本思想，对逼近误差不加任何限制；也就是说无论设计所得的误差有多大我们都接受，这种方法叫朴素设计法。

第二种方法则通过改变过渡带的样本值，努力使阻带中的误差极小化，以便产生一个较好的设计，这种方法叫最优设计法。

频率抽样法设计滤波器最大的优点是直接从频率域进行设计，比较直观，也适合于设计具有任意幅度特性的滤波器。

缺点是边缘频率不易控制。

如果增加采样点数N，对确定边缘频率有好处，但N加大会增加滤波器的成本。

因此它适合于窄带滤波器的设计。

MATLAB信号处理工具箱提供了一个频率抽样法的设计函数fir2，它的典型调用方法为:

h=fir2（M，f,A）。

其中M是FIR滤波器的阶数（滤波器的长度为N=M+1），长度为N的数组h为滤波器系数（或脉冲响应）。

数组f中包含各边缘频率，其单位为pi,0≦f≦1。

f=1对应于采样频率的一半，即奈奎斯特频率。

这些频率必须以递增顺序排列，从0开始，到1结束。

数组A为各指定频率上预期的幅度响应，f与A长度必须相等，plot（f,A）应该给出预期的滤波器幅频特性。

[4]

现在利用频率抽样直接法设计一FIR滤波器，低通滤波器技术指标为

ωp=0.2π,Rp=0.25dB

ωs=0.3π,As=50dB

选M=20，以使在ωp有一个频率样本，也即在k=2：

ωp=0.2π；下一个样本在ωs，也即在k=3：

ωs=0.3π。

这样在通带[0≤ω≤ωp]内有3个样本，在阻带[ωs≤ω≤π]内有7个样本。

（2-1）

由于M=20，

，并且这是一个

类线性相位滤波器，有

（2-2）

现在根据

（2-3）

将H（k）集合起来，并由

确定脉冲响应h（n），用MATLAB实现，得到图2.1。

图2.1直接频率抽样设计法

由图2.1可见，最小阻带衰减大约是16dB，这显然是不可接受的。

如果增加M，那么在过滤带内就一定有一些样本，而对这些我们又并不完全知道频率响应。

因此，在实际中很少采用直接设计法。

为了得到更大的衰减，就必须增大M，并让过渡带的样本作为自由样本；也就是说，改变它们的值以得到在已给定M下的最大衰减及其过渡带宽。

这个问题被认为是一个优化问题，可以用线性规划技术来解决。

现在我们用最优设计法来对上面的低通滤波器。

现选M=40以使有一个样本在过滤带

内。

过滤带内的样本是在k=5和在k=40-5=35。

现用T1表示这些样本值，且

，那么已采样的振幅响应是

（2-4）

由于

，相位响应的样本是

（2-5）

现在能够变化T1以得到最好的最小阻带衰减。

这将会引起过渡带加宽。

首先看看当T1＝0.5时得到图2.2。

图2.2　最优频率设计法：

T1=0.5

由图2.2可见，这个设计的最小阻带衰减现在是30dB，这就是比直接设计的衰减要好一些，但是仍然不在可以接受的50dB的水平上。

通过人为地改变T1值可以示得最佳T1值，并且找到接近最优解的值T1=0.39。

由图2.3可见，最优阻带衰减是43dB。

明显地要进一步增加衰减就必须在过渡带内改变更多的样本值。

图2.3最优频率设计法：

T1=0.39

很清楚，通过变化一个样本就能得到一种好得多的设计，从这点来看这种方法是很优越的。

实际上，过渡带宽一般是很小的，就只包括一个或两个样本。

因此需要优化最多也就是两个样本以获得最大的最小阻带衰减。

在绝对的意义上这也就是等效于使最大旁瓣幅度最小化，所以这类优化问题也称为最大最小化问题。

现在重新考虑上面的低通滤波器设计，用在过渡带内的两个样本来解决它，以求得到一个更好的阻带衰减。

先M=60使得在过渡带内有两上样本。

设这两个过渡带的样本值是T1和T2，那么，

给出为

（2-6）

利用T1=0.5925和T2=0.1099这两个值，可用MATLAB计算h（n），得图2.4。

图2.4　低通滤波器设计图

由图2.4可见，最小阻带衰减现在是大约是63dB,这就是可以接受的了。

[5]

3基于频率抽样法二维FIR滤波器的设计及图像滤波实现

3.1设计一维FIR滤波器

按频域采样定理，FIR数字滤波器的传输函数H（z）和单位脉冲响应h（n）可由它的N个频域采样值H（k）惟一确定，即

　　　　　　（3-1）

（3-2）

假定要设计的FIR数字滤波器频率响应为

，直接从频域出发，对理想频率响应在0~2π之间进行N点的等间隔休样，得到

为

（3-3）

取H（k）为理想滤波器在各频率采样点上的值，即令

（3-4）

然后按频域采样定理，由H（k）求H（z）或h（n）。

为了使设计的滤波器为线性相位滤波器，滤波器的频率响应

及频率采样值

必须满中一定的条件，同时，在有些情况下，给定的理想滤波器的幅度响应只给出了（0~π）范转频率采样的幅度值，因而上式中的k只能取0到N/2这（N/2+1）点，另外（N/2-1）点（从N/2+1到N-1）则需根据线性相位的奇偶对称条件求得。

对N为奇数，设计公式为

（3-5）

（3-6）

（3-7）

对N为偶数，设计公式为

（3-8）

（3-9）

（3-10）

（3-11）

3.2最小二乘优化设计

单位脉冲响应为{h（n1,n2）,n1=0,1,…,N1-1,n2=0,1,…,N2-1}的二维FIR滤波器,其频率响应为

（3-12）

二维FIR滤波器的约束最小二乘设计是指选择N1×N2个脉冲响应系数h（n1,n2）,使H（ejω1,ejω2）逼近某个期望响应,在逼近误差小于预定值的条件下,使平方逼近误差的积分最小.

一般情形下,N1×N2个脉冲响应是独立的。

当满足一定对称性时,独立脉冲响应系数有所减少,且滤波器具有大多数应用场合所要求的线性相位。

为描述简单,考虑N1,N2均为奇数、脉冲响应矩形对称的情形,即

h（n1,n2）=h（N1-1-n1,n2）=h（n1,N2-1-n2）=h（N1-1-n1,N2-1-n2）………………………………………………………………（3-13）

此时,共有数目为n=（N1+1）（N2+1）/4的独立脉冲响应,即

h（n1,n2）,n1=0,1,…,r1,n2=0,1,…，r2,其中r1=（N1-1）/2,r2=（N2-1）/2；频率响应

变为

（3-14）

其中,相频响应为-（r1ω1+r2ω2）,是频率ω1,ω2的线性函数（线性相位）,幅频响应P（ω1,ω2）等于

（3-15）

是独立脉冲响应的线性函数。

用{αi|i=1,2,…，n}来表示独立的脉冲响应,则幅频响应可表示为：

（3-16）

其中

是w1,w2的某个函数。

当αi=h（n1,n2）,0≤n1

（3-17）

当αi=h（n1,r2）,0≤n1

（3-18）

当αi=h（r1,n2）,0≤n2

当αi=h（r1,r2）时,

=1。

令

=[

，

，…，

]T，α=[α1,α2,…,αn]T,则幅频响应可进一步化为

（3-19）

当所设计的滤波器就是一个线性相位滤波器时,只需考虑幅频响应的逼近问题,问题简化为:

使P（ω1,ω2）按如下约束最小二乘准则逼近期望幅频响应D（ω1,ω2）：

（3-20）

（3-21）

其中Ω1及Ω2是频率区域{（ω1,ω2）-p≤ω1≤p,-p≤ω2≤p}的子集,ε（ω1,ω2）为预定的逼近误差上限函数。

容易看到,P（ω1,ω2）关于ω1,ω2矩形对称,因此要求D（ω1,ω2）也具有矩形对称性,并且只需考虑频率子区域Π≡{（ω1,ω2）0≤ω1≤p,0≤ω2≤p}上的逼近问题。

此时,逼近问题仍然表示为

（1）、

（2）,但其中的Ω1和Ω2是频率子区域Π的子集。

一般情况下,式（4-16）的积分需采用数值方法计算,式（4-17）的约束用稠密频率点集上的约束集代替。

为此,我们将频率子区域Π离散化为稠密的频率点集Ω={（ω1i,ω2j）=（ip/L,jp/L）|i,j=0,1,…,L}L是一个充分大的正整数。

将式（4-16）积分替换为Ω1∩Ω上的求和,将式（4-17）约束替换为Ω2∩Ω上的约束集。

二维线性相位FIR滤波器的约束最小二乘设计问题可化为:

[6]

（3-22）

（3-23）

3.3基于频率抽样法二维FIR滤波器的程序实现

Matlab信号处理工具箱中的firls函数用来设计具有最小二乘线性相位FIR滤波器。

该滤波器可使指定频段内的理想分段线性函数与滤波器频响应之间的误差平方和最小[7]。

b=firls（n,f,m）可得到n阶FIR滤波器，其幅频特性与由f和m指定的性能匹配。

滤波器系数矢量b共包含有n+1个系数，且满足对称关系

b（k）=b（n+2-k）k=1，…，n+1

对应滤波器即为Ⅰ型和Ⅱ型线性相位滤波器。

b=firls（n,f,m,w）可利用权值矢量w对各频率段进行加权拟合，w的长度为f和m长度的一半，用以指定各频率段的权值。

b=firls（n,f,m,`ftype`）和b=firls（f,m,w,`ftype`）可指定滤波器的类型：

当ftype=hilbert时，设计的滤波器为奇对称的线性相位滤波器（Ⅲ型和Ⅳ型），b的系数满足以b（k）=-b（n+2-k）,k=1,…,n+1。

这类滤波器包括Hilbert变换器。

当ftype=differentiator时，也可设计出Ⅲ型和Ⅳ型滤波器，但采用了特殊的加权技术，对非零幅值的频段，将误差乘以（1+f）2，这样可使低频段误差大大小于高频段误差。

[8]

fsamp2函数是使用频率抽样法设计二维FIR滤波器。

h=fsamp2（hd）,设计一个频率响应为hd的二维FIR滤波器，返回滤波器系数到矩阵h,滤波器h有一个频率响应经过点hd。

如果hd是在M、N之间，那么H也在M、N之间。

hd是一个包含理想频率响应采样点的矩阵，采样点在x和y频率轴-1.0和1.0之间的等间隔点上，1.0对应采样频率的1/2,或pi弧度。

为了结果的准确,使用FREQSPACE返回频率点得到hd。

h=fsamp2（f1,f2,hd,[MN]）产生一个从M到N的FIR滤波器，使滤波器响应hd在f1和f2之间。

得到的滤波器尽可能接近理想的频率响应。

最好的结果是必须有至少M*N个理想的频率点。

如果您指定的点少于M*N点,fsamp2将出现警告。

[3]

下面是利用该函数和二维FIR滤滤器的频率采样设计函数fsamp2进行设计：

设计要求：

利用频率抽样法设计一20阶的具有最小二乘线性相位的FIR带阻滤波器。

首先必须给出一维的相应FIR滤波器。

f=[00.10.20.30.40.60.70.80.91];

m=[1110000111]；

产生一个20阶的具有最小二乘线性相位的一维FIR带阻滤波器。

B=firls（20,f,m,`Hilbert`）;

[h0,w]=freqz（b,1,512,2）;

Plot（f,m,`--`,w,abs（h0）;

在一维滤波器函数取理想带阻滤波器时，采用频率抽样法产生一个20阶的具有最小二乘线性相位的二维FIR带阻滤波器h1。

h1=fsamp2（f,f,m,[55]）;

figure

（2）

freqz2（h1,[3232]）;

图3.1一维频率FIR响应

图3.1中虚线代表理想FIR频率响应曲线，实线代表实际的FIR频率响应曲线。

在一维波器函数取实际带阻滤波器时，采用频率抽样法产生一个20阶的具有最小二乘线性相位的二维FIR带阻滤波器h2。

ff=1/512:

1/512:

1;

h2=fsamp2（ff,ff,h0,[55]）;

figure（3）

freqz2（h2,[3232]）;

图3.2理想二维FIR频率响应

图3.3实际二维FIR频率响应

3.4基于频率抽样法图像滤波的实现

数字滤波器被广泛应用于语音、图像、无线电等领域，具有广阔的发展空间，特别是在图像处理中的应用非常广泛，因为本人所学知识有限，只能利用滤波器对图象进行一些简单的处理，但我坚信应用滤波器处理图象将是未来图象处理的一大趋势并占主导地位。

在本设计中,将对自己设计的滤波器进行一些简单的应用。

在现实生活当中滤波器越来越广泛地被应用于图象处理领域，因此，在本次设计中将把滤波思想应用于X线片图像处理，旨在通过FIR滤波器快速高效的对X线片图像进行滤波处理，使处理后的图片更具立体感。

图3.5图像滤波流程图

利用上面设计的滤波器进行图象分析，导入原始图像，对于彩色图片先要进行灰度化[10],图像滤波的流图如图3.5所示。

图3.4原始图像

f=imread（`lena256.bmp`）;

利用上面设计的FIR滤波器h1进行滤波

fh1=filter2（h1,f）;

figure（4）

imshow（fh1）;

利用上面设计的FIR滤波器h2进行滤波

fh2=filter2（h2,f）;

figure（5）

imshow（fh2）;

3.5结果及展望

3.5.1设计结果及所取得成绩

图3.6理想二维FIR滤波结果图3.7实际二维FIR滤波结果

图像滤波结果如图3.6，图3.7所示。

所完成的工作：

（1）研究了数字滤波的理论知识，为系统整体奠定了理论基础；

（2）了解FIR滤波器的一些设计方法，利用频率抽样法设计具有最小二乘线性相位的FIR带阻滤波器；（3）利用设计出的滤波器实现图像滤波；（4）研究了MATLAB软件在数字信号处理，尤其是数字滤波器处理中的应用。

3.5.2研究展望

在对医学知识的深入了解下,上述滤波器的快速处理X线片图象的特点可用于快速高效的去除人体手腕部X线片的噪声，从而方便医学人员准确评价骨龄。

下面将简要介绍骨龄评价，以便更好的理解这一应用背景。

人的骨骼年龄简称为骨龄，是一种与生活中人们常说的“时间年龄”不一致的、能准确反映人体生长发育信息的生物学年龄。

骨龄是出生后大多数正常儿童随年龄增长而出现的有规律的骨骺线解剖变化标志，是由于个体间营养发育的不同以及地区、种族、性别的不同而表现出的与时间年龄的不一致。

另外，通过研究大量X线片样本，比较各地区儿童或成年人骨骼年龄，找出其成长发育和身体状况规律，对于医学研究意义更大。

通过本次设计，自己从中取得了一些成绩，理论水平也得到了一定的提高，同时也暴露了一些问题。

首先，对一个课题必须要阅读大量的文献和书籍来获得一定的感性认识，然后才能有自己的想法，这是一条必经之路。

其次，理论基础知识很重要，论文设计了很多的算法，会用到很多基础知识，如果用的时候再去学会浪费时间。

最后，要有信心，遇到困难要向别人请教，这样可以大大加快研究进程。

以上是我做论文的一些心得体会，这些对我以后的工作会有很大的帮助。

附录A程序代码

%直接法设计低通滤波器

M=20;alpha=（M-1）/2;l=0:

M-1;wl=（2*pi/M）*l;

Hrs=[1,1,1,zeros（1,15）,1,1];

Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];

k1=0:

floor（（M-1）/2）;k2=floor（（M-1）/2）+1:

M-1;

angH=[-alpha*（2*pi）/M*k1,alpha*（2*pi）/M*（M-k2）];

H=Hrs.*exp（j*angH）;

h=real（iff