The synthesis algorithm for spatial filtering to maintain a constant level of the useful signal Текст научной статьи по специальности «Физика»

УДК 519.873

The Synthesis of Spatial Filtering Algorithm with a Constant Level of Useful Signal

Valeriy N. Tyapkin* Dmitry D. Dmitriev^ Yuri L. Fateev* Nikolay S. Kremez§

Military Engineering Institute Siberian Federal University

Svobodny, 79, Krasnoyarsk, 660041

Received 10.02.2016, received in revised form 26.02.2016, accepted 22.03.2016 The theoretical basis of two classical recursive algorithms is considered in the paper. They are beamform-ing algorithm based on the Kalman filter and the least squares algorithm based on QR decomposition. An algorithm that minimizes the noise at the output of an antenna array is proposed on the basis of these two algorithms. The algorithm allows one to maintain a constant level of the received useful signal.

Keywords: phased array, adaptive algorithm, Kalman filter, recursive least squares method, QR decomposition.

DOI: 10.17516/1997-1397-2016-9-2-258-267.

Introduction

Modern navigation receivers of global satellite systems have a significant drawback - low interference immunity [1]. To solve this problem, navigation receivers are based on phased array antennas that use adaptive beamforming techniques.

Pattern control is widely used in applications such as navigation, radar, wireless communication and others. The beamforming technique is one of the most popular methods to form digital antenna directional diagram. This technique allows one to form maximum in the antenna diagram in the direction of useful signal and to form null in the antenna diagram in the direction of interference. Then interference is suppressed and useful signal is accumulated from appropriate directions. The overall effect of this treatment is determined by the difference in useful signal and interference and statistical characteristics of useful signal and interference. Since beamforming takes place in digital form, it is necessary to use digital techniques that could be performed in real time. Therefore, the solution of this task depends on both hardware and software.

We consider theoretical foundations of digital beamforming algorithms such as the recursive Kalman filtering algorithm and the recursive least squares algorithm based on the QR decomposition (QR-RLS)[2].

On the basis of these two algorithms we propose an algorithm that minimizes the noise at the output of an antenna array and allows one to maintain a constant level of the desired signal.

* [email protected]

[email protected]

[email protected]

§ [email protected] © Siberian Federal University. All rights reserved

1. Kalman filtering algorithm

The Kalman filtering algorithm is an important adaptive digital signal processing algorithm. It is an efficient recursive filter that estimates the state vector of a dynamic system and uses a series of incomplete and noisy measurements. The Kalman filter for recursive estimation of the state vector of a dynamic system is known a priori. In order to calculate the current state of a system it is necessary to know the current measurement and the previous state of the filter. Thus, the Kalman filter, like other recursive filters, operates in time rather than in the frequency domain. The Kalman filter produces not only estimates of the current state but also uncertainty (distribution density) of the state vector that is based on Bayes' formula of conditional probability.

Traditional Kalman filter has the following parameters:

K(n) = K(n, n - 1) - F(n, n +1) • G(n) • C(n) • K(n, n - 1),

G(n) = F(n + 1,n) • K(n,n - 1) • CH • [C(n) • K(n,n - 1) • CH(n) + Q2(n)]_\ (1) a(n) = y(n) - C(n) • x(n|yn_i),

where C(n) is the observation matrix, F(n, n + 1) is the state transition matrix.

In the case of not induced noise dynamic model the process converges to white noise with zero mean and unit variance. Therefore, some parameters of the state space converge to a constant and they are

F(n +1,n) = VAI, C(n) = uH(n), Q2(n) = 1,

H

K(n - 1) g(n) • u(n) • K(n - 1)

K(n) = —---TA-

g(n) = K(n - 1) • u(n)

(2)

y/x (l + u(n)H • K(n - 1) • u(n)) ' a(n) = y(n) = uH(n) • x(u\yn-i).

Using expressions (2), we obtain parameters of the Kalman filter: ( ) K-1(n) • u(n)

g(n) = —ta—'

K-1(n) = A • K-1(n - 1) + u(n) • u(n)H, (3)

K-I(n) • x(n + 1\yn) = K-1(n - ^x(n|y"-l) + u(n) • y(n)

VA vA '

Then, by rewriting the expression K_1(n) = K_H/2(n) • K_1/2 and using the lemma on matrix inversion, we define the Kalman filter in the form of pre-array and post-array as follows

' VA • K-H (n - 1) xH(n\yn-1) • K-H (n - 1) 0T

a/A • u(n)

y*(n) 1

)=

K-H (n)

xH (n + 1\y„) • K-T (n)

0

a'(n)

Vr(n)

VA • uH (n) • K 2 (n) \J r(n)J

(4)

With the help of relation (4) coefficients of the Kalman filter can be recursively computed.

2. Recursive least squares algorithm based on the QR-decomposition (QR-RLS)

A recursive algorithm for adaptive filtering that uses QR-decomposition estimates the filter coefficients in the current step in terms of the filter coefficients calculated in the previous step. Due to its recursive nature the algorithm is called the recursive least squares algorithm (RLS). The basic idea of the QR-factorization is to convert a matrix into a product of an orthogonal matrix Q and an upper triangular matrix R.

Let us assume that m and n are any positive integers greater than 0. The QR-decomposition of matrix A of size m x n for every m > n has the form

A=QR. (5)

Let Q be a unitary matrix, i.e., Q • QH = I, where I is the identity matrix and R is m x n upper triangular matrix. Then equation (5) can be written in the divided form

A = QR = Q

Ri 0

[Qi Q2]

Ri 0

Q1R1,

(6)

where Ri is an n x n triangular matrix, Q1 is an m x n matrix and Q2 is an m x (m — n) matrix.

Matrix R can be obtained with the use of Givens rotations [3]. A Givens rotation is represented by a transformation matrix Gj of the form

G(i,j ) =

1

0

(7)

1

where c and s appear at the intersections of ith and jth rows and columns. These elements satisfy the following condition

2 , 2 s + c

1.

(8)

Matrix Gj is an orthogonal matrix and it represents the plane rotation because elements c and s can be represented as c = cos(^) and s = sin(^). Properly selecting elements c and s and multiplying matrix A by transformation matrices Gi, G2,..., Gn-1, one can reduce matrix A to an upper triangular matrix. Now we proceed to the QR-RLS algorithm.

Assume that the amplitude of navigation signal is significantly lower than the intensity of the interference. Then, the signal does not have significant influence on the estimate of the correlation matrix of the interference .

0

0

1

c

s

1

1

s

c

1

0

Suppose that at the inputs of M-element phased array there is an interference with the

instantaneous oscillations u1 (t),u2(t),.. .uM(t). Oscillation phases are shifted relative to each

2n

other due to the path difference AD from element to element by an amount Ap = — • AD. According to the Kotelnikov theorem, interfering signals at the output of the i-th element of a phased array can be represented as a set of discrete samples of instantaneous amplitude uii with sampling period TD = 1/2fm (fm is the maximum frequency in the spectrum of the received interference signal) (Fig. 1). Here, i = 1,2,..., M; l = 1,2,..., k are the numbers of discrete samples.

Fig. 1. Sampling signal interference

The cross-correlation between signals of the channels i and j is defined as the arithmetic

1 k

mean of the correlation points uiuj over the number of counts l from 1 to k — pij = T^^uuuji.

i=i

Upon introducing the complex amplitudes of instantaneous discrete values of Uj; and j, integrated mutual correlation time is defined as

1

V = tJ2uuu*i/2.

(9)

i=i

The weight vector at stage n is defined as

w(n) = [wo(n), wi(n), ..., wm—i(n)] . The objective function for the RLS algorithm is

(10)

(11)

where e(i) = d(i) — wH(n) ■ u(i). Let wH(n) be the optimum value such that the target function S(n) is minimal. This value satisfies the normal equation of the form [4]

$(n) ■ w(n) = z(n)

where $(n) = £ Xn—1 • u(i) • u(i)^, z(n) = £ Xn—1 • u(i) • d *(i). i=i i=i The last two equations can be written in a recursive form as follows

2

$(n) =A • $(n - 1) + u(n) • u(n)H, (13)

z(n) = A • z(n - 1) + u(n) • d*(n). ( )

Next we use the lemma on matrix inversion. The inverse of A + aBB,T has the form [5]

(A + — )_1 = A_1 - ■ (14)

Using (14), we obtain

P(n) = P(n - 1) - k(n) • u(n)H • P(n - 1), (15) AA

_, . _1, , , , , A_1 • P(n - 1) • u(n) N . . , ,

where P(n) = $ 1(n), k(n) =-—-h-= P(n) • u(n) = $ 1(n) • u(n).

1 + A • u(n) • P(n - 1) • u(n) Using relations given above, we can obtain expression for the weight vector in the constraint equation (4):

w(n) = $ 1(n) • z(n), w(n) = P(n) • z(n), w(n) = A • P(n) • z(n - 1) + P(n) • u(n) • d* (n),

w(n) = P(n - 1) • z(n - 1) - k(n) • u(n)H • P(n - 1) • z(n - 1) + P(n) • u(n) • d*(n), (16) w(n) = w(n - 1) - k(n) • u(n)H • w(n - 1) + P(n) • u(n) • d*(n), w(n) = w(n - 1) + k(n) • £*(n),

where £*(n) = d(n) - wH(n - 1) • u(n). From equations (15) and (3) one can observe correspondence between the RLS algorithm and the Kalman filter. Let us rewrite $(n) as

$(n) = $_1 (n) • $_T? (n), (17)

and introduce a new variable

p(n) = $(n) • w(n) = $_1 (n) • z(n). (18)

Then the RLS algorithm can be represented in the form of pre-array and post-array similar to the previously described Kalman algorithm:

a/A • $ ? (n - 1) u(n) VA • pH (n - 1) d(n) 0T 1

)=

$ 2 (n) pH (n)

0

a/y(")

uh (n) • $ 2 (n) vYm.

, wH(n) = pH(n) • $_2 (n). (19)

The QR-RLS algorithm is used to solve the problem of adaptive beamforming due to its recursive nature, computational efficiency and numerical stability [6].

3. An algorithm for minimizing the noise at the output of an antenna array and maintaining a constant level of the useful signal

One of the most promising methods for adaptation of the antenna arrays is the method of minimum variance distortionless response (MVDR). The MVDR algorithm allows one to increase

the signal to noise ratio in the desired direction while suppressing the interference of received signals.

The problem of minimizing the noise at the antenna array and maintaining a constant level of useful signal can be formulated as follows:

mm w(n)

J2Xn-i\e(i)\2,

(20)

where e(i) = wH (n) ■ u(i), with the constrain

wH (n) ■ s(6o) = 1,

(21)

where s(60) is the vector that determines the direction of the source signal. Taking into account constraints (21), a weight vector for the antenna array can be written as follows [7]:

w( n) =

$-1(n) ■ a(0o)

H(0o) ■ $-1(n) ■ s(9o)'

(22)

Let us define additional vector a(n) = $ 2 (n) * s(90). Then we rewrite expressions for the weight vector and error estimates as follows

w(n) = e(n) =

$ 2 (n) ■ a(n) \\a(n)\\2i ' aH(n) ■ $ 2 (n) ■ u(n)

(23)

e'(n) = aH (n) ■ $ 2 (n) ■ u(n).

When parameters of the MVDR algorithm are set, it is possible to solve the problem of maintaining a constant level of useful signal in the QR-RLS algorithm and present it in the form of pre-array and post-array as follows

y/X ■ $ 2 (n — 1) u(n) ■ aH (n — 1) d(n) 0T 1

9(n)

$ 2 (n) aH (n)

0

-e'(n)

a/yC")

uH (n) ■ $ 2 (n) ^Jy (n)_

(24)

After calculating parameters with the use of (24), one can estimate the error:

e(n)

— (—e'(n) VY(n))

\\a(n)\\2 ■

(25)

4. Representation of the MVDR algorithm as a systolic array

Expressions obtained above show that the calculation of the response of a phased array antenna with the use of MVDR algorithm or its predecessor the QR-RLS algorithm involves matrix multiplication. This requires a series of Givens rotations. It is known that the consistent implementation of the matrix multiplication is inefficient. Real time signal processing systems that use the QR-RLS and MVDR beamforming algorithms require a more efficient matrix multiplication algorithm.

s

In 1978, Kung and Leiserson have proposed systolic arrays for matrix calculations in very large-scale integration systems. The systolic array is based on the method of triangular complex rotations and it significantly increases productivity of calculations in comparison with the traditional method of complex Givens rotations.

A systolic array system contains individual processing cells arranged in a specific structure. Each individual cell has its own function and local memory. Furthermore, only adjacent cells are connected to each other.

Data are written in a systolic array of processing cells on the entrance. Then data are processed, stored in local memory and transmitted to adjacent cells. This processing and sending the processed data in each cell lasts as long as the data stream reaches the end of the system where final results of calculation are obtained. The propagation of data through a systolic array resembles the pulse of the human circulatory system. The name is derived from the medical term systole as an analogy to the regular pumping of blood by the heart. This approach allows one to significantly reduce the time required to perform QR-decomposition using the same computational resources. The second advantage of the proposed scheme is that the QR-decomposition is executed such that the upper triangular matrix R has only real diagonal elements. This greatly simplifies the subsequent inversion of the matrix R with the use of back-substitution algorithm.

Consider the MVDR algorithm. After applying Givens rotations to (24), we obtain lower triangular matrix in the right hand side of (24) and set input vector u(n) to zero. The number of elements in the input vector u(n) corresponds to the number of antennas.

To set the input vector to zero, it is necessary to apply a number of Givens rotations. The number of rotations is the number of elements in the input vector because each Givens rotation sets one element of the input vector to zero. In general, the MVDR beamforming algorithm with K antennas requires K Givens rotations.

Calculations of Givens rotations can be performed in parallel because there is no data dependency between elements.

Let us consider expressions (7) and (8). One can see that each Givens rotation changes only certain elements of the pre-array, other elements are remained unchanged. To illustrate this we present the MVDR algorithm in the case of the antenna array consisting of three elements.

where

$'1,1 = ($1,1 • cos) + (u1 • sin),

u'1 = ($1,1 • (- sin')) + (u1 • cos) = 0,

$'2,1 = ($2,1 • cos) + (u2 • sin),

u'2 = ($2,1 • (-si n')) + (u2 • cos), $'2,2 = $2,2,

$'3,1 = ($3,1 • cos) + (u3 • sin),

u'3 = ($3,1 • (-sin*)) + (u3 • cos),

$'3,2 = $3,2, $'3,3 = $3,3,

a' 1 = «1 • cos, a'2 = «2,

a'3 = «3, ¡' = «1 • (-sin*).

Step 1:

$1,1 0 0 u1

$2,1 $2,2 0 u2

$3,1 $3,2 $3,3 u3

a1 a2 a3 0

0 0 0 1

cos 0 0 - sin

0 10 0

0 0 1 0

sin 0 0 cos

(26)

Step 2:

$ 1,1 0 00

$ 2,1 $ 2,2 0 u '2

$' 3,1 $ 3,2 $'3,3 u'3

a 1 a 2 a' 3 /3'

. 3 ' 3' 3 3

where

$ ' '1,1 = $ ' 1,1, u 1 = 0,

$ ''2,1 = $'2,1, $ ''2,2 = ($ '2,2 • cos) + ( u'2 • sin),

u''2 = ($ '2,2 • (-sin')) + (u '2 cos) = 0,

$ ''3,1 = $ '3,1, $ ''3,2 = ($ '3,2 • cos) + ( u'3 • sin),

u''3 = ($ '3,2 • (-sin*)) + (u ' 3 cos) ,

$ ''3,3 = $ '3,3, a!' 1 = a! 1, a 2 = (a ' 2 cos) + (/'

a ''3 = a '3, /'' = ((a'2 • ( - sin* )) + (/ ' • cos)) .

Step 3:

r$ ' '1,1 0 0 0

$ 2,1 $ 2,2 0 0

$''3,1 $''3,2 $ 3,3 u 3

a 1 a 2 a 3 ¡''¡'' /3''

where

$ '''1,1 = $ ''1,1 , u ' ' '1 = 0,

$ '''2,1 = $ ''2,1 , u ' ' '2 = 0,

$ '''2,2 = $ ''2,2 , $'''3,1 = = $''3,1, $'''3,2 = $''3,2,

0

cos 0

sin

3'' 3''

$ ' ''3,3 = ($''3,3 • cos) + (u''3 • sin), u' ''3 = ($''3,3 • (-sin*)) + (u ''3 • cos) a'''1 = a''1, a!''2 = a''2, a'''3 = (a ''3 • cos) + (/ ' ' • sin), 3' '' = ((a''3 • (—sin*)) + (3' ' • cos)).

$ 1,1 0 0 0

$ 2,1 $'' '2,2 0 0

$ 3,1 $ 3,2 $ 3,3 0

a 1 a 2 a 3

/3''' 3''' 3''' 3'''

0

- sin 0

cos

(27)

"1 0 0 0

0 1 0 0

0 0 cos — sin '

0 0 sin cos

(28)

(29)

Fig. 2 shows the block diagram of the triangular systolic array that may be used to calculate the QR-decomposition of 3-by-3 matrix.

During systolic cycle processors execute one operation and transfer data to adjacent nodes of the triangular systolic array. Fig. 2 also shows the input and output data streams. The input vector is loaded after every systolic cycle. The square in Fig. 2 denotes the inside cell. It is responsible for the application of the elementary Givens rotations to the relevant matrix elements in accordance with expression (7).

— (—e'(n) • y— 2 (n) • y2 (n))

The error e(n) is calculated by the following formula: e(n) = ---.

IHn)H2

The circle with the index "x" corresponds to this operation on the diagram shown in Fig. 2.

Other circles shown on the diagram denote border cells. These cells are responsible for the calculation of diagonal elements of the transformed matrix.

Fig. 2. Block diagram of the triangular systolic array

Conclusion

The proposed structure of the triangular systolic array is optimized for implementation in large scale integrated circuits. It is based on triangular complex rotations and enables the efficient implementation of the QR-decomposition of complex matrices. The proposed structure can decrease the time of computation of QR-decomposition by 35% in comparison with the QR-RLS algorithm.

This work was supported by the Ministry of Education and Science of the Russian Federation (agreement No. 14.575.21.0081, project RFMEFI57514X0081).

References

[1] V.N.Tyapkin, D.D.Dmitriev, T.G.Moshkina, The potential noise immunity of consumer navigation equipment satellite navigation systems, Vestnik SibGAU, 43(2012), no. 3, 113-119 (in Russian).

[2] C.Dick, F.Harris, M.Pajic, D.Vuletic, Implementing a Real-Time Beamformer on an FPGA Platform: We designed a flexible QRD-based beamforming engine using Xilinx System Generator, XCell Journal, (2007), 36-40.

[3] D.Bindel, J.Demmel, W.Kahan, On Computing Givens Rotations Reliably and Efficiently., Acm Trans. Mathematical Software, (2002), 206-238.

[4] Ed. A.I.Perov, V.N.Kharisov., GLONASS. Principles of construction and operation. 4th ed., Revised, Moscow, Radio, 2010 (in Russian).

[5] Ya.D.Shearman, Electronic systems: Basics of construction and theory. Resource Book. Ed. 2nd, Revised. and add., Moscow, Radio, 2007 (in Russian).

[6] V.N.Tyapkin, D.D.Dmitriev, E.N.Garin, A.V.Sokolovsky, Management of amplitude-phase distributed adaptive phased array, Industrial Automatic Control Systems and Controllers, 4(2013), 66-71 (in Russian).

[7] V.N.Tyapkin, D.D.Dmitriev, V.G Konnov, A.N.Fomin, The method of determining the vector of spectral coefficients of the likelihood ratio, Vestnik SibGAU, 43(2012), no.3, 7679 (in Russian).

Синтез алгоритма пространственной фильтрации с сохранением постоянного уровня полезного сигнала

Валерий H. Тяпкин Дмитрий Д. Дмитриев Юрий Л. Фатеев Николай С. Кремез

Военно-инженерный институт Сибирский федеральный университет Академгородок, 13А, корпус S, 660036

Россия

В статье рассматриваются теоретические основы формирования диаграммы направленности на на основе двух классических рекурсивных алгоритмов — алгоритм фильтра Калмана и алгоритм наименьших квадратов на основе QR-разложения. На основе этих двух алгоритмов синтезирован алгоритм минимизации шума на выходе антенной решетки, что позволяет поддерживать постоянный уровень принимаемого полезного сигнала.

Ключевые слова: фазированная антенная решетка, адаптивные алгоритмы, фильтр Калмана, рекурсивный метод наименьших квадратов, QR-разложение.