function xhat=wienerfilter(reflectances,ys,A,sigma2)

% WIENERFILTER recovers a signal using the Wiener filter
%
%xhat=wFilter(reflectances,ys,A,sigma2)
%
%INPUTS:
%   reflectances: matrix giving the data (columns are for a single surface)
%           size(reflectances)=length(wave) x number of surfaces
%   ys:     matrix giving the noisy measurements (columns are for a single surface)
%           size(ys)=numfilters x number of surfaces
%   A:      matrix giving the measurement procedure (rows are the filters),
%           so y=A*reflectance + noise,    size(A)=numfilters x length(wave)
%   sigma2: variance of the noise added, scalar if global value is used,
%           matrix same size as ys if local estimates are used
%
%OUTPUTS:
%   xhat:        Matrix giving the estimates of the reflectances
%                size(xhat)=length(wave) x number of surfaces
%
%The Wiener filter is the optimal linear estimator under certain
%assumptions about the signal and noise distributions.  
%The primary assumption is that the signal and noise are independent.  
%This is not actually true but isn't a huge problem.

    
if length(sigma2(:))==1
    sigma2=repmat(sigma2,size(ys,1),size(ys,2));
end

xhat=zeros(size(reflectances,1),size(ys,2));
Rxx=reflectances*reflectances'/size(reflectances,2);     %autocorrelation
    
for inputnum=1:size(ys,2)   %iterate through each input vector
    y=ys(:,inputnum);   %measurements for current
    W = Rxx*A'*inv(A*Rxx*A' + diag(sigma2(:,inputnum)));    %Wiener filter for current noise level
    xhat(:,inputnum)=W*y;   %reflectance estimate is a linear combination of the measurements
end