function[deconv,psf]=blind_dconv_gausspsf(z,iterations,hesize,alpha,fe,k)

if(mod(hesize,2)~=1)
  error('Size of psf must be odd');
end

YSIZE=size(z,1);
XSIZE=size(z,2);
SIZEMAX=max([XSIZE YSIZE]);

%z is the linearly degraded image, fe is the image estimate, and k is the
%finite support. he is the PSF estimate.

%This is the PSF ESTIMATE (we estimate it as an gaussian PSF)
he=make_psf(hesize);

G=fft2(z,(2*SIZEMAX),(2*SIZEMAX)); %z is the degraded image

n=0;

while(n<iterations)
  Fe=fft2(fe,(2*SIZEMAX),(2*SIZEMAX));Hoe=fft2(he,(2*SIZEMAX),(2*SIZEMAX));

  %impose fourier constraints
  Hne=(G.*conj(Fe))./(abs(Fe).^2+(alpha./(abs(Hoe).^2)));
  hne=ifft2(Hne);

  %impose blur constraints
  hne=abs(hne(1:hesize,1:hesize));
  Hoe=fft2(hne,(2*SIZEMAX),(2*SIZEMAX));

  %Impose Fourier constraints
  Fne=(G.*conj(Hoe))./((abs(Hoe).^2+(alpha./abs(Fe).^2)));
  fne=ifft2(Fne);

  %impose image constraints
  fne=abs(fne(1:YSIZE,1:XSIZE));
  fne(~k)=0;
  fe=abs(fne);he=abs(hne);
  n=n+1;
end

deconv=fne;
psf=hne;
