% This part sets up the experiment according to the specifications laid out
% in the Expt_Setup.m file

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if Input_signal == 1
  load camera; X = camera;
elseif Input_signal == 2
  load lenna256; X = lenna;
elseif Input_signal == 3
  load theater; X = theater(50:305, 10:265);
elseif Input_signal == 4
  load boats; X = boats;
elseif Input_signal == 5
  load birthday; X = birthday;
else   
  error('You did not specify any valid Input_signal')
end

% The input signal is normalized between 0 and 1
X = X./255;
N = length(X);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Blurring function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if Blur == 1
  q = 9;
  h_row = [ ones(1,q) , zeros(1, N-q)];
  h_row = h_row(:);  
  h_col = h_row;
  Conv_mat = (h_col * h_row.') / q^2;
elseif Blur == 2
  q = floor(input('What is the blur size? Give a positive integer\n'));
  if q <= 0
    error('Enter a positive blur size')
  end
  h_row = [ ones(1,q) , zeros(1, N-q)];
  h_row = h_row(:);  
  h_col = h_row;
  Conv_mat = (h_col * h_row.') / q^2;
elseif Blur == 3
  Conv_mat = zeros(N,N);
  Conv_mat(1:4,1:3) = 1;
  Conv_mat(1:3,1:4) = 1;
  Conv_mat(5,1:2) = 1;
  Conv_mat(1:2,5) = 1;
  
  Conv_mat(N-3:N,1:2) = 1;
  Conv_mat(N-4:N,1:3) = 1;
  Conv_mat(N-1:N,1:4) = 1;
  Conv_mat(N,1:5) = 1;

  Conv_mat(1:2,N-3:N) = 1;
  Conv_mat(1:3,N-2:N) = 1;
  Conv_mat(1:4,N-1:N) = 1;
  Conv_mat(1:5,N) = 1;
  
  Conv_mat(N,N-3:N) = 1;
  Conv_mat(N-1,N-2:N) = 1;
  Conv_mat(N-2,N-1:N) = 1;
  Conv_mat(N-3,N) = 1;
  Conv_mat = Conv_mat/(pi*17);
elseif Blur == 4
  % The r used below is the same as that used in the Kalifa paper.
  r = 1;
  fft_row = [ones(1,N/4) (2^r) * (abs( 2*((N/4+1):(N/2))/(N)-1)).^r ] ; 
  sum1 = norm([ fft_row  fliplr(fft_row)]);
  fft_row = [ fft_row  fliplr(fft_row)];
  sum2 = norm(fft_row(:));
  fft_row = fft_row*sum1/sum2;  
  h_row   = real(ifft(fft_row));
  h_row = h_row(:);
  clear fft_col
  h_col = h_row;
  Conv_mat = h_col * h_row.';
else
  error('You did not specify any valid Blur')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Noise Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

temp_var = real(ifft2(fft2(X) .* fft2(Conv_mat)));
% org_sigma is the standard deviation of noise
%org_sigma= sqrt(norm(temp_var(:)-mean(temp_var(:)),2)^2 / (N^2*10^(BSNR/10)));
org_sigma = sqrt(norm(temp_var(:),2)^2 / (N^2*10^(BSNR/10)));
n = org_sigma*randn(N,N);
clear temp_var
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Blurred and noisy observation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y = real(ifft2(fft2(X) .* fft2(Conv_mat))) + n ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Defining the wavelet systems for WaRD. Redundant wavelet transforms are
% used.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[temp_w, temp_s] = mrdwt_cycle2D(Y,Wavelet_Filt_type1,Decomp_levels);
clear temp_s

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Noise variance estimation using median estimator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Variance == 'Estimate'
  temp = abs(temp_w(:,:,Decomp_levels,3));
  sigma = median(temp(:))/.67;
else
  sigma = org_sigma;
end
clear temp temp_w

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Power spectral density estimation
% This might not be robust for certain images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if PSD == 'Estimate'
  % Note: We add an offset to the PSD estimate to stabilize it at frequencies
  % where the convolution operator response goes to zero.
  Iterations = 15;
  delta = 8*10^-8;
  Pyy = abs(fft2(Y)).^2;
  Pxx = Pyy;
  G = fft2(Conv_mat);
  Phh = abs(G).^2;
  for i = 1:Iterations
    M = (conj(G) .* Pxx .* fft2(Y)) ./ (Phh .* Pxx + N^2*(sigma^2));
    Pxx_y = (Pxx .* N^2*sigma^2) ./ (Phh.*Pxx + N^2*(sigma^2));
    Pxx = Pxx_y + abs(M).^2;
    Norm_Pxx(i) = norm(Pxx(:));
  end
  Pxx = (Pxx+delta*max(Pyy(:))) *norm(Pxx(:))/norm(Pxx(:) + delta*max(Pyy(:)));
  temp = Pyy(2:256,2:256);  
  Weight = sqrt(Pxx);
else 
  Pxx = abs(fft2(X)).^2;
  Weight = sqrt(Pxx);
end
clear Pyy Pxx Phh










