%-------------------------------------------------------------------------------
% Wiener filter implmentation
%
%   Gonzalez and Woods. Digtial Image Processing 2nd Edition. 263 (eqn.
%   5.8-2).
%
% wienerFilter(I, PSF, UI, N, LAMBDA)
%   I               Degraded Image
%   PSF             Point Spread Function
%   UI              Undegraded Image
%   N               Noise Image
%   LAMBDA          Parametric Factor
%
%   Notes: This function assumes that inputs are properly zero-padded
%
%-------------------------------------------------------------------------------
function fHat = wienerFilter(I, PSF, UI, N, LAMBDA)

    IPS = abs(fft2(UI)).^2;
    NPS = abs(fft2(N)).^2;

% Set up equation variables
    H = fft2(PSF, size(I,1), size(I,2));
    HStar = conj(H);
    HSquared = abs(H).^2;
    G = fft2(I);
    
% Compute Noise to Signal Power Ratio
    NSPR = LAMBDA.*(NPS ./ IPS);
    
% Perform filtering
    den = HSquared + NSPR;
    num = HStar;
    FHat = num ./ den .* G;
    
    fHat = real(ifft2(FHat));
    
    % This line can be included to remove the circulat matrix shift done on
    % the image due to the circular convolution performed by the filtering.
    %
    %fHat = circshift(fHat,[1.*floor(length(PSF)/2) 1.*floor(length(PSF)/2)]);