%Adapted from Hussein Aly's program by Eric Dubois 2005-05-03
% Statistical parameters computation
%include mean, variance and autocorrelation
close all;
clear all;
N=64;
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u v]=meshgrid(ux,uy);
strin = input('Enter filename of the input grayscale image: ','s');
A = im2double(imread(strin));
[I J]=size(A);
prodsize=I*J;
mean=sum(sum(A))/prodsize;
Az=A-mean;
clear A;
var=((norm(Az,'fro'))^2)/prodsize;
%%%%%%%%%%Dubois' Design
n=0:1:N-1;
wd=zeros(N,1);
a0=0.40217;
a1=0.49703;
a2=0.09892;
a3=0.00183;
wd(n+1)=a0-a1*cos(2*pi*n/N)+a2*cos(2*pi*2*n/N)-a3*cos(2*pi*3*n/N);
W=wd*wd';
V=norm(W(:),2)^2;
PSD=zeros(N,N);
shx=N/4;     % overlap
shy=N/4;     % overlap
Mx=floor((I-N)/(N-shx));
My=floor((J-N)/(N-shy));

for i=0:1:Mx,
    for j=0:1:My,
        PSD=PSD+(abs(fft2(W.*Az(i*(N-shx)+1:i*(N-shx)+N,j*(N-shy)+1:j*(N-shy)+N))).^2);
    end;
end;

PSD=PSD/V/((Mx+1)*(My+1));
WelchPSD=[];
WelchPSD=[PSD(N/2+1:N,N/2+1:N) PSD(N/2+1:N,1:N/2);PSD(1:N/2,N/2+1:N) PSD(1:N/2,1:N/2)];
ux=-1/2:1/N:1/2-(1/N);
ux=ux';
uy=ux;
figure,mesh(ux,uy,WelchPSD);
maxplot = max(max(WelchPSD))*1.1;
axis ([-(1/2) (1/2) -(1/2) (1/2) 0  maxplot]);
xlabel('u (c/px)'), ylabel('v (c/px)');
title('Power Spectrum Density,(from Welch-periodgram): S_f(u,v) ')
set (gca, 'XTick',[-(1/2) -(1/4) 0 (1/4) (1/2)]);
set (gca, 'XTickLabel', '-1/2|-1/4|0|1/4|1/2');
%set(gca,'Fontname','Symbol');
set (gca, 'YTick',[-(1/2) -(1/4) 0 (1/4) (1/2)]);
set (gca, 'YTickLabel', '-1/2|-1/4|0|1/4|1/2');
set(gca,'ydir','reverse');  %make v axis point downward
colormap([0 0 0]);     % use black only;

dB = 20*log10(WelchPSD);
ground = -80;
dB (find(dB<ground)) = ground;
maxplot = ceil(max(max(dB))/10)*10;
figure,mesh(ux,uy,dB);
axis ([-(1/2) (1/2) -(1/2) (1/2) ground  maxplot]);
xlabel('u (c/px)'), ylabel('v (c/px)');zlabel('dB');
%title('Power Spectral Density (PSD) of F_{2}(u,v)')
set (gca, 'XTick',[-(1/2) -(1/4) 0 (1/4) (1/2)]);
set (gca, 'XTickLabel', '-1/2|-1/4|0|1/4|1/2');
%set(gca,'Fontname','Symbol');
set (gca, 'YTick',[-(1/2) -(1/4) 0 (1/4) (1/2)]);
set (gca, 'YTickLabel', '-1/2|-1/4|0|1/4|1/2');
set(gca,'ydir','reverse');  %make v axis point downward
colormap([0 0 0]);     % use black only

%Contour plot (dB)
vlev=ground:10:maxplot;             % contours will be in steps of 10 dB
figure, [Con,hval]=contour(ux,uy,dB,vlev); % generate the contour plot, including values to label contours
axis square;  %make the plot square
xlabel('Horizontal frequency (c/px)'), ylabel('Vertical frequency (c/px)')
clabel(Con,hval)             %label the contours
title('Contour plot of power spectral density')
set(gca,'ydir','reverse');  %make v axis point downward
colormap([0 0 0]);     % use black only


%   This section is to interpolate the spectrum for better resolution
%   display

% inter = 4;
% vv=-1/(2*X):(1/X)/N/inter:1/(2*X)-((1/X)/N/inter);
% uu=-1/(2*X):(1/X)/N/inter:1/(2*X)-((1/X)/N/inter);
% vv=vv';
% uu=uu';
% [Uy,Ux] = meshgrid(uy,ux);
% Sf = interp2(Uy,Ux,dB,vv',uu,'spline');
% figure,mesh(uu,vv,Sf);
% axis ([-(1/2/X) (1/2/X) -(1/2/X) (1/2/X) ground-5  maxplot]);
% xlabel('u (c/ph)'), ylabel('v (c/ph)');zlabel('dB');
% title('Power Spectral Density (PSD) of f_{2}(u,v)')
% set (gca, 'XTick',[-(1/2/X) -(1/4/X) 0 (1/4/X) (1/2/X)]);
% set (gca, 'XTickLabel', '-1/2X|-1/4X|0|1/4X|1/2X');
% %set(gca,'Fontname','Symbol');
% set (gca, 'YTick',[-(1/2/X) -(1/4/X) 0 (1/4/X) (1/2/X)]);
% set (gca, 'YTickLabel', '-1/2X|-1/4X|0|1/4X|1/2X');
% set(gca,'ydir','reverse');  %make v axis point downward
% colormap([0 0 0]);     % use black only