implementation for iterative wiener filter

This is my implementation in the iterative wiener filter in this paper : http://www.tsc.uc3m.es/~navia/LATDS07/IterativeWienerFilter.pdf

I wish it will help anyone

function wiener()

clc;
clear;

f=im2double(rgb2gray(imread('lena.jpg')));

imshow(f);
f=imresize(f,[32 32]);

figure,imshow(f)
[r c]=size(f);
h=fspecial('average');

g=imfilter(f,h,'circular');

s_avg = sum(sum(f))/(r*c);

SNR=90;

n_sigma=s_avg/(10^(SNR/20));

n=n_sigma*randn(size(f));

g=g+n;

[Nf,Mf]=size(g);

[Nh,Mh]=size(h);

L1=floor(Nh/2);

L2=floor(Mh/2);

H=zeros(Nf*Mf);

k=1;

for row=1:Mf,

      for col=1:Nf,

          hh=zeros(Nf,Mf);

          hh(1:Nh,1:Mh)=h;

          hh=circshift(hh,[col-1-L1,row-1-L2]);

          H(k,:)=hh(:)';

          k=k+1;
      end
  end

%%make vector of m^2*1 of the f,n,g

f=reshape(f',size(f,1)*size(f,2),1);

g=reshape(g',size(g,1)*size(g,2),1);

n=reshape(n',size(n,1)*size(n,2),1);

%%%calculate the autocorrelation matrix of f ,g,n

u=mean(g);

g1=autom(g-u);

Rg=toeplitz(g1);

n1=autom(n);

Rn=toeplitz(n1);%%%%%

Rf=Rg;

steps=10;

mse=zeros(1,steps);

for i=1:steps

      B=Rf*H'*inv( (H*Rf*H') +Rn);

      fHat=B*(g); 

      Rf=B*Rg*B';

      im=reshape(fHat,[32 32]);

      g=reshape(g,[32 32]);

      %figure,imshow(im',[]);

      mse(1,i) = sum(sum((im(:)-g(:))));

      g=reshape(g',size(g,1)*size(g,2),1);

end

t=1:steps;

mse

plot(t,mse);

end

function [Rxx]=autom(x)

N=length(x);

Rxx=zeros(1,N);

for m=1: N+1

      for n=1: N-m+1

          Rxx(m)=(Rxx(m)+x(n)*x(n+m-1))/N-m+1;

      end;

end

end

[Walter Roberson]

This should probably go into the File Exchange (FEX)

[Walter Roberson]

Please read the guide to tags and retag this; see http://www.mathworks.co.uk/matlabcentral/answers/43073-a-guide-to-tags