PolyVecNormalisation

u = PolyVecNormalisation(v,rho); 

  This function normalises a polynomial vector v to unit length, such that
  u^P(z).u(z) = 1. This is equivalent to normalising a constant vector to
  a Euclidean length of unity in the non-polynomial case.

  The result vector is trimmed to remove a fraction rho of the energy before
  returning u, in order to limit the order of this vector.

  Implemented in the DFT domain, the chosen DFT size exceeds 2^N times the 
  order of v, with a default of N=4. This order increase is for internal 
  precision, and the result is trimmed subsequently.
  
  Input parameter:
       v           Mx1xL matrix containing a polynomial vector of length L
       rho         fraction of energy to be trimmed for return variable
       N           internal order increase is 2^N (default: N=4)
       epsilon     small constant to avoid division by zero (default 2.2e-16)

  Output parameter:
       u           Mx1xJ matrix of normalised polynomial vector

0001 function u = PolyVecNormalisation(v,rho,N,epsilon);
0002 %u = PolyVecNormalisation(v,rho);
0003 %
0004 %  This function normalises a polynomial vector v to unit length, such that
0005 %  u^P(z).u(z) = 1. This is equivalent to normalising a constant vector to
0006 %  a Euclidean length of unity in the non-polynomial case.
0007 %
0008 %  The result vector is trimmed to remove a fraction rho of the energy before
0009 %  returning u, in order to limit the order of this vector.
0010 %
0011 %  Implemented in the DFT domain, the chosen DFT size exceeds 2^N times the
0012 %  order of v, with a default of N=4. This order increase is for internal
0013 %  precision, and the result is trimmed subsequently.
0014 %
0015 %  Input parameter:
0016 %       v           Mx1xL matrix containing a polynomial vector of length L
0017 %       rho         fraction of energy to be trimmed for return variable
0018 %       N           internal order increase is 2^N (default: N=4)
0019 %       epsilon     small constant to avoid division by zero (default 2.2e-16)
0020 %
0021 %  Output parameter:
0022 %       u           Mx1xJ matrix of normalised polynomial vector
0023 %
0024         
0025 %  S. Weiss, UoS, 27/7/2022
0026 
0027 if nargin<3,
0028    N = 4;
0029    epsilon = eps;
0030 elseif nargin<4,
0031    epsilon=eps;
0032 else
0033    % do nothing
0034 end;   
0035 %epsilon = 1e-16;
0036 %epsilon
0037 
0038 %------------------------------------------------------------------------------
0039 % determine normalisation function;
0040 %------------------------------------------------------------------------------
0041 rvv = squeeze(PolyMatConv(ParaHerm(v),v));
0042 Lrvv = length(rvv);
0043 Nfft = 2^(ceil(log2(Lrvv))+4);
0044 Rvv = sqrt(abs(fft(rvv,Nfft))); 
0045 %------------------------------------------------------------------------------
0046 % perform normalisation in the DFT domain
0047 %------------------------------------------------------------------------------
0048 U = fft(v,Nfft,3);
0049 for k = 1:Nfft,
0050    if Rvv(k) <= epsilon
0051       U(:,:,k) = 0;
0052    else   
0053       U(:,:,k) = U(:,:,k)/Rvv(k);
0054    end;   
0055 end;
0056 u = ifft(U,Nfft,3);
0057 
0058 %------------------------------------------------------------------------------
0059 % allow for non-causal solution
0060 %------------------------------------------------------------------------------
0061 dummy = u;
0062 u(:,:,1:Nfft/2) = dummy(:,:,Nfft/2+1:Nfft);
0063 u(:,:,Nfft/2+1:Nfft) = dummy(:,:,1:Nfft/2);
0064    
0065 %------------------------------------------------------------------------------
0066 % trimming
0067 %------------------------------------------------------------------------------
0068 E = zeros(Nfft,1);
0069 for k = 1:Nfft,
0070    E(k) = norm(u(:,:,k)).^2;
0071 end;         
0072 Etotal = sum(E);    
0073 % now take iteratively from the front or the end of u until a fraction rho of
0074 % the energy has been removed
0075 E = E/Etotal;
0076 Trim=1; StartIndex=1; EndIndex=Nfft; ClippedEnergyFraction=0;
0077 while (Trim==1),
0078   if E(StartIndex)<E(EndIndex),
0079      if (ClippedEnergyFraction+E(StartIndex))<=rho,
0080         ClippedEnergyFraction = ClippedEnergyFraction+E(StartIndex);
0081         StartIndex=StartIndex+1;
0082      else
0083         Trim=0;
0084      end;
0085   else
0086      if (ClippedEnergyFraction+E(EndIndex))<=rho,
0087         ClippedEnergyFraction = ClippedEnergyFraction+E(EndIndex);
0088         EndIndex=EndIndex-1;
0089      else
0090         Trim=0;
0091      end;
0092   end;
0093 end;
0094 u = u(:,:,StartIndex:EndIndex);                
0095 
0096