S = PHPolyMatInv(R,T,MaxIter,theta,epsilon);
PHPolyMatInv() calculates the inverse of parahermitian polynomial matrix as
described in [1].
S = PHPolyMatInv(R) returns the inverse of R. The input R represents an
MxMx(2L+1) parahermitian matrix R(z) of the form
R(z) = RL'z^L + ... + R1' z + R0 + R1 z^{-1} + ... + RLz^{-L}
whereby
R(:,:,1) = RL';
...
R(:,:,L) = R1';
R(:,:,L+1) = R0;
R(:,:,L+2) = R1;
...
R(:,:,2*L+1) = RL;
The function returns a paraunitary matrix S(z) such that approximately
S(z) R(z) = R(z) S(z) = I. The polynomial matrix S(z) is returned in the
same format at R, with the lag parameter forming a 3rd dimension.
The inversion is based in a polynomial EVD of R, and an inversion of the
polynomial eigenvalues. The PEVD uses the sequential best rotation (SMD, [2])
algorithm, and S(z) is determined from a minimum mean square error inverse
detailed in [1].
S = PHPolyMatInv(R,T) sets the order of the inverse polynomial eigenvalues
to 2T (default: T=10*L).
The SMD algorithm can be influenced by a number of optional inputs:
S = PHPolyMatInv(R,T,MaxIter) stops the SMD algorithm after MaxIter
iteration steps (default: MaxIter=100).
S = PHPolyMatInv(R,T,MaxIter,theta) stops the SMD algorithm either
after MaxIter iteration steps, or once the maximum off-diagonal element
falls below a threshold theta (default: theta=1e-5).
S = PHPolyMatInv(R,T,MaxIter,theta,epsilon) curtails the growth of
both the parahermitian and paraunitary matrices that arise from the PEVD
of R, by removing a proportion epsilon of the energy from these matrices.
Input parameters
R parahermitian matrix
T length of SISO inverses
MaxIter max. iterations for PEVD
(default = 200)
theta threshold for off-diagonal values
(default = 0.00001)
epsilon proportion of energy being trimmed to curtail polynomial order
Output parameter:
S inverse parahermitian matrix
References:
[1] S. Weiss, A. Millar, and R.W. Stewart: "Inversion of Parahermitian
Matrices,"European Signal Processing Conference, Aalborg, Denmark,
pp. 447-451, August 2010.
[2] S. Redif, S. Weiss, and J.G. McWhirter, "Sequential Matrix Diagonali-
sation Algorithms for Polynomial EVD of Parahermitian Matrices," IEEE
Transactions on Signal Processing, 63(1):81-89, January 2015.
0001 function S = PHPolyMatInv(R,T,MaxIter,theta,epsilon); 0002 %S = PHPolyMatInv(R,T,MaxIter,theta,epsilon); 0003 % 0004 % PHPolyMatInv() calculates the inverse of parahermitian polynomial matrix as 0005 % described in [1]. 0006 % 0007 % S = PHPolyMatInv(R) returns the inverse of R. The input R represents an 0008 % MxMx(2L+1) parahermitian matrix R(z) of the form 0009 % R(z) = RL'z^L + ... + R1' z + R0 + R1 z^{-1} + ... + RLz^{-L} 0010 % whereby 0011 % R(:,:,1) = RL'; 0012 % ... 0013 % R(:,:,L) = R1'; 0014 % R(:,:,L+1) = R0; 0015 % R(:,:,L+2) = R1; 0016 % ... 0017 % R(:,:,2*L+1) = RL; 0018 % The function returns a paraunitary matrix S(z) such that approximately 0019 % S(z) R(z) = R(z) S(z) = I. The polynomial matrix S(z) is returned in the 0020 % same format at R, with the lag parameter forming a 3rd dimension. 0021 % 0022 % The inversion is based in a polynomial EVD of R, and an inversion of the 0023 % polynomial eigenvalues. The PEVD uses the sequential best rotation (SMD, [2]) 0024 % algorithm, and S(z) is determined from a minimum mean square error inverse 0025 % detailed in [1]. 0026 % 0027 % S = PHPolyMatInv(R,T) sets the order of the inverse polynomial eigenvalues 0028 % to 2T (default: T=10*L). 0029 % 0030 % The SMD algorithm can be influenced by a number of optional inputs: 0031 % 0032 % S = PHPolyMatInv(R,T,MaxIter) stops the SMD algorithm after MaxIter 0033 % iteration steps (default: MaxIter=100). 0034 % 0035 % S = PHPolyMatInv(R,T,MaxIter,theta) stops the SMD algorithm either 0036 % after MaxIter iteration steps, or once the maximum off-diagonal element 0037 % falls below a threshold theta (default: theta=1e-5). 0038 % 0039 % S = PHPolyMatInv(R,T,MaxIter,theta,epsilon) curtails the growth of 0040 % both the parahermitian and paraunitary matrices that arise from the PEVD 0041 % of R, by removing a proportion epsilon of the energy from these matrices. 0042 % 0043 % Input parameters 0044 % R parahermitian matrix 0045 % T length of SISO inverses 0046 % MaxIter max. iterations for PEVD 0047 % (default = 200) 0048 % theta threshold for off-diagonal values 0049 % (default = 0.00001) 0050 % epsilon proportion of energy being trimmed to curtail polynomial order 0051 % 0052 % Output parameter: 0053 % S inverse parahermitian matrix 0054 % 0055 % References: 0056 % [1] S. Weiss, A. Millar, and R.W. Stewart: "Inversion of Parahermitian 0057 % Matrices,"European Signal Processing Conference, Aalborg, Denmark, 0058 % pp. 447-451, August 2010. 0059 % [2] S. Redif, S. Weiss, and J.G. McWhirter, "Sequential Matrix Diagonali- 0060 % sation Algorithms for Polynomial EVD of Parahermitian Matrices," IEEE 0061 % Transactions on Signal Processing, 63(1):81-89, January 2015. 0062 0063 % S.Weiss, UoS, updated 6/6/2015 0064 0065 %------------------------------------------------ 0066 % check for optional parameters 0067 %------------------------------------------------ 0068 if nargin<5, epsilon = 0.0; end; 0069 if nargin<4, theta = 0.00001; end; 0070 if nargin<3, MaxIter = 200; end; 0071 if nargin<2, T = 10*size(R,3); end; 0072 0073 %------------------------------------------------ 0074 % polynomial EVD 0075 %------------------------------------------------ 0076 [Q,Gamma] = SMD(R,MaxIter,theta,0.0); 0077 GammaInv = DiagInverse(Gamma,T,1e-5); 0078 S = PolyMatConv(ParaHerm(Q),PolyMatConv(GammaInv,Q)); 0079 0080 0081 0082 %------------------------------------------------ 0083 % subroutine DiagInverse() 0084 %------------------------------------------------ 0085 function Hinv = DiagInverse(D,T,gamma); 0086 % Inversion of a diagonal polynomial matrix D. D does not have 0087 % to be diagonal, but for the calculation of Dinv, only the on-diagonal 0088 % elements of D are considered. 0089 % 0090 % The inverse is calculated on the index interval [-T;+T]. 0091 % 0092 % The diagonal elements of D are truncated if there are tails with 0093 % elements smaller than the variable gamma. 0094 % 0095 % Input parameters: 0096 % D KxK polynomial matrix 0097 % T length of inverse (or lag of its autocorrelation sequence) 0098 % gamma element size for truncation 0099 % 0100 % Output parameters: 0101 % Dinv KxK diagonal matrix containing the inverse 0102 % polynomials of the on-diagonal elements of H. 0103 % 0104 % S. Weiss, 28/12/2006 0105 0106 [M,N,L] = size(D); 0107 Dinv = zeros(M,M,T); 0108 for m = 1:M, 0109 d = shiftdim(D(m,m,:)); 0110 ZeroPosition = find(abs(d)>gamma); % identify zeros 0111 d = d(ZeroPosition(1):ZeroPosition(end)); % truncate 0112 hinv = AcsMmseInv(d,T); 0113 Hinv(m,m,:) = hinv; 0114 end; 0115 0116 0117 %------------------------------------------------ 0118 % subroutine AcsMmseInv() 0119 %------------------------------------------------ 0120 function s = AcsMmseInv(r,T); 0121 % Inverts the response r which is assumed to have the 0122 % properties of an autocorrelation sequence. The inverse 0123 % is constrained to the same symmetry conditions as r, 0124 % with a length of 2T+1, i.e. corresponding to an inverse 0125 % autocorrelation sequence with maximum lag +/-T. 0126 % 0127 % Input parameters 0128 % r autocorrelation sequence 0129 % T maximum lag of inverse 0130 % 0131 % Output parameter 0132 % s inverse of r 0133 % 0134 % S. Weiss, univ of Strathclyde, 24/1/2010 0135 0136 %----------------------------------------------------- 0137 % setup matrix equation 0138 %----------------------------------------------------- 0139 T2 = (T+1)/2; 0140 R = toeplitz([r; zeros(T-1,1)],[r(1) zeros(1,T-1)]); 0141 R1 = R(:,1:T2) + fliplr(R(:,T2:end)); 0142 R2 = R(:,1:T2) - fliplr(R(:,T2:end)); 0143 R_real = [real(R1) -imag(R2); 0144 imag(R1) real(R2)]; 0145 d_real = zeros(2*(length(r)+T-1),1); 0146 d_real( (length(r)+T)/2 ) = 1; 0147 0148 %----------------------------------------------------- 0149 % pseudo-inverse 0150 %----------------------------------------------------- 0151 s_real = pinv(R_real)*d_real; 0152 0153 %----------------------------------------------------- 0154 % assemble inverse 0155 %----------------------------------------------------- 0156 s = [s_real(1:T2); zeros(T2-1,1)] + sqrt(-1)*[s_real(T2+1:2*T2); ... 0157 zeros(T2-1,1)]; 0158 s = s + conj(flipud(s)); 0159