六逆滤波和维纳滤波-Read

实验六 逆滤波和维纳滤波 SC07023017 陆金星一、实验内容1：先由原始图像（任选）产生待恢复的图像；(产生方法如下：冲激函数为240(,)h x y =，将原始图像与冲激函数卷积产生模糊，然后再迭加均值为0，方差为8，16，32的高斯随机噪声而得到一组待恢复的图像；2：分别用逆滤波和维纳滤波恢复上面图像。

实验图像：Cameraman.bmp二、实验原理 【退化模型】 图像恢复需要根据一定的图像退化模型来进行，一个简单的通用图像退化模型如下图所示，在这个模型中，图像退化过程被模型化为1个作用在输入图像(,)f x y 上的系统H ，它与1个加性噪声(,)n x y 的联合作用导致产生退化图像(,)g x y 。

根据这个模型恢复图像也就是要在给定(,)g x y 和H 的基础上得到对(,)f x y 的某个近似的过程（这里假设已知(,)n x y 的统计特性）。

一维情况：设()f x 中x 的取值范围是0,1,,1A -，()h x 中x 的取值范围是0,1,,1B -，为了避免卷积的各周期重叠(设每个采样函数的周期为M )，取1M A B ≥+-，并将函数用零扩展补齐。

用()e f x 和()e h x 表示扩展的函数，它们的卷积为：1()()()()M e e e e m g x f m h x m n x -==-+∑ 0,1,,M-1x =可以用矩阵形式表示为：(0)(0)(1)(1)(0)(0)(1)(1)(0)(2)(1)(1)(1)(1)(2)(0)(1)(1)e e e e e e e e e e e e e e e e e e g h h h M f n g h h h M f n g H f n g M h M h M h f M n M --+⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥=+==+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-----⎣⎦⎣⎦⎣⎦⎣⎦(,)f x y(,)g x y(,)n x y根据()e h x 的周期性可知()()e e h x h x M =+，所以上式中的H 可进一步写成：(0)(1)(1)(1)(0)(2)(1)(2)(0)e e e ee e e e e h h h h h h H h M h M h -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥--⎣⎦二维情况：由一维情况可以推广到二维情况，得到：(,)0 1 01(,)0 1 1e f x y x A y B f x y A x M B y N ≤≤-≤≤-⎧=⎨≤≤-≤≤-⎩和或 (,)0 1 01(,)0 1 1e h x y x C y D h x y C x M D y N ≤≤-≤≤-⎧=⎨≤≤-≤≤-⎩和或 又可以得到：11000,1,,1(,)(,)(,)(,)0,1,,1M N e e e e m n x M g x y f m n h x m y n n x y y N --===-=--+=-∑∑可以用矩阵形式表示为：011102120(0)(0)(1)(1)(1)(1)M e e e e ee M M H H Hf n f n H H Hg H f n f MN n MN H H H ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=+=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎣⎦⎣⎦⎣⎦ 其中每个i H 是由扩展函数(,)e h x y 的第i 行而来：i (,0)(,1)(,1)(,1)(,0)(,2)(,1)(,2)(,0)e e e ee e e e e h i h i N h i h i h i h i H h i N h i N h i -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥--⎣⎦【循环矩阵对角化】直接从g H f n =+的表达式中计算以求解f 的计算量会很大，可以通过对角化H 来简化。

对于一维循环矩阵，可以将H 的M 个本征矢量组成1个M M ⨯的矩阵W ，[(0),(1),,(1)]W w w w M =- 其中22(1)()[1,,,]jk jM kT MMw k eeππ-= 0,1,,1k M =-，这样可以将H 写成1H WDW -=，D 是对角矩阵，元素是H 的本征值。

对于二维块循环矩阵，定义1个MN MN ⨯的矩阵W ，其每个元素为：2(,)jim MN W i m eW π=，,0,1,,1i m M =-，其中N W 为1个N N ⨯的矩阵，其每个元素为：2(,)jkn MN W k n eπ=，,0,1,,1k n N =-，这样可以将H 写成1H W DW -=。

将循环矩阵对角化应用与退化模型计算，可以有n W f DW g W111---+=，令1,,1,0,1,,1,0-=-=N v M u ，则以下各式成立：1100110011001001(,)(,)exp 21(,)(,)exp 21(,)(,)exp 21(,)(,)M N e x y M N e x y M N e x y M N e x y ux vy G u v g x y j MN M N ux vy F u v f x y j MN M N ux vy N u v n x y j MNM N H u v h x y MNπππ--==--==--==--==⎡⎤⎛⎫=-+ ⎪⎢⎥⎝⎭⎣⎦⎡⎤⎛⎫=-+ ⎪⎢⎥⎝⎭⎣⎦⎡⎤⎛⎫=-+ ⎪⎢⎥⎝⎭⎣⎦=∑∑∑∑∑∑∑1exp 2ux vy j M N π⎡⎤⎛⎫-+ ⎪⎢⎥⎝⎭⎣⎦∑而D 的MN 个对角元素可表示为：(),mod , 0 k MN H k N i k D k i N i k ⎧⎛⎫⎢⎥⨯=⎪ ⎪⎢⎥=⎣⎦⎨⎝⎭⎪≠⎩如如，将MN 并入(,)H u v ，得到0,1,,1(,)(,)(,)(,)0,1,,1u M G u v H u v F u v N u v v N =-=+=-上式说明，求解退化模型的大系统方程，只需计算很少几个N M ⨯的傅里叶变换就可以了。

由Hf g n -=，在对n 没有先验知识的情况下，寻找1个f 的估计值∧f ，使n 的模或范数最小：⎪⎭⎫⎝⎛-⎪⎭⎫ ⎝⎛-=-==∧∧∧f H g f H g fH g n n nTT 22，这样可把恢复问题看作是对∧f 求式2∧∧-=⎪⎭⎫ ⎝⎛f H g f L 的最小值，将L 对∧f 求微分并设结果为零，再设N M =和1-H 存在，可得无约束恢复公式：()()g H g H H H g H H H f T T T T 1111----∧===当选取∧f 的1个线性操作符Q (变换矩阵)，使得2∧fQ 最小。

设l 为拉格朗日乘数，要找到最小化准则函数()222nHfg l fQ f L --+=⎪⎭⎫⎝⎛∧∧的∧f ，可得到有约束恢复公式(令l s 1=)：[]g H QsQ H H f T T T 1-∧+=。

【逆滤波】,逆滤波是一种无恢复约束，设M N =，则将1H WDW-=代入1f Hg ∧-=可以得到，11111()f H g WDW g WD W g ∧-----===，如将上式两边乘以1W -得到：111ˆW f D W g ---=。

式中各个元素能写成形式：(,)(,)(,)G u v F u v H u v ∧=，这种恢复方法称为逆滤波。

将(,)H u v 看作1个滤波函数，则它与(,)F u v 的乘积是退化图像(,)g u v 的傅里叶变换。

用(,)H u v 去除(,)G u v 就是1个逆滤波过程。

将(,)(,)(,)G u v F u v H u v ∧=的结果求反变换就得到恢复后的图像：11(,)ˆ(,)[(,)](,)G u v f x y F F u v F H u v ∧--⎡⎤==⎢⎥⎣⎦， 实际为了减小噪声的影响以及防止(,)H u v 取零或很小的影响，将恢复转移函数(.)M u v 取为 (.)(.)1/(.)k H u v dM u v H u v ≤⎧=⎨⎩如其它，其中k 和d 为小于1的常数，且d 选的较小为好。

【维纳滤波】维纳滤波是一种有约束恢复，是一种最小均方误差滤波器，设f R 和n R 分别是f 和n 的相关矩阵，则1f R WAW -=，1n R WBW -=，其中A 和B 中的元素对应f R 和n R 中的相关元素的变换，代入有约束恢复公式，可以得到：1111111()()T T f n f H H sR R H g WD DW sWA BW WD W g ∧-------=+=*+*也可以写成如下形式：22(,)1ˆ(,)(,)(,)(,)(,)/(,)n f H u v F u v G u v H u v H u v s S u v S u v ⎡⎤⎢⎥=⨯⎡⎤⎢⎥+⎣⎦⎣⎦ 其中，(,)n S u v 和(,)f S u v 分别是(,)e f x y 和(,)e n x y 的功率谱，可通过FFT 变换后模值的平方求得。

上式中，如果1=s ，则方括号中的项就是维纳滤波器；如果s 是变量，就称为参数滤波器；当没有噪声时，(,)n S u v =0，维纳滤波器退化为逆滤波器。

三、实验程序实验程序包含1个.m文件exp_6.m 执行文件，完成图像模糊与图像恢复功能，并给出实验结果具体程序内容在程序内有注释，实验时直接运行exp_6.m即可完成上述实验内容。

四、实验结果与分析【模糊噪声图】【逆滤波恢复结果】【维纳滤波恢复结果】【实验结果分析】（1）在带有噪声干扰的情况下，逆滤波的效果随着噪声的方差值的变大而变小，维纳滤波同样结果，同时恢复结果都带有不同程度的振铃效应。

（2）在带有噪声干扰的情况下，维纳滤波的效果要比逆滤波的效果好，且随着噪声的加强对比效果变得更明显。