逆滤波和维纳滤波
一、實驗目的
? ? ? ? 利用逆濾波和維納濾波，對Lena加噪運動模糊降質圖像進行復原，比較不同參數選擇對復原結果的影響。
二、實驗內容
1)?????輸入Lena圖像，對圖像進行運動降質；降質模型：
2)?????對圖像疊加高斯白噪聲；
3)?????尋找最佳逆濾波半徑r；
4)?????逆濾波；
5)? ? ?IFFT，展示結果；
6)?????再尋找最佳維納濾波K值；
7)?????維納濾波；
8)?????IFFT，展示結果。
三、實驗過程和數學原理
1、逆濾波
? ? ? ? 在無噪情況下，逆濾波是完美的：
????
? ? ? ? 然而實際情況都有噪聲：
? ? ? ? 如果H(u,v)存在零點，那么在H(u,v)零點附近進行復原，會導致第二項變得很大很大，復原效果很差。
? ? ? ? 理論上說，我們應該找出H(u,v)的所有零點，然后規避這些零點進行逆濾波。然而，降質模型零點十分分散（實際上是一條條斜線，實驗中會看到），并且要作無數個圓區域，編程非常麻煩。
? ? ? ? 實際編程中，我們采用如下思路：既然fft后的頻譜中，信號頻譜主要集中在低頻分量，那么，我們用fftshift，將頻譜移到中心；以頻譜中心為圓心，規定一個圓區域，在圓內正常逆濾波；而圓外是大量的較小的噪聲分量，不給任何機會，直接賦值為0。這樣，我們就不需要考慮H(u,v)的零點影響了。
? ? ? ? 此外，再通過PSNR最大準則，尋找最佳濾波半徑R。實驗結果如下：
? ? ? ? 實驗代碼：
function Imagerestoration %~~~~~~~~~~~~~~~~~~~~~~~~~~實驗說明~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% % 1、由于T=5,a=b=1效果太差太差，幾乎無法看到復原現象，因此本實驗采用T=1,a=b=0.02降質模型。 % 2、作為實驗講義的補充，本實驗加入均值為0、方差為1e-3的AGWN模型。 % 若方差過大，逆濾波效果也不理想。 % 3、一定要對fftshift后的頻譜進行運動模糊處理！！！否則現象都是錯誤的！！！ % 2018-5-19 by XING %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% clear;close all;clc; fprintf('-------------------------------逆濾波實驗-------------------------------\n'); I = im2double(imread('lena512.bmp'));% [0,1] [M,~] = size(I);% square% Display the original image.
figure;
subplot(1,3,1), imshow(I);
title(’\fontsize{20}原始圖像’);
%% Simulate a Motion Blur：H(u,v)
T=1;a=0.02;b=0.02;
v=[-M/2:M/2-1];u=v’;
A=repmat(a.u,1,M)+repmat(b.v,M,1);
H=T/pi./A.sin(pi.A).exp(-1ipi.*A);
H(A==0)=T;% replace NAN
%% Get the blurred Image
% Warning: fftshift should be written
F=fftshift(fft2(I));
FBlurred=F.*H;
% Display the blurred image
IBlurred =real(ifft2(ifftshift(FBlurred)));
subplot(1,3,2), imshow(uint8(255.*mat2gray(IBlurred)));
title(’\fontsize{20}運動模糊圖像’);
%% Deblur perfectly without Noise
FDeblurred=FBlurred./H;
IDeblurred=real(ifft2(ifftshift(FDeblurred)));
subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred)));
title(’\fontsize{20}無噪情況下直接逆濾波’);
%% Simulate Noise Model
noise_mean = 0;
noise_var = 1e-3;
noise=imnoise(zeros(M),‘gaussian’, noise_mean,noise_var);
FNoise=fftshift(fft2(noise));
%% Get the Blurred_Noised Image
FBlurred_Noised=FNoise+FBlurred;
% Display the blurred_noised image
IBlurred_Noised=real(ifft2(ifftshift(FBlurred_Noised)));
figure;
subplot(1,3,1), imshow(uint8(255.*mat2gray(IBlurred_Noised)));
title(’\fontsize{20}加噪運動模糊圖像’);
%% Deblur when Ignoring Noise
FDeblurred2=FBlurred_Noised./H;
FH1=abs(FDeblurred2);
IDeblurred2=real(ifft2(ifftshift(FDeblurred2)));
subplot(1,3,2), imshow(uint8(255.*mat2gray(IDeblurred2)));
title (’\fontsize{20}有噪情況下直接逆濾波’);
%% Find out the best Radius
maxPSNR=0;
bestRadius=0;
tic;
for Radius=33:1e-2:34 % 預實驗bestr約為33.8左右
FDeblurred2=zeros(M);
end
fprintf(’ 最佳濾波半徑: %.1f\n’, bestRadius);
fprintf(’ 最大PSNR： %d dB\n’, round(maxPSNR));
fprintf(’ 尋找最佳半徑耗時: %.1f s\n’, toc);
% Displace the best Restoration
FDeblurred2=zeros(M);
for a=1:M
for b=1:M
if sqrt((a-M/2).^2+(b-M/2).^2)<bestRadius
FDeblurred2(a,b)= FBlurred_Noised(a,b)./H(a,b);
end
end
end
FH2=abs(FDeblurred2);
IDeblurred2=real(ifft2(ifftshift(FDeblurred2)));
subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred2)));
title([’\fontsize{20}最佳半徑為 ', num2str(bestRadius),‘的圓內逆濾波’]);
fprintf(’\n 半徑逆濾波原理: \n’);
fprintf(’ 1：通過fft2，使有效信號集中在低頻區域，噪聲集中在高頻區域；\n’);
fprintf(’ 2：通過fftshift，把頻譜移到矩陣中心；\n’);
fprintf(’ 3：以矩陣中心為圓心，規定一個半徑為r的圓域；圓外直接賦0，抑制噪聲；圓內正常逆濾波。\n’);
fprintf(’ 這比尋找、避開H(u,v)的零點濾波，編程要簡單很多，效果也不錯！\n’);
fprintf(’\n 實驗結果說明: \n’);
fprintf(’ 1：對于無噪運動模糊圖像，逆濾波復原幾乎完美 ^ ^\n’);
fprintf(’ 2：對于有噪運動模糊圖像，直接逆濾波是災難，取半徑濾波效果尚可。\n’);
fprintf(’\n 同時說明: PSNR是多么不靠譜！\n’);
figure;
subplot(1,2,1),imshow(im2double(uint8(FH1)));
title (’\fontsize{20}有噪情況下直接逆濾波得到的圖像頻譜’);
subplot(1,2,2),imshow(im2double(uint8(FH2)));
title (’\fontsize{20}圓內逆濾波得到的圖像頻譜’);
fprintf(’\n~程序已暫停；按任意鍵進行維納濾波實驗~~\n’);
pause;
2、維納濾波
????? 理想維納濾波器為：
? ? ? 實際應用中，NSR難以被精確計算。因此，我們常常設為K，并尋找最佳K值。
????? 實驗結果：
fprintf(’\n-------------------------------維納濾波實驗-------------------------------\n’);
fprintf(’ 根據大量實驗，我發現PSNR無法作為尋找最佳K值的標準：\n’);
fprintf(’ 1：參數K將會停留在K=0處，即逆濾波；\n’);
fprintf(’ 2：最大PSNR達到70dB以上，實際圖像質量極差！\n’);
fprintf(’\n 因此，我通過觀察，選擇了最佳K=0.05~\n’);
fprintf(’ ps. 程序中保留了根據PSNR尋找K值的代碼，有興趣可以嘗試 ^ ^\n’);
% Display the blurred_noised image again
figure();
subplot(1,3,1);
imshow(uint8(255.*mat2gray(IBlurred_Noised)));
title(’\fontsize{20}加噪運動模糊圖像’);
% Deblur with theoretic NSR
buf=(abs(H)).^2; % Notice ‘.’ !!!
NSR=FNoise./F;
FDeblurred3=FBlurred_Noised./H.buf./(buf+NSR);
IDeblurred3=real(ifft2(ifftshift(FDeblurred3)));
subplot(1,3,2), imshow(uint8(255.mat2gray(IDeblurred3)));
title(’\fontsize{20}K=NSR的理論維納濾波’);
% Find out the best K
% tic;
% maxPSNR=0;
% beskK=0;
% for K=0:1e-2:1
% FDeblurred2=zeros(M);
% FDeblurred3=FBlurred_Noised./H.*buf./(buf+bestK);
% IDeblurred3=real(ifft2(ifftshift(FDeblurred3)));
%
% % Calculate PSNR and compare with the best
% PSNR=PSNRcal(IDeblurred3,I);
% if PSNR>maxPSNR
% maxPSNR=PSNR;
% bestK=K;
% end
% end
%
% fprintf(’ 最佳K值: %.2f\n’, bestK);
% fprintf(’ 最大PSNR: %d dB\n’, round(maxPSNR));
% fprintf(’ 尋找最佳K值耗時: %.1f s\n’, toc);
% Deblur with best K
bestK=0.05;
FDeblurred3=FBlurred_Noised./H.*buf./(buf+bestK);
IDeblurred3=real(ifft2(ifftshift(FDeblurred3)));
% Display the best restored Image
subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred3)));
title([’\fontsize{20}實際最佳K= ', num2str(bestK),‘的維納濾波’]);
fprintf(’ Written by XING\n’);
fprintf(’ 2018-5-19 Beijing\n’);
end
? ? ? ? 實驗中用到的PSNR計算子程序：
function PSNR=PSNRcal(I,I2)h=512;w=512;
B=8;% 編碼一個像素用8個二進制位
MAX=2^B-1;% 圖像有多少灰度級
MES=sum(sum((I-I2).^2))/(hw);% 均方差
PSNR=20log10(MAX/sqrt(MES));% 峰值信噪比
end
總結
- 上一篇： 小白能读懂的 《手把手教你学DSP(TM
- 下一篇： Package php5 have no