wiener_filter.m

Wiener filter in signal processing

% WienerFilter

function ex = wienerFilter(y,h,sigma,gamma,alpha);
%
% ex = wienerFilter(y,h,sigma,gamma,alpha);
%
% Generalized Wiener filter using parameter alpha. When
% alpha = 1,(The noise Power) it is the Wiener filter. It is also called
% Regularized inverse filter.
%

SIZE = size(y,1);
Yf = fft2(y);                     %Fourier transform of original image N by filter function
Hf = fft2(h,SIZE,SIZE);           % Fourier transform of filter function
Pyf = abs(Yf).^2/SIZE^2;

sHf = Hf.*(abs(Hf)>0)+1/gamma*(abs(Hf)==0);   %Fourier transform of filter function + its inverse
iHf = 1./sHf;         %inverse

iHf = iHf.*(abs(Hf)*gamma>1)+gamma*abs(sHf).*iHf.*(abs(sHf)*gamma<=1); 
%array element increases if bigger than noise variance, decreases if below noise variance. 
%iHf is the Power Spectral Frequency Fourier Transform


Pyf = Pyf.*(Pyf>sigma^2)+sigma^2*(Pyf<=sigma^2);    %Pyf is the noise power spectrum..increased for 
%pyf bigger than noise mean, if signal is small pyf will be made smaller.


Gf = iHf.*(Pyf-sigma^2)./(Pyf-(1-alpha)*sigma^2);   %Gf is the esitmated image Fourier Transform
%Gf is basically the filter functions effect on Pyf-the transformed signal N

% Restored image 
eXf = Gf.*Yf;   %The filter multiplies each pixel in the Fourier image(Yf) by this filter(Gf)

ex = real(ifft2(eXf));     
 %inverse fourier tranform of the filtered Fourier transform of the new image yields 
 %the hopefully noise free result.

return