维纳滤波及其应用

数学原理

维纳滤波(wiener filter)是一种以其输出与期望输出的差的平方最小为最优准则的线性滤波器。

信号退化模型可以用如下公式表示:

\vec{y}=\mathbf{H}\vec{x}+\vec{n}

其中\vec{x}是干净无噪信号,\vec{n}是噪声,对于噪声最常见的假设是均值为0,方差为\sigma^2的加性白高斯噪声。\mathbf{H}是退化算子,不同的退化形式对应不同的退化算子,如:在去噪问题中,\mathbf{H}=\mathbf{I};在去模糊问题中,\mathbf{H}是模糊卷积核h对应的矩阵。\vec{y}是观察到的退化信号。

下面为了简化书写,不再标出向量的箭头符号了。

以下推导基于添加均值为0,方差为\sigma^2的加性白高斯噪声的假设 

我们希望能利用观察到的退化信号y来获得干净信号x,这是一个病态问题,是无法实现的,我们只能获得去噪信号\hat{x}

\mathbf{G}来表示维纳滤波,有\hat{x}=\mathbf{G}y

误差:

\begin{aligned} e&=\mathbb{E}(|\hat{x}-x|^2)\\ &=\mathbb{E}(|\mathbf{G}y-x|^2)\\ &=\mathbb{E}(|\mathbf{G}(\mathbf{H}x+n)-x|^2)\\ &=\mathbb{E}(|(\mathbf{G}\mathbf{H}-\mathbf{I})x+\mathbf{G}n|^2) \end{aligned}

 展开来是:

(\mathbf{G}\mathbf{H}-\mathbf{I})^{\top}(\mathbf{G}\mathbf{H}-\mathbf{I})\mathbb{E}(x^2)+(\mathbf{G}\mathbf{H}-\mathbf{I})^{\top}\mathbf{G}\mathbb{E}(x^{\top}n)+(\mathbf{G}\mathbf{H}-\mathbf{I})\mathbf{G}^{\top}\mathbb{E}(xn^{\top})+\mathbf{G}^{\top}\mathbf{G}\mathbb{E}(n^2)

因为噪声和信号是相互独立的,且\mathbb{E}(n)=0,所以上述4项中的中间两项是0

 整理得e=(\mathbf{G}\mathbf{H}-\mathbf{I})^{\top}(\mathbf{G}\mathbf{H}-\mathbf{I})\mathbb{E}(x^2)+\mathbf{G}^{\top}\mathbf{G}\mathbb{E}(n^2)

由定义得:\mathbb{E}(x^2)=\int x^2 p_X(x)dx\mathbb{E}(n^2)=\int n^2 p_N(n)dn分别表示信号的能量s^x和噪声的能量s^n

根据维纳滤波的最优准则,我们的目标是误差e最小,即对e关于\mathbf{G}求导,导数为0。

 得:(\mathbf{H}^{\top}\mathbf{H}s^x+s^n)\mathbf{G}-\mathbf{H}^{\top}s^x=0

\therefore \mathbf{G} = (\mathbf{H}^{\top}\mathbf{H}+\frac{s^n}{s^x})^{-1}\mathbf{H}^{\top}

  • 如果我们在去噪任务中应用维纳滤波的话,因为\mathbf{H}=\mathbf{I},维纳滤波退化成一个缩放系数(Wiener shrinkage coefficient):\frac{s^x}{s^x+s^n}
  • 如果我们在去模糊任务中应用维纳滤波的话,那么退化算子\mathbf{H}必然存在对应的模糊卷积核h
  • 那么上式有另一种等价形式:G = \mathcal{F}((\mathcal{F}(\bar{h})\mathcal{F}(h)+\frac{s^n}{s^x})^{-1}\mathcal{F}(\bar{h})),其中\bar{h}表示翻转卷积核

 matlab代码示例

使用`deconvwnr`函数,用法示例:

J = deconvwnr(I,psf,nsr)
J = deconvwnr(I,psf,ncorr,icorr)
J = deconvwnr(I,psf)

其中,常用的是第一条和第三条:

  • I:值域为[0,1]的double数据类型图像;
  • psf:对I进行卷积的点扩散函数(point-spread function)[即卷积核];
  • nsr:加性噪声的噪信功率比。

读入图片:

I = im2double(imread('lena.png'));

 示例一:

图片被模糊&被噪声污染

% 构造运动模糊卷积核
LEN = 20;
THETA = 30;
PSF = fspecial('motion', LEN, THETA);
% 做卷积操作
blurred = imfilter(I, PSF, 'conv', 'circular');
% 被加性白高斯噪声污染
noise_mean = 0
noise_var = (15 / 255) ^ 2;
blurred_noisy = imnoise(blurred, 'gaussian', noise_mean, noise_var);

用维纳滤波做解卷积

% % 在去噪任务上应用维纳滤波,理论上需要被污染的图片的能量 和 噪声的能量
% 能量的计算方式为:\sum x[i]^2
% 其中,因为均值为0,所以噪声的能量就是noise_var * n
num = prod(size(I));
noise_energy = noise_var * num;
blurred_noisy_energy = norm(blurred_noisy(:)) ^ 2;     
nsr = noise_energy / blurred_noisy_energy;   
   

% % 在去模糊任务上应用维纳滤波
% 只需要模糊卷积核

% 综上,解卷积过程可以写成
wnrOnBlurred_noisy = deconvwnr(blurred_noisy, PSF, nsr); 
       

示例二:

图片被多个模糊核卷积

% 构造运动模糊卷积核
LEN = 20;
THETA = 30;
PSF1 = fspecial('motion', LEN, THETA);
% 做卷积操作
blurred1 = imfilter(I, PSF1, 'conv', 'circular');

% 构造高斯模糊卷积核
kernel_size = 3;
sigma = 25/255;
PSF2 = fspecial('gaussian', kernel_size, sigma);
% 做卷积操作
blurred2 = imfilter(blurred1, PSF2, 'conv', 'circular');

多层嵌套维纳滤波解卷积

wnrOnBlurred2 = deconvwnr(deconvwnr(blurred2, PSF2), PSF1);

  • 3
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值