CC BY
27
DOI: 10.17516/1997-1397-2020-13-2-160-169 УДК 519.873
A Recursive Algorithm for Estimating the Correlation Matrix of the Interference Based on the QR Decomposition
Valery N. Tyapkin* Dmitry D. Dmitriev^ Andrey B. Gladyshev* Peter Yu. Zverev§
Siberian Federal University Krasnoyarsk, Russian Federation
Received 02.11.2019, received in revised form 10.12.2019, accepted 20.01.2020 Abstract. Many tasks of digital signal processing require the implementation of matrix operations in real time. These are operations of matrix inversion or solving systems of linear algebraic or differential equations (Kalman filter). The transition to the implementation of digital signal processing on programmable logic device (FPGAs), as a rule, involves calculations based on the representation of numbers with a fixed point. This makes solving spatio-temporal processing problems practically impossible based on conventional computational methods. The article discusses the implementation of spatial-temporal signal processing algorithms in satellite broadband systems using QR decomposition. The technologies of CORDIC computations required for recurrent QR decomposition when used together in systolic algorithms are presented.
Keywords: phased antenna array, adaptive algorithms, Kalman filter, recursive least squares algorithm (RLS), QR decomposition, systolic algorithm.
Citation: V.N.Tyapkin, D.D.Dmitriev, A.B.Gladyshev, P.Yu.Zverev, A Recursive Algorithm for Estimating the Correlation Matrix of the Interference Based on the QR Decomposition, J. Sib. Fed. Univ. Math. Phys., 2020, 13(2), 160-169. DOI: 10.17516/1997-1397-2020-13-2-160-169.
Introduction
Modern satellite broadband radio communication systems have a significant drawback - low noise immunity. The solution to this problem is based on the use of phased array antennas. Such antennas are controlled using adaptive algorithms, the parameters of which can be quickly changed in accordance with the emerging interference environment.
To work effectively in real conditions of parametric a priori uncertainty and a dynamic change in the statistical characteristics of interference, adaptive systems are required, the parameters of which can be quickly changed in accordance with the interference environment. Currently, theoretically substantiated and tested in practice, the methods of the Markov theory of optimal filtration. This theory has been fully and strictly developed in a number of books and articles
[email protected] https://orcid.org/0000-0001-6438-6094
§ [email protected] © Siberian Federal University. All rights reserved
[1-12]. In [6,10,11], the solution to the problem of adaptive filtering of signals based on the use of the lemma on the inversion of the correlation matrix (MIL) of input signals of an adaptive filter is considered. This solution leads to a recursive least squares (RLS) algorithm.
The same problem can also be solved by reducing the matrix of input samples of the adaptive filter signals to a triangular form. In this case, the range of numbers involved in the calculations is reduced by comparing the solution to this problem with the estimation of the inverse correlation matrix of interference using MIL. This increases the stability of QRD RLS algorithms.
Assume that the signal and interference affect the input of a multichannel M-element phased array antenna. The set of signals from the outputs of the M-element antenna array is described by the time function y1(t),y2(t),... ,ym(t) and form a column vector y1(t) = = [yi(t), y2(t),..., ym(t)]T. Moreover, a single-channel reception (M = 1) is considered as a special case of multi-channel.
Discretization of the received useful and interference signals is performed at the radio frequency f0. A feature of this is the small sampling interval At, which is approximately half the period of the carrier frequency Td = 1/fo, At « 1/2fo = To/2. Discrete interference values obtained from the antenna are random numbers that are conveniently represented as a column vector yi = [yi(kAt)] = [yi(k)], k = 1,L, where L determines the duration of the observation interval T — L = T/At. In the case of multichannel reception, the vector of received oscillations will have the following form Y = (Y1Y2... Yk ... YT)T and dimension (T x M).
Most QR decomposition algorithms are based on Householder reflection and Givens rotations [13]. For the implementation of space-time processing, the most useful is the recursive version of the Givens method, which provides updates to the solutions at the rate of arrival of the input samples of the signal. High real-time performance provides a systolic version of the QR algorithm using pipelined implementation of Givens rotations on FPGA. High speed fixed-point number calculations on FPGAs are provided by the CORDIC processor. The principle of its operation differs significantly from the arithmetic-logical devices of existing processors. To implement Givens rotation, 10 shift-addition operations are sufficient. In this case, an accuracy sufficient to achieve an interference suppression ratio of more than 50 dB is ensured.
Thus, the development of this research area promises a significant improvement in the quality of reception and processing of broadband signals and noise immunity based on the existing element base and is relevant. Consider the implementation of spatial-temporal signal processing algorithms using QR decomposition. To solve the problem of recurrent QR decomposition, we will develop CORDIC computing technology in systolic algorithms.
1. Recursive adaptation algorithm using QR decomposition
A recursive adaptation algorithm using QR decomposition estimates the filter coefficients at the current time step through the calculated filter coefficient at the previous step. Due to its recursive nature, the algorithm is called QRD — a recursive QR decomposition algorithm. QR factorization consists in reducing a linear system to a triangular one. For this, the original matrix is represented as the product of the upper triangular matrix R and the orthogonal matrix Q.
Consider a system of linear equations
Ax = b, (1)
where A — (n x m) is the matrix, x is the vector of derivatives (for example, weights for the adaptive antenna array), b is the m-vector. The QR decomposition of matrix A of size (n x m)
for any n ^ m can be described as:
A = Q • R, (2)
Q is a unitary matrix: Q • QH = I, where I is the identity matrix, R is (n x m) the upper right triangular matrix. Equation (2) can be written in divided form:
A = Q • R = Q
Ri 0
[Qi Q2]
Ri 0
Qi • Ri,
(3)
where Ri is the (m x m) triangular matrix, Qi is the (m x n) matrix and Q2 is the ((n — m) x m) matrix.
In the spatio-temporal processing of broadband signals, the system of equations is usually redefined because n ^ m. Solution (2) minimizing the norm of the residual || Ax — by has the form x = (ATA)-1ATb.
Substitution of equality A = QTR gives the following form
A = (RT QT QR)-1RT QT b = R-1QT b.
Here the triangular (m x m)-matrix is subject to circulation, which requires (m2)/2 operations of addition and multiplication. Let us synthesize a recursive algorithm for estimating the correlation matrix of interference based on QR decomposition From the theory of matrices it is known that there exists a unitary matrix Qk(k), that for any AkM matrix can be obtained by decomposition
FM
Qk AkM
0(T-M )M_
(4)
where FM is the upper triangular square matrix, called the Cholesky decomposition, 0(T-M)M is the zero rectangular matrix. The superscript H means complex conjugation and transposition. For the unitary matrices Qk(k), the equalities Qk(k) • Qf (k) = Qf (k) • Qk(k) = I and
QH (k) = Q-'(k).
If designated
where
AkM =A0k5(k)YMk (k),
(5)
A
0.5
diag{\/\k-1,\/\k-2, • • • VA1, l|
'VAk-1
0 0 0
0 0 00
0 . y/x 0 0l
(6)
The parameter A is used to weight the signals and allows you to take into account changes in the statistics of signals if they are non-stationary and their statistical parameters change over time. The parameter A is also called the exponential weighting parameter or the "orgetting factor" parameter. Its value is usually limited by the limits (1 — 0.4K) < A < 1 [14,15]. For example, at A = 0.9, A0 = 1, A1 = 0.9, A2 = 0.81, A3 = 0.72, A4 = 0.66, A5 = 0.59,..., A10 = 0.35,..., A20 = 0.12,..., A30 = 0.04,..., A40 = 0.015,..., A50 = 0.005,..., A100 = 0.000027, that is, the parameter A determines the "memory" of the algorithm for solving the recursive problem.
For 0 < A < 1 and small values of the difference k — i, the summed terms are weighted with large weights, and for large values of this difference, with smaller weights. For A = 1 this
0
0
k
"memory" is equal to k samples of the observed signals. For 0 < A < 1, the contribution of the same i-th samples to the sum is different for different A. This contribution is greater for large A and smaller for smaller A. That is, with decreasing A, the effective memory decreases and vice versa. Substituting equation (5) in (4) we obtain
QH ' (k)YMk (k)
' Rm (k) _0(t-m )m_
where
1 Mk (k)K ' 5(k) kM Qk
Qk(k)
" Rm (k) ' h
0(T-M )M_
RM (k), 0(T-M)M
(7)
(8)
Reducing the matrix RM(k) to a triangular form using the observation matrix YHk(k) can be performed in various ways, the main of which is Givens rotation. The matrix RM (k) can be obtained recursively in time, performing calculations for the k-th iteration using the results from the previous, (k — 1)-th iteration. This is as follows.
Let us assume that at iteration k-1 there is a decomposition:
QH-^k - 1)Ak'-i(k - 1)YM(k-1)(k - 1) =
and conversion is required
Rm (k - 1)'
. 0(k-i-m)m .
Qk (mr (k)YMk (k)
Using the result of (9), we define the matrix
\qh- i(k -1)
Rm (k)
0(k-M )M_
Q H (k)
oLi
Ok-i 1
(9)
(10)
(11)
If the matrix Ak'5(k)YH (k) multiplied from left on the matrix (11), this operation modifies equation (9), adding it to the matrix in the right side of (k + 1) — w (bottom) row:
Q H (k)Ak'®i(k)Y H, (k) = Q h (k)
\0'5Al'-i(k - 1)YM(k-i)(k - 1)
QH-iA°'5Al-i(k - 1)YM(k-i)(k - 1)
yM (k)
YM (k)
rA0' 5flm(k - 1)'
0(k-i-M )M
H M
YM (k)
(12)
To perform the transformation (10), in equation (12) it is necessary to zero the last line. From equations (10) and (12) it follows that
QH Ak 5(k)YHk(k) = Q H (k)Q H (k)Ak'5 (k)YHk (k)
r\0 ' 5
= Q H (k) =
'A0' 5Rm(k - 1)'
0(k-i-M )M
M
YM (k)
R M (k)
= 0(k-i-M )M
- 0M -
■ R M (k)
0(k-M )M_
(13)
where the matrix Q^ (k) is the product of two matrices
QH (k) = Q H (k)Q H (k).
(14)
Thus, the reduction of matrix Ak'5(k)YMk (k) to a triangular form using the matrix Qf (k) (10) at iteration k can be done by zeroing the last row in the matrix (12) using the result of reducing the matrix Ak-^k — 1)Yf (k-1)(k — 1) to a triangular form obtained at iteration k — 1. This zeroing is carried out using matrix Qf (k), which is a product of matrices composed of Givens rotation matrices.
The recursive relationship between RM(k — 1) and RM(k)in a more compact form, i.e., when k matrices with a fixed number of elements (M +1) x (M + 1) are used at each iteration, can be represented using the following equation
G
M+1
(k)
A°'5R(k - 1)
. YH (k) .
R (k)
,oM.
(15)
where matrix GM+1(k) is unitary. This matrix can be formed using Givens spins. The structure of matrices GM+1 (k) is a "compressed" to size (M + 1) x (M + 1) matrix QH(k) with a variable number of elements (k) x (k).
The elements of Givens matrices are determined from equation (15), where
m
G
M+1
(k) = n GM +1(k).
(16)
¿=1
Givens GM (k) transformation is determined by the plane rotation matrices of the form:
G(i,j) =
0 1
1
(17)
Matrix Gi,j with fixed values of i,j € {1, 2,...,m — 1} differs from the identity n-matrix
E in that in it the 2 x 2-submatrix E occupying the cell formed by the intersection of the i-th
fc —s*\
and j-th rows and columns is replaced by the submatrix Gj =--, with elements c and
\s c J
s satisfying the condition
2 2 (18)
s2 + c2 = 1.
With this normalization condition, matrix Gj and matrix Gj are orthogonal. The elements c and s can be interpreted as the cosine and sine of a certain angle of rotation transformation.
Using a sequence of such orthogonal transformations, matrices G1, G2,..., Gm-1 of the form (17) can be reduced to the right triangular form by sequentially canceling the subdiagonal elements in the first, second, ..., (n — 1)-st columns. We consider the first of N steps leading to
1
0
0
c
s
1
1
s
c
1
0
transformation (15), for which matrix G^+^k) is used. Then
G
(1) M+1
(k) x
A0-5R n (k — 1)"
Y M)H (k)
CM,1(k) 0 0 1
00
SM,1(k) 0
0 -S
M,1
(k)
IN-3
10
0 CM,1(k)
'A0-5 RM,u(k — 1) A0-5 Am,12 (k — 1) 0 A0-5 Rm,22 (k — 1)
vM)! (k)
yA0,2 (k)
A0-5 Rm,1(m-1)(k — 1) A0-5 Rm,1m (k — 1) A0-5 Rm,2(m-1)(k — 1) A0-5 Rm,2M (k — 1)
Rm,(m-1)(M-1)(k — 1) Rm,(m-1)M (k — 1)
0 Rm,mm (k — 1)
yM,(M-1)(k) yM,M(k)
A0-5 Rm,11 (k) A0-5 i?M,12(k — 1) 0 A0-5 Rm,22 (k — 1)
y\M,2 (k)
a0-5 Rm,1(m-1)(k — 1) a0-5 Rm,1m (k — 1) A0-5 Rm,2(m-1)(k — 1) A0-5 Rm,2M (k — 1)
Rm,(m-1)(M-1)(k — 1) Rm,(m-1)M (k — 1)
0 Rm,mm (k — 1)
V((i,(M-1) (k) yM,M (k)
In equation (19), vector YN°)H (k) is defined as
(19)
Y (0)h
YN
(k)
VM)((k),yM)2((k),
•, yii) m (k)] = Yf (k)-
(20)
The superscript in parentheses means the number of the transformation performed on variable ytM*((k) with number i in the vector. Such a conversion over y(k) is performed once (after the first time the variable is reset), over y^^k) — twice (after the second time, the variable is
reset), etc. and over y(°)M(k) — M times (after the M-th time, the variable is reset). From equation (19) it follows that
CM,1(k)A0-5 R.M-u(k — 1) — s M,1(k)y(0)2(k)(k) = RM.u(k)
and
sM,1(k)^5 RM-u(k — 1) + CM,1(k)yM)2(k) = 0. Performing similar conversion of all i =1,2,..., M, can be established that
CM,i(k)A0-5 R,M-jj(k — 1) — sMij(k)y(Mj1)2(k) = RM.ii(k),
SM,i(k)A0-5 RM.jj(k — 1) + CM}i(k)y{M-p2(k) = 0. Then from equation (22) we can determine that
(21) (22)
(23)
sM,j(k) = — cm j(k)y(Mi)2(k)[A0-5 RM,u(k — 1)]
-l
0
0
x
x
0
0
0
and, given that c + ss* = 1, from the equation
c2M,i(k) + cM,i(k)y<i^,i1)*(k)[^3'5RM.ii(k — 1)] ^M,i(k)yM (k)[A0'5 RM.ii(k — 1)]
(i-i)*f
cMAk)[x0'5 RM.u(k — 1)]2 + cMAmM^mM-Pk]
M,i
= cM,i(k)[\RM,ii(k — 1) + M^'M^RMAk — 1)] * = 1
i
(25)
can determine that
cm, i(k) = \jARM.ii(k — 1) {^R2M.ii(k — 1) + yl^^viPk
= A0'5RM.ii(k — 1)RMii(k),
(26)
where
R-M.ii(k) =\[ARRMM~ii(k-l)+УM:J^k)УM^).
(27)
It is taken into account that the diagonal elements R-Nii(k) of the matrix RM(k) are real numbers. Then, using (26) in equation (24), we can determine
sN,i(k) = —A0'5RM.ii(k — 1)RM1u(k)+ vM-rk A0'5RM,u(k — 1)
i
—yM:iiy*(k)R-MLii(k)- (28)
i
Thus, equations (26)-(28) allow us to calculate the cosine and sine of a certain rotation angle. According to (19), they calculate the elements of the i-th row of matrix RM(k) , zero out the ¿-th element of vector YM 1)H (k), and modify the remaining nonzero elements of this vector as
0,0,..., 0,0,yMi+i(k),...,yMM (k)
(i)*
, i.e.
r(i- 1)H(i,\ _
M
(k)= 0, 00, yMf'tk), M+i (k),..., v MM (k)
(i-1)*/
~(i-1)*.
(29)
and
Y (i)H/
YM
(k)= 0, 0,..., 0, 0,yW*i+i(k),...,yM*M (k)
(i)*
(30)
These transformations for each value of i include calculations (26)-(28), and for all j = i + 1,..., M, calculations similar to (23) and (24), i.e.
RM'ij (k) = cM,i(k)A0'5RM'ij (k — 1) — sMi^yMJ*1* (k)
(31)
and
vM*j (k) = SM,i (k)A0'5RM.ij (k — 1) + cM,i(k)y(M1,j1)*(k).
(32)
Thus, transformation (13) can be performed either as QH(k)Ak'5(k)YHk)(k), applying a
k x k matrix QH (k) to a k x M matrix Ak (k)YHk)(k) at each iteration, or using a matrix
M
QH (k) = ]1 GMr(k) applied to a k x M matrix
i=i
A0'5R m (k — 1)'
(k-i-M)M MH (
YM (k)
i
2. Parallel implementation of the QR algorithm in a triangular systolic array
Givens transformation has good properties for use in a triangular systolic array. The architecture of the basic computations of the algorithm using such calculators is given in [16] and in Fig. 1.
Fig. 1. Block diagram of a triangular systolic array
The systolic array is based on the method of triangular complex rotations and allows to obtain a significant performance gain in comparison with the method of complex rotations of Givens.
In the system of a triangular systolic array, there are individual processing cells located in an ordered structure. Each individual cell of the system has its own processing functionality and local memory. Only neighboring cells are connected to each other and there is no direct connection between cells that are not adjacent. When data is fed into the systolic array system, the processing cells at the front end of the system will process the data, store them in local memory, and then forward them to adjacent cells. This processing and transfer of the processed data in each cell continues until the data stream reaches the end of the system, where the final calculation results are obtained. The proposed architecture provides a significant reduction in the time required to perform QR decomposition using the same computing resources (CORDIC computational cells). Another advantage of the proposed scheme is that during QR decomposition, the upper triangular matrix R has only real diagonal elements. This simplifies the subsequent inversion of the matrix R using the backward substitution algorithm, which requires division by the diagonal elements of the matrix R .
Algorithm (15) shows that the application of Givens rotations in a post-array to multiply matrices on the left side of the preliminary array allows one to obtain a triangular matrix and zero out the input vector Y(k). The number of elements in the input vector Y(k) corresponds to the number of antennas. The number of required Givens rotation operations is the same as the number of elements in the input. That is, each rotation of Givens will reset exactly one element of the input vector. Thus, the algorithm for generating a radiation pattern with K antennas requires the Givens rotation operation K in the calculation of the post-array.
In calculations, the Givens rotation operation can be performed in parallel, because there is no data dependence between the Givens rotation operation at one of the inputs and the Givens rotation operation at the same position in subsequent iterations. Thus, Givens rotation operations can be performed in parallel.
Conclusion
The proposed architecture of the triangular systolic array using the method of the triangular complex rotations optimized for implementation in large-scale integrated circuits, allowing you to effectively perform the operation QR-decomposition of complex matrices. Compared with the QR-RLS algorithm, the proposed architecture can provide a gain of up to 35% in the time of calculating the QR decomposition. The synthesized algorithm will make it possible to implement the methods of spatio-temporal processing of broadband signals of satellite communication systems.
This work was supported by the Ministry of Science and Higher Education of the Russian Federation in the framework of the Federal target program «Research and development on priority directions of development of the scientific-technological complex of Russia for 2014-2020» (agreement no. 05.605.21.0185, unique ID project RFMEFI60519X0185).
References
[1] Ya.D.Shirman, Radio-electronic systems: Fundamentals of construction and theory. Directory. Ed. 2nd, rev. and add., Moscow, Radio Engineering, 2007 (in Russian).
[2] V.N.Tyapkin, I.N.Kartsan, D.D.Dmitriev, S.V.Efremova, Algorithms for adaptive processing of signals in a flat phased antenna array, SIBCON, 2017.
[3] A.I.Perov, Statistical theory of radio systems, Moscow, Radio Engineering, 2003 (in Russian).
[4] I.N.Kartsan, V.N.Tyapkin, D.D.Dmitriev, A.E.Goncharov, P.V.Zelenkov, I.V.Kovalev, IOP Conf. Series: Materials Science and Engineering, 155(2016) no. 1.
DOI: 10.1088/1757-899X/155/1/012019
[5] I.N.Kartsan, V.N.Tyapkin, D.D.Dmitriev, A.E.Goncharov, I.V.Kovalev, IOP Conf. Series: Materials Science and Engineering, 255(2017). DOI: 10.1088/1757-899X/255/1/012009
[6] V.N.Tyapkin, D.D.Dmitriev, Yu.L.Fateev, N.S.Kremez, J. Sib. Fed. Univ. Math. & Phys, 9(2016), no. 2, 258-268. DOI: 10.17516/1997-1397-2016-9-2-258-267
[7] I.N.Kartsan, Y.L.Fateev, V.N.Tyapkin, D.D.Dmitriev, A.E.Goncharov, P.V.Zelenkov, I.V.Kovalev, IOP Conf. Series: Materials Science and Engineering, 155(2016), no. 1. DOI: 10.1088/1757-899X/155/1/012020
[8] R.A.Monzingo, T.W.Miller, Introduction to Adaptive Arrays, New York, John Wiley & Sons, 1980.
[9] Ya.D.Shirman et al., The first domestic studies of the adaptation of antenna systems to interfering influences, Moscow, Radio engineering, no. 11, 1989 (in Russian).
[10] V.N.Tyapkin, D.D.Dmitriev, V.G.Konnov, A.N.Fomin, A method for determining the vector of spectral coefficients by the likelihood ratio criterion, Bulletin of the Siberian state Aerospace University named after Acad. M. F. Reshetneva, 43(2012), no. 3, 76-79.
[11] V.I.Dzhigan, Equivalence conditions for recursive adaptive filtering algorithms by the least squares criterion, Telecommunications, 6(2006), 6-11 (in Russian).
[12] I.A.Lubkin, V.N.Tyapkin, The use of recurrent adaptive algorithms to solve the problem of suppressing active noise interference in satellite communication systems, Bull. Sib. state Aerospace Univ. named after Acad. M.F. Reshetneva, 28(2010), no. 2, 39-43 (in Russian).
[13] G.H.Golub, C.F.Van, Lone Matrix calculations, The Johns Hopkins University Press, Bal-timor, 1996.
[14] D.T.M.Slock, T.Kailath, Numerically stable fast transversal filters for recursive least squares adaptive filtering, IEEE Trans. Signal Processing, 39(1991), no. 1, 92-114.
[15] A.Benallal, A.Gilliore, A new method to stabilize fast RLS algorithm based on the firstorder model of the propagation of numerical errors, Proc. Int. Conf. on Acoustic, Speech and Signal Processing, Vol. 5, 1988, 1373—1376.
[16] J.G.McWhirter, Recursive least-squares minimization using a systolic array, Proceedings of the SPIE Intern. Sic. Opt. Eng. Vol. 43, 1983, 105-112.
Рекурсивный алгоритм оценивания корреляционной матрицы помех, основанный на QR-разложении
Валерий Н. Тяпкин Дмитрий Д. Дмитриев Андрей Б. Гладышев Пётр Ю. Зверев
Сибирский федеральный университет Красноярск, Российская Федерация
Аннотация. Многие задачи цифровой обработки сигналов требуют выполнения матричных операций в режиме реального времени. Это операции обращения матрицы или решения систем линейных алгебраических или дифференциальных уравнений (фильтр Калмана). Переход к реализации цифровой обработки сигналов на программируемых логических интегральных схемах (ПЛИС), как правило, предполагает расчеты, основанные на представлении чисел с фиксированной точкой. Это делает практически невозможным решение задач пространственно-временной обработки на основе традиционных вычислительных методов. В статье рассматривается реализация алгоритмов пространственно-временной обработки сигналов в широкополосных спутниковых системах с использованием QR-разложения. Представлены технологии вычислений CORDIC, необходимые для повторного QR-разложения при совместном использовании в систолических алгоритмах.
Ключевые слова: фазированная антенная решетка, адаптивные алгоритмы, фильтр Калмана, рекурсивный алгоритм по критерию наименьших квадратов, QR-разложение, систолические алгоритмы.
