PUPolyMatTrim

[H2,rho] = PUPolyMatTrim(H1,gamma)
 
  H2=PUPolyMatTrim(H1) trims the time dimension of a paraunitary matrix by
  removing outer matrix coefficients. The trimming is performed from both
  ends, and will remove the matrix coefficient with the smallest Frobenius
  norm, such that the total energy of the trimmed parts does not exceed  
  one permille of the total energy in H1.
 
  H2=PUPolyMatTrim(H1,gamma) trims such that the ratio between the removed 
  and total energy of R1 is less than gamma, with 0<=gamma<1.

  [H2,rho]=PUPolyMatTrim(H1) or [H2,rho]=PUPolyMatTrim(H1,r)
  additionally returns the error in paraunitarity --- the norm of 
  (H2*~H2-I) --- after trimming in rho.
  
  This function is based on the 'trim' operation described in [1,2].

  Input parameters:
     H1      paraunitary matrix
     gamma   maximum ratio between the maximum energy removed at outer lags 
             and the total energy 
             default: 1/1000
  
  Output parameters:
     H2      trimmed parahermitian matrix
     rho     error in paraunitarity

  References:
  [1] J.G. McWhirter, P.D. Baxter, T. Cooper, S. Redif, and J. Foster, "An EVD 
      Algorithm for Para-Hermitian Polynomial Matrices," IEEE Transactions on 
      Signal Processing, vol. 55, no. 5, pp. 2158-2169, May 2007.
  [2] C.H. Ta and S. Weiss, "Shortening the Order of Paraunitary Matrices in  
      SBR2 Algorithm", 6th International Conference on Information, Communi-
      cations & Signal Processing, Singapore, pp. 1-5, Dec. 2007.

0001 function [H2,rho] = PUPolyMatTrim(H1,gamma)
0002 %[H2,rho] = PUPolyMatTrim(H1,gamma)
0003 %
0004 %  H2=PUPolyMatTrim(H1) trims the time dimension of a paraunitary matrix by
0005 %  removing outer matrix coefficients. The trimming is performed from both
0006 %  ends, and will remove the matrix coefficient with the smallest Frobenius
0007 %  norm, such that the total energy of the trimmed parts does not exceed
0008 %  one permille of the total energy in H1.
0009 %
0010 %  H2=PUPolyMatTrim(H1,gamma) trims such that the ratio between the removed
0011 %  and total energy of R1 is less than gamma, with 0<=gamma<1.
0012 %
0013 %  [H2,rho]=PUPolyMatTrim(H1) or [H2,rho]=PUPolyMatTrim(H1,r)
0014 %  additionally returns the error in paraunitarity --- the norm of
0015 %  (H2*~H2-I) --- after trimming in rho.
0016 %
0017 %  This function is based on the 'trim' operation described in [1,2].
0018 %
0019 %  Input parameters:
0020 %     H1      paraunitary matrix
0021 %     gamma   maximum ratio between the maximum energy removed at outer lags
0022 %             and the total energy
0023 %             default: 1/1000
0024 %
0025 %  Output parameters:
0026 %     H2      trimmed parahermitian matrix
0027 %     rho     error in paraunitarity
0028 %
0029 %  References:
0030 %  [1] J.G. McWhirter, P.D. Baxter, T. Cooper, S. Redif, and J. Foster, "An EVD
0031 %      Algorithm for Para-Hermitian Polynomial Matrices," IEEE Transactions on
0032 %      Signal Processing, vol. 55, no. 5, pp. 2158-2169, May 2007.
0033 %  [2] C.H. Ta and S. Weiss, "Shortening the Order of Paraunitary Matrices in
0034 %      SBR2 Algorithm", 6th International Conference on Information, Communi-
0035 %      cations & Signal Processing, Singapore, pp. 1-5, Dec. 2007.
0036 
0037 % S. Weiss and J. Corr, University of Strathclyde, 14/11/2014
0038 
0039 % check input parameters
0040 if nargin<2,
0041    gamma = 0.001;
0042 end;
0043 
0044 % check limit
0045 E = gamma*PolyMatNorm(H1);
0046 [M,N,L] = size(H1);
0047 
0048 % trim from sequentially from front or back, where ever the matrix with the
0049 % smallest norm can be found, until the desired amount of energy has been
0050 % removed
0051 Etrim = 0;
0052 Index_front = 1; Index_back = L;
0053 Efront = norm(H1(:,:,Index_front),'fro').^2;
0054 Eback = norm(H1(:,:,Index_back),'fro').^2;
0055 while (Etrim + min([Efront Eback])) <= E,
0056    if Eback > Efront,
0057       % trim at the front
0058       Etrim = Etrim + Efront;
0059       Index_front = Index_front + 1;
0060       Efront = norm(H1(:,:,Index_front),'fro').^2;
0061    else
0062       % trim at the back
0063       Etrim = Etrim + Eback;
0064       Index_back = Index_back - 1;
0065       Eback = norm(H1(:,:,Index_back),'fro').^2;
0066    end;
0067 end;    
0068 
0069 % output
0070 H2 = H1(:,:,Index_front:Index_back);
0071 
0072 % prepare additional output parameter if error in paraunitarity is requested
0073 if nargout>1,
0074    HH2 = PolyMatConv(H2,ParaHerm(H2));
0075    L2 = (size(HH2,3)+1)/2;
0076    HH2(:,:,L2) = HH2(:,:,L2)-eye(N);
0077    rho = PolyMatNorm(HH2);
0078 end;