function [xn fu conver]=ACT_byu(xp,xs,M,N,iter,iflag)
%
%  [xn fu conver]=act_byu(xp,xs,M,N,iter,iflag);
%
% outputs:
%  xn = reconstructed signal of length N
%  fu = fft of xn
%  conver = convergence sequence (blank if iflag==0)
%
% inputs:
%  xp = sample locations within range (0..N-1)
%  xs = sample values
%  M  = signal bandlimit
%  N  = signal length
%  iter = <optional> restrict number of iterations (iter <= 2*M+1)
%  iflag = <optional> output conver if set to 1 and print condition
%          number of matrix. usage requires specification of iter
%

% code is based on the original ACT.m code by T. Strohmer , Nov. '92 
%
% modified by DGL at BYU  2 Feb 2013
%  T. Strohmer's original code is modified to use Long and Franz's
%  M-bandlimit and to output the fft of the reconstructed signal and 
%  (optionally) the sequence of convergence errors.  Initialization
%  and reconstruction steps slightly modified to improve numerical accuracy

xp=xp(:).'+1; xs=xs(:).'; % insure input in desired form and range
yp1=[N-M+1:N 1:M+1]; % frequecy bin order

if exist('iter','var')==0 % optional arg
  iter=2*M+1;
end
if exist('iflag','var')==0 % optional arg
  iflag=0;
end

lxp=length(xp); xa=zeros(1,N); xe=zeros(1,N); t=zeros(1,N);
lyp=length(yp1); yn=zeros(1,lyp); bb=zeros([1 lyp]);
xp1=[xp(2:lxp), N+xp(1)];
xp2=[xp(lxp)-N, xp(1:(lxp-1))];
w=(xp1-xp2)/2;
xe(xp)=w; 
ye=fft(xe);
t=ye(1:lyp)/N; % first row of the Toeplitz matrix
zeropad=2^(nextpow2(2*lyp-1))-2*lyp+1;
c=[t, zeros(1,zeropad), conj(t(lyp:-1:2))]; % first row of circulant
                                            % embedding of Toeplitz matrix
fc=fft(c); % preparation for fast multiplication

if 1==0 %% original Strohmer initialization
  xe(xp)=xs.*w; ye=fft(xe);
  r=ye(yp1);  % start value for CG algorithm
else    %% BYU initialization
  for k=1:lyp
    bb(k)=sum(xs.*w.*exp(-2*pi*j*(k-M-1)*(xp-1)/N));
  end
  r=bb; % start value for CG algorithm

  %%check
  %xe2=zeros(1,N);xe2(xp)=xs.*w; ye2=fft(xe2); r2=ye2(yp1);
  %[r-r2]
  
end

if iflag % optionally keep convergence error sequence
  conver=zeros([1 iter]);
else
  conver=[];
end

p=r;
rrnew=r*r';
for jj=1:iter;
  pz=[p zeros(1,lyp-1+zeropad)];
  fp=fft(pz);
  dummy=ifft(fc.*fp);
  Ap=dummy(1:lyp);
  alpha=rrnew/(p*Ap');
  if iflag % explicity compute convergence error
    aap=alpha*p;
    yn=yn+aap;
    conver(jj)=(aap)*(aap)';
  else     % faster not to compute conver
    yn=yn+alpha*p;
  end
  r=r-alpha*Ap;
  rrold=rrnew;
  rrnew=r*r';
  beta1=rrnew/rrold;
  p=r+beta1*p;
end;

if 1==0 %% original Strohmer reconstruction
  ye1=zeros(1,N); ye1(yp1)=yn; xn=ifft(ye1).'; 
  fu=ifftshift(yn);
else % alternate reconstruction approach
  ye=zeros([1 N]); 
  fu1=ifftshift(yn);
  ye(1:M+1)=fu1(1:M+1);
  ye(N-M+1:end)=fu1(M+2:end);
  xn=ifft(ye).';
  fu=fu1;
end

if iflag % create full matrix from first row and compute and 
         % print condition number
  Tw=toeplitz(t);
  disp(sprintf('ACT Cond:  %f',cond(Tw)))
end
