function [f varargout]=irregreconAp(fr,xr,Ap,N,M)
%
% Reconstruct a periodic signal from irregular samples with variable
% aperture functions  Each sample is from a (possibly different)
% aperture-filtered version of the signal. 
%
%  Fy=irregreconAp(fr,xr,Ap,N,M)
%
%  [Fy fu]=irregreconAp...
%  [Fy fu xu]=irregreconAp...
%  [Fy fu xu Do]=irregreconAp...
%  [Fy fu xu Do D]=irregreconAp...
%  [Fy fu xu Do D K]=irregreconAp...
%
% Discrete irregular reconstruction (with variable apertures)
%
% inputs:
%  fr = column vector of irregular samples of length R
%  xr = column vector of (center) locations of irregular samples (0..N-1)
%  Ap = R x 2*L+1 matrix of aperture functions.  Each row is the aperture
%       function of the corresponding sample in fr and is centered at
%       L+1. L can be 0, i.e. a single point.  Each aperture function
%       should be "sampled" at the high rate given N and d
%       Thus, fr(i)=sum_{n=1}^{2*L+1} Ap(i,n)*Fy(n-L+xr(k))
%       where Fy is the true signal to be reconstructed expressed at
%       the high rate given by N and d.  The sum is computed
%       modulo N, the true period.  Each row of Ap should be normalized 
%       so that sum_{n=1}^{2*L+1} Ap(i,n)=1
%  N  = signal period (N>L>R)
%  M =  signal bandlimit (2*M+1<=R & 2*M+1<=N)
%
% outputs:
%  Fy = N values of one period of reconstructed signal
%  fu = (optional) 2*M+1 samples of equivalent uniform samples
%  xu = (optional) 2*M+1 locations of equivalent uniform samples (0..N-1)
%  Do = (optional) matrix relating irregular samples to equiv uniform
%        samples fr=Do*fu
%  D  = (optional) matrix relating equivent uniform samples to output
%        signal Fy=D*fu
%  K  = (optional) condition number of Do
%
% Basic theory:  Given periodic discrete signal Fy(n) with period N
% which is bandlimited to FFT(Fy)(k)=0 for |k|>M, Fy(n) can be
% reconstructed using the formula
%
% Fy(n)=sum_{j=-M}^M fu(j) d D_M(n - j d,N)
%
% where d=N/(2M+1) and D_M(n,N) = sin(pi (2M+1) x/N) / sin(pi x/N) is
% the Dirichlet function.  Note that fu(j)=Fy(j d) if d is an integer
% (d>1) with N=d*(2M+1).
%
% Given R (R>2*M+1) samples of Fy(n) at locations xr (xr in 0...N),
% then fu=pinv(Do)*fr and Fy=D*fu where Do and D are matrices of
% sampled Dirichlet functions.  The k,l entry of the R by 2*M+1 matrix
% Do is d D_M(xr(k)-(k-1)*d,N) while the k,l entry of the 
% N by 2*M+1 matrix D is d D_M((l-1)-(k-1)*d,N).  The matrix D, in 
% effect does sinc (actually Dirichlet) interpolation. Note that when d
% is not an integer, the locations of the equivalent uniform samples
% are not integers.  When d is an integer and the sampling is regular
% Do is an identity matrix.
%
% see Long and Franz, Band-limited signal reconstruction from irregular
% samples with variable apertures, IEEE Trans. Geosci. Remote Sens., 2015
%
% Note: While theory is exact, the matrix Do can be poorly conditioned 
% for large N,M and poor irregular sampling, resulting in reconstruction
% errors.
%

% written by D.G. Long at BYU  30 Dec 2006
% revised by D.G. Long at BYU  31 Dec 2014

if length(fr) ~= length(xr)
  error('Number of samples and their location must be the same length in irregrecon');
end

[R, L2]=size(Ap);     % number of samples (R) and aperture length
L=floor((L2-1)/2);
if L<1
  L=1;
end

d=N/(2*M+1);          % equivalent uniform sample spacing, must be >= 0
if d<1
  error(sprintf('d=N/(2*M+1) must be greater than one but N=%d M=%d so d=%f',N,M,d));
end
if R<2*M+1
  error(sprintf('R=%d is less than 2*M+1 where M=%d, reconstruction impossible',R,M));
end

%fprintf('Equivalent samples: d=%f  N=%d  M=%d  R=%d  minR=%d\n',d,N,M,R,2*M+1);

x=[0:N-1]';
xu=d*[0:2*M]';  % locations of estimated uniform samples
                % will not be integers unless d is an integer!

% generate Do matrix with variable aperture (converts irregular
% aperture-filtered samples to uniform samples w/o aperture)
Do=zeros([R 2*M+1]);
for c=1:R
  Do(c,:)=d*ApDm(-xr(c)+xu,M,N,Ap(c,:));
end
% alternately, here is a growing array method to make Do
%Do=[];
%for c=1:R
%  Do=[Do; d*ApDm(-xr(c)+xu,M,N,Ap(c,:))];
%end

% compute equivalent uniform samples
fu=Do\fr;           % same as pinv(Do)*fr but uses matlab's solver which
                    % is faster and more accurate

% get number of output arguments
Nout=nargout-1;

if Nout<=3 % D matrix is not output, so don't compute it
  
  % compute f=D*fu without explicitly computing D
  % Note that apertures are not involved with this step.
  % Note: potentially this could be done more efficiently using FFTs
  f=zeros([N 1]);
  for c=-M:M
    f=f+d*Dm(x-xu(c+M+1),M,N)*fu(c+M+1);
  end

else       % D matrix is included in output 

  % generate the general D matrix which reconstructs signal from
  % uniform samples.  Note that apertures are not involved with this step.
  % generate the general D matrix, then compute f=D*fu
  D=zeros([N 2*M+1]);
  for c=-M:M
    D(:,c+M+1)=d*Dm(x-xu(c+M+1),M,N)';
  end
  
  % explicit definition of general D matrx
  %D=zeros([N 2*M+1]);
  %for k=1:N
  %  for l=1:2*M+1
  %    D(k,l)=d*Dm((k-1)-(l-1)*d,M,N);
  %  end
  %end
  % Here is a growing array method for computing D
  %D=[];
  %for c=-M:M
  %  D=[D d*Dm(x-xu(c+M+1),M,N)'];
  %end

  % reconstruction full signal from equivalent uniform samples
  f=D*fu;     % reconstructed signal
end

% output optional arguments
if Nout>0
  varargout(1)={fu};
  if Nout>1
    varargout(2)={xu};
    if Nout>2
      varargout(3)={Do};
      if Nout>3
	varargout(4)={D};
	if Nout>4
	  varargout(5)={cond(Do)};
	end
      end
    end
  end
end

