% Reference: Coding Prony’s method in MATLAB and applying it to biomedical
% signal filtering, BMC Bioinformatics, 19, 451, 2018
%
% syntax: MatrixOut = LP1D(FID, LPalgorithm, LPdirection, NPP, NPU, NCOEFS)
%
% FID = matrix containing the data
% LPalgorithm=[1] linear predicition using SVD
%             [2] itmpm linear predicition using matrix pencil method
% LPdirection=[1] forward
%             [0] backward
% NPP = number of points to predict
% NPU = number of points of FID to use typically 0.75*N
% NCOEFS = number of signals in the spectrum / order of polynomial

function newFID = LP1D(FID, Options)

fig=uifigure;
  d = uiprogressdlg(fig,'Title','Linear Prediction in Progress',...
        'Indeterminate','on');
    drawnow;

 NCOEFS=Options(5);
 NPU=Options(4);
 NPP=Options(3);
 LPdirection=Options(2);
 LPalgorithm=Options(1);
 
 Size=length(FID);
 LPC = [];

  if (LPalgorithm == 1)	%lpsvd method
    if (NCOEFS < 1)
      disp('LP1D WARNING: number of frequencies cannot be smaller than 1 for lpsvd! Aborting ...')
      return;
    end

    LPC = lpsvd2(FID(1:NPU), NCOEFS);

  elseif (LPalgorithm == 2)			%itmpm method
    LPC = itcmp2(FID(1:NPU), NCOEFS);
  end


% Check whether the solution is sensible, return empty if not.
  if (length(LPC) < 1)
    disp('LP1D ERROR: algorithm did not find a proper result! Aborting ...')
    return
  else  
%% predict points and add to the spectrum
    if (LPdirection == 0)			%backward prediction
      newFID = compexp( (-NPP):(-1) , LPC);
      newFID = [newFID; FID];
    
    
    else				%forward prediction
      newFID = compexp( (Size):(Size+NPP-1) , LPC);%,1000);
      newFID = [FID; newFID];
    end
  end
  close(d)
  close(fig)
end 

% UTILITY FUNCTIONS 
function y=compexp(lp,para)
% y the predicted FID data points
% para LP coefficents from LPSVD or MPM
  lp=lp(:);
[M,~]=size(para);
ypure=zeros(size(lp));
ypure=ypure(:);
for ii=1:M
   ss=-para(ii,1)+sqrt(-1)*2*pi*para(ii,2);
   ypure=ypure+(para(ii,3)*exp(ss*lp+sqrt(-1)*para(ii,4))); % FID w/o noise     
end %ii

y=ypure;
end
function y=compexp2(lp,para,sn)

lp=lp(:);
[M,temp]=size(para);
delta=sqrt((sum(para(:,3).^2))/(2*10^(sn/10))); % delta: noise standard deviation
ypure=zeros(size(lp));
ypure=ypure(:);
for ii=1:M
   ss=-para(ii,1)+sqrt(-1)*2*pi*para(ii,2);
   ypure=ypure+(para(ii,3)*exp(ss*lp+sqrt(-1)*para(ii,4))); % FID w/o noise     
end %ii
randn('seed',100*randn(1));
w=delta*randn(size(ypure))*sqrt(-1); % noise
w=w+delta*randn(size(ypure));
y=ypure+w; % y: FID with noise    
varw=std(w).^2;

end % In case we need to include signal-to-noise parameter
function [para]=lpsvd2(y,M)
% linear prediction with singular value decomposition
% reference: R. Kumaresan, D. W. Tufts IEEE Trans. Acoust. Speech Signal Processing 
% vol. ASSP-30, 837-840, 1982. 

ScalingFactor = max(real(y)); 			

y=y(:) / ScalingFactor;

N=length(y);						
L=floor(N*3/4);						
A=hankel(conj(y(2:N-L+1)),conj(y(N-L+1:N)));		
h=conj(y(1:N-L));					

try
    
   [U, S, V] = svdecon(A); % Singular value decomposition 

    S = diag(S);
    bias = mean(S(M+1:min([N-L, L]))); 
    b = -V(:, 1:M) * (diag(1./(S(1:M)-bias)) * (U(:, 1:M)' * h)); 
    s = conj(log(roots([b(end:-1:1); 1]))); 
    s = s(real(s) < 0); 
  Z=zeros(N,length(s));   
  for k=1:length(s) 
    Z(:,k)=exp(s(k)).^[0:N-1].'; 
  end 
  a=Z\y;							
  para=[-real(s) imag(s)/2/pi abs(a) imag(log(a./abs(a)))]; 
 
 
  %% Check the values are reaonable
 % Check the values are reasonable
    
    % Identify rows with NaN values
    nanRows = any(isnan(para), 2);

    % Remove rows with NaN values
    para = para(~nanRows, :);

    % Protect against NaN in amplitude
    ampThreshold = sum(abs(para(:, 3))) / 1e4;
    para = para(para(:, 3) > ampThreshold & para(:, 3) > 1e-6, :);

    % Protect against inconsistent sign of real and imaginary parts
    consistentSign = para(:, 1) .* abs(sign(para(:, 2))) >= 0;
    para = para(consistentSign, :);

    para(:, 3) = para(:, 3) * ScalingFactor;


catch
  disp('SVD needed for linear prediction failed to converge!');
  para = [];
end
end 
function [para,M,itc]=itcmp2(y,M) 

ScalingFactor = max(real(y)); 			

y=y(:) / ScalingFactor;
N=length(y); 
L=floor(N/3);					
Y=toeplitz(y(L+1:N),y(L+1:-1:1));		

try
  [U,S,V]=svdecon(Y(:,2:L+1));	% takes a long time might be a faster way		

  S=diag(S);
  itc=zeros(1,L); 
  if M==-1					
    for k=0:L-1 
      itc(k+1)=-2*N*sum(log(S(k+1:L))) ... 
               + 2*N*(L-k)*log((sum(S(k+1:L))/(L-k))) + 2*k*(2*L-k); 
    end 
    [~, tempI]=min(itc); 
    M=tempI-1; 
  end
  
  if M==-2					
    for k=0:L-1 
      itc(k+1)=-N*sum(log(S(k+1:L))) + N*(L-k)*log((sum(S(k+1:L))/(L-k))) + k*(2*L-k)*log(N)/2;
    end
    [~, tempI]=min(itc); 
    M=tempI-1;
  end
  
  s=log(eig(diag(1./S(1:M))* ((U(:,1:M)'*Y(:,1:L))*V(:,1:M)))); 
  Z=zeros(N,M); 
  for k=1:M
    Z(:,k)=exp(s(k)).^[0:N-1].';
  end
  
  a=Z\y;	
  
  para=[-real(s) imag(s)/2/pi abs(a) imag(log(a./abs(a)))];
  
  
  % check the values are reasonable
  
  % remove rows that contain NaN
    nanRows = any(isnan(para), 2);
    para = para(~nanRows, :);

    % Protect against very small amplitudes
    ampThreshold = sum(abs(para(:, 3))) / 1e6;
    para = para(para(:, 3) > ampThreshold & para(:, 3) > 1e-6, :);

    % Refuse to accept negative dampening factors
    consistentSign = para(:, 1) .* abs(sign(para(:, 2))) >= 0;
    para = para(consistentSign, :);

    % Rescale the amplitude back
    para(:, 3) = para(:, 3) * ScalingFactor;

catch
  disp('SVD (needed for linear prediction) failed to converge!');
  para = [];
  M = [];
  itc = [];
end
end

function [U,S,V] = svdecon(X)
% Input:
% X : m x n matrix
%
% Output:
% X = U*S*V'
%
% Description:
% Does equivalent to svd(X,'econ') but faster
%
% Vipin Vijayan (2014)
%X = bsxfun(@minus,X,mean(X,2));
[m,n] = size(X);
if  m <= n
    C = X*X';
    [U,D] = eig(C);
    clear C;
    
    [d,ix] = sort(abs(diag(D)),'descend');
    U = U(:,ix);    
    
    if nargout > 2
        V = X'*U;
        s = sqrt(d);
        V = bsxfun(@(x,c)x./c, V, s');
        S = diag(s);
    end
else
    C = X'*X; 
    [V,D] = eig(C);
    clear C;
    
    [d,ix] = sort(abs(diag(D)),'descend');
    V = V(:,ix);    
    
    U = X*V; % convert evecs from X'*X to X*X'. the evals are the same.
    %s = sqrt(sum(U.^2,1))';
    s = sqrt(d);
    U = bsxfun(@(x,c)x./c, U, s');
    S = diag(s);
end
end 
