PUPolyMatRand

H = PUPolyMatRand(M,L,N,mode);

  H = PUPolyMatRand(M,L) generates a random MxM paraunitary matrix of order L.
  The matrix is assembled from L first order elementary paraunitary matrices
  described in [1],

  V(z) = I - u u' + u u' z^{-1}
 
  where u is an arbitrary unit norm vector. The resulting paraunitary matrix 
  H(z) is the product of L such elementary paraunitary matrices, s.t.
    H(z) = H0 + H1 z^{-1} + H2 z^{-2} + ... + HL z^{-L}
  is represented in a 3-dimensional matrix as
    H(:,:,1) = H0;
    H(:,:,2) = H1;
    H(:,:,3) = H2;
      ...
    H(:,:,L) = HL;
  Paraunitarity means that both PolyMatConv(H,ParaHerm(H)) and 
  PolyMatConv(ParaHerm(H),H) will result in an identity matrix.

  H = PUPolyMatRand(M,L,N) generates a random paraunitary matrix with a seed
  value of N. Repeated calls with the same N will results in the same para-
  unitary matrices.

  H = PUPolyMatRand(M,L,N,'real') generates a real-valyed random paraunitary 
  while matrix for H = PUPolyMatRand(M,L,N,'complex'), the output will be
  complex-valued.

  Input parameters:
     M       spatial dimension of paraunitary matrix
     L       polynomial order of paraunitary matrix
     N       seed value for random number generator
             (optional); default is random;
     mode    real- or complex valued operation (default is real)

  Output parameter:
     H       MxMx(L+1) paraunitary matrix

  Reference:
  [1] P.P. Vaidyanathan: "Multirate Systems and Filter Banks" Prentice
      Hall, 1993.

0001 function H = PUPolyMatRand(M,L,N,mode);
0002 %H = PUPolyMatRand(M,L,N,mode);
0003 %
0004 %  H = PUPolyMatRand(M,L) generates a random MxM paraunitary matrix of order L.
0005 %  The matrix is assembled from L first order elementary paraunitary matrices
0006 %  described in [1],
0007 %
0008 %  V(z) = I - u u' + u u' z^{-1}
0009 %
0010 %  where u is an arbitrary unit norm vector. The resulting paraunitary matrix
0011 %  H(z) is the product of L such elementary paraunitary matrices, s.t.
0012 %    H(z) = H0 + H1 z^{-1} + H2 z^{-2} + ... + HL z^{-L}
0013 %  is represented in a 3-dimensional matrix as
0014 %    H(:,:,1) = H0;
0015 %    H(:,:,2) = H1;
0016 %    H(:,:,3) = H2;
0017 %      ...
0018 %    H(:,:,L) = HL;
0019 %  Paraunitarity means that both PolyMatConv(H,ParaHerm(H)) and
0020 %  PolyMatConv(ParaHerm(H),H) will result in an identity matrix.
0021 %
0022 %  H = PUPolyMatRand(M,L,N) generates a random paraunitary matrix with a seed
0023 %  value of N. Repeated calls with the same N will results in the same para-
0024 %  unitary matrices.
0025 %
0026 %  H = PUPolyMatRand(M,L,N,'real') generates a real-valyed random paraunitary
0027 %  while matrix for H = PUPolyMatRand(M,L,N,'complex'), the output will be
0028 %  complex-valued.
0029 %
0030 %  Input parameters:
0031 %     M       spatial dimension of paraunitary matrix
0032 %     L       polynomial order of paraunitary matrix
0033 %     N       seed value for random number generator
0034 %             (optional); default is random;
0035 %     mode    real- or complex valued operation (default is real)
0036 %
0037 %  Output parameter:
0038 %     H       MxMx(L+1) paraunitary matrix
0039 %
0040 %  Reference:
0041 %  [1] P.P. Vaidyanathan: "Multirate Systems and Filter Banks" Prentice
0042 %      Hall, 1993.
0043 
0044 %  S. Weiss, University of Strathclyde, 14/12/14
0045 %       amended 31/3/23
0046 
0047 %-----------------------------------------------------
0048 %   check for optional seed value
0049 %-----------------------------------------------------
0050 if nargin>2,
0051    randn('seed',N);
0052 end;
0053 CmplxVld=0;
0054 if nargin>3,
0055    if strcmp(mode,'complex')==1,
0056       CmplxVld=1;
0057    end;
0058 end;  
0059       
0060 %-----------------------------------------------------
0061 %   generate L-th order paraunitary matrix
0062 %-----------------------------------------------------
0063 if CmplxVld==1,
0064    H=randn(M,M) + sqrt(-1)*randn(M,M);
0065 else,   
0066    H=randn(M,M);
0067 end;   
0068 [u,s,v] = svd(H);
0069 H = u*v';
0070 for i = 1:L,
0071    % generate elementary paraunitary matrix
0072    if CmplxVld==1,
0073       u = randn(M,2)*[1; sqrt(-1)];
0074    else
0075       u = randn(M,1);
0076    end;
0077    u = u/norm(u,2);
0078    U(:,:,1)=eye(M)-u*u';
0079    U(:,:,2)=u*u';
0080    % apply ith elementary PU matrix
0081    H = PolyMatConv(H,U);
0082  end;
0083