1317

Optimal Cooperative Jamming to Enhance Physical Layer Security Using Relays Gan Zheng, Li-Chia Choo, and Kai-Kit Wong

Abstract—This correspondence studies cooperative jamming (CJ) to increase the physical layer security of a wiretap fading channel via distributed relays. We first provide the feasible conditions on the positiveness of the secrecy rate and then show that the optimization problem can be solved using a combination of convex optimization and a one-dimensional search. Distributed implementation to realize the CJ solution and extension to deal with per group relays’ power constraints are discussed. Index Terms—Communication system security, cooperative systems, relays, wireless communication.

I. INTRODUCTION Wireless communication networks are widely deployed today, and for many applications it is required that the data transmitted be secure. Technically, the system design goal is to prevent eavesdroppers from decoding any message exchanged by legitimate users. At the physical layer, using an information-theoretic approach, the seminal works by Wyner [1] and shortly afterwards by Csiszär and Körner [2] showed that by using channel codes and signal processing, secure communication is in fact possible without using key encryption in the presence of the eavesdropper. The key feature was that the channel for the legitimate users is better than the eavesdropper’s channel. Therefore, under unfavorable channel conditions (from the point of view of the legitimate users), the secrecy rate would be very low or even zero. To improve the secrecy rate under unfavorable conditions, a possible approach is to use multiple-input multiple-output (MIMO) systems where multiple antennas at the transmitter and receiver allow the use of beamforming to steer a null to the eavesdropper [3], or the generation of artificial noise in the direction of the range space of the eavesdropper’s channel [4]. Another, arguably more flexible, technique is to use cooperation via relaying, as proposed in [5]–[7], where the source to destination transmission is helped by a bank of relays, in the presence of one or more eavesdroppers. In these works, [5] studied the use of decode-and-forward (DF) relays with multiple eavesdroppers. Dong et al. in [6] looked at both the use of DF and amplify-and-forward (AF) relays; in the case of DF with a single eavesdropper and a total relay power constraint, the optimal collaborative relay weights have been obtained in closed form. The most recent work [7] addressed the optimization of the collaborative DF and AF relays with more practical individual relay power constraints and imperfect channel state information (CSI). In providing security, a sensible approach is to utilize the relays for cooperative jamming (CJ). In this case, each relay node transmits a Manuscript received June 30, 2010; revised October 18, 2010, November 08, 2010; accepted November 09, 2010. Date of publication November 15, 2010; date of current version February 09, 2011. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Mathini Sellathurai. G. Zheng is with the Interdisciplinary Centre for Security, Reliability and Trust (SnT), The University of Luxembourg, L-1359 Luxembourg (e-mail: gan. [email protected]). L.-C. Choo and K.-K. Wong are with the Department of Electronic and Electrical Engineering, University College London, WC1E 7JE, U.K. (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org Digital Object Identifier 10.1109/TSP.2010.2092774

Fig. 1. The system model.

weighted jamming signal to degrade the eavesdropper’s signal. The works in [5] and [6] were the first few to investigate the use of CJ. However, their main limitation is that a total relay power constraint was considered and the optimal CJ solution was unknown. This motivates our work which addresses the optimization of collaborative relay weights for CJ in maximizing the secrecy rate with individual relay power constraints. In particular, our contributions are as follows: • we first give conditions under which positive secrecy rate is possible; • given these conditions, we propose an algorithm to obtain the optimal CJ relay beamforming solution using a combination of convex optimization and a one-dimensional search; • the proposed algorithm is also extended to cope with the per group relays’ power constraints; • we further develop a distributed implementation algorithm which permits each individual relay to derive its own weight based on local CSI for achieving a near-optimal secrecy rate. II. SYSTEM MODEL AND PROBLEM FORMULATION Consider a wireless communication system in Fig. 1 with one source S, one destination D; M trusted relays, labeled as fR1 ; . . . ; RM g, and one eavesdropper E. All nodes are assumed to have a single antenna. There is a direct link between S and D, and all relays work synchronously in half-duplex mode. The message from S is uniformly distributed over W = f1; . . . ; 2nR g, for n channel uses and the message rate R0 . The message is mapped to the length-n source codewords xn s, and the codewords are transmitted using n time units in a single transmission slot in a time division system. The source codewords are assumed to be independent zero-mean Gaussian to enable evaluation of the achievable secrecy rate. The channels are assumed to undergo flat fading with CSI perfectly known at S; D, and also E.1 Let hD denote the channel between S and (R) D; hm denote the channel between S and Rm , and hE denotes the channel between S and E. In addition, the channel between Rm to D (D) (E) and that between Rm to E are denoted, respectively, by gm and gm . For ease of exposition, we define the channel vectors (D)

gD

g1

gE

g1

(E)

(D)

; g2

(E)

; g2

(D) T M (E) T M

;...;g ;...;g

(1)

where the superscript ( 1 )T denotes the transpose operation. Both S and fRm g transmit simultaneously to both D and E. The intention of relays is to send jamming signals to interfere with E. To describe this further, let us focus on a source symbol xs (with Efjxs j2 g = PS ) which appears within one time unit of the transmission slot of n 1This assumption is reasonable when the eavesdropper is an active user; and for the high SNR case when the eavesdropper is passive [8].

1053-587X/$26.00 © 2010 IEEE

1318

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 3, MARCH 2011

time units. In the same time unit, the relays transmit the CJ codewords fx(cm) g that are assumed to be independent zero-mean Gaussian sig(m) nals, with Efjxc j2 g = 1. The CJ codeword at Rm is weighted by wm , before sending at the same time unit as S transmitting to D via the di[w1 ; . . . ; wM ]T as the CJ beamforming rect link. We shall define vector. Thus, only a 1-stage transmission protocol is considered here as opposed to the conventional AF and DF relay protocols. The received signals at D and E can be, respectively, written as

w

yD yE

= hD xs + = hE xs +

M

m=1 M m=1

(D) gm wm x(cm) + nD (E) gm wm x(cm) + nE

(3)

(4)

2

D

(5)

E2

(6)

2

0E = jwyPgSEjhj E+j 2

2

2

where ( 1 )y denotes the conjugate transposition. The secrecy rate expression in (4) is recognized as the difference of the capacity at D and the capacity at E. This simple form is obtained under the following assumptions [5]: First, the received signals at D and E at time i only depend on the transmitted codewords at the relays at time i (referred to as the memoryless relay channel assumption); second, the relays use independent, zero-mean Gaussian codewords, similar to the source, to send the jamming signal. Existing results for the MIMO wiretap channel (see [9]) can be applied to (4) under the first assumption. Our aim is to maximize the secrecy rate via the design of the beamforming vector . That is,

w

max w R s:t: jwm j

2

pm 8m

(7)

where pm is the transmit power constraint of Rm and

1 + jwyPgsDjhj D+j 2

D2 :

1 + jwyPgsEjhj E+j 2

(8)

E

2

2

Ps jhD j

Ps jhE j

2

wy gE j + E > jhE j wy gD j + D jhD j

(10)

2

2

or equivalently

j j

2

2

2

2

2

2

(11)

(12)

t s:t: p 8mg w jmax

2 2 tjwy gD j2 + (tD 0 E ) wy gE :

(13)

wg

wg

2 2 It can be easily seen that t = maxw (j y E j2 + E =j y D j2 + D ) 2 2 2 2 > (E =D ) (when = ) or tD 0 E > 0. Hence, (13) belongs to a second-order cone programming (SOCP) [10] which can be solved using a bi-section search over t. Note that a similar technique was used in [7] for AF relay beamforming but [7] did not show that all terms within the square root are nonnegative to enable the formulation (13). In the remaining of this paper, we shall study (7) assuming a positive secrecy rate can be achieved.

w 0

wy gD AND ONE-DIMENSIONAL SEARCH

IV. METHODOLOGY: FIXED

Solving (7) is challenging because it is non-convex and R is a complicated function of . Our approach is to first study a subproblem with j y D j2 fixed and then use one dimension search to find the solution to (7) which is guaranteed to be optimal by its analytical properties. The philosophy behind is similar to that in [11] but the system, problem formulation and techniques used have no overlap as those in [11].

wg

w

wy gD j

A. Subproblem With Fixed j

Consider the subproblem where the interference received at D, i.e., j y D j2 , is fixed to some scalar t 0. As a result, the objective is reduced to maximizing the interference to the eavesdropper E, i.e.,

wg

y max w jw gE j s:t: 2

wg

j y D j2 = t; jwm j2 pm 8m:

(14)

w (t) and the corresponding op-

Suppose that the optimal solution is timal objective value is f0 (t). Define

0

1 + jw Pt jgh jj 1 + Pt jh j R(t) = : (15) 1 + fP tjh j 1 + fP tjh j Therefore, we can focus on finding the maximum of R(t) over t 0. Unfortunately, R(t) is difficult to evaluate due to the equality constraint ( )

+

+

( )+

in (14). To overcome this, we consider the modified problem

wg

j y D j2 t jwm j2 pm 8m:

(16)

Denote the optimal solution to (16) as (t) and the corresponding optimal objective value as f (t). Then, we can define

R1 (t)

1 + Pt jh j 1 + fP tjh j +

(17)

( )+

(9)

2

jhE j2 : jhD j2

w

wy gD j + D > jwy gE j + E ; 9w 2

>

2 t s:t: jwy gE j2 + E2 t(jwy gD j2 + D ) p 8mg w jmax

2

The optimization (7) makes sense only when a positive secrecy rate is possible. Therefore, we find it essential to first derive the conditions for that. In order to obtain a positive secrecy rate, we need

j

2

fj

y max w jw gE j s:t:

III. CONDITIONS FOR POSITIVE SECRECY RATE

2

2

2

( )+

2

R

2

We rewrite the left-hand-side problem in (11) as

fj

where 0D and 0E are the signal-to-noise ratios (SNRs) at D and E, respectively, and are given by

0D = jwyPgSDjhj D+j

jwy gE j + E max fjw j p 8mg jwy gD j + D

which can further be re-expressed as (2)

where nD and nE are the zero-mean Gaussian noises at D and E, respec2 2 and E[jnE j2 ] = E . The jamming signal tively, with E[jnD j2 ] = D introduces interference at both D and E, and it has recently been known that the achievable secrecy rate for this channel can be expressed as

Rs = log2 (1 + 0D ) 0 log2 (1 + 0E )

whose feasibility can be checked by

and have the following result stated in Theorem 1. Theorem 1: R1 (t) and R(t) have the same maximizer and the same maximum function value. Proof: Suppose that the maximizer of R(t) is t3 and the associated beamforming vector solution to (14) is 0 (t3 ). Obviously, 0 (t3 )

w

w

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 3, MARCH 2011

is also a feasible solution to (14), meaning f (t3 ) f0 (t3 ) and as a result, we have

max R1 (t) t

jh j + R1 (t ) = P jh j 1 + f (t )+

jh j + P jh j 1 + f (t )+

1+

3

max

w

wy g

s:t:

E

0

(18)

t

g

jwy g j2 t jw j2 p 8m D

m

(19)

m

wg

where we have used the fact that y E can be made positive (see [10] for a similar example) without loss of optimality and the advantage is that (19) becomes convex and can be optimally solved. Denote the optimal solution to (19) by 3 (t) for a given t and f (t) j 3 (t)y E j2 . We have the following result. Theorem 2: f (t) is a concave function of t. Proof: The Lagrangian of (19) is given by

w

w

w

g

L( ; ; ) =

=

=

0wy g

+ (

0wy g

+

E

E

jwy g j2 0 t) + D

wy g gy D

D

M

D

wy g gy

D

()

+

m=1

m (jwm j2

()

+

gg 0w yg E

y

E E

w0

0p

m)

(20)

M

t +

m pm

w0

t +

m pm m=1

w

=

w (t) 0

g g

D

y

D

m pm m=1

E

E

E

w (t) g

E

M

min t + ;>0

0

g gy

m pm s:t: m=1

D

D

()

+

0 g gy 0: E

E

(25) It is easy to see that the optimal solutions (; ) to (25) are just a scaled version of those to (24) by a factor of 3 (t)y E , so they give the optimal value of (25) since the objective function is linear. Due to the convexity of (16), strong duality holds and the optimal objective value of (24) is (t)3y E , same as (16). Multiplied by another 3y E , the optimal objective value of (25) becomes j 3y E j2 which is exactly f (t) as defined before. It is easily checked that f (t) is a point-wise minimum of a family of affine functions and as a result concave for t 0 [12, p. 80]. This completes the proof.

w

w

g

g

w g

w g

B. Search for the Optimal Solution Even if f (t) is available and it is concave, R1 (t) still appears to be too complicated for optimization. Our next result studies the analytical properties of R1 (t) which will permit us to design an efficient algorithm that finds the optimal solution. The main result is summarized in the following theorem. Theorem 3: R1 (t) is quasi-concave in t and its maximum can be found via a one-dimensional search. Ps jhD j2 ; b Proof: For convenience, we define a 2 2 Ps jhE j2 ; c D ; d E . Then, we have a

R1 (t) =

1 + t+c 1+

b f (t)+d

:

The proof follows the result in [12] that if R0 (t) = 0 implies that R00 (t) < 0 for any t 0, then R1 (t) is quasi-concave in t 0. This result was also used in [11]. The first-order derivative of R1 (t) can be derived as

R10 (t) = 0a(t + c)02 1 + 1+

a t+c

b

01

f (t) + d

bf 0 (t) 1 +

b f (t) + d

02

(f (t) + d)

02

:

E

Setting it to zero, we have

E

:

D

y 0 g3 gy 0:

w

M

t +

()

(26)

g gy w3 (t) () 0 w3 (t)y g

+

D

It should be noted that this is not the usual way of deriving the dual problem but it will prove to be useful later when the property of f (t) is studied. Obviously, (24) does not lead to a solution to the dual problem due to the involvement of the optimal value of primal variable 3 (t) which is unknown at the moment. However, it will result in the same objective value of the original dual problem and therefore that of the primal problem as well. Now, let us consider the following modified dual problem:

+

y

m=1

M

0

3

g gy +

(21)

where and are Lagrange multipliers, and () is a square matrix with the diagonal elements being . The dual objective is G(; ) = minw L( ; ; ) which reaches the maximum at the optimal value of the primal variables 3 (t), i.e.,

G(; )

0

m pm s:t:

m=1

(22)

w

g

M

= max R(t):

w w

g

gg w

min t + ;>0

On the other hand, suppose that the maximizer of R(t) is t1 with the corresponding solution to (16) as (t1 ). Then, we must have j (t1 )y D j2 = t1 , which can be easily seen by the method of contradiction as follows. Assuming that j (t1 )y D j2 t2 < t1 , then f (t1 ) = f (t2 ). By definition of R1 in (17), we get R1 (t1 ) < R1 (t2 ), which contradicts that t1 is the maximizer. The above observation means that at the optimum of R1 (t), the first inequality in (16) becomes an equality, and thus this optimum can be attained by R(t). Combining this and (18) completes the proof. As (14) is not solvable due to the equality constraint, Theorem 1 is very useful in that we can instead solve (16) (see below) to get R1 (t) in order to find the optimal solution of the original problem. The remainder of this section focuses on maximizing R1 (t). To proceed, we first rewrite (16) as

w

g g

y To ensure that G(; ) is lower bounded, we must have D D + y y 3 () 0 ( E E = (t) E ) (positive-semidefiniteness). Then, its dual problem can be expressed as

(24)

P t

1+

P t

1319

(23)

a (t + c)02 1+

b f (t) + d

2 (f + d) =

b 1+

a t+c

f 0 (t): (27)

1320

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 3, MARCH 2011

R1 (t)

Further, the second-order derivative of 3 4 (1 + (b=f + d)) (f + d) gives

R1 (t) 00

b

multiplied by

3

4 (28a) f (t) + d (f (t) + d) 2 b 3 4 = 2a(t + c) 1+ (28b) f (t) + d (f (t) + d) 0 2abf (t)(t + c) 2 1 + f (t)b+ d (f (t) + d)2 (28c) a 2 2 1+ + 2b (28d) t + c (f (t)) 0 2b 1 + t +a c 1 + f (t)b+ d (f (t) + d)(f (t))2 (28e) a b 2 + 1+ t + c b 1 + f (t) + d (f (t) + d) f (t): (28f)

1+

0

0

0

11: Solve problem (16) with t + 1t for very small 1t > 0. 12: Evaluate R10 (t) using the above two solutions and (17). 13: if R10 (t) > 0, tmin = t. 14: 15: else tmax = t. 16: 17: end 18: end 19: end 20: Output: .

w

0

0

00

Combining (28d) and (28e) leads to

02b(f (t))2 0

1+

a t + c (f (t) + d):

02b(f t

2 ( ))

1+

a

t + c (f (t) + d + b):

b2 (t + c) a

1+

a

t+c

2

f (t))2 :

(

0

2 w gE s:t: jw gD j jwmtj2 pl 8l: y

y

max

w

0;>0 t +

L

min

(31)

(33)

m2N

The only difference is that (25) becomes

(30)

Squaring (27), (28b) becomes 2

Suppose there are L groups of relays fN1 ; . . . ; Nl ; . . . ; NL g, which form a partition of f1; 2; . . . ; M g. It is natural to consider per group relays power constraints, which can be viewed as a generalization to individual relay power constraints, and can also reflect the scenario where each group has multiple antennas. In this case, the optimization problem is given by

(29)

Substituting (27) into (28c) and combining it with (29), we then have 0

V. GENERALIZATION TO PER GROUP RELAYS’ POWER CONSTRAINTS

l=1

l pl s:t: gD gD +

0 gE gE 0

y

y

( )

(34)

where (1 ; . . . ; 1 ; 2 ; . . . ; 2 ; . . . ; L ; . . . ; L ). It can be seen that all analysis is still valid, so the proposed algorithm can be easily generalized to handle the per group relays constraints.

On the other hand, 0(30)=(31) gives

a(f (t) + d + b) b(t + c) 1 + t+a c a(f (t) + d + b) = a(1 + f (tb)+d ) > a(1 + t+a c ) = 1 = b(a + t + c) (a + t + c) (a + t + c)

VI. DISTRIBUTED IMPLEMENTATION

(32)

where we have used the assumption that the secrecy rate is positive, i.e., (t + c=a) < (f (t) + d)=(b). This indicates that the sum of (30) and (31), and the sum of first four terms (28b)–(28e), are negative. Also, (28f) is negative since f (t) is concave. Thus, we have proved that R100 < 0, which means that R1 (t) is a quasi-concave function and its maximum can be efficiently found by a one-dimensional search [12, p. 101]. The complete algorithm is summarized in Algorithm 1.

The overall optimization can be performed at either the source S or the destination D, which needs to know all the necessary system 2 2 parameters. In practice, D or S may learn E ; E , and D ; D from is obtained, D the relays fRm g. After the optimal relay weights or S informs each individual relay about its own amplification coefficient, wm , which however will take some bandwidth to do so. It would be much preferred if each individual relay can derive its own beamforming weight based on its local CSI. Here, we shall discuss such distributed implementation algorithm, with the assumption that Rm knows (E) (D) perfectly gm and gm . To facilitate the design, we rewrite the dual problem (25) as

g

0;>0 t+

Algorithm 1: Proposed Algorithm

w

m pm s:t: gD gD +

m=1

g

0gE gE 0:

y

min

2 2 1: Input: Ps ; hD ; hE ; gD ; gE ; D , and E . 2: begin 3: Use the conditions in Section III to check whether a positive secrecy rate is possible. 4: if a positive secrecy rate is not possible, then 5: secrecy rate is zero, exit. 6: end 7: Initialize tmin and tmax . 8: While tmax 0 tmin where is a preset threshold, t = (tmin + tmax )=2. 9: 10: Solve problem (16) with the above t and get solution .

M

w

( )

y

(35)

Due to one of the complementary conditions, we have

gD gD +

0 gE gE )w = 0

y

(

y

( )

(36)

which gives

w = (gE w)gE 0 (gD w)gD : y

( )

y

(37)

g w > 0, it leads to

y Since E

w = gE 0 (gD w) gD :

( )

g w) y

( E

y

g w) y

( E

(38)

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 3, MARCH 2011

1321

Therefore, if the coefficient m > 0, we have

g w gmE 0 g w gmD p wm = pm g w gmE 0 g w gmD ( )

(

(

( )

(

(

)

( )

)

)

(39)

( )

)

and the reason why Rm should use its full transmit power pm is due to another complementary condition of (19), i.e., m ( wm 2 pm ) = 0 and m > 0 implies that wm 2 = pm . From (39), as long as m > 0; Rm can learn its own weight based (E) (D) y = Ey ) on local CSI gm and gm while the common scalar (gD can be broadcast to all relays. In the ideal case where m > 0 m, y only one positive scalar ( D = Ey ) needs to be broadcast to allow that relays themselves can find the optimal weights. While if m = 0, the individual relay gets no information about its beamforming weight from the common information, and in this case, (39) may not be the optimal solution. Our simulation results in Section VII will reveal that the distributed algorithm can achieve secrecy rate very close to the optimal one.

j j

g wgw

j j 0

wg w 8

Fig. 2. Secrecy rate versus source-destination distance.

VII. SIMULATION RESULTS Computer simulations are performed to evaluate the achievable secrecy rate of the proposed algorithm. For the simulations, we assumed a one-dimension system model, and place the source, relays, destination and eavesdropper along a line. Furthermore, the distance between relays is assumed to be short compared to the source-relay, relay-destination, relay-eavesdropper distances. Channels are modeled by a line-ofsight channel model including the path loss, to take into account the effects of distances. Therefore, the channels can be expressed in general as h = d0(c=2) ej , where d is the distance, and c is the path loss exponent chosen as 3.5, and the random phase is uniformly distributed over [0; 2 ). Also, the path loss for nodes to and from the relays can be assumed to be the same since the distance between relays is small. The source-eavesdropper distance is fixed at 50 m. The noise power is –100 dBm while the source power and individual relay power constraints are chosen to be 40 dBm. We first let the source-relays distance fixed at 25 m. We change the position of the destination so that the source-destination distance varies from 10 m to 100 m. Fig. 2 shows the secrecy rate comparison between the proposed algorithms and the direct transmission (without the aid of relays). We use 10 relays in the simulations, i.e., M = 10. The zero-forcing (ZF) solution (labeled as ZF-CJ in the figure) corresponds to forcing the jamming signal at the destination to zero (i.e., t = 0 in (14)). Note that it is different from that in [6] where only total power constraint was considered. We can see that the secrecy rate is significantly increased using our optimal CJ solution, as compared to the secrecy rate using direct-only transmission. We also observe that by using CJ, there is a positive secrecy rate even when the channel for the direct transmission from source to destination is not favorable, when compared to the source to eavesdropper channel. This happens when the source-destination distance is 50 m onwards, corresponding to the eavesdropper masking the direct transmission. Results also illustrate that ZF-CJ can achieve nearly optimal performance in most cases, but when the eavesdropper is close to the destination (e.g., source-destination distance is around 50 m), a significant gap is observed, and this can be explained by the fact that the ZF-CJ solution fails to manage the interference properly at the destination while the optimal solution is able to use a more intelligent strategy to balance the interference to the eavesdropper and the destination. It can also be seen that the secrecy rate achieved by the distributed algorithm is very close to the optimal rate and this justifies the proposed distributed implementation. We then fix the source-destination distance at 60 m. The position of the relays is changed so that the source-relays distance varies from 5 m

Fig. 3. Secrecy rate versus source-relay distance.

to 45 m. Fig. 3 shows a similar secrecy rate comparison as has already been observed in Fig. 1 with regard to the source–relays distance. It is observed that the secrecy rate increases as the relays are closer to the eavesdropper. However, the performance gain is not significant, which can be explained by the fact that as the relays are closer to the eavesdropper, they also generate more interference to the destination. VIII. CONCLUSION This correspondence has studied the CJ via distributed relays to increase the physical layer security. The conditions for positive secrecy rate have been derived and we have shown that the optimal CJ solution can be obtained by a combination of convex optimization and a one-dimensional search. Extension to per group relay power constraints and distributed implementation have also been considered.

REFERENCES [1] A. D. Wyner, “The wire-tap channel,” Bell Syst. Tech. J., vol. 54, no. 8, pp. 1355–1387, 1975. [2] I. Csiszár and J. Körner, “Broadcast channels with confidential messages,” IEEE Trans. Inf. Theory, vol. 24, no. 3, pp. 339–348, May 1978.

1322

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 3, MARCH 2011

[3] A. Khisti and G. W. Wornell, “Secure transmission with multiple antennas: The MISOME wiretap channel,” IEEE Trans. Inf. Theory [Online]. Available: http://allegro.mit.edu/pubs/posted/journal/2007khisti-wornell-it.pdf, to be published [4] S. Goel and R. Negi, “Guaranteeing secrecy using artificial noise,” IEEE Trans. Wireless Commun., vol. 7, no. 6, p. 2180, Jun. 2008. [5] J. Li, A. P. Petropulu, and S. Weber, “Optimal cooperative relaying schemes for improving wireless physical layer security,” Jan. 2010 [Online]. Available: http://arxiv.org/abs/1001.1389, submitted for publication [6] L. Dong, Z. Han, A. Petropulu, and H. V. Poor, “Improving wireless physical layer security via cooperating relays,” IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1875–1888, Mar. 2010. [7] J. Zhang and M. C. Gursoy, “Collaborative relay beamforming for secrecy,” Dec. 2009 [Online]. Available: http://arxiv.org/abs/0910.4132, submitted for publication [8] P. K. Gopala, L. Lai, and H. El Gamal, “On the secrecy capacity of fading channels,” IEEE Trans. Inf. Theory, vol. 54, no. 10, pp. 4687–4698, Oct. 2008. [9] F. Oggier and B. Hassibi, “The secrecy capacity of the MIMO wiretap channel,” in Proc. IEEE Int. Sym. Inf. Theory, Toronto, ON, Canada, Jul. 2008, pp. 524–528. [10] M. Bengtsson and B. Ottersten, , L. C. Godara, Ed., “Optimal and suboptimal transmit beamforming,” in Handbook of Antennas in Wireless Commun.. Boca Raton, FL: CRC Press, Aug. 2001. [11] A. Wiesel, Y. C. Eldar, and A. Beck, “Maximum likelihood estimation in linear models with a Gaussian model matrix,” IEEE Signal Process. Lett., vol. 13, no. 5, pp. 292–295, May 2006. [12] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.: Cambridge Univ. Press, 2004.

Rician Distributed FMRI: Asymptotic Power Analysis and Cramér–Rao Lower Bounds Joonki Noh and Victor Solo, Fellow, IEEE

Abstract—Traditional functional MRI detection statistics for activation and hemodynamic response modeling assume Gaussian data, which is true only for high signal-to-noise ratio (SNR). The correct distribution is Rician. In this correspondence, we provide two new developments in the Rician case. First, a derivation of the theoretical asymptotic power function (as the number of samples goes to infinity). Second, the derivation of a Cramér–Rao lower bound. This allows a correct assessment of the impact of various signal and noise levels on detection power for activation and/or hemodynamic response parameter estimation accuracy. Based on our analysis, we are able to extend existing definitions of SNR by considering variation not only in baseline but also in drifts. Index Terms—Activation detection, fMRI, hemodynamic response modeling, Rician distribution, SNR.

I. INTRODUCTION From the very beginning of functional magnetic resonance imaging (fMRI) data analysis, most attention has focused on analyzing the magnitude voxel time courses under a Gaussian assumption. However, since the magnitude voxel time series are mostly obtained from complex-valued raw MR data, they in fact are Rician distributed. The prevailing Gaussian model has been supported by the claim that a Rician probability density function (PDF) is well approximated by a Gaussian PDF when the signal-to-noise ratio (SNR) is sufficiently high [1]. Under this high SNR assumption, a dynamic linear model (erroneously called the general linear model (GLM) [2]) has been the dominant statistical model used to analyze sequence of functional images collected from the human brain by fMRI [3], [4]. As pointed out by Rowe and Logan [5], however, when the high SNR assumption is in doubt, e.g., in studies with higher spatial resolution, voxels with high spatial dropout and studies with inherently low SNR such as pain studies, then the correct detection of activated regions requires the use of a more accurate Rician model. To model magnitude voxel time courses more accurately and to construct the associated activation statistics that hold for all SNRs, a few methods based on the Rician distribution have been proposed in the fMRI literature. Dekker and Sijbers [6] developed a detection statistic for activation with a Rician distribution through a generalized likelihood ratio test (GLRT). For a simple signal model involving one baseline and one activation amplitude, a nonstandard numerical optimization to estimate parameters under the null and alternative hypotheses

Manuscript received November 14, 2009; revised July 13, 2010, October 25, 2010; accepted November 18, 2010. Date of publication December 10, 2010; date of current version February 09, 2011. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Yufei Huang. The material of this correspondence was presented in part at the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2010. J. Noh is with the Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI48109 USA (e-mail: [email protected]). V. Solo is with the School of Electrical Engineering and Telecommunications, University of New South Wales, Sydney, NSW 2052, Australia (e-mail: [email protected]). Color versions of one or more of the figures in this correspondence are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TSP.2010.2098400

1053-587X/$26.00 © 2010 IEEE