function Fy=recon(varargin)
%
% General routine to reconstruct a periodic signal from irregular samples with 
% possibly variable aperture functions  Each sample is from a (possibly 
% different) aperture-filtered version of the signal. 
%
%  Fy=recon(fr,xr,N,M)
%  Fy=recon(fr,xr,N,M,<method>,<Na,xa,Ap>,<alpha>,<ctrl>)
%
% Discrete irregular reconstruction (with variable apertures)
%
% required inputs:
%  fr = column vector of irregular samples of length R
%  xr = column vector of locations of length R of the
%        locations of the irregular samples within one period (0..N-1)
%  N  = signal period (N>L>R)
%  M  = signal bandlimit index (2*M+1<=R & 2*M+1<=N)
% optional inputs:
%  method = choice of method
%       0: matrix approach
%       1: pseudo-inverse matrix approach
%       2: ACT implemenation of FGS (only useable for Na<=1)
%       3: Long and Franz demonstration code
%  Na = Number of separate apertures (Na>=0)
%        Na=0 indicates no aperture (ideal delta function)
%  xa = index of aperture to sample number (vector of length R)
%  Ap = Na x 2*L+1 matrix of aperture functions.  Each row is an aperture
%       function centered at L+1. L can be 0, i.e. a single point.  
%       Each aperture function must be "sampled" at the rate given
%       by N and d where d=N/(2*M+1)
%       Thus, fr(i)=sum_{n=1}^{2*L+1} Ap(xa(i),n)*F(n-L+xr(k))
%       where F is the true signal to be reconstructed.  The sum is 
%       computed modulo N, the true period.  Each row of Ap should
%       be normalized so that 1=sum_{n=1}^{2*L+1} Ap(i,n)
% alph = regularization value (default=0) for methods 0 and 1
%        iteration conrol vlaue(default=2*M+1) for ACT for method 2
% ctrl = control flag
%        when non-zero, print additional information 
%
% outputs:
%  Fy = N values of one period of reconstructed signal
%
% Basic theory:  Given periodic discrete signal F(n) with period N
% which is bandlimited to FFT(F)(k)=0 for |k|>M, Fy(n) can be
% reconstructed using the formula
%
% Fy(n)=sum_{k=-M}^M a(j) exp(j 2 pi n k/N)
%
% where a(j) are the Fourier coefficients.
%
% Given R (R>=2*M+1) samples of Fy(n) at locations xr (X in 1..N),
% then fu=pinv(D)*fr and Fy=D*fu where D is a matrix of sampled
% sine waves, D(n,j)= exp(j 2 pi n k/N).  If each sample is actually
% a sample of the aperture filtered signal, i.e. 
% Fy_ap(n)=sum_{l=-L}^L Fy(n+l) Ap(n,l+L+1) then D is replaced
% by the element-wise product of D and the bandlimited Fourier
% tansform of Ap.
%
% Note: While theory is exact, the matrices can be poorly conditioned 
% for large N,M and poor irregular sampling, resulting in reconstruction
% errors.  If the Fourier transform of any of the rows of Ap have 
% zeros within the signal spectrum, the inverse may not exist.
%
% see Long and Franz, Band-limited signal reconstruction from irregular
% samples with variable apertures, IEEE Trans. Geosci. Remote Sens., 2015
%
% The following code implements the matrix formulation exactly, and
% has options to use a reduced-memory pseudo-matrix inverse.  When
% one or less apertures is used, FGS's ACT approach can be used. The
% ACT iterative approach uses less memory and can be numerically more
% efficient than the other methods, but is less acruate and cannot
% be used for multiple apertures.
%

% written by D.G. Long at BYU  15 Mar 2013
% written by D.G. Long at BYU  30 Dec 2014

Narg=nargin;
if Narg<4
  error('Minimum of 4 arguments required for recon');
end

% decode minimum argument list
fr=varargin{1};
xr=varargin{2};
N=varargin{3};
M=varargin{4};

% set defaults for other optional arguments
method=0;   % default method
Na=0;       % default: no apertures
Ap=[];      % empty Ap array
xa=[];      % empty pointer to Ap array
iter=2*M+1; % default number of iterations
ctrl=0;     % no extra printing
alpha=0;    % no regularization

% set internal variables and check for valid values
R=length(fr);
MM=2*M+1;
if R<MM
  error(sprintf('R=%d is less than 2*M+1 where M=%d, reconstruction impossible',R,M));
end

if Narg>4 % method specified
  method=varargin{5};
  if method<0 | method>3
    error('Invalid method specified');
  end
end

if Narg>5 % Na specified
  Na=varargin{6};  % number of apertures
  if Na<0
    error('Number of apertures must be 0 or more in recon')
  end
end

if Narg>6 % xa specified
  xa=varargin{7};  % xa input
  xa=xa(:); % force to column vector
  if length(xa)<R
    error('Number of aperture locations must be R or more');
  end
end

if Narg>7 % Ap specified
  Ap=varargin{8};  % Ap input
  if size(Ap,1)<Na
    error('Number of rows of Ap must be equal to Na');
  end
end

if Narg>8 % regularization value specified
  alpha=varargin{9};  % alpha input
end

if Narg>9 % ctrl specified
  ctrl=varargin{10};  % ctrl input
end


% compute useful constants
L2=size(Ap,2); % aperture input length
L=floor((L2-1)/2);
if L<1
  L=0;
  Na=0;
end

% check method specification
method_used=method;
if method==2 & Na>1 % override specified method
  fprintf('*** Method=2 specified, with Na>1, overriding method to use 0.  Regularization switched to default.\n');
  method_used=0;
elseif method==2 & alpha > 0
  iter=alpha;
end

% processing depends on method specification
switch method_used
 case 0, % A more efficient implementation of Long and Franz's matrix 
  % formulation structured to reduce computation
  % this method supports regularization and multiple apertures
  % Note: for efficiency for a single aperture, no-aperture
  % reconstruction is employed, followed by deconvolution of the aperture

  % frequency domain matrix formation
  w0=2*pi*j/N;
  Tw=zeros([R MM]);
  YP=[N-M+1:N 1:M+1]; % frequency ordering
  for k=1:R
    % row of matrix
    Tw(k,:)=exp(w0*(YP-1)*xr(k));
    % explicit definition
    %for jj=1:2*M+1
    %  Tw(k,jj)=exp(w0*(YP(jj)-1)*xr(k));
    %end
  end
  bw=fr;
  if ctrl>0
    fprintf('Condition number w/o apertures: %f\n',cond(Tw))  
  end
  
  if Na>1 % apply variable apertures to reconstruction matrix as a Hadamard product
    yp2=[1:L+1 N-L+1:N]; % frequency mapping
    yp3=[L+1:2*L+1 1:L]; % frequency mapping
    for k=1:R
      % compute zero-padded Fourier transform of each aperture function
      FAp1=zeros([1 N]);  % row vector    
      FAp1(yp2)=Ap(xa(k),yp3);
      FAp1=fft(FAp1);
      % apply aperture function
      Tw(k,:)=Tw(k,:).*FAp1(YP);
    end
    if ctrl>0
      fprintf('Condition number after apertures: %f\n',cond(Tw))
    end
  end
  
  if alpha(1)>0 % apply regularization (quite computationally expensive)
    if length(alpha)>=MM % weighted regularization
      T=zeros([R R]);
      for k=1:R
	T(k,k)=alpha(k);
      end
    else                % single-valued regularization
      T=alpha*eye(R);
    end
    Tw=(Tw*pinv(Tw)+T)*Tw;  % regularize Tw matrix
    clear T
    if ctrl>0
      fprintf('Regularization parameter %f\n',alpha)  
      fprintf('Cond after regularization %f\n',cond(Tw))
    end
  end
  
  % solve linear system to compute Fourier spectrum coefficients
  % (may generate a warning message is matrix is poorly conditioned)
  fu=Tw\bw;
  % could use instead
  % fu=pinv(Tw)*bw;
  
  if Na==1 % in single aperture case, apply inverse aperture
    if ctrl>0
      fprintf('Apply single aperture correction\n');
    end
    % compute Fourier transform of aperture
    FAp=zeros([N 1]);
    FAp(1:L+1)=Ap(1,L+1:2*L+1);
    FAp(N-L+1:end)=Ap(1,1:L);
    FAp=fft(FAp);
    fu=fu./FAp(YP); 
  end
  
  % reconstruct signal from Fourier spectrum
  FF=zeros([1 N]);
  fu1=ifftshift(fu);
  FF(1:M+1)=fu1(1:M+1);
  FF(N-M+1:end)=fu1(M+2:end);
  Fy=ifft(FF).'*N;
  
  % if input signal real, force real output signal to be real
  if isreal(fr)==1
    Fy=real(Fy);
  end  
  
  
 case 1, % a pseudo inverse approach: uses less memory, but slower
    % this method supports regularization and multiple apertures
    % Note: for efficiency for a single aperture, no-aperture
    % reconstruction is employed, followed by deconvolution of aperture
    
  if Na>1
    % pre-compute the Fourier transforms for the aperture functions
    YP=[N-M+1:N 1:M+1];
    FAp=zeros([size(Ap,1) N]);
    FAp(:,1:L+1)=Ap(:,L+1:2*L+1);
    FAp(:,N-L+1:end)=Ap(:,1:L);
    FAp=fft(FAp,[],2);
    FAp1=FAp(:,YP);
    clear FAp
  end
  if ctrl>0
    fprintf('apertures: %d\n',Na);
  end

  TT=zeros([MM MM]);
  Tbw=zeros([MM 1]); 
  w0=2*pi*j/N;
  YP=[N-M+1:N 1:M+1]; % frequency ordering

  if Na>1 % multiple apertures
   for jj=1:MM
     Tbw(jj)=fr.'*(exp(-w0*(YP(jj)-1)*xr).*FAp1(xa,jj));
   end
  else    % one or fewer
   for jj=1:MM
     Tbw(jj)=fr.'*exp(-w0*(YP(jj)-1)*xr);
   end
  end
    
  if Na>1
    for jj=1:MM
      for kk=1:MM
	TT(jj,kk)=(exp(w0*(YP(jj)-1)*xr).*FAp1(xa,jj))'*(exp(w0*(YP(kk)-1)*xr).*FAp1(xa,kk));
      end
    end    
  else
    for jj=1:MM
      for kk=1:MM
	TT(jj,kk)=exp(w0*(YP(jj)-1)*xr)'*exp(w0*(YP(kk)-1)*xr);
      end
    end
  end
  if ctrl>0
    fprintf('Condition number: %f\n',cond(TT))  
  end

  if alpha(1)>0 % apply regularization      
    if length(alpha)>=MM % weighted regularization
      for k=1:MM
	TT(k,k)=TT(k,k)+alpha(k);
      end
    else
      TT=TT+alpha*eye(MM);
    end
    if ctrl>0
      fprintf('Regularization parameter %f\n',alpha)  
      fprintf('Cond after regularization %f\n',cond(TT))  
    end
  end
  
  % solve linear system to compute Fourier spectrum coefficients
  % (may generate a warning message is matrix is poorly conditioned)
  fu=TT\Tbw;

  if Na==1 % in single aperture case, apply inverse aperture
    if ctrl>0
      fprintf('Apply single aperture correction\n');
    end
    % compute Fourier transform of aperture
    FAp=zeros([N 1]);
    FAp(1:L+1)=Ap(1,L+1:2*L+1);
    FAp(N-L+1:end)=Ap(1,1:L);
    FAp=fft(FAp);
    fu=fu./FAp(YP);
    clear FAp
  end
  
  % reconstruct signal from Fourier spectrum
  FF=zeros([1 N]);
  fu=ifftshift(fu);
  FF(1:M+1)=fu(1:M+1);
  FF(N-M+1:end)=fu(M+2:end);
  Fy=ifft(FF).'*N;
  clear FF

  % if input signal real, force real output signal
  if isreal(fr)==1
    Fy=real(Fy);
  end

 case 2, % ACT option
  % suports one or fewer apertures, does not support explicit regularization
  % (regularization can be done by fixing iter)
  % Note: for a single aperture, no-aperture
  % reconstruction is employed, followed by deconvolution of aperture

  if Na>1, error('ACT method does not support multiple apertures'); end;

  % use ACT method of FGS for no aperture irregular reconstruction
  %fprintf('call ACT_byu with: %d %d  %d %d  %d %d\n',M,N,length(xr),length(fr),iter,ctrl)
  [Fy Fu conver]=ACT_byu(xr,fr,M,N,iter,ctrl);
  if ctrl>0
    fprintf('ACT convergence:\n');
    conver
  end
  
  if Na==1   % single aperture compensation if aperture
    if ctrl>0
      fprintf('Apply single_aperture correction\n');
    end
    % compute Fourier transform of aperture
    FAp=zeros([N 1]);
    FAp(1:L+1)=Ap(1,L+1:2*L+1);
    FAp(N-L+1:end)=Ap(1,1:L);
    FAp=fft(FAp);
    tmp=fft(Fy);
    YP=[1:M+1 N-M+1:N];
    tmp(YP)=tmp(YP)./FAp(YP); 
    tmp(M+2:end-M)=0; % to reduce noise, set terms not needed to zero
    Fy=ifft(tmp); 
  end  

  % if input signal real, force real output signal
  if isreal(fr)==1
    Fy=real(Fy);
  end

case 3, % Long and Franz demonstration code
  % this method supports regularization and multiple apertures

  if Na>0   % one or more apertures
    if ctrl>0
      fprintf('use irregreconAp\n');
    end
    Fy=irregreconAp(fr,xr,Ap,N,M);
  else      % no aperture (i.e., ideal aperture)
    if ctrl>0
      fprintf('use irregrecon\n');
    end
    Fy=irregrecon(fr,xr,N,M);
  end
  
end
