% wiener1_template.m
% wiener filter problem

rng(0) % initialize random number generator so everyone has same data
nx = 62; ny = 50;

p = 0.004;
h1 = zeros(27,27); h1(10:18, 11:17) = 1;
yy = ((rand(nx,ny) < p) - p) / sqrt(p*(1-p)); % normalized bernoulli r.p.
ff = conv2(yy, h1, 'same');	% random rect signal

gg = ff + 5 * randn(size(ff));		% add gaussian noise - what variance?

%Pf = ?			% power spectrum of f[m,n]
%H = ?			% wiener filter H[k,l]
% h

%f_hat = ?( ? gg ? )?	% wiener filter estimate from g[m,n]

%rf = ?			% empirical autocorrelation estimate


%{
% Below here are the plotting commands I used, provided as a reference.
% You need not use these exact commands.

clf, pl = @(i) subplot(330+i); colormap(1-gray(256)),

pl(1), imagesc(1:nx, 1:ny, ff'), axis xy, axis image
xlabel m, ylabel n, title 'True: f[m,n]', colorbar horiz

% n=?
% m=?
%pl(2), imagesc(n, m, ?'), axis xy, axis image, axis([-1 1 -1 1]*20)
%xlabel m, ylabel n, title 'R_f[m,n] true', colorbar horiz

% n=?
% m=?
%pl(3), imagesc(n, m, rf'), axis xy, axis image, axis([-1 1 -1 1]*20)
%xlabel m, ylabel n, title 'R_f[m,n] empirical', colorbar horiz

data_rmse = sqrt(mean2((gg-ff).^2));
pl(4), imagesc(1:nx, 1:ny, gg'), axis xy, axis image
xlabel m, ylabel n, colorbar horiz
title(sprintf('Noisy: g[m,n], RMSE=%.2f', data_rmse))

%k=?
%l=?
%pl(5), imagesc(k, l, ?'), axis xy, axis image
%xlabel k, ylabel l, colorbar horiz
%title 'Signal power spectrum: P_f[k,l]'

%pl(6), imagesc(k, l, ?'), axis xy, axis image
%xlabel k, ylabel l, title 'Wiener filter: H[k,l]', colorbar horiz

%pl(8), imagesc(1:nx, 1:ny, ?'), axis xy, axis image
%xlabel m, ylabel n, colorbar horiz
%wiener_mse = mean2((f_hat-ff).^2);
%title(sprintf('Wiener output: RMSE=%.2f', sqrt(wiener_mse)))

fhat2 = npls_sps(gg);

pl(9)
imagesc(1:nx, 1:ny, fhat2'), axis xy, axis image
xlabel m, ylabel n, colorbar horiz
npls_mse = mean2((fhat2-ff).^2);
title(sprintf('NPLS output: RMSE=%.2f', sqrt(npls_mse)))

%orient tall, print('fig', '-deps')
%}
