%Function to compute a power density spectrum estimate of a single grayscale image
%using the Welch method of averaging modified periodograms.
%Eric Dubois 2017-12-11
function WelchPDS = WelchPDS_est(A,N)
% A is a grayscale image in double format
% N x N is the block size. N should be a power of 2. Usually N=64
% WelchPDS is the N x N power density spectrum estimate on a linear scale
% It is computed at normalized frequencies ux=-1/2:1/N:1/2-(1/N); vx=ux;
%
% Remove the mean from the image
[I, J]=size(A);
prodsize=I*J;
mean=sum(sum(A))/prodsize;
Az=A-mean;
% Compute the 2D separable window (4-term Harris-Nutall Window)
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 = sum(sum(W.*W));  % total power of the window
%
% Compute the spectral estimate using an x and y overlap of N/4
PDS=zeros(N,N);
shx=N/4;     % x overlap
shy=N/4;     % y overlap
% Number of blocks in x and y direction
Mx=floor((I-N)/(N-shx));
My=floor((J-N)/(N-shy));

for i=0:1:Mx
    for j=0:1:My
        PDS=PDS+(abs(fft2(W.*Az(i*(N-shx)+1:i*(N-shx)+N,j*(N-shy)+1:j*(N-shy)+N))).^2);
    end
end
%Compute the average
PDS=PDS/V/((Mx+1)*(My+1)); 
WelchPDS=[];
%Shuffle the FFT blocks to put the origin at the center
WelchPDS=[PDS(N/2+1:N,N/2+1:N) PDS(N/2+1:N,1:N/2);PDS(1:N/2,N/2+1:N) PDS(1:N/2,1:N/2)];