function [Q, Qm, Qn]=psd2(X, DIM1, DIM2, DIM3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Q, Qm, Qn]=psd2(X, DIM1, DIM2, DIM3)
% Power Spectrum Density 2 (function)
%	Using Bartlett's method of averaging periodograms, 
%       this function averages several samples of a 2D signal X
%       to estimate the power spectrum.  These samples are 
%       extracted from the 2D signal according to the rectangular
%       window function specified by DIM1, DIM2 and DIM3.
%
% INPUT ARGUMENTS:
%	X    -> 2D signal
%       DIM1 -> [M N] where M is the number of rows and N the 
%               number of columns of the sample window.  The
%               default value is [256 256], and can be empty.
%       DIM2 -> [O P] where O is the number of rows and P the
%               number of columns used by FFT2 to transform the
%               sample window to the spectral domain.  The
%               default value is [256 256], and can be empty.
%       DIM3 -> [R S] where R is the number of rows and S the
%               number of columns of overlap between samples of
%               the original signal X.  The default value is 
%               [0 0], and can be empty.
%
% OUTPUT ARGUMENT:
%	P    -> The estimate power spectrum density of the 2D
%               signal X.
%       Pm   -> The spectral frequencies corresponding to the 
%               rows of P.
%       Pn   -> The spectral frequencies corresponding to the 
%               columns of P.
%
% EXAMPLE:
%	Y = halftone(ones(600, 600), 'error', [0 0 0; 0 0 7; 1 5 3]/16);
%       [P, Pm, Pn]=psd2(Y, [256 256], [256 256], [128 128]);
%       P(1)=0;
%       imagesc(Pm, Pn, P);
%       axis image;
%
% Daniel Leo Lau
% lau@ee.udel.edu
%
% October 12, 1998
% Copyright 1998 Daniel Leo Lau
%
% Last modified on October 12, 1998
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin<2) DIM1=[256 256]; end;
if (nargin<3) DIM2=[256 256]; end;
if (nargin<4) DIM3=[0 0];     end;
if (isempty(DIM1)) DIM1=[256 256]; end;
if (isempty(DIM2)) DIM2=[256 256]; end;
if (isempty(DIM3)) DIM3=[0 0];     end;

if (length(DIM1)==1)
     M=DIM1;
     N=DIM1;
else
     M=DIM1(1);
     N=DIM1(2);
end;

if (length(DIM2)==1)
     O=DIM2;
     P=DIM2;
else
     O=DIM2(1);
     P=DIM2(2);
end;

if (length(DIM3)==1)
     R=DIM3;
     S=DIM3;
else
     R=DIM3(1);
     S=DIM3(2);
end;

[Xm, Xn]=size(X);

dn=N-S;
dm=M-R;
Dn=floor((Xn-N)/dn)+1;
Dm=floor((Xm-M)/dm)+1;
Q=zeros(O,P);
j=0;
for m=1:Dm
     for n=1:Dn
          hm=(m-1)*dm;
          hn=(n-1)*dn;
          H=double(X(hm+1:hm+M,hn+1:hn+N));
          H=abs(fft2(H, O, P)).^2/(M*N);
          Q=Q+H;
          j=j+1;
     end;
end;
Q(1)=0;
Q=fftshift(Q/j);
if (nargout>1)
     [Qm, Qn]=freqspace(DIM2);
end;
return;