On optimizations with magnitude constraints on frequency or angular responses

Accepted Manuscript

On Optimizations with Magnitude Constraints on Frequency or Angular Responses Junli Liang, H.C. So, Jian Li, Alfonso Farina, Deyun Zhou PII: DOI: Reference:

S0165-1684(17)30422-X 10.1016/j.sigpro.2017.12.009 SIGPRO 6678

To appear in:

Signal Processing

Received date: Revised date: Accepted date:

13 August 2017 2 December 2017 6 December 2017

Please cite this article as: Junli Liang, H.C. So, Jian Li, Alfonso Farina, Deyun Zhou, On Optimizations with Magnitude Constraints on Frequency or Angular Responses, Signal Processing (2017), doi: 10.1016/j.sigpro.2017.12.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

1

On Optimizations with Magnitude Constraints on Frequency or Angular Responses Junli Liang*1 , H. C. So2 , Jian Li3 , Alfonso Farina4 , and Deyun Zhou1 School of Electronics and Information, Northwestern Polytechnical University, Xi’an, China (e-mail: heery [email protected]). 2

Department of Electronic Engineering, City University of Hong Kong, Hong Kong, China (e-mail: [email protected]).

Department of Electrical and Computer Engineering, University of Florida, USA (e-mail: [email protected]). Leonardo Company Consultant, Rome (Italy)(e-mail: [email protected]).

M

4

AN US

3

CR IP T

1

ED

Abstract

We consider optimization problems with nonconvex magnitude constraints on frequency or angular responses encountered in many signal processing problems. To provide a solution to these nonconvex problems, an alternating direction method of

PT

multipliers based solution framework is developed in this paper. We apply variable splitting to introduce auxiliary two-element vector variables for each frequency or angular response corresponding to each magnitude constraint. In doing so, this class of

CE

problems can be solved in an alternate manner, namely, iteratively tackling two subproblems with respect to the original and introduced parameters. Each new constraint may act on only the corresponding two-element vector variables rather than all of

AC

them. The subproblem is also simplified such that it is a function of only one parameter, and we prove the convexity of the resultant single-variable piecewise optimization problem. Finally, the effectiveness of the proposed approach is demonstrated via its successful applications to array pattern synthesis, waveform design and robust beamforming.

Index Terms Array signal processing, Waveform design, Array pattern synthesis, Beamforming, Filter design. *Corresponding author.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

2

1. I NTRODUCTION In many signal processing problems, including adaptive beamforming [1], [2], [3], [4], [5], array pattern synthesis [6], [7], waveform design [8], [9], [10], [11], and finite impulse response (FIR) filter design [13], nonconvex magnitude constraints [14], [15] on frequency or angular responses are often imposed for performance improvement: 1) A beamformer [1], [2] is used for spatial filtering and interference rejection to detect and estimate the signal of interest at

CR IP T

the output of a sensor array. Since it plays an important role in radar, sonar, audio and speech processing, medical imaging, and communications, considerable attention has been received. However, in practical applications the provided steering vector of the desired signal may be imprecise due to look direction errors, imperfect array calibration, and distorted antenna shape, etc. [1]. As a result, the performance of the Capon or minimum power distortionless response (MPDR)

AN US

beamformer [3], which assumes perfect a priori knowledge of the steering vector, may degrade severely. The main reason is that the desired signal is suppressed due to the mismatch between the presumed and actual array responses [1]. To be robust against the imprecise steering vector, nonconvex magnitude constraints on the potential responses of target incoming angles can be imposed to improve the beamformer performance, including exceeding-unit [17] or two-sided magnitude constraints [18], [19] on the angular region of interest.

M

2) Compared to a single antenna, an array of multiple elements can produce different beampatterns by modifying the amplitudes and phases of the so-called weights at the receiver array [6], [7]. Thus, the array can detect and process

ED

signals arriving from different directions. The array pattern synthesis problem is to find weights that satisfy a set of specifications on the beampattern. Although Lebret and Boyd have addressed this topic via convex optimization [7],

PT

such as minimization of the beampattern level over a given sidelobe zone, there exist a large number of more difficult

CE

nonconvex array synthesis problems in practice [20]. 3) Through transmission of probing waveforms toward an inspected region and analysis of the reflected signals, target or

AC

propagation medium estimation can be achieved. The probing sequence plays an important role in such active sensing applications [8], [9], [10], [11], [12], [21], [22], [23], [24]. A well-designed waveform set should possess impulse-like auto-correlation functions so that correlations among the echoes with different time delays are minimized to improve the target detection performance. Based on the fact that the auto-correlation function and power spectral density form a Fourier transform pair, a flat-spectrum sequence is desired [10], [8]. To enable the spectrum of the designed sequence to be as flat as possible, we need to minimize the ripple term in the frequency domain. Also, in multiple-input multipleoutput (MIMO) radar systems [9], waveforms with flexible transmit beampatterns are required [25], [26]. As stated in [27], important issues such as ripples within the energy focusing section and attenuation of the sidelobes should be

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

3

considered. 4) FIR filter has always been one of the prominent building blocks in digital signal processing. To ensure the passband to be as flat as possible while the stopband as low as possible [28], [34], upper and lower bounds are applied on the corresponding magnitudes of the frequency responses [29]. Although the spectral masks have been employed to enable the designed filter to have the desired property, there is room for further improvement via ripple or stopband level control,

CR IP T

resulting in nonconvex optimization problems. To deal with the aforementioned nonconvex problems, the semidefinite relaxation (SDR) [15], [30] technique typically transforms the original vector-variable problem to a high-dimensional matrix-variable problem without the rank-1 constraint, followed by a postprocessing procedure with randomization or eigenvalue decomposition strategy, which possibly violates the nonconvex constraints. To improve the rank-1 approximation accuracy, matrix rank minimization in the objective function can

AN US

be employed [39]. The spectral factorization method [38], [42] is a well-developed technique for uniformly sampled sequences or uniform linear arrays only, where a Toeplitz matrix is utilized to formulate an equivalent problem on the autocorrelation function. The Toeplitz property is also applied in [43] and [44]. In order to reduce the number of optimization variables, avoid possibly inaccurate rank-1 approximation under the nonconvex bounded constraints and rank minimization objective functions,

M

and be applicable to arbitrary array configurations or non-uniformly sampled sequences, this paper focuses on deriving a new solution based on the alternating direction method of multipliers (ADMM) [31], [32], which blends the decomposability of dual

ED

ascent with the superior convergence property of the multiplier method. Our approach has wider applicability as it can solve signal processing problems with bounded magnitude constraints on frequency or angular responses, including the nonconvex

PT

constrained sidelobe, ripple, and output power controls. Note that we have recently developed a solution to the unimodular sequence design problem via linear and quadratic objective subfunction separation [24]. Aspects of the proposed approach are

CE

highlighted as follows:

1) We apply the variable splitting trick to introduce two-element auxiliary variables for each frequency or angular response

AC

corresponding to each magnitude constraint in order to enable the original variables to get out of the complicated nonconvex constraints but maintain the equivalence between the former and transformed optimization problems. Then the newly-formed optimization problem is solved in an alternating manner [31], i.e., addressing two subproblems with respect to (w.r.t.) the original and auxiliary optimization variables.

2) More specifically, each new constraint may act on only the corresponding two-element vector variables rather than all of them. We then derive the relationship between the optimal objective function value of the nonconvex subproblem [14] and the norm of the introduced vector variables, and simplify the subproblem as a function of a single variable.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

4

3) According to the turning points for the single variable, we design special piecewise functions and develop feasible methods to determine the unknown parameter. 4) Since the proposed method does not rely on the uniform array (or uniform sampling interval for filters) or Toeplitz matrix assumption, it is applicable to arbitrary array configurations and non-uniformly sampled sequences. It is also more computationally attractive than the SDR approach because fewer variables are employed.

in three aspects:

CR IP T

It is worth pointing out that our methodology is different from the consensus-ADMM (CADMM) optimization method [45]

1) Regarding problem formulation: In our work, the single sidelobe level or ripple term control variable in the objective function is coupled with the quadratic or unimodular constraints, whereas both the objective function and constraints in

AN US

[45] are restricted to be quadratic terms (see (1) in [45]);

2) Regarding auxiliary variables: We introduce two-element auxiliary variables for each frequency or angular response, whereas [45] uses auxiliary vector variables with the same high dimension as the original optimization problem for each quadratic constraint;

3) Regarding optimization strategies: We transform the challenging subproblem of the ripple or sidelobe level variable with

M

complicated constraints into an optimization problem of only a single variable without constraints by introducing the step functions. Furthermore, we determine the unknown parameter via finding the global minimum from the local minima

ED

of specially designed piecewise functions. On the other hand, [45] utilizes the strategies of the consensus or distributed optimization [14] and nonlinear equation root finding for Lagrange multiplier [4], [16].

PT

The rest of this paper is organized as follows. The general problem with magnitude response constraints is formulated and its

CE

different signal processing applications are illustrated in Section 2. The proposed method is developed in Section 3. Numerical examples are included in Section 4. Finally, conclusions are drawn in Section 5.

and

T

,

AC

Notation: Vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively. The E{ }, k k, H

, stand for the expectation, Frobenius norm, transpose, and Hermitian transpose, respectively. The In denotes the n × n

identity matrix while <{} and ={} mean the real and imaginary parts, respectively. 2. P ROBLEM F ORMULATION AND M OTIVATION In this section, different magnitude response constraints are first introduced, followed by a general constrained optimization problem. In order to motivate the importance of this problem formulation, we then relate it to representative signal processing applications.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

5

Generally, there are three types of magnitude response constraints: left-sided (LS), right-sided (RS) and two-sided (TS), denoted by ≤ •, • ≤, and ≤ • ≤, respectively. RS constraints are convex while the LS and TS constraints are nonconvex [18]. Since the approach developed in this work is suitable for arbitrary magnitude response constraints or their mixtures, we address the following general problem: min f (w, L(θn ), U (θn )) w

CR IP T

s.t. L(θn ) E |wH a(θn )|2 E U (θn ), n = 1, · · · , N,

(1)

where w ∈ CM ×1 is the parameter or weight vector to be estimated, L(θn ) E |wH a(θn )|2 E U (θn ) are the so-called bounded magnitude (inequality) constraints, a(θn ) is the steering vector or discrete Fourier transform (DFT) coefficient vector, θn are the grid points of the region of interest (such as frequency and angle), and L(θn ) ≥ 0 and U (θn ) > 0 stand for the lower and upper

AN US

bounds of the nonconvex constraints, n = 1, · · · , N . Note that w ∈ RM ×1 is considered as a special case of (1). In addition, L(θn ) E |wH a(θn )|2 E U (θn ) refers to one of the following constraints, namely, |wH a(θn )|2 ≤ U (θn ), L(θn ) ≤ |wH a(θn )|2 , or L(θn ) ≤ |wH a(θn )|2 ≤ U (θn ).

The real-valued objective function f (w, L(θn ), U (θn )) can be a function of w only, and a representative example is [1]: (2)

M

f (w) = wH Rw,

ED

which denotes the output power of an array with R ∈ CM ×M being the corresponding sample covariance matrix. On the other hand, the objective can be functions of L(θn ) and/or U (θn ), such as in ripple control [35]:

PT

f (L(θn ), U (θn )) = U (θn ) − L(θn ), θn ∈ Mainbeam or Passband,

(3)

AC

CE

or sidelobe level control [28]:

f (U (θn )) = U (θn ), θn ∈ Sidelobe or Stopband.

(4)

Note that L(θn ) and U (θn ) are functions of θn in (1), but we will mainly discuss its simplified forms L and U , which are constants, in the rest of this paper. Nevertheless, the developed methodology can be straightforwardly extended to the more complicated case with L(θn ) and U (θn ). In what follows, we present several representative applications which are specific forms of (1).

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

6

A. Robust Narrowband Beamforming The Capon beamformer [3] cannot work properly under various mismatch conditions, and state-of-the-art proposals include the exceeding unity [17], [37] and two-sided [18], [19] magnitude response constraints to improve the robustness in the main beam to prevent target signal cancellation. For the former, the array output power is minimized under I exceeding unity constraints, i.e.,

CR IP T

min wH Rw w

s.t.

|wH a(θi )|2 ≥ 1,

i = 1, · · · , I;

(5)

whereas the latter introduces two-sided magnitude response constraints to the I discrete grids {θi }Ii=1 of the region of interest

min wH Rw w

s.t.

AN US

[θl , θu ] where the signal arrives with a high probability:

L(θi ) ≤ |wH a(θi )|2 ≤ U (θi ),

B. Array Beampattern Synthesis

i = 1, · · · , I.

(6)

M

With I magnitude constraints {L ≤ |wH a(θ˜i )|2 ≤ U }Ii=1 on the I angle grids in the shaped beam region and K constraints {|wH a(θk )|2 ≤ ε}K k=1 on the sidelobe region [6], [7], the array beampattern synthesis problem is also nonconvex due to the

ED

LS constraints L ≤ |wH a(θ˜i )|2 . It is worth pointing out that since the rank minimization of the matrix in the SDR formulation [20] is used as the objective function [39], other merits, e.g., sidelobe level minimization and mainband ripple minimization,

PT

cannot be used as the objective functions, and thus such a problem formulation is unhelpful to find the minimal sidelobe or ripple level. In fact, there may exist much adjustable room for the sidelobe level or mainlobe ripple level control. Therefore, in

CE

order to maximally suppress interferences, the sidelobe level minimization problem with the two-sided shaped beam constraints

AC

can be formulated as:

min ε w,ε

s.t.

L ≤ |wH a(θ˜i )|2 ≤ U, |wH a(θk )|2 ≤ ε,

i = 1, · · · , I

k = 1, · · · , K.

(7)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

7

The mainlobe ripple minimization problem with the sideband level constraints can be formulated as: min η w,η

s.t.

U +L U +L − η ≤ |wH a(θ˜i )|2 ≤ + η, 2 2 |wH a(θk )|2 ≤ ε,

i = 1, · · · , I (8)

k = 1, · · · , K,

CR IP T

to achieve robustness against the imprecise incoming target angles in radar systems [18] or estimate the channel gain after transmit/receive beam alignment in millimeter wave systems [33], where η is the ripple level.

C. Waveform Design

Denote x = [x0 · · · xN −1 ]T as the discrete-time periodic probing sequence in each cycle where N is the fundamental

AN US

period. When the designed waveform x has an impulse-like autocorrelation function [8], [10], [11], the spectrum is flat. To enable the spectrum of the designed waveform to be as flat as possible, we consider the following problem with ripple control: min η x,η

s.t. N − η ≤ xH fnH fn x ≤ N + η, n = 0, · · · , N − 1, n = 0, · · · , N − 1,

M

|xn | = 1,

(9)

2π

ED

where fn ∈ C1×N is the DFT coefficient vector with the kth element fn (k) = e−j N kn , k = 0, · · · , N − 1. Additionally, η here denotes the maximum deviation from N . Note that in (9) there are also nonconvex unimodular constraints for the best

CE

D. FIR Filter Design

PT

power efficiency, namely, |xn | = 1, n = 0, · · · , N − 1. Note that the model in (9) is different from (11) of [24].

With a set of desired constraints on the frequency response of the FIR filter coefficient vector w ∈ RN ×1 [28], [29], [34], p = 1, · · · , P on the squared magnitude in the passband and S RS constraints

AC

i.e., P TS constraints L ≤ wT ˜fpH ˜fp w ≤ U, wT fsH fs w ≤ ε,

s = 1, · · · , S, on the stopband, minimization of the stopband energy, stopband level, and ripple level can

be achieved, respectively, using the following formulations:

min w

s.t.

S X

wT fsH fs w

s=1

L ≤ wT ˜fpH ˜fp w ≤ U, wT fsH fs w ≤ ε,

p = 1, · · · , P

s = 1, · · · , S,

(10)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

8

min

ε

w,ε

s.t.

L ≤ wT ˜fpH ˜fp w ≤ U, wT fsH fs w ≤ ε,

p = 1, · · · , P (11)

s = 1, · · · , S,

min w,η

s.t.

η U +L U +L − η ≤ wT ˜fpH ˜fp w ≤ + η, 2 2 wT fsH fs w ≤ ε,

s = 1, · · · , S.

CR IP T

and

p = 1, · · · , P

(12)

AN US

Other relevant applications of (1) include MIMO radar transmit beampattern design with ripple control [27], radar code design [40], transmit B1 shimming [15], and physical layer multicasting [41].

3. P ROPOSED A PPROACH

In this section, we present an ADMM solution framework for the problem with magnitude response constraints as shown

M

in (1). We consider three kinds of objective functions: stopband level ((7),(11)), ripple level ((8), (9), (12)), and output power

ED

level ((5), (6), (10)). For each of them, a different strategy is devised to solve the corresponding subproblem.

A. Objective Function for Sidelobe Control

PT

The problem of array beampattern synthesis in (7) is considered as an illustration. To solve it, we exploit the idea that multiple inequality constraints should be associated with different optimization variables so that they can be tackled separately.

AC

CE

Introducing real-valued two-element vector variables ˘ ∈ R2×1 , yi = A(θ˜i )w

i = 1, · · · , I,

(13)

˘ ∈ R2×1 , zk = B(θk )w

k = 1, · · · , K,

(14)

where







 <{w}  , ˜ = w   ={w}

 ˜ ˜  <{a(θi )} ={a(θi )}   , i = 1, · · · , I, A(θ˜i ) =    ˜ ˜ −={a(θi )} <{a(θi )}

(15)

(16)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

9

and B(θk ), k = 1, · · · , K, are defined in a similar manner, (7) becomes ε

˘ s.t. yi = A(θ˜i )w, ˘ zk = B(θk )w, L ≤ kyi k2 ≤ U, kzk k2 ≤ ε,

i = 1, · · · , I, k = 1, · · · , K, i = 1, · · · , I,

k = 1, · · · , K.

CR IP T

min

˘ w,ε,y i ,zk

(17)

˘ 2 ≤ U and kB(θk )wk ˘ 2 ≤ ε become L ≤ kyi k2 ≤ U and kzk k2 ≤ ε, respectively. Thus, the Note that L ≤ kA(θ˜i )wk constraints L ≤ kyi k2 ≤ U and kzk k2 ≤ ε play their roles only in determining yi and zk , respectively. For this reason, we

AN US

define the augmented Lagrangian: ˘ ε, yi , zk , λ1,i , λ2,k ) L1,ρ (w, =ε+

I X i=1

+

K X


ρ ˘ + kzk − B(θk )wk ˘ 2 λT2,k (zk − B(θk )w) 2

M

k=1


ρ ˘ 2 ˘ + kyi − A(θ˜i )wk λT1,i (yi − A(θ˜i )w) 2

(18)

ED

where λ1,i ∈ R2×1 and λ2,k ∈ R2×1 are the Lagrange multiplier vectors, while ρ > 0 is a user-defined parameter representing the primal step size. We then determine all unknowns in (17) as well as λ1,i and λ2,k using the iterative manner of the ADMM

PT

˘ ε, yi , zk , λ1,i , λ2,k ): from L1,ρ (w,

˘ with given {yi (t), zk (t), λ1,i (t), λ2,k (t), ε(t)} from Step 1: Solve w

AC

CE

min ε(t) ˘ w

+

I X i=1

+

K X

k=1


ρ ˘ + kyi (t) − A(θ˜i )wk ˘ 2 λT1,i (t)(yi (t) − A(θ˜i )w) 2


ρ ˘ + kzk (t) − B(θk )wk ˘ 2 λT2,k (t)(zk (t) − B(θk )w) 2

(19)

˘ (19) is rewritten as: After overlooking the terms independent of w, ˘ T R1 w ˘ + dT1 (t)w, ˘ min w ˘ w

(20)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

10

where I

R1 =

K

ρX T ˜ ρX T A (θi )A(θ˜i ) + B (θk )B(θk ), 2 i=1 2

(21)

k=1

−

I X i=1

K X

k=1


AT (θ˜i ) λ1,i (t) + ρyi (t)


BT (θk ) λ2,k (t) + ρzk (t) .

˘ is simply computed as: From (20), w 1 ˘ + 1) = − R−1 d1 (t). w(t 2 1 ˘ + 1), λ1,i (t), λ2,k (t)} from: Step 2: Solve {zk , yi , ε} with given {w(t ε+

zk ,yi ,ε

i=1

I X ρ

2

i=1

+ +

K X

k=1 K X

˘ + 1)) λT1,i (t)(yi − A(θ˜i )w(t

˘ + 1)k2 kyi − A(θ˜i )w(t

˘ + 1)) λT2,k (t)(zk − B(θk )w(t

k=1

ρ ˘ + 1)k2 kzk − B(θk )w(t 2

M

+

I X

(23)

AN US

min

(22)

CR IP T

d1 (t) = −

L ≤ kyi k2 ≤ U, i = 1, · · · , I,

s.t.

ED

kzk k2 ≤ ε, k = 1, · · · , K.

(24)

Similarly, after removing the irrelevant terms, (24) is rewritten as: I

PT min

s.t.

K

ρX ρX ˜ i (t)k2 + ˜k (t)k2 kyi − y kzk − z 2 i=1 2 k=1

L ≤ kyi k2 ≤ U, i = 1, · · · , I, kzk k2 ≤ ε, k = 1, · · · , K,

AC

CE

zk ,yi ,ε

ε+

(25)

˜ i (t) = A(θ˜i )w(t ˘ + 1) − ρ1 λT1,i (t) and z ˜k (t) = B(θk )w(t ˘ + 1) − ρ1 λT2,k (t). where y Note that yi is independent of both zk and ε in (25). Therefore, (25) can be separated into two optimization subproblems: I

min yi

s.t.

ρX ˜ i (t)k2 kyi − y 2 i=1 L ≤ kyi k2 ≤ U, i = 1, · · · , I,

(26)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

11

and K

min ε + zk ,ε

s.t.

ρX ˜k (t)k2 kzk − z 2 k=1

kzk k2 ≤ ε, k = 1, · · · , K.

(27)

More importantly, (26) can be decomposed into I subproblems: ˜ i (t)k2 min kyi − y s.t.

yi (t + 1) =

              

L ≤ kyi k2 ≤ U,

√ U ˜ k˜ yi (t)k yi (t), √ L ˜ k˜ yi (t)k yi (t),

k˜ yi (t)k >

√

k˜ yi (t)k <

(28)

U

√

.

L

(29)

AN US

for i = 1, · · · , I, and the solution is given by:

CR IP T

yi

˜ i (t), y

√

L ≤ k˜ yi (t)k ≤

√

U

˜k (t)k2 However, in (27), ε is coupled with zk for k = 1, · · · , K. Similar to (28), we can draw a conclusion that kzk − z ˜k (t). depends on the relationship between ε and z

ED

M

For ease of exposition, we define the step function as follows:   
 0, ˜k (t) = S1 ε, z    1,

k˜ zk (t)k ≤ k˜ zk (t)k >

√ √

ε

(30)

ε.

˜k (t)k2 has the following form with given ε: According to (29), it can be readily shown that kzk − z

CE

PT

˜k (t)k2 kzk − z   √   0, k˜ zk (t)k ≤ ε = .  √ √   ( ε − k˜ zk (t)k)2 , k˜ zk (t)k > ε

(31)

AC

Therefore, (27) can be rewritten in another form as: min ε + ε

K
√ ρX ˜k (t) ( ε − k˜ S1 ε, z zk (t)k)2 , 2 k=1

which becomes a function of a single variable, namely, ε.

(32)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

12

Define [kx1 (t)k · · · kxK (t)k] as the ascending order of [k˜ z1 (t)k · · · k˜ zK (t)k], (32) can be expressed in a piecewise function as: K
√ ρX ˜k (t) k ε − k˜ S1 ε, z zk (t)k)2 2  k=1  PK √ √   0 ≤ ε ≤ kx1 (t)k ε + ρ2 k=1 ( ε − kxk (t)k)2 ,       P √ √  2  ε+ ρ K ε ≤ kx2 (t)k k=2 ( ε − kxk (t)k) , kx1 (t) ≤ 2

=

.. .

          

Furthermore, the kth function (i.e., kxk (t)k ≤

ε, √

√

1, Bk = −ρ

PK

l=k

√

ε

ε ≤ kxk+1 (t)k) of the K + 1 piecewise functions can be written in a

AN US

that

l=k

(34)

kxK (t)k ≤

√ Ak ε + Bk ε + Ck , PK

(33)

.. .

common form as:

where Ak = 1 +

.

CR IP T

ε+

kxk (t)k, and Ck =

ε is within the kth subregion (i.e., kxk (t)k ≤

(kx1 (t)k, kx2 (t)k], · · · , (kxK (t)k, +∞).

√

ρ 2

PK

l=k

kxk (t)k2 for k = 1, · · · , K, and “k” denotes

ε ≤ kxk+1 (t)k) of all the K + 1 subregions [0, kx1 (t)k],

PT C−

B2 4A ,

Akxk (t)k2 + Bkxk (t)k + C,        Akxk−1 (t)k2 + Bkxk−1 (t)k + C,

CE

vˆk =

       

ED

M

Therefore, for the kth segment, the optimal value of ε and the corresponding minimal value are:    B2 B  − 2A ∈ [kxk−1 (t)k, kxk (t)k]  4A2 ,    εˆk = kxk (t)k2 , 2Akxk (t − 1)k + B < 0 and 2Akxk (t)k + B < 0 ,        kxk−1 (t)k2 , 2Akxk (t − 1)k + B > 0 and 2Akxk (t)k + B > 0

(35)

B − 2A ∈ [kxk−1 (t)k, kxk (t)k]

2Akxk (t − 1)k + B < 0

and

2Akxk (t)k + B < 0 .

2Akxk (t − 1)k + B > 0

and

2Akxk (t)k + B > 0

(36)

AC

By selecting the minimum from all of the K + 1 minimal values of the K + 1 segments, denoted as vˆi = mink (ˆ v1 , · · · , vˆk , · · · , vˆK+1 ), we determine the optimal ε(t + 1) as: ε(t + 1) = εˆi .

(37)

zk (t + 1)  √ p  ε(t+1)   ˜ zk (t)k ≥ ε(t + 1) k˜ zk (t)k zk (t), k˜ = , p    ˜k (t), z k˜ zk (t)k < ε(t + 1)

(38)

Once ε(t + 1) is obtained, the optimal zk is then:

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

13

for k = 1, · · · , K. Step 3: Update λ1,i and λ2,k as:


˘ + 1) , λ1,i (t + 1) = λ1,i (t) + ρ yi (t + 1) − A(θ˜i )w(t

(39)

and

CR IP T


˘ + 1) . λ2,k (t + 1) = λ2,k (t) + ρ zk (t + 1) − B(θk )w(t

(40)

Via iterating Steps 1 to 3, the proposed solution for the problem of array pattern synthesis with sidelobe control is summarized in Algorithm 1: Algorithm 1: Sidelobe Control for Array Pattern Synthesis

Initialization: {yi (0)},{zk (0)}, {λ1,i (0)}, and {λ2,k (0)}; for t = 0, . . . , T ,where T is predefined maximum iteration number ˘ + 1) using (23); Obtain w(t Determine yi (t + 1) using (29); Determine ε(t + 1) and zk (t + 1) using (35) and (38); Update λ1,i (t + 1) and λ2,k (t + 1) using (39) and (40);

AN US

———————————————————————————————————–

M

˘ + 1)k < δ where δ > 0 is tolerance parameter. end for t = T or maxn kyn (t + 1) − A(θn )w(t

ED

˘ Obtain weight vector w from w.

PT

B. Objective Function for Ripple Control

For the ripple control problems (8),(9) and (12), we use waveform design in (9) as a representative example. We first

AC

CE

transform (9) into the following real-valued optimization problem: min η ¯ ,η x

s.t. x ¯2n + x ¯2n+N = 1,

n = 0, · · · , N − 1,

¯ T FTn Fn x ¯ ≤ N + η, n = 0, · · · , N − 1, N −η ≤x η ∈ [0, ηU ],

(41)

¯ and Fn are the real-valued vector and matrix of x and fn as in (15) and (16) respectively, and x where x ¯n represents the nth ¯ , n = 0, · · · , 2N − 1. Note that here we also include the constraint η ∈ [0, ηU ] where ηU is the given upper bound element of x to ensure that the optimal η is determined from the feasible region.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

14

Again, based on the variable splitting concept, we define real-valued two-element vector variables: ¯ ∈ R2×1 , v n = Fn x

(42)

n = 0, · · · , N − 1.

Then, (41) can be rewritten in an equivalent form as: η

s.t. x ¯2n + x ¯2n+N = 1,

n = 0, · · · , N − 1,

¯ , n = 0, · · · , N − 1, vn = Fn x N − η ≤ kvn k2 ≤ N + η,

CR IP T

min

¯ ,η,vn x

n = 0, · · · , N − 1,

η ∈ [0, ηU ].

(43)

AN US

¯2n+N = 1, N − η ≤ kvn k2 ≤ N + η, and η ∈ [0, ηU ] play their roles only in the determination Similarly, the constraints x ¯2n + x ¯ , vn , and η. From (43), the augmented Lagrangian is: of variables x ¯ , η, vn , λ3,n L2,ρ x =η+

N −1 X n=0





ρ ¯ ) + kvn − Fn x ¯ k2 , λT3,n (vn − Fn x 2

(44)

M

where λ3,n ∈ R2×1 is the Lagrange multiplier vector. Then, based on ADMM, we determine {¯ x, η, vn , λ3,n } via the following

ED

alternating iteration steps:

PT

¯ with given {η(t), {vn (t)}, {λ3,n (t)}} from: Step 1:Determine x

AC

CE

¯ , η(t), {vn (t)}, {λ3,n (t)} ¯ (t + 1) = arg min L2,ρ x x = arg min η(t) +

N −1 X n=0





ρ ¯ ) + kvn (t) − Fn x ¯ k2 λT3,n (t)(vn (t) − Fn x 2

s.t. x ¯2n + x ¯2n+N = 1,

n = 0, · · · , N − 1.

(45)

¯ , (45) can be simplified as: Ignoring the irrelevant terms to x ¯ T R2 x ¯ + dT2 (t + 1)¯ min x x ¯ x

s.t. x ¯2n + x ¯2n+N = 1, where d2 (t + 1) = −

PN −1 n=0


FTn λ3,n (t) + ρFTn vn (t) .

n = 0, · · · , N − 1,

(46)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

Since R2 =

PN −1 n=0

15

¯ T R2 x ¯ = N 2 , (46) becomes: FTn Fn = N IN and x min dT2 (t + 1)¯ x ¯ x

s.t. x ¯2n + x ¯2n+N = 1,

n = 0, · · · , N − 1,

(47)

pair {¯ xn , x ¯n+N } is easily derived as:    



x ¯n x ¯n+N

   

d2,n (t + 1)

d2,n+N (t + 1)   = −  

 d2,n (t + 1)





d

2,n+N (t + 1)

  

 ,

  

(48)

AN US

where d2,n (t + 1) is the nth element of d2 (t + 1), n = 0, · · · , 2N − 1.



CR IP T

which can be separated into N subproblems, for the variable pair {¯ xn , x ¯n+N }, n = 0, · · · , N − 1. The solution to the variable

Step 2: Solve {η(t + 1), vn (t + 1)} with given {¯ x(t + 1), λ3,n (t)} from: {η(t + 1), vn (t + 1)}


¯ (t + 1), η, vn , λ3,n (t) = arg min L2,ρ x

+

N −1 X


ρ ¯ (t + 1)) + kvn − Fn x ¯ (t + 1)k2 λT3,n (t)(vn − Fn x 2

ED

n=0

M

= arg min η

s.t. N − η ≤ kvn k2 ≤ N + η,

n = 0, · · · , N − 1, (49)

PT

η ∈ [0, ηU ].

CE

Let

1 ˜ n (t + 1) = Fn x ¯ (t + 1) − λ3,n (t), v ρ

n = 0, · · · , N − 1,

(50)

AC

(49) can be simplified as:

min η + η,vn

N −1 X n=0

ρ ˜ n (t + 1)k2 kvn − v 2

s.t. N − η ≤ kvn k2 ≤ N + η, η ∈ [0, ηU ].

n = 0, · · · , N − 1, (51)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

16

Apparently, η is coupled with vn due to the constraints N − η ≤ kvn k2 ≤ N + η. Hence once η is provided, the optimal vn (t + 1) is given by:

k˜ vn (t + 1)k ≥ k˜ vn (t + 1)k ≤

√ √

N +η N −η

(52)

otherwise.

CR IP T

vn (t + 1)  √   N +η  ˜  k˜ vn (t+1)k vn (t + 1),    √ N −η = ˜ k˜ vn (t+1)k vn (t + 1),        ˜ n (t + 1), v

˜ n (t + 1)k2 can be represented as a function of η: Additionally, the term kvn − v

(53)

S2 (˜ vn (t + 1), N + η)   √   1, k˜ vn (t + 1)k ≥ N + η =    0, otherwise,

(54)

S3 (˜ vn (t + 1), N − η)   √   1, k˜ vn (t + 1)k ≤ N − η =    0, otherwise.

(55)

AN US

˜ n (t + 1)k2 kvn − v   √ √   ( N + η − k˜ vn (t + 1)k)2 , k˜ vn (t + 1)k ≥ N + η     √ √ = vn (t + 1)k)2 , k˜ vn (t + 1)k ≤ N − η ( N − η − k˜        0, otherwise.

Based on the idea above, we convert (51) into an optimization problem with a single variable η. For simplicity, we define

CE

PT

ED

M

two step functions as follows:

AC

With the use of the step functions, (51) is rewritten as: min η η

+ + s.t.

N −1 p
2 ρ X S2 (˜ vn (t + 1), N + η) N + η − k˜ vn (t + 1)k 2 n=0

N −1 p
2 ρ X S3 (˜ vn (t + 1), N − η) N − η − k˜ vn (t + 1)k 2 n=0

η ∈ [0, ηU ],

(56)

which implies that the objective function in (56) is actually a piecewise function that depends on the values of the step functions characterized by η.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

17

To solve (56), the reasonable and nonredundant K turning points (for simplicity, K denotes the number of turning points) {k˜ vn (t+1)k2 −N, N −k˜ vn (t+1)k2 }, which belong to (0, ηU ), are determined in the ascending order as [η1 (t+1) · · · ηK (t+1)]. Then, we can express the objective function in (56) as K + 1 subfunctions. Especially, when η ∈ [ηk−1 (t + 1), ηk (t + 1)], the subfunction has the same form as that of (56) but with the constraint η ∈ [ηk−1 (t + 1), ηk (t + 1)]. We now analyze the convexity of subfunction in the region [ηk−1 (t + 1), ηk (t + 1)]. The first and second derivatives of (56) w.r.t. η are given by: N −1 ρ X k˜ vn (t + 1)k ) S2 (˜ vn (t + 1), N + η)(1 − √ 2 n=0 N +η

+

N −1 ρ X k˜ vn (t + 1)k ) S3 (˜ vn (t + 1), N − η)(−1 + √ 2 n=0 N −η

and

AN US

N −1 ρ X k˜ vn (t + 1)k S2 (˜ vn (t + 1), N + η) 3 4 n=0 (N + η) 2

+

CR IP T

1+

N −1 k˜ vn (t + 1)k ρ X S3 (˜ vn (t + 1), N − η) > 0. 3 4 n=1 (N − η) 2

(57)

(58)

Obviously, each subfunction is convex, and we obtain the optimal value of the kth subfunction for the following cases: 1) If the first derivative of (56) in η ∈ [ηk−1 (t+1), ηk (t+1)] w.r.t. η is negative (or positive), then the function is decreasing

M

(or increasing) in the region, and ηk (t + 1) (or ηk−1 (t + 1)) is the optimal value of η, denoted as ηˆ(k); 2) Otherwise, the first derivatives of (56) in η ∈ [ηk−1 (t + 1), ηk (t + 1)] w.r.t. η at ηk−1 (t + 1) and ηk (t + 1) are negative

ED

and positive, respectively. Let ∆ > 0 be a small tolerance parameter (e.g., 10−8 ), we can obtain the root ηˆ(k) for the first derivative of (56) w.r.t. η by applying a simple bisection method in the region η ∈ [ηk−1 (t + 1), ηk (t + 1)].

PT

Inserting ηˆ(k) into the objective function in (56) with η ∈ [ηk−1 (t + 1), ηk (t + 1)], we obtain the potential minimal value of the kth subfunction, denoted as w(k). ˆ By selecting the minimal value from all of the K + 1 subfunction values of the K + 1

CE

segments, denoted by w ˆi = min {w ˆ1 , · · · , w ˆk , · · · , w ˆK+1 }, we determine the optimal η(t + 1): (59)

AC

η(t + 1) = ηˆi .

Once η(t + 1) is determined, the optimal vn (t + 1) is computed from (52). Step 3: Update λ3,n as:
¯ (t + 1) , λ3,n (t + 1) = λ3,n (t) + ρ vn (t + 1) − Fn x n = 0, · · · , N − 1.

(60)

Steps 1-3 are repeated until a stopping criterion is reached. The proposed procedure for waveform design with magnitude response constraints is summarized in Algorithm 2.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

18

Algorithm 2: Ripple Control for Waveform Design ———————————————————————————————————– ¯ (0), {vn (0)}, and {λ3,n (0)}; Initialization: x for t = 0, . . . , T ,where T is predefined maximum iteration number ¯ (t + 1) using (48); Obtain x Determine η(t + 1) using (59); Determine vn (t + 1) according to the obtained η(t + 1) and (52);

¯ (t + 1)k < δ where δ > 0 is tolerance parameter. end for t = T or maxn kvn (t + 1) − Fn x ¯. Obtain waveform sequence x from x

C. Objective Function for Output Power Minimization

CR IP T

Update λ3,n (t + 1) using (60);

We now address the problem in (1) where the objective function has the form of wH Rw as in robust adaptive beamforming

applied to solve the resultant optimization problem.

AN US

and FIR filter design, as shown in (6) and (10). We first utilize the variable splitting idea to simplify (1). Then, ADMM is

˜ R3 , and E(θn ) as the real-valued vector or matrix forms of w, R, and Consider real-valued implementation and define w, a(θn ), respectively, as in (15)–(16). Then we have

˜ T R3 w, ˜ wH Rw = w

M

(61)

ED

and |wH a(θn )|2

PT

= <2 {wH a(θn )} + =2 {wH a(θn )}

CE

˜ = [<(a(θn )) =(a(θn ))]w ˜ 2. = kE(θn )wk


2

˜ + [−=(a(θn )) <(a(θn ))]w


2

(62)

AC

In doing so, (6) can be represented as a real-valued optimization problem: ˜ T R3 w ˜ min w ˜ w

˜ 2 E U, n = 1, · · · , N. s.t. L E kE(θn )wk

(63)

Define N two-element vector variables, i.e., ˜ ∈ R2×1 , un = E(θn )w

n = 1, · · · , N.

(64)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

19

With more optimization variables, (63) becomes: ˜ T R3 w ˜ w

min

˜ w,{u n}

s.t.

n = 1, · · · , N,

˜ un = E(θn )w,

L E kun k2 E U, n = 1, · · · , N.

(65)

Thus, the augmented Lagrangian of (65) is: ˜ un , λ4,n ) = w ˜ T R3 w ˜+ L3,ρ (w,

N X

n=1

CR IP T

˜ for n = 1, · · · , N . We introduce the Lagrange multiplier vector λ4,n ∈ R2×1 corresponding to the constraint un = E(θn )w

˜ λT4,n un − E(θn )w

N X ρ ˜ 2. + kun − E(θn )wk 2 n=1




(66)

AN US

˜ and un alternately from (66): With initializations {un (0), λ4,n (0)}, we determine w

˜ + 1) is determined via solving the following optimization subproblem: Step 1: With given un (t) and λ4,n (t), w(t ˜ T R3 w ˜ min w ˜ w

+

N X

n=1


ρ ˜ 2 . ˜ + kun (t) − E(θn )wk λT4,n (t)(un (t) − E(θn )w) 2

M

Removing the independent terms, (67) is rewritten as:

˜ T R4 w ˜ + 2dT3 (t)w, ˜ min w

(68)

ρ 2

PN

ED

˜ w

where R4 = R3 +

T n=1 E (θn )E(θn ) and d3 (t) = −

CE

PT

The solution to (68) is simply:

PN

n=1

h

1 T 2 λ4,n (t)E(θn )

(67)

iT + ρ2 uTn (t)E(θn ) .

˜ + 1) = −R−1 w(t 4 d3 (t).

(69)

AC

˜ + 1) and λ4,n (t), un (t + 1) is determined via solving the following optimization subproblem: Step 2: With given w(t ˜ T (t + 1)R3 w(t ˜ + 1) min w

{un }

N X

+ +

n=1 N X

n=1

˜ + 1) λT4,n (t) un − E(θn )w(t ρ ˜ + 1)k2 . kun − E(θn )w(t 2


(70)

Similarly, considering the inequality constraint L E kun k2 E U, n = 1, · · · , N , (70) is rewritten as: min

{un }

s.t.

N X

n=1

˜ n (t + 1)k2 , kun − u

L E kun k2 E U, n = 1, · · · , N,

(71)

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

20

˜ n (t + 1) = E(θn )w(t ˜ + 1) − ρ1 λ4,n (t). where u This problem can be divided into N subproblems, i.e., ˜ n (t + 1)k2 , min kun − u un

s.t.

L E kun k2 E U,

(72)

for n = 1, · · · , N .

CR IP T

When the constraint L E kun k2 E U corresponds to L ≤ kun k2 ≤ U , L ≤ kun k2 , or kun k2 ≤ U , the solution to (72) is, respectively, given by:

(73)

un (t + 1)  √ √   L  ˜ k˜ un (t + 1)k ≤ L k˜ un (t+1)k un (t + 1), = √    ˜ n (t + 1), u k˜ un (t + 1)k > L,

(74)

un (t + 1)  √ √   U  ˜ k˜ un (t + 1)k ≥ U k˜ un (t+1)k un (t + 1), = √    ˜ n (t + 1), u k˜ un (t + 1)k < U .

(75)

M

AN US

un (t + 1)  √ √   U  ˜ k˜ un (t + 1)k ≥ U  k˜ un (t+1)k un (t + 1),    √ √ L = ˜ u (t + 1), k˜ u (t + 1)k ≤ L n n k˜ u (t+1)k  n    √ √    ˜ n (t + 1), L < k˜ un (t + 1)k < U , u

CE

PT

ED

or

AC

˜ + 1), λ4,n is updated as: Step 3: With given un (t + 1) and w(t


˜ + 1) . λ4,n (t + 1) = λ4,n (t) + ρ un (t + 1) − E(θn )w(t

(76)

Steps 1-3 are repeated until a stopping criterion is reached. The corresponding solution is summarized in Algorithm 3:

4. N UMERICAL E XAMPLES

MATLAB based computer simulations are conducted to evaluate the performance of the proposed approach. The problems of array pattern synthesis, waveform design, and robust adaptive beamforming are considered. Among the three examples, the first and third correspond to the angle response whereas the second is a frequency response problem. In our study, the values of the primal step size, maximum iteration number and tolerance are ρ = 0.1, T = 20000 and ∆ = 10−8 .

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

21

Algorithm 3: Output Power Minimization for Robust Beamforming ———————————————————————————————————– Initialization: {un (0)} and {λ4,n (0)}; for t = 0, . . . , T ,where T is predefined maximum iteration number ˜ + 1) using (69); Obtain w(t Determine un (t + 1) by (73) or (74) or (75); Update λ4,n (t + 1) using (76); ˜ + 1)k < δ where δ > 0 is tolerance parameter. end for t = T or maxn kun (t + 1) − E(θn )w(t

CR IP T

˜ . Obtain weight vector w

Note that the CADMM method [45] can solve the output power minimization problem, but it is not applicable for sidelobe and ripple level control problems due to the coupled level variable and weight variables in (7)–(12) as well as the nonconvex

AN US

unimodular constraint.

Example 1: Array Pattern Synthesis

In the first test, we consider the array pattern synthesis problem using a 10-element uniform linear array (ULA) with halfwavelength element spacing, where the element positions are [0, 1, · · · , 9] with unit of half-wavelength, and we need to find

M

w ∈ C10×1 . A flat-magnitude response in the mainlobe region [−15◦ , 15◦ ] is desired, namely, the maximum ripple level is 0.1 or U = 1.1 and L = 0.9. The magnitudes in the sidelobe regions [−90◦ , −26◦ ] and [26◦ , 90◦ ] should be no larger than 0.01.

ED

Here we discretize the mainlobe and sidelobe regions into 31 and 130 angular grids, respectively, with successive separation of 1◦ . The results of the proposed approach to synthesize the patterns in (7) and (8) are plotted in Figs. 1 and 2, where the

PT

latter shows the enlarged version of the mainlobe region, and the found minimal ripple level with the given sidelobe level and minimal sidelobe level with given ripple level are 0.0539 and −21.6362 dB, respectively.

CE

We have applied the SDR method [15] to solve the problems in (7) (sidelobe control problem) and (8) (ripple control problem). After the SDR procedure, the randomization strategy [15] is employed to generate scaled vectors as the weight w.

AC

However, even though the number of generated vectors is equal to 20000, none of them can satisfy all of the 31 two-sided magnitude response constraints, which implies that all are infeasible. The eigenvalue decomposition (EVD) strategy is also utilized for the rank-one approximation. However, the resultant vector still cannot satisfy all the two-sided constraints. To improve it, we use the rank minimization of the corresponding matrix variable as the objective function rather than the ripple minimization or sidelobe level minimization, and the corresponding result is marked with “SDR (Rank minimization))” [39] in Figs. 1 and 2. Note that there may exist some adjustable room for ripple term or sidelobe levels, which could not be found by “SDR (Rank minimization))” due to the objective function formulation [20].

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

22

0 −10

−30

CR IP T

−40 −50 −60

Proposed (Ripple control) Proposed (Sidelobe control) SDR (Rank minimization) SF (Ripple control) SF (Sidelobe control) Mask

−70 −80 −90 −100

−80

−60

−40

−20

AN US

Power pattern (dB)

−20

0

20

Angle (degree)

40

60

80

Fig. 1. Array pattern synthesis for uniform linear array: ripple control problems are subject to the sidelobe level constraints; whereas sidelobe control problems

M

are subject to the ripple term constraints.

Since the current array is a ULA, the corresponding matrices {a(θn )aH (θn )} and {a(θ˜i )aH (θ˜i )} are Toeplitz and the SF

ED

methods [38,42] are modified to solve the problems in (7) and (8). Figs. 1 and 2 indicate that the SF approach has similar performance as the proposed method in the ULA configuration case and they are the same for the sidelobe control problem.

PT

In the second test, we run these methods for array pattern synthesis with a nonuniform linear array (NLA) configuration, which is generated by modifying the last 9 element positions with a uniform distribution [−0.4, 0.4]. The simulation results are

CE

plotted in Fig. 3. The minimal ripple level (under the sidelobe level constraint) and minimal sidelobe level (under the ripple term constraint) measured by the proposed method are 0.1134 and −14.8985 dB, respectively. Note that the SDR method with

AC

the rank minimization as the corresponding objective function [39] still works for the infeasible problem (beyond the mainlobe mask) and we cannot set the minimal ripple level or sidelobe level as an objective function for SDR. Note that in the SDR implementation the number of the auxiliary matrix variables (i.e., M 2 ) is the square of that (i.e., M ) of the complex vector variables in the original problem. Whereas in our method, the number of the introduced auxiliary variables is just twice of the number of magnitude constraints. Apparently, the SDR approach is not suitable for large antenna arrays, which may result in “Out of Memory” problem in MATLAB. On the other hand, the SF approach fails in the NLA case due to the fact that the resulting matrices {a(θn )aH (θn )} and {a(θ˜i )aH (θ˜i )} are no longer Toeplitz.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

23

4 Proposed (Ripple control) Proposed (Sidelobe control) SDR (Rank minimization) SF (Ripple control) SF (Sidelobe control) Mask

2 1

CR IP T

Power pattern (dB)

3

0 −1

−3 −20

−15

−10

−5

M 0.5 0

PT

−10

−1

−20

0

5

10

15

20

Proposed (Ripple control) Proposed (Sidelobe control) SDR (Rank minimization) Mask

20

CE

Power pattern (dB)

−0.5

ED

20

0

0

Angle (degree)

Fig. 2. Zoomed mainlobe of Fig. 1

10

AN US

−2

AC

−20

−30

−40

−80

−60

−40

−20

0

20

Angle (degree)

40

60

80

Fig. 3. Array pattern synthesis for nonuniform linear array: SF cannot be applied; whereas SDR with rank minimization as the objective function cannot satisfy mainlobe constraints.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

24

Example 2: Periodic Waveform Design In the second example, we consider generating a periodic unimodular sequence of period N = 128 with nearly flat spectrum, as described in (9). A random phase sequence with length 128 is used for initialization. Comparison with several state-ofthe-art methods, namely, PEriodic Cyclic Algorithm New (PeCAN) [8]-[11], Lagrange Programming Neural Network (LPNN) [23], Monotonic Integrated Sidelobe Level (MISL) [12], and LInear and QUAdratic Subfunction Separation (LIQUASS) [24],

CR IP T

is made. It is worth pointing out that the SDR method [15] has also been tried to tackle the ripple control problem in (9). After the SDR procedure, the randomization and EVD strategies [15] are used for generating scaled vectors as the designed sequence. However, they cannot satisfy all the 128 two-sided constraints.

The autocorrelation performance is plotted in Fig. 4, and it is seen that the the maximal peak of autocorrelation sidelobe levels for the PeCAN, MISL, LPNN, LIQUASS, and proposed ripple control methods are −188.894, −129.787, −263.367,

AN US

−285.202, and −287.938 dB, respectively. To further assess the flatness of the spectrum of the designed sequence, we define the following error measure:

e(n) = |xH fnH fn x − N |, , n = 0, · · · , N − 1,

(77)

M

which is the absolute difference between the actual and ideal spectra. The corresponding results are shown in Fig. 5, where the maximal errors of the PeCAN, MISL, LPNN, LIQUASS, and proposed ripple control methods are 7.729 × 10−7 , 6.113 × 10−4 ,

ED

6.148 × 10−4 ,1.423 × 10−11 , and 9.109 × 10−12 , respectively. Again we see the superiority of the proposed ripple control

PT

scheme in terms of spectrum flatness.

Example 3: Robust Beamforming

CE

We now address the robust beamforming problem with a 10-element linear array. In the first test, the same ULA as in Example 1 is employed while we consider another nonuniform linear array of Example 1 in the second test. The DOA of the

AC

desired signal is 0◦ while those of two interferences are −40◦ and 60◦ . Both the signal-to-noise ratio and interference-to-noise ratio are 10 dB.

First, we consider robust beamforming with two-sided magnitude response constraints of (6) in the ULA configuration, where the upper and lower magnitudes are 0.95 and 1.05, respectively. Here, θ ∈ [−10◦ , 10◦ ] is the constrained mainlobe region. For comparison purposes, we implement the SDR [15], SF [18][44][45], and CADMM [43] methods. The result yielded from SDR [15] cannot satisfy the nonconvex constraints. Therefore, only the results from the SF, CADMM, and proposed methods are shown in Fig. 6. To see the mainlobe and null notches clearly, we enlarge the corresponding regions in Fig. 7. It is observed

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

25

0

PeCAN MISL LPNN LIQUASS Proposed

−100

CR IP T

Autocorrelation (dB)

−50

−150 −200

AN US

−250 −300 −100

−50

0 Time lag

50

100

ED

M

Fig. 4. Autocorrelation function of synthesized waveform in second example: the lower the sidelobe is, the better the designed waveform is.

PT

PeCAN MISL LPNN LIQUASS Proposed

AC

CE

Autocorrelation (dB)

−5

10

−10

10

0

20

40

60 80 Time lag

100

120

Fig. 5. Spectrum difference of synthesized waveform in second example: the smaller the error is, the better the designed waveform is.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

26

10

Proposed CADMM SF Mask

0

−20 −30

CR IP T

Beampattern (dB)

−10

−40 −50

−70

−80

−60

−40

−20

AN US

−60 0

20

Angle (degree)

40

60

80

Fig. 6. Robust beamforming with mainlobe width control for uniformly linear array: the SDR method fails in this case with 21 two-sided constraints.

M

that the proposed approach possesses lower responses at the interferences than [18], indicating that the former achieves better

ED

interference suppression in the ULA configuration.

Finally, we study the sensitivity of the methods to nonuniform linear arrays. However, the SDR method [15] cannot satisfy

PT

the constraints. Since the problem cannot be reformulated as that with respect to the auto-correlation of variables, the SF-type methods cannot work [38][42][43][44] due to the same reason as for Example 1. Therefore, only the results from the proposed

CE

and CADMM methods are plotted in Fig. 8. Comparing Figs. 6 and 8, it is seen that: i) the SF-type methods fail in suppressing the interferences in a nonuniform linear array configuration because it is designed for exact ULAs [44]; ii) the proposed scheme

AC

still works well for such a nonuniform linear array since its development is independent of the array configuration; iii) the beamforming task is to minimize the output power and thus the SDR method with the rank minimization as the objective function cannot work in this case; and (iv) although the CADMM and proposed method achieve almost the same performance, the former needs to introduce 420 (21 × 10 × 2) complex-valued weight vector and corresponding Lagrange multiplier vector variables resulted from 21 two-sided constraints. While the proposed method requires only 42 real-valued two-element vector variables, indicating its computational efficiency over [45].

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

27

0.1 0 −0.1

−62.7

−55

−62.75

−56

−62.8

−57

−62.85

−58 −59 −60 −61

−62.9

CR IP T

0.2

Beampattern (dB)

Beampattern (dB)

0.3

−54

Beampattern (dB)

Proposed CADMM SF Mask

−62.95 −63

−63.05

−62 −0.2

−63 0

−64 −41 −40 −39

10

Angle (degree)

60.1

M Proposed CADMM Mask

CE

PT

0

Beampattern (dB)

60

Angle (degree)

ED

10

−20

−63.2 59.9

Angle (degree)

Fig. 7. Zoomed mainlobe and notches of Fig. 6

−10

−63.15

AN US

−10

−63.1

0.2

−30

AC

0

−40

−0.2 −10

0

10

−50 −60

−80

−60

−40

−20

0

20

Angle (degree)

40

60

Fig. 8. Robust beamforming for nonuniformly linear array: both SF and SDR methods fail in this array configuration

80

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

28

5. C ONCLUSION In this paper, we have devised an ADMM framework for solving optimization problems with nonconvex magnitude constraints on frequency or angular responses. Comparing with existing solutions, our approach reduces the number of optimization variables, avoids rank-1 approximation under the nonconvex bounded constraints as well as rank minimization objective functions, and is applicable to arbitrary array and non-uniformly sample scenarios, indicating its computational efficiency

pattern synthesis, waveform design and robust beamforming.

ACKNOWLEDGEMENT

CR IP T

and wider applicability. The high performance of the proposed methodology is demonstrated via the applications in array

This work was supported in part by Natural Science Foundation of China (NSFC 61471295), and Central University Funds

AN US

(G2016KY0308,G2016KY0002,17GH030144), and Key Research Program of Frontier Sciences, CAS (QYZDY-SSW-JSC035).

R EFERENCES [1] J. Li, and P. Stoica. Robust Adaptive Beamforming. Hoboken, NJ: Wiley, 2006.

[2] A. Farina, Antenna Based Signal Processing Techniques for Radar Systems. Artech House, Inc., Norwood (MA), USA, 1992.

M

[3] J. Capon, “High-resolution frequency wavenumber spectrum analysis,” Proceedings of the IEEE, vol. 57, no. 8, pp. 1408-1418, Aug. 1969. [4] J. Li, P. Stoica, and Z. Wang, “On robust Capon beamforming and diagonal loading,” IEEE Transactions on Signal Processing, vol. 51, no. 7, pp.

ED

1702-1715, Jul. 2003.

[5] J. Li, P. Stoica, and Z. Wang, “Doubly constrained robust Capon beamforming,” IEEE Transactions on Signal Processing, vol. 52, no. 9, pp. 2407-2423, Sep. 2004.

PT

[6] F. Wang, V. Balakrishnan, P. Y. Zhou, J. J. Chen, R. Yang, and C. Frank, “Optimal array pattern synthesis using semidefinite programming,” IEEE Trans. Signal Processing, vol. 51, no. 5, pp. 1172-1183, 2003. [7] H. Lebret, and S. Boyd, “Antenna array pattern synthesis via convex optimization,” IEEE Trans. Signal Processing, vol. 45, no.3, pp. 526-531, 1997.

CE

[8] H. He, J. Li, and P. Stoica. Waveform Design for Active Sensing Systems: A Computational Approach. Cambridge, U. K.: Cambridge Univ. Press, 2012. [9] J. Li, and P. Stoica. MIMO Radar Signal Processing, Hoboken, NJ: Wiley, 2009.

AC

[10] P. Stoica, H. He, and J. Li, “On designing sequences with impulse-like periodic correlation,” IEEE Signal Processing Letters, vol. 16, no. 8, pp.703-706, 2009.

[11] P. Stoica, H. He, and J. Li, “New algorithms for designing unimodular sequences with good correaltion properties,” IEEE Trans. Signal Processing, vol. 57, no. 4, pp. 1415-1425, 2009.

[12] J. Song, P. Babu, and D. P. Palomar, “Optimization methods for designing sequences with low autocorrelation sidelobes,” IEEE Trans. Signal Processing, vol. 63, no. 15, pp. 3998-4009, Aug. 2015. [13] L. R. Rabiner, J. H. McClellan, and T. W. Parks, “FIR digital filter design techniques using weighted Chebyshev approximation,” Proceedings of the IEEE, vol. 63, no. 4, pp. 595-610, 1975. [14] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

29

[15] Z.-Q. Luo, W.-K. Ma, A. M.-C. So and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” IEEE Signal Processing Magazine, vol. 27, no. 3, pp. 20-34, May 2010. [16] R. G. Lorenz, and S. P. Boyd,“Robust minimum variance beamforming,” IEEE Trans. Signal Processing, vol. 53, no. 5, pp. 1684-1696, May 2005. [17] J. Xu, G. Liao, S. Zhu, and L. Huang, “Response vector constrained robust LCMV beamforming based on semidefinite programming,” IEEE Trans. Signal Processing, vol. 63, no. 21, pp. 5720-5732, Nov. 2015. [18] Z. L. Yu, M. Er, and W. Ser. “A novel adaptive beamformer based on semidefinite programming with magnitude response constraints,” IEEE Trans. Antennas and Propagation, vol. 56, no. 5, pp. 1297-1307, Nov. 2008.

IEEE Trans. Signal Processing, vol. 57, no. 7, pp. 2615-2628, Jul. 2009.

CR IP T

[19] Z. L. Yu, W. Ser, M. Er, Z. Gu, and Y. Li,“Robust adaptive beamformers based on worst-case optimization and constraints on magnitude response,”

[20] B. Fuchs, “Application of convex relaxation to array synthesis problems,” IEEE Trans. Antenans and Propagation, vol. 62, no. 2, pp. 634-640, 2014. [21] A. Aubry, V. Carotenuto, A. De Maio, A. Farina, and L. Pallotta,“Optimization theory-based radar waveform design for spectrally dense environments,” IEEE Aerospace and Electronic Systems Magazine, vol. 31, no. 12, pp. 14-25, 2016.

[22] A. Aubry, A. De Maio, M. M. Naghsh, “Optimizing Radar waveform and Doppler filter bank via generalized fractional programming,” IEEE Journal of

AN US

Selected Topics in Signal Processing, vol. 9, no. 8, pp. 1387-1399, 2015.

[23] J. Liang, H. C. So, C. S. Leung, J. Li, and A. Farina, “Waveform design with unit modulus and spectral shape constraints via lagrange programming neural network,” IEEE Journal of Selected Topics in Signal Processing, vol. 9, no. 8, pp. 1377-1386, Dec. 2015. [24] J. Liang, H. C. So, J. Li, and A. Farina, “Unimodular sequence design based on alternating direction method of multipliers,” IEEE Trans. Signal Processing, vol. 64, no. 20, pp. 5367-5381, 2016.

[25] P. Stoica, J. Li, and Y. Xie, “On probing signal design for MIMO radar,” IEEE Transactions on Signal Processing, vol. 55, no. 8, pp. 4151-4161, Aug.

M

2007.

[26] X. Zhang, Z. He, L. Rayman-Bacchus, and J. Yan, “MIMO radar transmit beampattern matching design,” IEEE Trans. Signal Processing, vol. 63, no. 8,

ED

pp. 2049-2056, 2015.

[27] G. Hua, and S. S. Abeysekera,“MIMO radar transmit beampattern design with ripple and transition band control,” IEEE Trans. Signal Processing, vol. 61, no. 11, pp. 2963-2974, 2013.

PT

[28] A. Mutapcic, S. J. Kim, and S. Boyd,“Robust Chebyshev FIR equalization,” IEEE Global Telecommunications Conference, pp. 3074-3079, 2007. [29] K. M. Tsui, S. C. Chan, and K. S. Yeung, “Design of FIR digital filters with prescribed flatness and peak error constraints using second-order cone

CE

programming,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol. 52, no. 9, pp. 601-605, 2005. [30] S. Zhang, “Quadratic maximization and semidefinite relaxation,” Math. Program., vol. 87, no.3, pp. 453-465, 2000. [31] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of

AC

multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1-122, 2011.

[32] D. Gabay,“Applications of the method of multipliers to variational inequalities,” In Augmented Lagrangina Methods: Applications to the solution of Boundary-Value problems, North-Holland: Amsterdam, 1983.

[33] S. Hur, T. Kim, D. J. Love, J. V. Krogmeier, T. A. Thomas and A. Ghosh, “Millimeter wave beamforming for wireless backhaul and access in small cell networks,” IEEE Trans. Communications, vol. 61, no. 10, pp. 4391-4403, 2013. [34] M. Dedeoglu, Y. Kemal Alp, and O. Arlkan, “FIR filter design by convex optimization using directed iterative rank refinement algorithm ,” IEEE Trans. Signal Processing, vol. 64, no. 3, pp. 2209-2219, 2016. [35] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, “Applications of second-order cone programming,” Linear Algebra Appl., vol. 284, no. 1, pp. 193-228, Jan. 1998. [36] O. L. Forst, “An algorithm for linearly constrained adaptive processing,” Proceedings of the IEEE, vol. 60, no. 8, pp. 926-935, Aug. 1972.

ACCEPTED MANUSCRIPT ELSEVIER: SIGNAL PROCESSING

30

[37] C. Chen, and P. P. Vaidyanathan, “Quadratically constrained beamforming robust against direction-of-arrival mismatch,” IEEE Trans. Signal Processing, vol. 55, no. 8, pp. 4139-4150, Aug. 2008. [38] A. H. Sayed, and T. Kailath, “A survey of spectral factorization methods,” Numer. Linear Algebra Appl., vol. 8, no. 6-7, pp. 467-496, 2001. [39] M. Fazel, H. Hindi, and S. Boyd, “Rank minimization and applications in system theory,” Proceedings of the 2004 American Control Conference, vol. 4, pp. 3273-3278, Boston, 2004. [40] A. De Maio, Y. Huang, and M. Piezzo,“A Doppler robust max-min approach to radar code design,” IEEE Trans. Signal Processing, vol. 58, no. 9, pp. 4943-4947, Sep. 2010.

CR IP T

[41] L. N. Tran, M. F. Hanif, and M. Juntti, “A conic quadratic programming approach to physical layer multicasting for large-scale antenna arrays,” IEEE Signal Processing Letters, vol. 21, no. 1, pp. 114-117, Jan. 2014.

[42] S. P. Wu, S. Boyd, and L. Vandenberghe, “FIR filter design via spectral factorization and convex optimization,” Applied and Computational Control, Signals, and Circuits, pp. 215-245, Birkhauser Boston, 1999.

[43] A. Konar and N. D. Sidiropoulos, “Hidden convexity in QCQP with Toeplitz-Hermitian quadratics,” IEEE Signal Processing Letters, vol. 22, no. 10, pp. 1623-1627, 2015.

AN US

[44] K. Huang, Y. C. Eldar, and N. D. Sidiropoulos, “Phase retrieval from 1D Fourier measurements: Convexity, uniqueness, and algorithms,” IEEE Trans. Signal Processing, vol. 64, no. 23, pp. 6105-6117, 2016.

[45] K. Huang and N. D. Sidiropoulos, “Consensus-ADMM for general quadratically constrained quadratic programming,” IEEE Trans. Signal Processing,

AC

CE

PT

ED

M

vol. 64, no. 20, pp. 5297-5310, Oct. 2016.