优秀毕业设计—维纳自适应滤波器设计及Matlab实现

优秀毕业 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 —维纳自适应滤波器设计及Matlab实现 2005 Matlab 2001 2001 级 电子信息 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 2004年05月24日 I II 本文从随机噪声的特性出发，分析了传统滤波和自适应滤波基本工作原理和 性能，以及滤波技术的现状和发展前景。然后系统阐述了基本维纳滤波原理和自 适应滤波器的基本结构模型，接着在此基础上结合最陡下降法引出LMS算法。在MSE准则下，设计了一个定长的自适应最小均方横向滤波器，并通过MATLAB编程实现。接着用图像复原来验证该滤波器的性能，结果表明图像的质量在MSE准则下得到了明显的改善。最后分析比较了自适应LMS滤波和频域维纳递归滤波 之间的性能。本文还对MATLAB里面的自适应维纳滤波函数wiener2进行了简单分析。 ：退化图像 维纳滤波 自适应滤波 最陡下降法 LMS Abstract This paper analyses the basic work theory, performance of traditional filter and adaptive filter based on the property of random noise, and introduce the status quo and the foreground of filter technology. Then we explain basic theory of wiener filter and basic structure model of adaptive filter, and combine the method of steepest descent to deduce the LMS. Afterward according to the MSE rule, we design a limited length transversal filter, and implement by MATLAB. And then we validate performance of adaptive LMS filter by restoring images, Test result show that the quality of the degrade images were improved under the rule of MSE. Finally, we compare the performance of adaptive LMS filter and iterative wiener filter. We also simply analyses the wiener2 () which is a adaptive filter in MATLAB. Keywords: degrade image；wiener filter；adaptive filter；ADF；LMS algorithm 1 绪论…………………………………………………………………………………1 1. 1 引言…………………………………………………………………………1 1. 2 研究目标及现状……………………………………………………………1 1. 2 .1 图像复原技术的目标……………………………………………1 1. 2 .2 图像复原技术的研究现状………………………………………1 2 理论基础 …………………………………………………………………………3 2. 1 基本自适应滤波器的模块结构……………………………………………3 2. 2 基本维纳滤波原理…………………………………………………………4 3 自适应滤波原理及算法 ………………………………………………………6 3.1 横向滤波结构的最陡下降算法……………………………………………7 3.1.1 最陡下降算法的原理……………………………………………7 3.1.2 最陡下降算法稳定性……………………………………………10 3.2 LMS滤波原理及算法……………………………………………………11 3.2.1 从最陡下降算法导出LMS算法 ………………………………11 3.2.2 基本LMS算法的实现步骤 ……………………………………11 3.2.3 基本LMS算法的实现流程图 …………………………………12 3.2.4 LMS算法的Matlab实现 ………………………………………12 3.2.5 wiener2()的原理 ……………………………………………12 3.2.6 LMS性能分析——自适应收敛性……………………………13 4 Matlab 实验结果 …………………………………………………………14 4.1.LMS滤波器的收敛性 ………………………………………………14 4.2.LMS滤波器和频域迭代维纳滤波器的性能比较 ……………………16 5 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf ………………………………………………………………………………18 致谢 …………………………………………………………………………………19 参考文献 ……………………………………………………………………………20 附录A ………………………………………………………………………………21 附录B ………………………………………………………………………………22 附录C ………………………………………………………………………………27 1 1.1 人类传递信息的主要媒介是语言和图像。据统计，在人类接受的信息中，听 觉信息占20％，视觉信息占60％，其它如味觉、触觉、嗅觉总的加起来不过占 [1]20％，所以图像信息是十分重要的信息。然而，在图像的获取和图像信号的传 输过程中，图像信号中不可避免的混入各种各样的随机噪声，造成图像失真（图 像退化）。造成人类所获取的信息和实际是有偏差的，成为人类从外界获取准确 信息的障碍。因此，对图像信号中的随机噪声的抑制处理是图像处理中非常重要 的一项工作。 在图像的获取和传输过程中所混入的噪声，主要来源于通信系统中的各种各 样的噪声，根据通信原理及统计方面的知识，可以知道在通信系统中所遇到的信 号和噪声，大多数均可视为平稳的随机过程[15]。又有“高斯过程又称正态随机过 程，它是一种普遍存在和重要的随机过程，在通信信道中的噪声，通常是一种高 斯过程，故又称高斯噪声。因此，在大多数的情况下，我们可以把造成图像失真 的噪声可视为广义平稳高斯过程。 本文针对图像信号中混入的随机噪声，在怎样把现有的滤波算法应用到实际 的图像复原中去的问题上提出了解决方法，并且应用Matlab 软件编程对图像进 行处理。 1.2 1.2.1 为了从含有噪声的数据中提取我们所感兴趣的、接近规定质量的图像，我们 需要设计一个系统满足：当信号与噪声同时输入时，在输出端能将信号尽可能精 确地重现出来，而噪声却受到最大抑制，即最佳滤波器。 1.2.2 目前的图像复原技术，即去噪的滤波技术可以分为两大类：传统滤波和现代 滤波。传统滤波技术是建立在已知有用信号和干扰噪声的统计特性（自相关函数 或功率谱）的基础上的噪声去除；现代滤波技术则是不需要知道图像的先验知识， 只是根据观测数据，即可对噪声进行有效滤除。 早在20世纪40年代，就对平稳随机信号建立了维纳滤波理论。根据有用信 号和干扰噪声的统计特性（自相关函数或功率谱），以线性最小均方误差（MSE） 估计准则所设计的最佳滤波器，称为维纳滤波器。这种滤波器能最大程度的滤除 干扰噪声，提取有用信号。但是，当输入信号的统计特性偏离设计条件，则它就 不再是最佳的了，这在实际应用中受到了限制。到60年代初，由于空间技术的 发展，出现了卡尔曼滤波理论，即利用状态变量模型对非平稳、多输入多输出随 机序列作最优估计。卡尔曼滤波器既可以对平稳的和平稳的随机信号作线性最佳 滤波，也可以作为非线性滤波[2]。 然而只有在对信号和噪声的统计特性已知的情况下，这两种滤波器才能获得 最优解。在实际的应用中，往往无法得到这些统计特性的先验知识，或者统计特 性是随时间变化的，因此，这两种滤波器就实现不了真正的最佳滤波。 Widrow B.和Hoff于1967年提出的自适应滤波理论，可使在设计自适应滤 波器时不需要事先知道关于输入信号和噪声的统计特性的知识，它能够在自己的 工作过程中逐渐估计出所需的统计特性，并以此为依据自动调整自己的参数，以 达到最佳滤波效果。一旦输入信号的统计特性发生变化，它又能够跟踪这种变化， 自动调整参数，使滤波器性能重新达到最佳。 自适应滤波器自动调节参数可以通过各种不同的递推算法来实现，由于它采 用的是逼近的算法，使得实际估计值和理论值之间必然存在差距，也就造成了自 适应滤波问题没有唯一的解。依照各种递推算法的特点，我们把它应用于不同的 场合。现在广为应用的自适应滤波方法主要是基于以下几种基本理论，再融合递 推算法导出来的： （1） 基于维纳滤波理论的方法 维纳滤波是在最小均方误差准则下通过求解维纳—霍夫方程来解决线性最 优滤波问题的。基于维纳滤波原理，我们利用相关的瞬时值通过在工作过程中的 逐步调整参数逼近信号的统计特性，实现最优滤波。由此，我们得到一种最常用 的算法——最小均方算法，简称LMS算法。 （2） 基于卡尔曼滤波理论的方法 卡尔曼滤波是线性无偏最小方差滤波递推滤波，它能使滤波器工作在平稳的 或非平稳的环境，得到最优解。利用卡尔曼滤波理论的递推求解法导出自适应滤 波器更新权矢量得不同递推算法。比LMS算法有极快的收敛速率，可是计算复 杂度也增大了，它需要计算卡尔曼矩阵。 （3） 基于最小二乘准则的方法 维纳滤波和卡尔曼滤波推导的算法是基于统计概念的，而最小二乘估计算法 是以最小误差平方和为优化目标的。根据滤波器的实现结构，有以下3种不同的最小二乘自适应滤波算法：自适应递归最小二乘法（RLS），自适应最小二乘格型算法，QR分解最小二乘算法。 （4） 基于神经网络理论的方法 神经网络是有大量的神经元相互连接而成的网络系统，实质上它是一个高度 非线性的动力学网络系统，这个系统具有很强的自适应、自学习、自组织能力， 以及巨量并行性、容错性和坚韧性，因而，它可以做很多传统的信号和信息处理 技术所不能做的事情。因其超强的自动调节能力，使得它在自适应信号处理方面 有着广阔的前景[2]。 在一系列的自适应算法中，虽然基于后面3种基本理论的方法在收敛速率和 稳定、坚韧性方面有着更好的性能，但是， 基于维纳滤波理论的LMS算法因其算法简单，而且能达到满意的性能，得到了青睐，成为了应用最广泛的自适应算 法。 为此，本文主要研究LMS自适应滤波器在图像去噪方面的应用。 2. 2.1 自适应滤波器通常由两部分构成，其一是滤波子系统，根据它所要处理的功 能而往往有不同的结构形式。另一是自适应算法部分，用来调整滤波子系统结构 的参数，或滤波系数。在自适应调整滤波系数的过程中，有不同的准则和算法。 算法是指调整自适应滤波系数的步骤，以达到在所描述的准则下的误差最小化。 自适应滤波器含有两个过程，即自适应过程和滤波过程。前一过程的基本目标是 调节滤波系数w(k),(.)i，使得有意义的目标函数或代价函数最小化，滤波器输 y(k)e(k)d(k)出信号逐步逼近所期望的参考信号，由两者之间的误差信号驱动某 种算法对滤波系数进行调整，使得滤波器处于最佳工作状态以实现滤波过程。所 以自适应过程是一个闭合的反馈环，算法决定了这个闭合环路的自适应过程所需 x(k)d(k),(.)要的时间。但是，由于目标函数是输入信号,参考信号及输出信号 (.),,[x(k),,d(k),y(k)] y(k)的函数，即,因此目标函数必须具有以下两个性质： （1） 非负性 ,(.),,[x(k),d(k),y(k)] ,0 ,,x(k),d(k),y(k) ( 2.1 ) （2） 最佳性 ,(.),,[x(k),d(k),y(k)] , 0 , when y(k),d(k) ( 2.2 ) ,(.)y(k)在自适应过程中，自适应算法逐步使目标函数最小化，最终使逼近 www(k)d(k)optopti于，滤波参数或权系数收敛于，这里是自适应滤波系数的最优 解即维纳解。因此，自适应过程也是自适应滤波器的最佳线性估计的过程，既要 d(k)估计滤波器能实现期望信号的整个过程，又要估计滤波权系数以进行有利于 主要目标方向的调整。这些估计过程是以连续的时变形式进行的，这就是自适应 滤波器需要有的自适应收敛过程。如何缩短自适应收敛过程所需要的收敛时间， 这个与算法和结构有关的问题时人们一直重视研究的问题之一[2]。 当然滤波子系统在整个自适应滤波器的设计中也占有很重要的地位，因为它 对最终的滤波性能有很大的影响。本文要研究的是基于维纳滤波原理的LMS算法，那么下面我们需要介绍一下基本维纳滤波原理。 2.2 基本维纳滤波就是用来解决从噪声中提取信号问题的一种过滤(或滤波)方法。它基于平稳随机过程模型，且假设退化模型为线性空间不变系统的。实际上 这种线性滤波问题，可以看成是一种估计问题或一种线性估计问题。基本的维纳 滤波是根据全部过去的和当前的观察数据来估计信号的当前值，它的解是以均方 误差最小条件下所得到的系统的传递函数H(z)h(k)或单位样本响应的形式给出 的，因此更常称这种系统为最佳线性过滤器或滤波器。设计维纳滤波器的过程就 H(z)h(k)是寻求在最小均方误差下滤波器的单位样本响应或传递函数的表达 式，其实质是解维纳-霍夫(Wiener-Hopf)方程。 基本维纳滤波器是这样的，有两个信号x(k)和y(k)同时加在滤波器上。典型 地y(k)包含一个与x(k)相关地分量和另一个与x(k)不相关地分量。维纳滤波器则 产生y(k)中与x(k)相关分量地最优估计，再从y(k)中减去它就得到e(k)。 假定一个N个系数（权值）的FIR滤波器的结构，维纳滤波和原始信号y(k)之间的差信号e(k)为： N,1,Te,y,n,y,wx,y,w(i)xk ( 2.3 ) kkkkkk,i,i,0 xk其中和w分别为输入信号矢量和权矢量，由下式确定： x，，w(0)，，k,,,,xw(1)k,1,,,,w,x,k ( 2.4 ) ,,,,？？ ,,,,xw(N,1),,k,(N,1)，，，， 误差平方为： 22TTTe,y,2yxw，wxxw ( 2.5 ) kkkkkk ,对（3）式两边取期望得到均方误差(MSE) ，若输入x(k)与输出y(k)是联合平稳的，则： 2,,E[e]k 2TTT,E[y],2E[yxw]，E[wxxw] ( 2.6 ) kkkkk 2TT,,，2Pw，wRw 22P,E[yx],,E[y]y(k)E,,kkk其中代表期望，是的方差，是长度为N地互 TR,E[xx]k相关矢量，k是N×N的自相关矩阵。一个MSE滤波系数的图形是碗形地，且只有唯一地底部，这个图称为性能曲面，它是非负的。性能曲面地梯度 可由下式给出： ,d ,,,,2P，2Rw ( 2.7 ) dw 每组系数w(i)(i=1,2,…N-1)对应曲面是一点，在曲面是地最小点梯度为0， wopt滤波权矢量达到最优, ,1w,RP ( 2.8 ) opt 即著名的维纳—霍夫方程的解。自适应滤波地任务是采用合适的算法来调节 ww(0),w(1),？,w(N-1)opt滤波权重iii，从而找到性能曲面地最优点。 维纳滤波的实际用途有限，因为： （1） 它需要已知自相关矩阵R和互相关矢量P，这两个量通常是未知的。 （2） 它包含了矩阵的求逆，非常的耗时。 （3） 若信号为非平稳的，则R和P是时变的，导致必需重复计算wopt。 wopt对于实际的应用需要一种能够依次加入地抽样点而得到的算法。自适应算法就就是用于达到这个目的，而且不需显式计算R和P或进行矩阵求逆[3]。 3 在实际应用中常常会遇到这样的情况：随机信号的统计特性是未知的，或者 信号的统计特性是缓慢的变化着的（非平稳信号），这就促使人们去研究一类特 殊的滤波器，这类滤波器具有以下特点：当输入过程的统计特性未知时，或者输 入过程的统计特性变化时，能够相应的调整自身的参数，以满足某种准则的 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 ， 由于这类滤波器能变动自身的参数以“适应”输入过程统计特性的估计或变化，因此，就把这类滤波器称为自适应滤波器[4] 。在本文中我们研究的是退化图像复原的问题，由于图像自身的多样性和所混入的噪声的随机性和多样性，我们选择 自适应滤波取出图像中混入的噪声。 3.1 3.1.1 首先考虑如下图所示的横向FTR自适应滤波器 它的输入序列以向量的形式记为： T,,X(k),x(k)x(k,1)？x(k,M，1) ( 3.1 ) X(k)R假设取自一均值为零，自相关矩阵为的广义平稳随机过程，而滤波 器的系数矢量（加权矢量）为： T,,W(k),w(k)w(k)？w(k) ( 3.2 ) M12 X(k)W(k)kk以上二式中括号内的为时间指数，因此，和分别表示时刻的 y(k)滤波器输入序列和加权值，滤波器的输出为： M ( 3.3 ) y(k),w(n)x(n,i，1)i,i,1 式中M为滤波器的长度。 d(k)图3.1 中的称为“期望理想响应信号”，有时也可称为“训练信号”，它决 W(k)定了设计最佳滤波器加权向量的取值方向。在实际应用中，通常用一路参 d(k)e(k)y(k)考信号来作为期望响应信号。是滤波器输出相对于的误差，即 (3.4) e(k)显然，自适应滤波控制机理是用误差序列按照某种准则和算法对其系数 ,,wi)n),i,1,2,？,M进行调节的，最终使自适应滤波的目标（代价）函数最小化， 达到最佳滤波状态。 按照均方误差（MSE）准则所定义得目标函数是 2,(k),E[e(k)] ( 3.5 ) 22,E[d(k),2d(k)y(k)，y(k)] 将式（3.4）代入式（3.5），目标函数可以化为 2,(k),E[e(k)] ( 3.6 ) 2TTT,E[d(k)],2E[d(k)W(k)X(k)]，E[W(k)X(k)X(k)W(k)]当滤波系数固定时，目标函数又可以写为 2TT,(k),[d(k)],2W(k)P，W(k)RW(k)] ( 3.7 ) P,E[yx]kk其中，是长度为N的期望信号与输入信号的互相关矢量， TR,E[xx]kk是N×N的输入向量得自相关矩阵。 ,(k)由式（3.7）可见，自适应滤波器的目标函数是延迟线抽头系数（加权 W(k)或滤波系数）的二次函数。当矩阵R和矢量P已知时，可以由权矢量直接求其解。现在我们将式（3.7）对W求倒数，并令其等于零，同时假设R是非奇 wopt异的，由此可以得到目标函数最小的最佳滤波系数为 ,1w,RP ( 3.8 ) opt W(k)这个解就是维纳解，即最佳滤波系数值。因为均方误差函数是滤波系数 的二次方程，由此形成一个形如图（2.2）的超抛物面，当滤波器工作在平稳随 机过程的环境下，这个误差性能曲面就具有固定边缘的恒定形状。自适应滤波系 ,,w(0),i,1,2,？,M数的起始值i是位于性能曲面上的某一点，结果自适应调节过 程，使对应于滤波系数变化的点移动，朝碗底最小点方向移动，最终到达碗底的 最小点，实现了最佳维纳滤波。 最陡下降法就是实现上述搜索最佳值的一种优化技术，它利用梯度信息分析 自适应滤波性能和追踪最佳滤波状态。梯度矢量是由均方误差,(,)的梯度来定义 ,(,)的，在多维超抛物面上任意一点的梯度矢量是对应于均方误差对滤波系数 w(k)i的一阶倒数，由起始点变化到下一点的滤波系数变化量正好是梯度矢量的 负数。换句话说，自适应过程是在梯度矢量的负方向接连的校正滤波系数的，即 在误差性能曲面的最陡下降方向移动和逐步校正滤波系数，最终到达均方误差为 最小的碗底最小点，获得最佳滤波或准优工作状态。 令,(k)M，1代表k时刻的维梯度矢量，这里M等于滤波器滤波系数的数 w(k)目，为自适应滤波器在k时刻的滤波系数和权矢量。按照最陡下降法调节 w(k，1)滤波系数，则在k+1时刻的滤波系数或权矢量可以用下列简单递归关系 来计算： 1w(k，1),w(k)，,[,,(k)] ( 3.9 ) 2 ,式中，为自适应收敛系数或步长，是一个正实常数。根据梯度矢量定义， ,(k)可写成 2,E[e(k)],(k),,w(k) ，，,,,,(k),(k),(k)？, ( 3.10 ) ,,,w(k),w(k),w(k)12M，， ,e(k),E[2e(k)],,E[2e(k)x(k)],w(k) ,(k)当滤波系数为最佳值，即是维纳解时，梯度矢量应等于零。将式（3.7） 代入（3.10）得到 ( 3.11 ) ,(k),,2P，2Rw(k) 因此，在最陡下降算法中，当相关矩阵R和互相关矢量P已知时，由滤波 w(k),(k)系数矢量可以计算梯度矢量，把式（3.11）代入（3.9）中，可以计算出滤波系数的更新值 ( 3.12 ) w(k，1),w(k)，,[P，Rw(k)],k,1,2,？,M 式（3.12）所描述的即是最陡下降算法自适应迭代的基本 公式 小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载 ，且由该式我 们可以不用再直接求R的逆。（3.12）式所示的迭代算法是一个反馈模型，因此 算法的收敛性（稳定性）就非常重要。下面我们简单讨论一下最陡下降法的收敛 性。 3.1.2 首先我们定义k时刻的加权误差矢量为 ,w(k),w(k),w ( 3.13 ) opt 则最陡下降算法式（3.12）可以写成另一种形式 ,w(k，1),,w(k)，,(P,Rw(k)) ( 3.14 ) ,,w(k),,R,w(k),(I,,R),w(k) ,这样，我们得到最陡下降法的稳定性取决于两个因素：自适应步长参数和 X(k)输入信号矢量R的自相关矩阵。根据线性代数里的酉相似变换原理，用酉 Q矩阵R将相关矩阵对角线化，即 HR,QDQ ( 3.15 ) QDR式中为对角线矩阵，它的元素是的特征值，矩阵的列矢量是相关矩阵 HQQR的特征值对应的特征向量的正交集，是的共轭转置。酉矩阵的性质是HHQQ,QQ,I。把式（3.15）代入式（3.14），得到 H,w(k，1),(I,,QDQ),w(k) ( 3.16 ) HQ两边乘以得到， HHHQ,w(k，1),(I,,QDQ)Q,w(k) (3.17a) 或写为 (3.17b) V(k，1),(I,,D)V(k) HV(k),Q,w(k)V(k)其中，为旋转参数矢量或旋转滤波系数矢量误差，的 起始值为 H ( 3.18) V(1),Q[w(1),w]opt 由（3.17b）式，可以推算出 kV(k，1),(I,,D)V(1) ( 3.19 ) DI把单位矩阵和对角线矩阵展开，上式可以写为 k，，(1,)00,,？1,,k01,)0,,？2,, ( 3.20 ) V(k，1),,V(1),,？？？？,,k001,)？,,,,，，M 上式表明，为了保证最陡下降算法的稳定性（收敛性），矩阵中的每一个元 1,,,,k,1,2,？,Mk素的绝对值必须小于1，由此可以得到算法稳定性的收敛条件为 2,,, ( 3.21 ) 0,max ,maxR式中是相关矩阵的最大特征值。在此条件下，对角线矩阵中全部元素 V(k，1),0k,,k当而趋近于零，结果使得。当很大时，意味着自适应滤波 系数矢量趋近于最佳维纳解。 3.2 LMS 3.2.1 LMS 如上节所述，最陡下降算法不需要知道误差特性曲面的先验知识，其算法就 可以收敛到最佳维纳解，且与起始条件无关。但是最陡下降算法的主要限制是它 需要准确测得每次迭代的梯度矢量，这妨碍了它的应用。为了减少计算复杂度和 缩短自适应收敛时间，1960年，美国斯坦福大学的Widrow等提出了最小均方（LMS）算法，这是一种用瞬时值估计梯度矢量的方法，即 2,,E[e(k)],(k), ( 3.22 ) ,w(k) ,,[2e(k)x(k)] , E[,(k)]可见，这种瞬时估计法是无偏的，因为它的期望值确实等于式（3.10） ,(k)的梯度矢量。所以，按照自适应滤波器滤波系数矢量的变化与梯度矢量估 计的方向之间的关系，可以写出LMS算法的公式如下： ,,,1w(k，1),w(k)，,[,,(k)]2 ( 3.23 ) , ,w(k)，,X(k)d(k) 3.2.2 LMS w(k)i（1） 初始化 令所有权重为任一固定值或为0，对每一个接下来地抽样时刻（k=1,2,…,N）执行（2）到（4）。 （2） 计算滤波输出 M y(k),w(i)x(k,i，1), ,1 ( 3.24 ) i （3） 计算估计误差 e(k),d(k),y(k) ( 3.25 ) （4） 更新下一时刻的权 w(k，1),w(k),2,e(k)x(k) ( 3.26 ) 从上面看出，LMS算法具有简洁和易于实现地特点使它成为许多实时系统 的首选算法，LMS算法对每组输入和输出抽样大约需2N－1次乘法和2N－1次加法。太多数信号处理器陡适宜进行乘法和累加运算，使直接实现LMS算法更具有吸引力。 3.2.3 LMS 3.2.4 LMSMatlabwiener2 LMS算法实现的代码可见附录A 3.2.5 wiener2 我们下面介绍一下wiener2()的基本处理方法和过程： 根据图像的局部统计性质，一个像素与它周围的局部区域内的像素相关，因 此去噪图像中某一点像素值可由退化图像中相应于该点的局部区域内的像素值 予以恢复，我们做一个S×S（S=2×L+1,L为窗口半宽）的矩形窗口，对退化图像 中的每个像素在以该像素维中心的窗口内进行相应的滤波处理，可得去噪图像 为： LL r(i,j),w(m，L，1,n，L，1)g(m，i,n，j) ( 3.27 ) ,,,,,,mLnL w(m,n)g(i,j)式中为所要求的权值，为退化图像，写成矩阵形式有： TR,wG ( 3.28 ) 其中 w(1,1)g(i,L,j,L)，，，， ,,,,w(1,2)g(i,L,j,L，1),,,, ， ( 3.29 ) w,G,,,,,？？,,,,w(s,s)g(i，l,j，L)，，，， w(m,n)于是问题转换成求，我们选取一块有较多信号的区域作为特征区域， SA令其为，并对其做自适应滤波处理，以得到该区域的近似原始图像 I(i,j),(i,j),SA，称为理想图像，滤波步骤如下： （1） 先求特征区域中的局部图像均值和方差 LL mean(i,j),g(m，i,n，j) ( 3.30 ) ,,mLnL,,,, LL2var(i,j),[g(i，m,j，n),mean(i,j)] ( 3.31 ) ,,,,,,mLnL （2） 再根据加权最小二乘法，求得自适应滤波后的理想图像 ( 3.32 ) I(i,j),mean(i,j)，c(i,j)[g(i,j),mean(i,j)] 2,,var(i,j)2n,,其中，这里表示噪声方差。 c(i,j)nvar(i,j) Wiener2() 函数的解析可见附录B 3.2.6 LMS—— 自适应滤波器系数矢量的起始值是任意的常数，应用LMS算法调节滤w(1) 波系数具有随机性而是系数矢量w(k)带来非平稳过程。通常为了简化LMS算法的统计分析，往往假设算法连续迭代之间存在以下的充分条件： （1） 每个输入信号样本矢量x(k)与其过去全部样本矢量是统计独立的，不 相关的。 （2） 每个输入信号样本矢量与全部过去的期望响应信号d(k)也是统计独立的。 （3） 期望响应信号依赖于输入样本矢量，但全部过去的期望信号d(k)x(k)样本是统计独立的。 （4） 滤波器抽头输入信号矢量于期望信号包含着全部k的共同的d(k)x(k) 高斯分布随机变量。 通常将基于上述基本假设的LMS算法的统计分析称为独立理论。我们发现 （1）和（2）的观点，假设滤波系数矢量与输入信号矢量和期望信号的独立无关 是很有用的。 但是，在实际中，许多问题对输入过程和期望信号并不满足上述基本假设。 尽管如此，LMS算法的实践经验证明，在有足够的关于自适应过程结构信息的 条件下，基于这些假设所分析的结果仍可用作可靠的设计指导准则，即使某些问 题带有依赖的数据样本。 经过推导，我们得到LMS算法和它的基础——最陡下降算法有系统的精确 数学表达形式。因此，可以得出它们的收敛条件也是一致的，即 2,,, ( 3.33 ) 0,max 当迭代次数k趋近于无穷的时候，自适应滤波系数矢量近似等于最佳维纳 解，但是事实上我们知道，LMS算法的滤波系数总是在最佳维纳解附近波动， 而不能精确等于最佳维纳解。 在满足收敛条件的情况下，自适应步长,决定了收敛速率，当迭代次数一定 ,时，的值越大，则收敛的越快，滤波器的均方误差（MSE）也越来越小；当自 ,适应步长一定时，随着迭代次数的增大，滤波器的均方误差（MSE）越来越小。 4 Matlab 4.1 LMS 下面我们通过实验来查看LMS算法的滤波性能。在实验中我们通过改变自 ,,适应步长的值，来观测的变化对LMS滤波器的性能的影响。以下为实验结果： 图4.1 图4.2 图4.3 2,,,0,,max通过上面的三幅图，我们可以看出当值超出收敛条件时，LMS 滤波器不能收敛，呈发散形式，所得到的复原图像的质量很差，甚至比原图像更 ,差。而仔细观测图4.2和图4.3，我们又发现当满足收敛条件时，的值小，反 ,而效果要差，或者说是，MSE的下降较小，事实上，这是因为当值越小，下 降越平滑，效果本应该更好，但是由于我们设计的是定长的滤波器，也就是滤波 ,器的迭代次数是固定的，在这样的情况下，当然值越大，收敛越快效果越好了。 4.2 LMS 在这里我们把所设计的自适应LMS滤波器跟频域迭代维纳滤波器进行比 较，涉及到了基本频域迭代和加权频域迭代两种形式的滤波器。在实验中我们对 图像人为的加入20db的高斯白噪声。在性能比较方面，频域加权的维纳滤波我 们计算的是复原图像和原图像的MSE。而在LMS滤波器方面，我们是按照目标 函数来计算的MSE。 图4.4 图4.5 通过以上两幅图片，我们可以看出统一的迭代次数的情况下，频域滤波器的 性能要比LMS算法要好。这主要是因为频域迭代维纳滤波器还属于传统滤波器， 它是通过已知的统计知识——功率谱密度来进行迭代的。是一种比较精确的迭代 算法，而我们的LMS滤波器，则是仅仅根据退化图像，进行迭代滤波的，其中 采用的是用瞬时值来估计统计特性，所以在系统的迭代次数下，所取得的效果是 要稍微差一些的。 5 在整个毕业设计的过程中，我们从基本的维纳滤波到最陡下降法再到自适应 LMS算法，进行了系统的分析和实验验证。知道在相关统计特性已知的情况下， 传统滤波器能取得最佳滤波，但是在没有相关的先验知识的情况下，传统滤波器 就不能满足我们的质量要求，这就需要我们的自适应滤波器来实现了。但是，自 适应滤波技术是一种迭代的运算，采用的是逼近的策略，所以这有限次数的迭代 下，我们还是不能精确恢复图像。可是，实验证明，在有足够的关于自适应过程 结构信息的情况下，我们所设计的自适应滤波器还是能取得很好的效果，并且远 远优于传统滤波器。 在目前的自适应技术中怎样解决滤波器中的参考输入是一个很关键的问题， 由自适应滤波算法的原理可知，参考信号必须是与实际噪声信号相关的信号，且 与有用信号不相关。参考信号与噪声信号的相关性越强，则估计出来的噪声才会 越接近真实噪声。在实际中我们往往无法得到符合理想要求的参考信号，但是我 们只要采用于噪声类型一致，统计特性相似的信号就可以取得较好的效果的。在 本课题实验中我们采用的是随机信号作为参考信号，理论上大多数的随机噪声都 是广义平稳的高斯噪声或瑞利噪声，我们所采用的randn()函数生成的参考信号 应该能取得与理论要求很相符的值，但是因为计算机中的随机信号是在内存中随 机获得的，虽然统计特性可以符合要求，但是可能会出现大能量信号，导致图像 质量受到影响。不过，实验证明在随机信号的数值取得较少的情况下，还是能满 足广义平稳的，也能取得较好的效果。 在实验中，我们验证了LMS的收敛性能，比较了它与频域迭代维纳的性能，可以知道，我们所采用的基本LMS算法还是有待改进的。我们需要在收敛速度和处理速度方面，以及稳定性方面进行改善。 本文的顺利完成受到了我的导师周铁工程师，以及海南大学信息科学院的邓 家先老师的大力支持和精心指导。在毕业设计进行阶段，周老师总是严格的要求， 及时为我指引方向并不断给予督促，他的鼓励使我能有信心去克服毕业设计遇到 的困难。邓老师则是在设计的末尾阶段就性能评估方面给了我理论上的指导，使 得毕业设计得以顺利完成。两位老师对学生平易近人、精心负责的态度使我受益 非浅，在此，我谨向两位老师表示衷心的感谢！ 感谢呼志刚同学和罗正华同学多次同我进行算法和编程实现上的交流和讨 论，给予我技术上的指导。另外，还有感谢我所有的朋友和亲人，感谢他们在生 活上和学习上给予我的支持和帮助！ 最后，感谢所有参加论文评审和对本文提出宝贵意见的专家和老师！ [1] 阮秋琦 著.数字图像处理学.北京：电子工业出版社，2001 [2] 何振亚 著.自适应信号处理.北京：科学出版社，2002 [3] [英] Emmanuel C. Ifeacher , Barrie W . Jervis 著，罗鹏飞，杨世海，朱国富，谭全之 译， 数字信号处理实践方法（第二版）。电子工业出版社，2004 [4] 沈凤麟 陈和晏．《生物医学随机信号处理》．合肥：中国科学技术大学出版社．1999 [5] [美] Simon Haykin 著.郑宝玉 等译.自适应滤波器原理（第四版）.北京：电子工业出版 社.2003 [6] Dimitris G. Manokakis,Vinay K. Ingle,Stephen M.Kogon著, 周正 等译.统计与自适应滤波. 北京：电子工业出版社，2003 [7] 施晓红，周佳编著，精通GUI 图形界面编程。北京大学出版社，2003 [8] 董长虹主编，赖志国，余啸海编著，MATLAB 图像处理与应用。国防工业出版社，2004 [9] [加]Simon Haykin,[美]Barry Van Veen 著，信号与系统（第二版）。电子工业出版社，2004 [10] Miroslaw D. Lutovac , Dejan V . Tosic , Brian L . Evans 著，朱义胜，董辉等译，信号处理 滤波器设计——基于MATLAB和Mathematica 设计方法。电子工业出版社，2004 [12] Paulo S.R. Diniz , Eduardo A.B.da Silva , Sergio L. Netto 著，门爱东，杨波，全子一译， 数字信号处理系统分析与设计。北京：电子工业出版社，2004 [13] 胡广书编著，数字信号处理——理论、算法与实现（第二版）。北京：清华大学出版社， 2003 [14] 张玲华，郑宝玉编著，随机信号处理。北京：清华大学出版社，2003 [15] 通信原理，樊昌信，徐炳祥，吴成柯，国防工业出版社，北京，2001 A clear clc F=checkerboard(2); %原图像 D=imnoise(F,'gaussian',0,0.02); %退化图像=参考信号 nhood=[3 3]; % Estimate the local mean of f. localMean = filter2(ones(nhood), D) / prod(nhood); H=D-localMean;％仿照wiener2函数里的求取类似的“均值” [h k]=size(D); L=h*k; D=reshape(D,L,1); %将图像矩阵变为矢量形式 f=zeros(size(D)); %系统输出＝误差信号 y=f; %噪声逼近信号初始化 E=f; %声明误差矢量 %设置输入噪声信号＝基本输入信号 %X=randn(h,k); X=zeros(h,k); X=imnoise(X,'gaussian',0,0.02); X=reshape(X,L,1); W=zeros(L,L); %设置权矢量初值 lmsMSE=f; a=0; %核心算法 u=0.000005; for i=1:L for n=1:L for m=1:i if n-m+1>0 y(n)=y(n)+W(m,n)*X(n-m+1);％滤波器输出 end end end E=D-y;％计算误差 a=a+1; lmsMSE(a)=mean(E.^2);％根据目标函数计算MSE for n=1:L-1 W(i,n+1)=W(i,n)+u*E(n)*X(n);％更新权值 end end a=linspace(1,L,L);%设置画MSE变化曲线的横坐标 f=reshape(E,h,k)+H;％重构图像矩阵 figure,subplot(2,2,1),imshow(F); title('原图像') D=reshape(D,h,k); subplot(2,2,2),imshow(D); title('退化图像') myMSE=mean2((F-f).^2) subplot(2,2,3),imshow(f); title('复原图像') subplot(2,2,4) plot(a,lmsMSE,'r-'); title('当u＝0.000001时的 MSE 变化曲线') xlabel('迭代次数'),ylabel('MSE'); legend('lmsMSE(a)',0); B function [f,noise] = wiener2(varargin) %WIENER2 Perform 2-D adaptive noise-removal filtering. % WIENER2 lowpass filters an intensity image that has been % degraded by constant power additive noise. WIENER2 uses a % pixel-wise adaptive Wiener method based on statistics % estimated from a local neighborhood of each pixel. % % J = WIENER2(I,[M N],NOISE) filters the image I using % pixel-wise adaptive Wiener filtering, using neighborhoods of % size M-by-N to estimate the local image mean and standard % deviation. If you omit the [M N] argument, M and N default to % 3. The additive noise (Gaussian white noise) power is assumed % to be NOISE. % % [J,NOISE] = WIENER2(I,[M N]) also estimates the additive % noise power before doing the filtering. WIENER2 returns this % estimate as NOISE. % % Class Support % ------------- % The input image I can be of class uint8, uint16, or double. % The output image J is of the same class as I. % % Example % ------- % I = imread('saturn.tif'); % J = imnoise(I,'gaussian',0,0.005); % K = wiener2(J,[5 5]); % imshow(J), figure, imshow(K) % % See also FILTER2, MEDFILT2. % The following syntax is grandfathered: % % J = WIENER2(I,[M N],[MBLOCK NBLOCK],NOISE) or [J,NOISE] = % WIENER2(I,[M N],[MBLOCK NBLOCK]) processes the intensity % image I as above but in blocks of size MBLOCK-by-NBLOCK. Use % J = WIENER2(I,[M N],SIZE(I),NOISE) to process the matrix all % at once. % Copyright 1993-2002 The MathWorks, Inc. % $Revision: 5.18 $ $Date: 2002/03/15 15:29:28 $ % Uses algorithm developed by Lee (1980). % Reference: "Two-Dimensional Signal and Image Processing" by % Jae S. Lim, pp.536-540. [g, nhood, block, noise, msg] = ParseInputs(varargin{:});%调用ParseInput()处理不同的输入参 数的情况,且将输入参数以数组的形式表示 if (~isempty(msg)) error(msg); end classin = class(g); classChanged = 0; if ~isa(g, 'double') classChanged = 1; g = im2double(g); end % Estimate the local mean of f. localMean = filter2(ones(nhood), g) / prod(nhood); %用一个所有的值全相等的小矩阵中的二维FIR滤波器对二维数据g进行滤波 %其实质是利用小矩阵对g进行二维卷积运算，最终结果的维数大小 %等于g的维数，取的是卷积结果的中心部分。类似于求均值 % Estimate of the local variance of f. localVar = filter2(ones(nhood), g.^2) / prod(nhood) - localMean.^2; %求取类似于方差的localVar。这样的算法正好体现了图像矩阵中 %元素间的 相关性 % Estimate the noise power if necessary. if (isempty(noise)) noise = mean2(localVar); %若没有输入noise，则取localVar的均值作为noise end % Compute result % f = localMean + (max(0, localVar - noise) ./ ... % max(localVar, noise)) .* (g - localMean); % % Computation is split up to minimize use of memory % for temp arrays. f = g - localMean; %去除图像信号中的直流分量，使得图像信号成为零均值的信号 g = localVar - noise; g = max(g, 0); localVar = max(localVar, noise); f = f ./ localVar; f = f .* g; f = f + localMean; %对处理后的信号再重新加上直流分量 if classChanged, f = changeclass(classin, f); %使图像矩阵的数据类型与输入时相同 end %%% %%% Subfunction ParseInputs %%% function [g, nhood, block, noise, msg] = ParseInputs(varargin) g = []; nhood = [3 3]; block = []; noise = []; msg = ''; switch nargin%表示输入参数个数的系统参数 case 0 msg = 'Too few input arguments.'; return; case 1 % wiener2(I) g = varargin{1};%把输入的第一个参数赋值给g case 2 g = varargin{1}; switch prod(size(varargin{2})) case 1 % wiener2(I,noise) noise = varargin{2}; case 2 % wiener2(I,[m n]) nhood = varargin{2}; otherwise msg = 'Invalid input syntax'; return; end case 3 g = varargin{1}; if (prod(size(varargin{3})) == 2) % wiener2(I,[m n],[mblock nblock]) OBSOLETE warning(['WIENER2(I,[m n],[mblock nblock]) is an obsolete syntax.',... 'Omit the block size, the image matrix is processed all at once.']); nhood = varargin{2}; block = varargin{3}; else % wiener2(I,[m n],noise) nhood = varargin{2}; noise = varargin{3}; end case 4 % wiener2(I,[m n],[mblock nblock],noise) OBSOLETE warning(['WIENER2(I,[m n],[mblock nblock],noise) is an obsolete syntax.',... 'Omit the block size, the image matrix is processed all at once.']); g = varargin{1}; nhood = varargin{2}; block = varargin{3}; noise = varargin{4}; otherwise msg = 'Too many input arguments.'; return; end % checking if input image is a truecolor image-not supported by WIENER2 if (ndims(g) == 3)%判断是否为三维图像 msg = 'WIENER2 does not support 3D truecolor images as an input.'; return; end; if (isempty(block)) block = bestblk(size(g)); C % iterative wiener filter clear all; tic % orignal image I=checkerboard(2); % num of iterations num_iter=128; for snr =20% 20:10:40 %db res = size(I,1) % image size res2 = res^2 % image size square %%% degradation %%% h=[6;1;1;1;1;1;zeros(res2-11,1);1;1;1;1;1]; % blur h=h/sum(h); % normalize Dh= fft(h); % diagonal for c=1:1 % process RGB components individually f=I; % orignal f_vec=reshape(f, res2,1); % vector image mn_f=mean(f_vec); % mean amp_f = norm(f_vec)/res; % amplitude noise_vec = randn(res2, 1)*amp_f/10^(snr/20); % Guassian noise given SNR Dn = abs(fft(noise_vec )).^2/res2; % psd of noise %%% observation %%% g_vec = noise_vec; for i=1:res2 % degraded image g_vec(i) = g_vec(i) + [h(res2-i+2: res2)', h(1:res2-i+1)']*f_vec; % H %is circulant matrix end mn_g = mean(g_vec); g_vec = g_vec -mn_g; % forced to be zero-mean Dg = abs(fft(g_vec )).^2/res2; % psd of degraded image %%%%%%%%%%%% Wiener %%%%%%%%%%%%%%% Df = Dg; % basic iterative Df_add = Dg; % additive Dh_sq = abs(Dh).*abs(Dh); for i=1:num_iter i B = real(ifft((Df.*conj(Dh))./(Dh_sq.*Df + Dn))) ; B_add = real(ifft((Df_add.*conj(Dh))./(Dh_sq.*Df_add + Dn) )); for k=1:res2 fhat_vec(k,1) = [B(res2-k+2: res2)', B(1:res2-k+1)']*(g_vec); fhat_add_vec(k,1) = [B_add(res2-k+2: res2)', B_add(1:res2-k+1)']*(g_vec); end Df = Dg.*(Df.*Df.*Dh_sq)./(( Df.*Dh_sq + Dn).*( Df.*Dh_sq + Dn)); Df_add= ( (Dg-Df_add.*Dh_sq-Dn) .*(Df_add.*Df_add.*Dh_sq) )./(( Df_add.*Dh_sq + Dn).*( Df_add.*Dh_sq + Dn)) + Df_add; if i==1 m=1; Fhat_ft = reshape(fhat_vec,res, res) +mn_g; elseif i==num_iter m=2; Fhat_nd = reshape(fhat_vec,res, res)+ mn_g; Fhat_add_nd = reshape(fhat_add_vec, res, res) +mn_g; end dis(i)=norm(fhat_vec-(f_vec-mn_f) )^2/res2; % MSE dis_add(i)=norm(fhat_add_vec - (f_vec-mn_f))^2/res2; % MSE end % num_iter dis_g(c) = norm(g_vec -(f_vec-mn_f))^2/res2; %MSE G(:,:) = reshape(g_vec +mn_g, res, res) ; end % color %%%%%%% my LMS filter %%%%%%%%%%%%%%%%%%%%%% my_g=g_vec; my_mean_g=mn_g; my_f=zeros(size(my_g)); %系统输出＝误差信号 my_y=my_f; %噪声逼近信号初始化 my_E=my_f; %声明误差矢量 %设置输入噪声信号＝基本输入信号 %X=randn(h,k); my_X=zeros(res); my_X=imnoise(my_X,'gaussian',0,0.02); my_X=reshape(my_X,res2,1); my_W=zeros(num_iter,res2); %设置权矢量初值 lmsMSE=zeros(num_iter,1); a=0; %核心算法 u=0.00001; for i=1:num_iter for n=1:res2 for m=1:i if n-m+1>0 my_y(n)=my_y(n)+my_W(m,n)*my_X(n-m+1); end end end my_E=my_g-my_y; a=a+1; lmsMSE(a)=mean(my_E.^2); for n=1:res2-1 my_W(i,n+1)=my_W(i,n)+u*my_E(n)*my_X(n); end end a=linspace(1,res2,res2); my_f=reshape(my_E,res,res)+my_mean_g; %%% end of my LMS filter%%%%%%%%%%%%%%%%%%%%%%%% if snr==20 % save results_fruits_snr20_iter25_data dis dis_g dis_add snr res num_iter G Fhat_ft Fhat_nd Fhat_nd Fhat_add_nd; elseif snr==30 % save results_fruits_snr30_iter25_data dis dis_g dis_add snr res num_iter G Fhat_ft Fhat_nd Fhat_add_nd; elseif snr==40 % save results_fruits_snr40_iter25_data dis dis_g dis_add snr res num_iter G Fhat_ft Fhat_nd Fhat_add_nd; elseif snr==50 save results_jing2c128_snr50_iter25_data dis dis_g dis_add snr res num_iter G Fhat_ft Fhat_nd Fhat_add_nd; end toc figure subplot(2,3,1) imshow(mat2gray(I),[]) title('Orignal') subplot(2,3,2) imshow(mat2gray(G),[]) %str=['Gaussian noise (SNR=',num2str(snr), 'dB) ', 'and LSI blurring ']; str=['Degraded', ' (blur filter, ', 'SNR=',num2str(snr), 'dB)']; title(str); %title('degradated') strx=['MSE=',num2str(dis_g)]; xlabel(strx) %title('degradated') subplot(2,3,4) imshow(my_f) title('LMS filter') strx=['MSE=',num2str(lmsMSE(num_iter))]; xlabel(strx) subplot(2,3,5) imshow(mat2gray(Fhat_nd),[]) title('Iterative 25th') strx=['MSE=',num2str(dis(num_iter))]; xlabel(strx) subplot(2,3,6) imshow(mat2gray(Fhat_add_nd),[]) title('Additive 25th') strx=['MSE=',num2str(dis_add(num_iter))]; xlabel(strx) figure plot(dis,'- r') hold on plot(dis_add,'. b') hold on plot(lmsMSE,'g:') legend('Iterative','Additive','LMS') %title('MSE v.s. num of iterations') %str2=['degradated image MSE:',num2str( sum(dis_g)/3 )]; %xlabel(str2); ylabel('MSE of Restored Image') xlabel('Num of iterations') str2=['Degraded image: ', 'SNR=',num2str(snr), ' MSE=', num2str(dis_g) ]; title(str2) axis([-1 num_iter+1 -0.001 0.05]); hold off end %snr