The source code for the FFT Wiener Filter was created in Matlab as this is a strong mathematical programming language with fast inbuilt Fourier functions which will increase the efficiency of the program. My experience in Matlab has been fairly limited in the past and hence studying Wiener Filtering in Matlab gives me an opportunity to learn a new language.
The basic premise of the code below is taking the Fourier transform of the image and multiplying it by the transform of the filter function(h). The process elimates low frequency signals in the Fourier domain which is the noise and emphasizes the high frequency signals which is the line. The output image is then created by taking the inverse Fourier transform of the product. The filter multiplies each pixel in the Fourier image by this filter to give the final filtered image.
The original model for the FFT Wiener Filter code was taken from a website http://www.owlnet.rice.edu/~elec539/Projects99/BACH/proj2/wiener.html although many alterations have been made to the original Matlab functions and the program restructured. The original code was about twice as long and did not operate as efficiently or as effectively as it does now. Originally the code could only just make out a line at an intensity of 1.5 above the noise mean. Now it clearly shows a line at an intensity of only 0.5 above the noise mean. I have also commented all steps in the code so the Wiener Filtering Process may be followed and understood in relation to the code.
clear; clf; %generate random noise array SIZE = 128; N = randn(SIZE,SIZE); %generates a random Gaussian matrix with SIZExSIZE components N(50,:)=N(50,:) + 2; %line is x=50 and raised by 2 above the noise mean. figure(1) colormap(jet) image(N*50) %multiplies each pixel in the image by 50 h = ones(10,10)/16; %Filter function sigma=1; % should always be 1 for a gaussian noise distribution as stated in the Matlab RANDN function Xf = fft2(N); %Fourier Transform of supplied image Hf = fft2(h,SIZE,SIZE); %Fourier Transform of blurring filter y = real(ifft2(Hf.*Xf)); %y is the real inverse fourier transform of the Fourier transform of %the initial image multiplied by the Fourier Transform of the filter %function (h). This creates the output image y % restoration using generalized Wiener filtering gamma = 1; alpha = 1; ewx = wienerFilter(y,h,sigma,gamma,alpha); figure(2) colormap(jet) image(ewx*50) %magnifies each pixel in the image by 50 return
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