function [y0,y1,th0,th1] = autophas(x);
% Autophase performs automatic zero and first order phase correction
% on the complex spectrum in x.  The criterion is to maximize the sum of the real 
% components of x.   y0 is the spectrum with zero order correction applied, y1 is the                 
% corresponding phase angle and phase coefficient.   
%

% Roger J. McNichols
% Texas A&M University
% 12-19-97

n = length(x);
x = x(:);

% zero order correction...
a0 = sum((real(x)));
b0 = sum((imag(x)));
th0 = atan(-b0/a0);

y0 = x*exp(j*th0);


% First order correction...
w = [linspace(0,2*pi,n)]';
a1 = sum(real(y0).*cos(w));
b1 = sum(real(y0).*sin(w));
c1 = sum(imag(y0).*sin(w));
d1 = sum(imag(y0).*cos(w));
th1 = atan((a1-d1)/(b1+c1));

y1 = y0.*exp(j*w*th1);

y = y1;

% That's all folks!