数学原理
维纳滤波(wiener filter)是一种以其输出与期望输出的差的平方最小为最优准则的线性滤波器。
信号退化模型可以用如下公式表示:
其中是干净无噪信号,
是噪声,对于噪声最常见的假设是均值为0,方差为
的加性白高斯噪声。
是退化算子,不同的退化形式对应不同的退化算子,如:在去噪问题中,
;在去模糊问题中,
是模糊卷积核
对应的矩阵。
是观察到的退化信号。
下面为了简化书写,不再标出向量的箭头符号了。
以下推导基于添加均值为0,方差为的加性白高斯噪声的假设
我们希望能利用观察到的退化信号来获得干净信号
,这是一个病态问题,是无法实现的,我们只能获得去噪信号
。
用来表示维纳滤波,有
。
误差:
展开来是:
因为噪声和信号是相互独立的,且,所以上述4项中的中间两项是0
整理得
由定义得:,
分别表示信号的能量
和噪声的能量
。
根据维纳滤波的最优准则,我们的目标是误差最小,即对e关于
求导,导数为0。
得:
- 如果我们在去噪任务中应用维纳滤波的话,因为
,维纳滤波退化成一个缩放系数(Wiener shrinkage coefficient):
- 如果我们在去模糊任务中应用维纳滤波的话,那么退化算子
必然存在对应的模糊卷积核
- 那么上式有另一种等价形式:
,其中
表示翻转卷积核
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);