Computation of the Normalized Detection Threshold for the FFT Filter ...

Report 2 Downloads 14 Views
JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

35

Computation of the Normalized Detection Threshold for the FFT Filter Bank-Based Summation CFAR Detector Sichun Wang Defence Research and Development Canada-Ottawa Email: [email protected] Franc ¸ ois Patenaude Communications Research Centre Canada Email: [email protected] Robert Inkol Defence Research and Development Canada-Ottawa Email: [email protected]

Abstract—The FFT lter bank-based summation CFAR

procedures such as the Newton-Ralphson algorithm. As

detector is widely used for the detection of narrowband

practical implementations of these numerical procedures

signals embedded in wideband noise. The simulation and implementation of this detector involves some problems concerning the reliable computation of the normalized detection

require good lower and upper bounds of

T

for reliable

initialization, several theoretical lower and upper bounds

T

for

presents a comprehensive theoretical treatment of major

However, a self-contained comprehensive treatment of

aspects of the numerical computation of the normalized

the FFT lter bank-based summation CFAR detector has

detection threshold for an AWGN channel model. Equations are derived for the probability of false alarm,

Pf a ,

for

both non-overlapped and overlapped input data and then

have also been obtained [8].

not appeared in the literature and details of relevant mathematical derivations have not been published. This paper presents a comprehensive theoretical treat-

used to compute theoretical upper and lower bounds for

T.

for a given

Pf a

threshold for a given probability of false alarm. This paper

A very useful transformation is

ment of major aspects of the numerical computation of

introduced that guarantees the global quadratic convergence

the normalized detection threshold for an AWGN channel

the detection threshold

of the Newton-Ralphson algorithm in the computation of

T

for overlapped data with an overlap ratio not exceeding

50%.

It is shown that if the product of the number of FFT

model. Formulas for computing

Pf a

as a function of

T

are derived from rst principles. Several theoretical lower

T

bins assigned to a channel for signal power estimation and

and upper bounds for

the number of input data blocks is relatively small, e.g.,

transformation is introduced that guarantees the global

less than

60,

the theoretical normalized detection threshold

can be accurately computed without numerical problems. To handle other cases, good approximations are derived. Index Terms—Digital FFT lter bank, Detection and estimation, Constant false alarm rate (CFAR) detection.

are derived and a very useful

quadratic convergence of the Newton-Ralphson algorithm in the computation of the normalized detection threshold for overlapped data. Finally, accurate approximations for

T

are presented. The results presented in this paper are

very useful for the simulation and implementation of the FFT lter bank-based summation CFAR detector.

I. I NTRODUCTION The FFT lter bank-based summation CFAR detector

This paper is organized as follows. Section II introduces the FFT lter bank-based summation CFAR detector.

is an efcient technique for the detection of narrowband

Section III derives the formulas for computing

signals embedded in wideband noise and has important

a given

T

Pf a

for

for overlapped and non-overlapped signal data.

T

applications in civilian spectrum monitoring, electronic

Section IV proves several lower and upper bounds of

warfare, radio astronomy and instrumentation [1]-[7].

for overlapped and non-overlapped input data. Section V

This detector operates by adding estimates of spectral

introduces a very useful transformation that guarantees

power of FFT bins in channels, each of which corresponds

the global quadratic convergence of the Newton-Ralphson

to a group of one or more contiguous FFT bins, and

algorithm while Section VI derives two approximations

comparing the sum of the channel power estimates over

for

multiple input data blocks against a detection threshold,

Throughout this paper, typical numerical examples are

T.

given to illustrate the results.

In its general form, this detector has been the subject

T.

Finally Section VII presents concluding remarks.

of many recent studies [7]-[14]. In particular, closed-form algebraic formulas for computing the probability of false alarm,

Pf a , as a function of T

have been derived [5], [7],

[13]. These results can be used to compute the normalized detection threshold (to be dened later) using numerical

© 2007 ACADEMY PUBLISHER

II. T HE FFT F ILTER BANK -B ASED S UMMATION CFAR D ETECTOR Consider a uniformly sampled band-limited signal embedded in additive white Gaussian noise. Assume that

M

36

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

channels are uniformly distributed across the frequency

detector or simply the summation detector if there is no

range contained within the Nyquist bandwidth and that

risk of confusion. A block diagram is given in Figure 1

K N

for the special case where

FFT bins are assigned to each channel with the FFT bins (N

≤ K)

N

and

K

are equal.

centered within each channel

used to estimate the power contained within the channel.

MK

Consequently, an FFT of length the power levels for the

M

Input data blocks

Detection threshold

is used to compute

channels. For notational

convenience and without loss of generality, assume

2

L stage delay

| |

K−N

w = [w0 , · · · , wM K−1 ]t be a symmetric window of length M K , where the superscript t denotes matrix transposition. Let L ≥ 1 be a positive integer and consider L consecutive overlapping sample vectors, Rl , constructed as follows:

is an even integer. Let

2

| |

S +

S

. . .

Channel 1 detection decision

2

Rl = [rl(1−γ)M K+M K−1 , · · · , rl(1−γ)M K ]t ,

| |

(1)

IFFT

sample of the input data stream and

2

| |

L stage delay

| |2

S

. . .

+

0 ≤ l ≤ L − 1, 0 ≤ γ ≤ 1/2, rn is the n-th γM K is an integer. In practice, the overlap ratio, γ , is often selected to be either 0 or 1/2. These two cases correspond to zero and 50% overlap, respectively. For each l, the two input vectors, Rl and Rl+1 , have γM K samples in common. The vectors, Rl , are windowed by w, resulting in the

where

Channel 2 detection decision

windowed sample vectors:

Xl = [w0 rl(1−γ)M K+M K−1 , · · · , wM K−1 rl(1−γ)M K ]t . The vectors,

Xl ,

Block diagram for the

summation CFAR detector.

F of dimensions M K ×

to yield the FFT lter bank output sample vectors:

III. T HE P ROBABILITY OF FALSE A LARM FOR THE S UMMATION CFAR D ETECTOR

Yl = FXl = [sl,0 , sl,1 , · · · , sl,M K−1 ]t ,

rn

Assume the input data stream

where

complex-valued

F=  1  ···    1   ··· 1

(2)

1 ··· 2πjl e MK ··· e

of length

M

··· 1 ··· ··· 2πj(M K−1)l MK ··· e ··· ··· 2πj(M K−1)(M K−1) MK ··· e

2πj(M K−1) MK

From each vector,

N

L-block

are then transformed by the inverse

discrete Fourier transform matrix

MK

Fig. 1.

Yl ,

a vector



   .  

is formed by summing the power from the

zl,k =

N −1 X

|sl,Ik +m | ,

N

and the implementation details of the detector. For brevity, the FFT lter bank-based

L-block

sum-

mation CFAR detector shall be referred to as the

L-

block summation CFAR detector or the summation CFAR

with

(L−1 X

) zl,k ≥ T

.

(4)

l=0 The following three theorems provide the mathematical

T > 0, Pf a

(3)

Ik , Ik + 1, · · ·, Ik + (N − 1) are summed to form the power estimate for the k -th channel for the data block Rl where Ik = kK + K−N 2 . The detection criterion is dened as follows: if for a given detection threshold, T , PL−1 l=0 zl,k ≥ T , a signal is declared to exist in the k th channel. The choice of T is dependent on the desired probability of false alarm, Pf a , the noise spectral density,

© 2007 ACADEMY PUBLISHER

Pf a = Pr

given

FFT bins with indices

sequence

detector is then given by:

Theorem 1. Assume

0 ≤ l ≤ L − 1, 0 ≤ k ≤ M − 1,

is a zero-mean

noise

foundation for computing the threshold

2

m=0

i.e., the power estimates from the

Gaussian

E(rp rq∗ ) = σ 2 δpq , where σ 2 > 0 is the noise variance (noise oor), δpq = 1 if p = q and δpq = 0 if p 6= q . For a given threshold, T , the probability of false alarm, Pf a , for the k -th channel of the L-block summation CFAR

zl = [zl,0 , · · · , zl,M −1 ]t

FFT bins centered within each channel:

white

Pf a =

LN X

L ≥ 2

and

T

for a given

0 < γ ≤ 1/2.

Pf a .

For a

is given by:

Y

−1 λLN l

l=1

(λl − λm )



e

T σ 2 λl

,

(5)

1≤m≤LN,m6=l

λl , 1 ≤ l ≤ LN , are the LN distinct positive L × L block matrix H:   A B 0 ··· 0 0  BH A B 0 ··· 0    H  0 B A B ··· 0   H= (6)  ··· ··· ··· ··· ··· ··· .    0 ··· 0 BH A B  H 0 ··· 0 ··· B A

where

eigenvalues of the

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

In (6),

H

LN × LN , 0 is the N × N N × N matrices dened

is of dimensions

zero matrix and

A

and

B

are

respectively by:



τ11  ···  A=  τp1  ··· τN 1 τpq

τ12 ··· τp2 ··· τN 2

· · · τ1q ··· ··· · · · τpq ··· ··· · · · τN q

MX K−1

=

wl2

l=0



γ11 ··· γp1 ··· γN 1

  B=  

γ12 ··· γp2 ··· γN 2

37

··· ··· ··· ··· ···

· · · τ1N ··· ··· · · · τpN ··· ··· · · · τN N

In practical systems, the noise oor,

   ,  

For a given

(7)

l=0

L ≥ 1, N > 1 ³

Amp

m=1 p=1

                      

p−1 X

(8)

T /σ 2 ,

which

is computed using one of the equations (5), (11) or (13).

threshold in this paper. This can be loosely interpreted as single data block.

  ,  

IV. L OWER AND U PPER B OUNDS FOR THE

(9)

D ETECTION T HRESHOLD In practice, the computation of

T

T /σ 2

involves the

numerical solution of one of equations (5), (11) or (13) algorithm. To avoid potential numerical problems and

(10)

and

γ = 0.

For a

ensure fast convergence, a good lower or upper bound for

T

is required to initialize the numerical procedures. This



T σ 2 µm

,

(11) A. Overlapped Input Data

Amp , 1 ≤ m ≤ N , 1 ≤ p ≤ L,

Y

The cases for overlapped and non-overlapped signal

the appendix.

e

s!

T.

data are treated separately, with full proofs provided in

are

L

 AmL =  

−1 µN m

(µm − µl )

  , 

1≤l≤N,l6=m

L ≥ 2, 0 < γ ≤ 1/2, and let λ1 > λ2 > · · · > λLN be the LN distinct positive eigenvalues of the LN × LN Hermitian matrix H dened by (6). For a given threshold T > 0, let the probability of false 2 alarm be denoted by Pf a . Let T1 (Pf a ) = −λ1 σ ln Pf a and denote the solution of the following equation in x by Tn (Pf a ): Theorem 4. Assume

Amp = AmL × X

Pf a =

(12)

1≤l≤N,l6=m

µl − µm

(kl )!(L − 1)!

γ = 0.

Pf a =

¡ T ¢s λ

s!

s=0

T

e− λ ,

MX K−1 l=0

© 2007 ACADEMY PUBLISHER

,

(15)

For a given

n ≥ 1, let the unique solution of z be denoted by Bn (Pf a ): ¶ µ 2 z n−1 z −z = Pf a . (17) + ··· + e 1+z+ (n − 1)! 2!

Then for

0 < Pf a ≤ 2/e,

(13)

and for

(14)

(16)

the following equation in

Bn (Pf a ) ≤

wl2 .

x σ 2 λm

For any positive integer

where

λ = σ2



T1 (Pf a ) ≤ T2 (Pf a ) ≤ · · · ≤ Tl (Pf a ) ≤ Tl+1 (Pf a ) ≤ · · · ≤ TLN (Pf a ) = T.

.

is given by:

L−1 X

(λm − λl )

e

Then

kq , 1 ≤ q ≤ N , are non-negative integers and µq , 1 ≤ q ≤ N , are the N distinct positive eigenvalues of the Hermitian matrix A dened by (7). and

λn−1 m

2 ≤ n ≤ LN.

Note that

L ≥ 1, N = 1

Y

1≤l≤n,l6=m

Γm (k1 , · · · , km−1 , km+1 , · · · , kN ), Γm (k1 , · · · , km−1 , km+1 , · · · , kN ) = ¸kl · Y µl (L + kl − 1)!

n X m=1

k1 +···+km−1 +km+1 +···+kN =L−p

T > 0, Pf a

is obtained by

by

using numerical procedures such as the Newton-Ralphson



Theorem 3. Let

T

the detection threshold for one FFT bin computed from a

´s

T

σ 2 µm

s=0

where the coefcients

                       

σ2

T /(LN σ 2 ), which will be called the normalized detection



· · · γ1N ··· ··· · · · γpN ··· ··· · · · γN N

is given by:

N X L X

dened by:

the detection threshold

section presents theoretical lower and upper bounds for

T > 0, Pf a

Pf a =

Pf a ,

multiplying the estimated noise oor

found that it is more appropriate to focus on the quantity,

2πjl(p − q) . wl wl+(1−γ)M K exp MK

Theorem 2. Assume given

needs to be

In the analysis of the summation CFAR detector, it is

γpq = e−2πj(1−γ)(Ik +q−1) × γM K−1 X

σ2 ,

estimated from the output samples of the FFT lter bank.

2πjl(p − q) , exp MK

γ1q ··· γpq ··· γN q

The proofs of these theorems are relatively lengthy and will be given in the appendix.

n+



¸ · 2 2πn , ln ePf a 2

(18)

0 < Pf a < 1, T



λ1 σ 2 BLN (Pf a ).

(19)

38

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

L ≥ 1, N ≥ 2, γ = 0 and let µ1 > N distinct positive eigenvalues of A dened by (7). We have

Theorem 5. Assume

µ 2 > · · · > µN

be the

the Hermitian matrix

2

2

µ1 σ BL (Pf a ) ≤ T ≤ µ1 σ BLN (Pf a ). T1

Moreover, if

T2

and

g 00 (p) < 0

Blackman Window (N=L=5, MK=1024) 1

x: 0.9 0.8

1

0.7

(21)

0.6 Pfa

−1 µN − x m Q e σ2 µm 1≤l≤N,l6=m (µm − µl )

m=1

for

(20)

= (Pf a ) L , N X

g(p) is strictly increasing and (0, 1). Specically, g 0 (p) > 0 and 0 < p < 1.

Theorem 6. The function concave on the interval

−1 µN − x m Q e σ2 µm 1≤l≤N,l6=m (µm − µl )

m=1

is rather

non-trivial and will be given in the appendix:

are respectively the solutions of

the equations (21) and (22) in

N X

g(p)

computed. The proof of the concavity of

B. Non-Overlapped Input Data

0.5

1

= 1 − (1 − Pf a ) L ,

(22)

0.4

then

0.3

with

LT1 ≤ T ≤ LT2 ,

(23)

´ ³ 1 T2 ≤ µ1 σ 2 BN 1 − (1 − Pf a ) L .

(24)

0.2 0.1 0

V.

0

0.002

0.004

0.006

0.008

0.01

p

N UMERICAL C OMPUTATION OF THE N ORMALIZED

D ETECTION T HRESHOLD FOR THE S UMMATION CFAR D ETECTOR

Fig. 2.

Using (5), (11) and (13) and in conjuction with the

e−T /(σ

2

Pf a , as a function of p = N = L = 5, M K = 1024, γ = 0.5, Blackman window.

The probability of false alarm,

λ1 ) ,

lower and upper bounds derived in the previous section, the normalized detection threshold

T /(LN σ 2 )

can be

computed numerically. It has been observed that in gen-

Blackman Window (N=L=5, MK=1024) 1

eral, equation (5) is numerically more unstable and hence more difcult to apply in practice than equation (11). To

0.9

avoid numerical problems with (5), the transformation

function

has proved to be very useful. Dene the

g(p)

0.8

by:

g(p) =

LN X

Y

−1 λLN l

l=1

(λl − λm )

0.7

p

λ1 λl

,

0.6

1≤m≤LN,m6=l

0 < p < 1.

(25)

0.5 0.4

Then the formula (5) can be rewritten as:

µ −T ¶ g e σ2 λ1 = Pf a .

Pfa

p = e

−T σ 2 λ1

0.3 (26)

0.2

T can be obtained by rst solving the equation g(p) = Pf a for p and then using the 2 relationship T = λ1 σ ln(1/p). It can be shown that g(p) is concave on the interval (0, 1), as is demonstrated by Figure 2. For comparison, Pf a is plotted as a function 2 of T /σ using (5) in Figure 3. The concavity of g(p) The detection threshold

is a very desirable property to have in the numerical computation of

T /(LN σ 2 );

it guarantees the fast global

0.1 0

0

20

40

60

80

100

T/σ2 Fig. 3.

The probability of false alarm,

N = L = 5, M K = 1024, γ = 0.5,

Pf a ,

as a function of

T /σ 2 ,

Blackman window.

quadratic convergence of the Newton-Ralphson algorithm provided that the functions

g(p) and g 0 (p) can be reliably

© 2007 ACADEMY PUBLISHER

For

LN

not too large, we have successfully computed

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

T /(LN σ 2 ),

the normalized detection threshold,

L-block

39

for the

3

summation CFAR detector. Typical results are

γ = 0.

plotted in Figures 4-5 for

Note that the windows

used to produce these results have been normalized in the

PM K−1

2.8 Hanning Blackman Blackman−Harris Hamming Nuttall Flat−Top Rectangular

Normalized Threshold T/(LNσ2)

2.6 2.4 2.2

2.6 Normalized Threshold T/(LNσ2)

wl2 = 1. The results for γ = 1/2 are l=0 very similar, but are not shown due to space limitations. sense that

Hanning Blackman Blackman−Harris Hamming Nuttall Flat−Top Rectangular

2.8

2.4 2.2 2 1.8 1.6 1.4

2 1.2

1.8 1 −10 10

1.6

−8

−6

10

−4

10

−2

10

10

0

10

PFA

1.4

Fig. 5.

The normalized detection threshold,

of the probability of false alarm

L = 10, γ = 0.

1.2 1 −10 10

−8

−6

10

−4

10

−2

10

0

10

10

with

T /(LN σ 2 ), as a function N = 8, M K = 1024,

The window functions are normalized.

This implies that the normalized detection threshold,

T /(LN σ 2 ),

PFA

Pf a ,

can be approximately computed via a chi-

square distribution. Using this approach, the following Fig. 4.

T /(LN σ 2 ), as a function Pf a , N = 5, M K = 1024, L = 18,

The normalized detection threshold,

of the probability of false alarm

γ = 0.

The window functions are normalized.

two theorems have been derived for normalized windows

PM K−1

wl2 = 1). l=0 Theorem 7. Assume

(i.e.,

Pf a , T VI. ACCURATE A PPROXIMATIONS TO THE D ETECTION T HRESHOLD For very large

LN ,

there are serious numerical dif-

culties in computing the normalized detection threshold

T /(LN σ 2 )

using (5) or (11). In such cases accurate

approximations to

T /(LN σ 2 )

are necessary. As will be

demonstrated in the appendix in the proofs of Theorems 1-3, the detection decision statistic,

PL−1 l=0

zl,k ,

of the

k-

th channel is distributed as a weighted sum of chi-square random variables. More specically,

L−1 X

  

zl,k ∼

l=0

PLN

  PN

m=1

0 < γ ≤ 1/2,

σ 2 µm 2 2 χ2L (m),

γ = 0, (27)

where

λl , 1 ≤ l ≤ LN ,

are the

eigenvalues of the Hermitian matrix

µm , 1 ≤ m ≤ N ,

are the

of the Hermitian matrix

N A

LN distinct positive H dened by (6) and

distinct positive eigenvalues dened by (7). As pointed

out in [15], one can approximate a weighted sum of chi-square random variables by a random variable of the form

c + dχ2h

where

c

and

d

are constants and

chi-square random variable with

© 2007 ACADEMY PUBLISHER

h

and

γ = 0.

For a given

T ≈

³P

N m=1

  σ 2 L N − PN

µ2m

´2 

3 m=1 µm  ³ ´3    PN 2   m=1 µm   h = L ³P ´2  , N 3 m=1 µm

 +

PN



µ3m Bh (Pf a ) Pm=1 , N 2 m=1 µm

(28)

µm , 1 ≤ m ≤ N , are the N distinct positive A dened by (7) and bxc denotes the greatest integer less than or equal to x. Theorem 8. Assume L ≥ 2 and 0 < γ ≤ 1/2. For a given Pf a , T is approximately computed by:   ´2 ³P LN 2 PLN 3 λ l=1 l λl   , + Bh (Pf a ) Pl=1 T ≈ σ 2 LN − PLN LN 2  3 l=1 λl l=1 λl  ³  PLN 2 ´3    l=1 λl   (29) h =  ³P ´2  , LN 3 λ l=1 l

where

eigenvalues of

σ 2 λl 2 2 χ2 (l),

l=1

L ≥ 2

is approximately computed by:

χh

is a

degrees of freedom.

where

λl , 1 ≤ l ≤ LN , are the LN H dened by (6).

eigenvalues of

distinct positive

40

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

The excellent accuracy of approximations (28) and

A PPENDIX

(29) is apparent in Figure 6. This gure shows the true normalized detection thresholds for the cases of

γ = 1/2

γ = 0 and

and the corresponding approximate normalized

detection thresholds plotted as a function of

Pf a

for a

K−N 2 . For

typical example. In this gure's legend, “NO” and “O” denote non-overlapped (γ

= 0) and overlapped (γ = 1/2),

zl,k

N −1 X

=

|sl,Ik +m |2 =

respectively. It follows that, the approximations (28 ) and

m=0

(29) can be used to avoid numerical problems in the direct

= ZH l,k Zl,k ,

computation of the normalized detection threshold.

Blackman Window

2

Normalized Threshold T/(LNσ )

12

10

8

N −1 X

and let

sl,Ik +m (sl,Ik +m )∗

m=0 (30)

= [sl,Ik , sl,Ik +1 , · · · , sl,Ik +N −1 ]t = Fk WRl ,

Zl,k

14 True Threshold (NO) True Threshold (O) Upper Bound: Eq.(19) Upper Bound: Eq. (23) Lower Bound: Eq. (20) Lower Bound: Eq. (23) Approximation: Eq. (28) Approximation: Eq. (29)

0 ≤ γ ≤ 1/2 0 ≤ l ≤ L − 1, we have

Proof of Theorems 1-3. Assume

Ik = kK +

(31)

Rl is the l-th sample vector dened by (1), M K × M K diagonal matrix with its diagonal dened by w:

where in (31),

W

is the

elements



w0 0  0 w1 W=  ··· ··· 0 ···

6

 ··· ··· 0  ··· ··· 0 ,  ··· ··· ··· · · · · · · wM K−1

(32)

Fk is the N × M K matrix consisting of N rows of M K × M K inverse discrete Fourier transform matrix F with row indices Ik , Ik +1, · · ·, Ik +m, · · ·, Ik +N −1:

and

4

the

2

0 −10 10

−8

−6

10

−4

10

−2

10

Fk =  1  ···   1   ··· 1

0

10

10

PFA Fig. 6.

Upper and lower bounds and true and approximate normalized

detection thresholds as a function of the probability of false alarm

N = L = 5, M K = 1024, γ = 0

or

γ = 1/2,

Pf a ,

··· αIk +N −1

L−1 X



  ,  

The problem of computing the normalized detection detector has been comprehensively treated. Formulas for computing the probability of false alarm,

Pf a ,

as a

Hence

zl,k =

l=0

threshold for the FFT lter bank-based summation CFAR where

Z

L−1 X

H ZH l,k Zl,k = Z Z,

(34)

l=0

is the length-LN zero mean Gaussian random

vector dened by:

were derived from

£ ¤t Z = Zt0,k , Zt1,k , · · · , ZtL−1,k .

rst principles. In addition, accurate approximations for

T

Ik +m

2πj α = exp M K.

VII. C ONCLUSIONS

T,

α

··· α(M K−1)Ik ··· ··· (M K−1)(Ik +m) ··· α ··· ··· · · · α(M K−1)(Ik +N −1)

Blackman window.

where

function of the detection threshold,

(33)

α Ik ···

(35)

that avoid potential numerical problems are presented.

These results are very useful for the simulation and implementation of the FFT lter bank-based summation CFAR detector.

From the identity (4a) of [16], the characteristic function of the quadratic form

φ(t),

are

grateful

to

the

reviewer

for

φ(t) = very

l=0

zl,k = ZH Z,

1 , det |I − jtE (ZZH )|

(36)

helpful

comments and suggestions which greatly improved the

where

readability of this paper. We would also like to thank

LN .

I

denotes the identity matrix of dimensions

LN ×

We shall derive the equations (5), (11) and (13)

φ(t)¡. First,¢it is necessary E ZZH explicitly. In

our colleague Dr. Sreeraman Rajan at DRDC Ottawa for

using the characteristic function

carefully reviewing an earlier version of this paper.

to compute the Hermitian matrix

© 2007 ACADEMY PUBLISHER

denoted by

can be shown to be given by:

VIII. ACKNOWLEDGEMENTS We

PL−1

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

fact,

¡ ¢ E ZZH

   ¡ H¢  E ZZ = E    

       

¡ ¢ E ZZH

is given by:





41

³ ´ E Z0,k ZH 0,k · · · ³ ´

Z0,k Z1,k ··· Zp,k ··· ZL−1,k ··· ···

E Zp,k ZH ··· 0,k · · · ³ ´ ··· H E ZL−1,k Z0,k ···

       

Z0,k Z1,k ··· Zp,k ··· ZL−1,k

H        

   =   

³ ´ E Z0,k ZH L−1,k · · · ³ ´

E Zp,k ZH L−1,k · · · ³

E ZL−1,k ZH L−1,k



    .   ´  (37)

For any

l, m, 0 ≤ l ≤ m ≤ L − 1, ¡ ¢ E Zl,k ZH m,k ³ ´ H = E (Fk WRl ) (Fk WRm ) ¡ ¢ H H = E Fk WRl RH m W Fk ¡ ¢ H = Fk WE Rl RH m WFk ,

simplies to the following tridiagonal

(38)

¡ ¢ E ZZH =  A0 A1 0 ···  AH A A 0 0 1  1  ··· ··· ··· ···   0 ··· 0 AH 1 0 ··· 0 ··· ³ ´ H where A0 = E Z0,k Z0,k and A1 of dimensions N × N .

where



δJ−M K+1  ···  δJ−M K+n  ,  ··· δJ

J = (m − l)(1 − γ)M K . δq = 0, if q¡ 6= 0 and ¢ δq = 1 if q = 0. H It is clear that the matrix E Rl Rm is shift invariant in the sense that it depends only on the value of m − l. Let ¡ ¢ Ap = E Zl,k ZH (40) m,k ¡ ¢ H H = Fk WE Rl Rm WFk , p = m − l, m ≥ l. ¡ H¢ The covariance matrix E ZZ can then be put in the with

Here in (39),

following block Toeplitz form:

¢ ¡ E ZZH =  A0 A1  AH A0 1  H  A2 AH 1   ··· ···   ··· ···  H  A ··· L−3  H  AL−2 · · · AH AH L−2 L−1

(41)

A2 A1 A0 ··· ··· AH 1 AH 2 ···

··· A2 A1 ··· ··· A0 AH 1 AH 2

AL−2 ··· ··· ··· ··· A1 A0 AH 1

AL−1 AL−2 AL−3 ··· ··· A2 A1 A0



     .     

m−l ≥ 2, then |(m−l)(1−γ)M K| ≥ 2(1−γ)M ¡ K ≥¢ 2(1 − 1/2)M K = M K and hence entries in E Rl RH m are all equal to zero. Thus Ap = 0 if p ≥ 2 and H =

0 0 ··· 0   ··· ···  , A0 A1  AH A0 1 ³ ´ = E Z0,k ZH 1,k

= σ 2 Fk W2 FH k . Ik = kK + 2πj M K , we obtain:

With

exp

are

K−N 2 ,

(44)

J = MK

and

α =

2 A0 = σ 2 Fk W2 FH k =σ ×  2 (M K−1)Ik · · · wM w12 αIk w02 K−1 α  ··· ··· ··· ···   2 Ik +l   w2 ··· ··· w1 α   0  ··· ··· ··· ··· 2 2 Ik +N −1 ··· ··· w0 w1 α   1 ··· 1   ··· ··· ···   −lIk −l(Ik +N −1)  α · · · α ×     ··· ··· ··· α−(M K−1)Ik · · · α−(M K−1)(Ik +N −1)   τ11 τ12 · · · τ1q · · · τ1N  ··· ··· ··· ··· ··· ···     = σ2   τp1 τp2 · · · τpq · · · τpN   ··· ··· ··· ··· ··· ···  τN 1 τN 2 · · · τN q · · · τN N



= σ 2 A, where, for

τpq

(45)

1 ≤ p ≤ q ≤ N, = =

If

© 2007 ACADEMY PUBLISHER

(42)

¡ ¢ A0 = E Z0,k ZH 0,k ¡ ¢ H = Fk WE R0 RH 0 WFk

(39)

··· ··· ··· ··· ···



We now consider the two separate cases γ = 0 and 0 < γ ≤ 1/2. ¡ H¢ 1) γ = 0. In this case, A1 = 0 and E ZZ becomes the following diagonal L × L block matrix ¢ ¡ E ZZH = (43)   A0 0 0 ··· 0 0  0 A0 0 0 · · · 0     ··· ··· ··· ··· ··· ··· ,    0 ··· 0 0 A0 0  0 · · · 0 · · · 0 A0

where it can be veried that

¡ ¢ 2 E Rl RH m =σ ×  δJ ··· δJ−m  ··· ··· ···   δJ+n−1 ··· δJ+n−1−m   ··· ··· ··· δJ+M K−1 · · · δJ+M K−1−m

L×L

block matrix:

MX K−1 l=0 MX K−1

wl2 αl(p−q) wl2 exp

l=0 and

A

is dened by (7).

2πjl(p − q) , MK

(46)

42

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

2)

0 < γ ≤ 1/2.

In this case

and it remains to compute

A0 is computed A1 . In fact,

¡ ¢ E R0 RH = 1 · 0γM K×(1−γ)M K σ2 0(1−γ)M K×(1−γ)M K

by (45)

and

F1k =  JI α k ···   αJ(Ik +p−1)  ··· αJ(Ik +N −1)

(47)

¸

IγM K×γM K 0(1−γ)M K×γM K

,

0γM K×(1−γ)M K , 0(1−γ)M K×(1−γ)M K and 0(1−γ)M K×γM K are zero matrices of dimensions γM K ×(1−γ)M K, (1−γ)M K ×(1−γ)M K and (1 − γ)M K × γM K respectively and IγM K×γM K is the identity matrix of dimensions γM K ×γM K .

where

Let



w0  0 W1 =   ··· 0  w0  0  W2 =  ··· 0

0 w1 ··· 0

··· 0  0 ··· ,  ··· ··· · · · wγM K−1  ··· 0  0 ··· ,  ··· ··· · · · w(1−γ)M K−1

0 w1 ··· 0

wγM K  0   ··· 0

wγM K+1 ··· 0

W4 = w(1−γ)M K  0   ··· 0

0 w(1−γ)M K+1 ··· 0

¡

··· ··· ··· ···

0 0 ···

where

(49)

·

0 W3

γpq =

  ,  

γM K−1 X

αl(Ik +p−1) ×

(57)

l=0

wl wl+(1−γ)M K α−(l+(1−γ)M K)(Ik +q−1)

 . 

(52)

¸

=·σ 2 Fk × 0γM K×(1−γ)M K 0(1−γ)M K×(1−γ)M K

= e−2πj(1−γ)(Ik +q−1) × γM K−1 X 2πjl(p − q) . wl wl+(1−γ)M K exp MK

W1 W4 0(1−γ)M K×γM K

Summarizing the preceding calculations, we see that

¡ ¢ E ZZH = σ 2 H where H is dened by (6) with B = 0 if γ = 0. Theorems 1-3 can now be proved using the following identity:

1 2π =

¸

Z

+∞

T p−1 X s=0

where

×FH k .

·Z

¡ T ¢s λ

s!

T > 0, λ > 0

+∞

−∞

¸ e−2πjtx dt dx (1 − jtλ)p

T

e− λ ,

and

p

(58)

is a positive integer. This

identity can be proved using an elaborate but standard

2πj α = exp( M K ), J = (1 − γ)M K

F0k =  1 α Ik ··· ···   1 αIk +p−1  ··· ··· 1 αIk +N −1

γ1N ··· γpN ··· γN N

(56)

with



= σ Fk × ¸ · 0γM K×(1−γ)M K IγM K×γM K 0(1−γ)M K×(1−γ)M K 0(1−γ)M K×γM K · ¸ W2 0 × FH k 0 W4

Let

··· ··· ··· ··· ···



l=0

¢

W1 0

¡ 1 ¢H Fk · · · γ1q ··· ··· · · · γpq ··· ··· · · · γN q

C = F0k W1 W4  γ11 γ12  ··· ···  γ γp2 = p1   ··· ··· γN 1 γN 2

wM K−1

A1 = E Z0,k ZH 1,k ¡ ¢ H = Fk WE R0 RH 1 WFk 2

(55)

= σ 2 C,

(51)



α  ···  (M K−1)(Ik +p−1)  . α   ··· (M K−1)(Ik +N −1) α

¡ ¢ = E Z0,k ZH 1,k · ¸ 0 W1 W4 = σ 2 Fk FH k 0 0 ¡ ¢ H = σ 2 F0k W1 W4 F1k

(48)

 ··· 0 0  0 ··· 0 ,  ··· ··· ··· · · · 0 wM K−1

0

··· ··· ··· ··· ···

(54)

can be further simplied to yield:

A1

(50)



A1



W3 =

We obtain

Then



(M K−1)Ik

contour integral argument from the theory of functions of

and set

one complex variable. Details are omitted due to space

(53)

(γM K−1)Ik



··· α  ··· ···  (γM K−1)(Ik +p−1)  , ··· α   ··· ··· (γM K−1)(Ik +N −1) ··· α

© 2007 ACADEMY PUBLISHER

constraints.

0 < γ ≤ 1/2. H dened by (6) denoted by λ1 >

We rst prove Theorem 1. Assume Let the distinct positive eigenvalues of be arranged in decreasing order and

λ2 > · · · > λLN > 0. From (36) it follows that the PL−1 H characteristic function φ(t) of l=0 zl,k = Z Z is given

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

43

by:

φ(t)

1 (59) = det |I − jtE (ZZH )| 1 1 = QLN = 2 det |I − jtσ 2 H| l=1 (1 − jtσ λl ) LN X

=

Y

λlLN −1

1 . 2λ 1 − jtσ l (λl − λm )

l=1 1≤m≤LN,m6=l This implies that

PL−1 l=0

zl,k = ZH Z

=

LN X σ 2 λl

zl,k = ZH Z ∼

l=0 where

2

l=1

χ22 (l), 1 ≤ l ≤ LN, ∼

m=1

Z

+∞

= =

T LN X

1 2π

Z

L−1 X

χ22 (l),

(60)

Y

T

=

2L

of

degrees of freedom. For a given

given by:

Pf a =

) zl,k > T

Z

is

+∞

=

p(x)dx

(65)

T

l=0

¸ Z +∞ 1 −2πjtx φ(t)e dt dx = 2π −∞ T ¸ ·Z +∞ Z N X L X e−2πjtx Amp +∞ dt dx = 2 p 2π T −∞ (1 − jtσ µm ) m=1 p=1 ´s ³ T p−1 N X L X X 2 σ µm − T e σ2 µm . = Amp s! m=1 p=1 s=0

Z

(61)

+∞

·

This completes the proof of Theorem 2.

¸

−2πjtx

(L−1 X

T > 0, Pf a

It remains to prove Theorem the covariance matrix

dt dx

¡ ¢ E ZZH

3. Since

N = 1, L×L

reduces to a

diagonal matrix with the diagonal elements all equal to

−1 λLN l

σ2

(λl − λm )

−∞

Y

are independent and

−∞

1≤m≤LN,m6=l Z +∞ ·Z +∞ T

(64)

identically distributed chi-square random variables each

φ(t)e−2πjtx dt, x > 0,

l=1

1 2π LN X

is distributed

N X σ 2 µm 2 χ2L (m), 2 m=1

χ22L (m), 1 ≤ m ≤ N,

where

means having identical prob-

φ(t)e

zl,k = ZH Z ∼

l=0

are independent and identi-

+∞

zl,k

variables:

T > 0, Pf a is then computed by: (L−1 ) Z +∞ X = zl,k > T = p(x)dx l=0

l=0

as the following weighted sum of chi-square random

−∞

Z +∞ ·

PL−1

the detection decision statistic

and for a given

Pf a

1 (1−jtσ 2 µm )L

are omitted due to space constraints. It then follows that

σ 2 λl 2 2 χ2 (l) is given by 1 . Let p(x) denote the probability density function 2 1−jtσ λl PL−1 H of the detection decision statistic l=0 zl,k = Z Z. We

1 2π

are dened by

QN

decomposition of the rational function

that the characteristic function of

p(x) =

(63)

(12). The routine but tedious derivations of the fractional

ability density functions. Note that here we used the fact

have

1 , (1 − jtσ 2 µm )p

Amp , 1 ≤ m ≤ N , 1 ≤ p ≤ L,

where

is distributed as the

cally distributed chi-square random variables, each of two degrees of freedom, and

Amp

m=1 p=1

following weighted sum of chi-square random variables:

L−1 X

N X L X

e 1 − jtσ 2 λl

−1 λLN l

l=1

2 from (36) that the characteristic l=0 wl . It follows PL−1 H function φ(t) of l=0 zl,k = Z Z is given by:

×

−2πjtx

(λl − λm )

e



φ(t) =

¸ dt dx

T σ 2 λl

,

= (62)

p = 1

was used. This

L−1 X

completes the proof of Theorem 1. To prove Theorem 2, let the distinct positive eigenvalues of the positive denite Hermitian matrix

A

dened

by (7) be arranged in decreasing order and denoted by

µ1 > µ2 > · · · > µN . It follows PL−1 from (36) that the φ(t) of l=0 zl,k = ZH Z is given

characteristic function by:

1 det |I − jtE (ZZH )| 1 = det |I − jtσ 2 H| 1 = QN 2 L m=1 (1 − jtσ µm )

φ(t) =

© 2007 ACADEMY PUBLISHER

1 det |I − jtE (ZZH )| 1 ³ ´L . P MK 1 − jtσ 2 l=0 wl2

(66)

PL−1

H l=0 zl,k = Z Z is distributed as a weighted chi-square random variable:

This implies that

1≤m≤LN,m6=l where the formula (58) with

PM K

zl,k = ZH Z ∼

σ2

PM K−1 l=0

2

l=0

χ22L

where

χ22L ,

(67)

is a chi-square random variable with

degrees of freedom. For a given computed by:

Pf a =

wl2

(L−1 X l=0

) zl,k > T

Z

T > 0, Pf a

2L

is then

+∞

=

p(x)dx T

¸ Z +∞ 1 −2πjtx φ(t)e dt dx = 2π −∞ T # " Z +∞ Z +∞ e−2πjtx 1 = PM K−1 2 L dt dx 2π T wl ) −∞ (1 − jtσ 2 l=0

Z

+∞

·

44

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

¶s

µ =

L−1 X

σ

PMTK−1 2 l=0



wl2

e

s!

s=0

σ2

PMTK−1 l=0

The function w2 l

.

(68)

ψn0 (x) = −

This completes the proof of Theorem 3. Proof of Theorem 4. For

fn (x)

by:

fn (x) = Pr

( n X σ 2 λl 2

l=1

n

where

x > 0,

dene the function and

) χ22 (l) > x ,

is an integer in the range

1 ≤ n ≤ LN .

fn+1 (x) = Pr ≥ Pr

( n X σ 2 λl 2

l=1 Also we have

f1 (x) = e

fn (x) =

n X

χ22 (l)

x σ 2 λ1



2

l=1

Y

) χ22 (l)

>x

)

>x

= fn (x).

λn−1 m (λm − λl )

e



x σ 2 λm

(71)

Let

µ0

be the solution in

gn (x) as follows: ¶ µ xn x2 . + ··· + gn (x) = e−x 1 + x + n! 2!

(72)

xn e−x , gn0 (x) = − n!

(73)

and

xn−1 (x − n)e−x . n!

(74)

gn is a decreasing function on the interval (0, ∞), gn0 (x) < 0 for x > 0. gn is concave on (0, n) 00 and convex on (n, ∞) since gn (x) < 0 on (0, n) and gn00 (x) > 0 on (n, ∞). It can be proved that gn (n) > 1/2 for all n ≥ 1 and gn (n) monotonically decreases to 1/2 as n approaches innity. Let the solution in x to the equation Clearly since

Dene

ψn (x)

µ.

0 < p < 1,

(75)

as the natural logarithm of

gn (x),

n X xl l=0

© 2007 ACADEMY PUBLISHER

ψn (x), it follows

x

(80)

to the equation

This linear equation in

x

is easily solved and

(81)

µ0

is

obtained as:

µ0

= =

We claim that

µ0

ln p − ψn (x0 ) ψn0 (x0 ) ln p − ψn (n) . n+ ψn0 (n)

x0 +

is an upper bound for the solution

(82)

µ

to

ln p ≥ ψn (µ0 ) = ln(gn (µ0 )). Hence p ≥ gn (µ0 ). Since gn (µ) = p, we have gn (µ) ≥ gn (µ0 ). Hence µ ≤ µ0 , since gn (x) is a decreasing function of x. Thus µ0 is an upper bound for µ. We can simplify the expression (82). Assume 0 < p ≤ 2/e. Since gn (n) is a monotone decreasing sequence, gn (n) ≤ g1 (1) = 2/e and therefore ψn (n) ≤ ln(2/e). It follows that

µ0

(83)

ln p − ln(2/e) ln p − ψn (n) ≤n+ =n+ ψn0 (n) ψn0 (n) # µ ¶ " n n3 n2 1 + n + 2! + 3! + · · · + nn! 2 . ln ≤n+ nn pe n! This upper bound can be further simplied. Using the

We now derive an upper bound for

ψn (x) = ln(gn (x)) = −x + ln

(79)

the equation (75). In fact, from (80) and (81), we see that

It can be veried that

be denoted by

ψn00 (x) < 0

ln p = (x − x0 )ψn0 (x0 ) + ψn (x0 ).

dene

= p,

(77)

(x − x0 )ψn0 (x0 ) + ψn (x0 ) ≥ ψn (x), x ∈ (0, ∞).

,

We next prove the inequality (18). For each positive

l!

,

y = (x − x0 )ψn0 (x0 ) + ψn (x0 ).

This proves the chain of inequalities in (16).

l=0

xn n!

that

fn (Tn (Pf a )) = Pf a , 1 ≤ n ≤ LN . Since fn (x) ≤ fn+1 (x), 1 ≤ n ≤ LN − 1, and since fn (x) are decreasing functions of x, we see immediately that Tn (Pf a ) ≤ Tn+1 (Pf a ), 1 ≤ n ≤ LN − 1. Also T1 (Pf a ) = −σ 2 λ1 ln(Pf a ) and TLN (Pf a ) = T .

n X xl

+ ··· +

for all x ∈ (0, ∞). Hence ψn (x) is concave on the positive real axis (0, ∞) and is always below its tangent lines. Let x0 = n. The equation of the tangent line of the curve of ψn (x) at the point (x0 , ψn (x0 )) is given by:

From the denition (15), we see that

gn (x) = e−x

x2 2!

As this tangent line is above the curve of

2 ≤ n ≤ LN.

gn00 (x) =

1+x+

The above expression shows that

1≤l≤n,l6=m

n,

xn n!

(78)

(70)

and

m=1

integer

fact,

Since

probability density functions are supported on the positive

(n+1 X σ 2 λl

shown to be concave on the

i h Pn−1 n 1 )xk − (k−1)! xn−1 n + k=1 ( k! . ψn00 (x) = − ¡ n ¢2 2 n! 1 + x + x2! + · · · + xn!

(69)

chi-square random variables are non-negative (i.e., their real axis), we have:

ψn (x) can be (0, ∞). In

positive real line

l!

.

µ.

that is, (76)

following result for large

n,

gn (n) = (84) ¶ µ n 3 2 n n n ∼ + ··· + + e−n 1 + n + = 1/2, n! 3! 2!

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

45

Proof of Theorem 5. From (64) we see that

the expression

1+n+

n2 2!

n3 3!

+

+ ··· +

nn n!

nn n!

Pf a = Pr

(85)

√ can be shown to be bounded from above by



2πn

³ n ´n e

(86)

we see that

Ã

1+n+ ³



exp(n) 2 nn n!

√ ≈

2πn

n2 2!

+

n3 3!

+ ··· +

nn

nn n!

!

= ¡ n ¢n e

n! exp(n) 2nn

T ≥ BL (Pf a ) σ 2 µ1

exp(n)

2nn

=

2πn . 2

Pf a = Pr

(87)

µ≤n+

2πn ln 2

µ

2 pe

¶ , 0 < p ≤ 2/e.

(88)

The inequality (18) then follows directly from (88). Finally we prove (19). Since

Pf a = Pr

(L−1 X

) zl,k > T

where

l=0

= Pr

( LN X σ 2 λl 2

l=1 ( LN X

)

where

χ22LN

³ Pf a ≤

s=0

T σ 2 λ1

s!

χ22LN

e

.

T ≤ BLN (Pf a ), σ 2 λ1

T ≤ σ 2 µ1 BLN (Pf a ).

(95)

in increasing order as follows: (96)

Then we have

(89)

Ly1 ≤

2LN

(90)

L−1 X

zl,k ≤ LyL ,

(97)

l=0 and it follows that

(L−1 X

) zl,k > T

(98)

l=0

≥ Pr {Ly1 > T } = Pr {y1 > T /L} = Pr {zl,k > T /L, 1 ≤ l ≤ L − 1} L−1 Y = Pr {zl,k > T /L} l=0

(91)

" =

which completes the proof of (19).

or

y1 ≤ y2 ≤ y3 ≤ · · · ≤ yL .

Pf a = Pr T σ 2 λ1

2LN

It remains to prove (23). Rearrange the power levels

This implies

© 2007 ACADEMY PUBLISHER

is a chi-square random variable with

zl,k , 0 ≤ l ≤ L − 1,

´s −

(94)

Combining (93) and (95) yields (20).

degrees of freedom, we have

LN −1 X

zl,k > T

l=0

T ≤ BLN (Pf a ) σ 2 µ1

)

is a chi-square random variable of

)

degrees of freedom. This implies that

χ22 (l) > T

σ 2 λ1 2 χ2 (l) > T 2 l=1 ¾ ½ 2 σ λ1 2 χ2LN > T = Pr 2 ´s ³ T LN −1 X σ 2 λ1 − T e σ2 λ1 , = s! s=0

≤ Pr

(93)

) N X σ 2 µm 2 χ2L (m) > T = Pr 2 m=1 ) ( N X σ 2 µ1 χ22L (m) > T ≤ Pr 2 ¾ ½ m=1 σ 2 µ1 2 χ2LN > T = Pr 2 ´s ³ T LN −1 X σ 2 µ1 − T e σ2 µ1 , = s! s=0

µ: √

(92)

σ 2 µ1 BL (Pf a ) ≤ T.

(L−1 X

(

Substituting (87) into (83), we obtain the following upper bound for

or

Similarly from (64) we also see that



)

This implies that

n!

´

zl,k > T

N X σ 2 µm 2 χ2L (m) > T = Pr 2 m=1 ¾ ½ 2 σ µ1 2 χ2L (1) > T ≥ Pr 2 ´s ³ L−1 X σ2Tµ1 − T e σ 2 µ1 . = s! s=0

2πn 2 . In

,

)

l=0

(

fact, using Stirling's formula

n! ∼ =

(L−1 X

N X m=1

T

−1 µN − L m Q e σ2 µm 1≤l≤N,l6=m (µm − µl )

#L .

46

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

or

This implies

N X m=1

−1 µN m Q e 1≤l≤N,l6=m (µm − µl )

T − 2L σ µm



1/L Pf a

k6=i (αk − α) αim Q k6=i (αk − αi ) i=1

(99)

T ≥ T1 , L

LT1 ≤ T.

or

n X

(100)

(−1)n−1

i=1

Similarly, we have

(L−1 X

)

If

zl,k > T

m = n − 1,

(101)

n X

l=0

≤ Pr {LyL > T } = Pr {yL > T /L} = 1 − Pr {yL ≤ T /L} = 1 − Pr {zl,k ≤ T /L, 1 ≤ l ≤ L − 1} L−1 Y =1− Pr {zl,k ≤ T /L} =1−

l=0 L−1 Y

=1− 1−

Letting

α

n X

αin−1−m Q

n−1

(−1)

(−1)

βi

T

#L

α)

= 1.

k6=i (αk − αn−1

α→+∞

Q lim

−1 µN − L m Q e σ2 µm 1≤l≤N,l6=m (µm − µl )

k6=i (αk − αn−1

βi lim

i=1

α→+∞

(109)

k6=i (αk − αn−1

α)

α)

= 1.

= (−1)n−1 .

(110)

(111)

Substituting this result into (110) yields:

.

n X

βi = 1.

(112)

i=1

N X

T

Q

m=1

−1 µN − L m e σ2 µm 1≤l≤N,l6=m (µm − µl ) 1/L

≥ 1 − (1 − Pf a )

,

This proves the rst identity in (105). To prove the

T ≤ T2 , L

or

α=0

(102)

T ≤ LT2 .

(103)

in (108). We then obtain the identity

n X

Y

βi αk n−1−m α i=1 i k6=i

Dividing both sides of (113) by

Combining (100) and (103) yields the inequalities of

n X

(23). The inequality (24) is obtained by using (102) and applying the argument in the proof of (19) in Theorem 4. The proof of Theorem 5 is thus completed. To prove Theorem 6, we need to prove two lemmas. Lemma 1. Let

n ≥ 2.

Assume

0 < α1 < α2 < · · ·
0.

(116)

Proof of Lemma 2. We prove this lemma by mathemat-

1 ≤ m ≤ n−1 be an integer and m dene the polynomial Pm by Pm (α) = α . According to Lagrange's Interpolation Theorem, the following identity holds for all

n X

© 2007 ACADEMY PUBLISHER

i=1

where

Proof of Lemma 1. Let

α ∈ (−∞, +∞): Q k6=i (αk − α) = Pm (α), Pm (αi ) Q k6=i (αk − αi ) i=1

1 ≤ m ≤ n−1

remaining identities in (105), assume and let

which in turn yields

Then

(108)

k6=i

Q

This implies

αn

(αk − α) = αm .

go to innity yields the identity

n−1

[1 − Pr {zl,k > T /L}]

m=1

(107)

(108) can be rewritten as

But

N X

Y

βi

i=1

l=0

"

= αm .

This identity can be rewritten as

which in turn yields

Pf a = Pr

Q

n X

(106)

ical induction. First consider the case

n = 2.

We have

C(p, α , α ) = (117) · ¸ ·1 2 ¸ α2 α2 − 1 pα2 /α1 β1 α1 α1 ¸ ¸ · · ¸· α1 α2 α2 < 0, − 1 pα2 /α1 = α1 − α2 α1 α1

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

since

α2 > α1

and

p ∈ (0, 1).

47

∂E ∂p (p, α1 , α2 , · · · , αn ), viewed as a function of p, is strictly on the interval [0, 1]. Using the Pn decreasing βi identity, i=1 α = 0, proved in Lemma 1, we see that

Hence,

On the other hand,

D(p, α1 , α2 ) · ¸h · ¸h i i α2 α2 α2 /α1 pα2 /α2 β2 p β1 + = α2 α ¸ ¸ · · 1¸· α2 α1 α2 p pα2 /α1 − = α1 − α2 α1 − α2 α1 ¸h · i α2 pα2 /α1 − p > 0, (118) = α1 − α2 α2 > α1 and p ∈ (0, 1). This proves Lemma 2 for the case n = 2. Next assume Lemma 2 holds for positive integer n − 1. We show that Lemma 2 must also hold for the integer n. To do this, we reformulate the expression C(p, α1 , · · · , αn ) as follows:

i

∂E (p, α1 , α2 , · · · , αn )|p=1 ∂p ¡ −1 ¢ = p D(p, α1 , α2 , · · · , αn ) p=1 = D(1, α1 , α2 , · · · , αn ) n X βi = 0. = αn α i=1 i

again using the fact that

C(p, α1 , α2 , · · · , αn ) (119) ¸ ¸ · · n−1 X αn αn − 1 pαn /αi βi = α α i i i=1 n−1 X · αn ¸ · αn − αi ¸ αin−1 Q pαn /αi = αi αi k6=i (αi − αk ) i=1 = −αn = −

n−1 X

αin−2 1 Q pαn /αi αi 1≤k≤n−1,k6=i (αi − αk )

i=1 n−1 αn X

αn−1

i=1

³

αn−1 αn /αn−1 p αi

´αn−1 /αi

βi∗ ,

It follows that for all

(127)

and consequently

D(p, α1 , α2 , · · · , αn ) ¶ µ ∂E (p, α1 , α2 , · · · , αn ) > 0. =p ∂p

(128)

By the principle of mathematical induction, Lemma 2

n ≥ 2.

must hold for all integers

Proof of Theorem 6. Applying Lemma 2 to the

λLN < λLN −1 < · · · < λ2 < λ1 , we g 0 (p) > 0 and g 00 (p) for 0 < p < 1.

see immediately that

αin−2 . 1≤k≤n−1,k6=i (αi − αk )

βi∗ = Q

Proof of Theorem 7. Since (c.f. (27))

(120)

L−1 X

Hence

zl,k ∼

l=0

C(p, α1 , α2 , · · · , αn ) = αn αn D(p αn−1 , α1 , α2 , · · · , αn−1 ). − αn−1 D(p, α1 , α2 , · · · , αn−1 ) > 0,

(121)

To prove

Pf a = Pr

we have

D(p, α1 , α2 , · · · , αn ) > 0,

E(p, α1 , α2 , · · · , αn ) =

n X

( = Pr

where

(123)

i=1 It can be veried that

      

∂E ∂p (p, α1 , α2 , · · · , αn ) = p−1 D(p, α1 , α2 , · · · , αn ), ∂2E ∂ 2 p (p, α1 , α2 , · · · , αn ) = p−2 C(p, α1 , α2 , · · · , αn ).

C(p, α1 , α2 , · · · , αn ) < 0,

) zl,k > T

l=0 N X

)

µm χ22L (m) > 2T /σ 2 (129)

¡ ¢ 1/2 h0 = c32 /c23 , y = 2T /σ 2 − c1 (h0 /c2 ) +h0 and  PN   c1 = 2L Pm=1 µm = 2LN, N (130) c2 = 2L m=1 µ2m ,  P  c = 2L N µ3 . 3 m=1 m

Here we used the fact that

PN m=1

µm = N , which follows

from the assumption that the window is normalized. We have

³

it follows that

∂2E (p, α1 , α2 , · · · , αn ) ∂2p = p−2 C(p, α1 , α2 , · · · , αn ) < 0. © 2007 ACADEMY PUBLISHER

(124)

(L−1 X

© m=1 ª ≈ Pr χ2h0 > y ,

let

pαn /αi βi .

N X σ 2 µm 2 χ2L (m), 2 m=1

we see from Eq. (4.1) of [15] that

C(p, α1 , α2 , · · · , αn ) = (122) αn αn αn−1 , α1 , α2 , · · · , αn−1 ) < 0. D(p − αn−1

Since

0 < p < 1,

∂E (p, α1 , α2 , · · · , αn ) ∂p ¶ µ ∂E (p, α1 , α2 , · · · , αn ) |p=1 = 0, > ∂p

increasing sequence

where

Since

(126)

h0

=

2L ³ 2L

(125)

³P

PN

´3

PN

´2 = 2L ³P

2 m=1 µm 3 m=1 µm

N m=1

µ2m

N m=1

µ3m

´3 ´2 , (131)

48

JOURNAL OF COMPUTERS, VOL. 2, NO. 6, AUGUST 2007

and

[8] S. Wang, F. Patenaude and R. Inkol, “Upper and Lower Bounds for the Threshold of the FFT Filter Bank-Based Summation CFAR

1/2

0

(h /c2 ) ¢3 1/2  ¡PN µ2 m=1 m ¡ ¢2  P 2L PN N  µ3 µ2m   m=1 m . = Pm=1 =  PN N 3  2L m=1 µ2m  m=1 µm

Detector,” Proceedings of ICASSP 2006, Toulouse, France, vol.3, pp. 289-301, May 2006. [9] S. Wang, R. Inkol, S. Rajan and F. Patenaude, “Theory of the FFT Filter Bank-Based Majority and Median CFAR Detectors,” DRDC

(132)

Ottawa TR 2007-088, May 2007. [10] S. Wang and R. Inkol, “Operating Characteristics of the Wideband FFT Filter Bank

J -out-of-L

CFAR Detector,” DRDC Ottawa TR

2003-235, December 2003. [11] R. Inkol and S. Wang, “A Comparative Study of FFT-Summation

It follows that

and Polyphase-FFT CFAR Detectors,” Canadian Conference on

¡ ¢ 1/2 y = 2T /σ 2 − c1 (h0 /c2 ) + h0 ¡

= 2T /σ 2 − 2LN

¢

³P

PN

2 m=1 µm PN 3 m=1 µm

(133)

N m=1

µ2m

N m=1

µ3m

+ 2L ³P

´3

´2 .

© ª © ª Pf a ≈ Pr χ2h0 > y ≈ Pr χ22h > y ª © ≈ Pr χ22h /2 > y/2 ¶ µ z h−1 z2 , + ··· + = e−z 1 + z + (n − 1)! 2! $ ¡P ¢3 % N

z = y/2

Summation CFAR Detectors: A Comparative Study,” Canadian Conference on Electrical and Computer Engineering, pp. 10391044, May 2004. [13] S. Wang and R. Inkol, “Theoretical Performance of the FFT Filter Bank Based Summation Detector,” DRDC Ottawa TR 2005-153, November 2005.

Consequently

where

Electrical and Computer Engineering, pp. 1175-1178, May 2004. [12] S. Wang and R. Inkol, “FFT Filter Bank-Based Majority and

and

h = L ¡P

m=1 N m=1

µ2m

µ3m

¢2

[14] S. Wang, R. Inkol and S. Rajan, “A Formula for the Probability of False Alarm for the FFT Filter Bank-Based August 2005.

(134)

¡ ¢ T /σ 2 − LN

2 m=1 µm PN 3 m=1 µm

in Normal Variables,” Biometrika, vol.48, pp. 419-426, December 1961. [16] G.L. Turin, “The Characteristic Function of Hermitian Quadratic Forms in Complex Normal Variables,” Biometrika, vol.47, pp. 199201, June 1960.

. Hence

[17] H. Solomon and M. A. Stephens, “Distribution of a Sum of Weighted Chi-Square Variables,” Journal of the American Statistical Association, vol.72, no.360, pp. 881-885, December 1977.

³P

PN

CFAR

[15] J. P. Imhof, “Computing the Distribution of Quadratic Forms

z = y/2 ≈ Bh (Pf a ), or

J -out-of-L

Detector,” Proceedings of MWSCAS 2005, Cincinnati, Ohio, USA,

N m=1

µ2m

N m=1

µ3m

+ L ³P

≈ Bh (Pf a ).

´3 Sichun Wang (DRDC Ottawa and CRC) received the B.Sc. and M.Sc.

´2

degrees, both in Mathematics, from Nankai University in 1983 and 1989, respectively. He received the PhD degree in Mathematics from McMaster University in 1996. He was a postdoctoral fellow at the

(135)

Communications Research Laboratory of McMaster University from February 1996 to August 1996 and a research scientist at Telexis

The approximation (28) then follows immediately from

and Intrinsix Corporations from September 1997 to May 2002. From

(135).

June 2002 to March 2007 he worked as a research consultant for

The proof of Theorem 8, as it is almost identical to that of Theorem 7, is omitted here.

Defence Research and Development Canada. Currently he is a research scientist at the Communications Research Centre (CRC) on secondment to Defence Research and Development Canada. His recent research has focused on the FFT lter bank-based CFAR detection algorithms and

R EFERENCES

turbo codes. Dr. Wang is a senior member of the IEEE.

[1] E. O. Brigham, The Fast Fourier Transform and Its Applications, Prentice-Hall, Inc., NJ, USA, 1988. [2] G. A. Zimmerman and S. Gulkis, “Polyphase-Discrete Fourier Transform Spectrum Analysis for the Search for the Extraterrestrial Intelligence Sky Survey,” TDA progress report 42-107, Jet Propulsion Laboratory, Pasadena, CA., pp. 141-154, July-September, 1991. [3] G. Lanyi and R. Khan, “Tone Detection via Incoherent Averaging of Fourier Transforms to Support the Automated SpacecraftMonitoring Concept,” TDA Progress Report 42-129, Jet Propulsion Laboratory, Pasadena, CA., pp. 1-22, May 1997. [4] R. S. Walker, “The Detection Performance of FFT Processors for Narrowband Signals,” D.R.E.A. Techincal Memorandum 82/A, Defence Research Establishment Atlantic, Research and Development

Franc ¸ ois Patenaude (CRC) received the B.A.Sc. degree from the University of Sherbrooke in 1986 and the M.A.Sc. and Ph.D. degrees from the University of Ottawa in 1990 and 1996, all in Electrical Engineering. His graduate work was conducted in collaboration with the Mobile Satellite Group of the Communications Research Centre (CRC), Ottawa, Canada. In 1995, he joined CRC to work on signal processing applications for communications and for spectrum monitoring. His main research interests include automatic signal recognition, detection and estimation in the spectrum monitoring context, and real-time signal processing.

Branch, Department of National Defence, Canada, February 1982. [5] B. H. Maranda, “On the False Alarm Probability for an Overlapped FFT Processor,” IEEE Transactions on Aerospace and Electronic Systems, vol.32, no.4, pp. 1452-1456, October 1996. [6] H. C. So, Y. T. Chan, Q. Ma, P. C. Ching, “Comparisons of Various

Robert Inkol (DRDC Ottawa) received the B.Sc. and M.A.Sc. degrees

Periodograms for Sinusoidal Detection and Frequency Estimation,”

in Applied Physics and Electrical Engineering from the University of

IEEE Transactions on Aerospace and Electronic Systems, vol.35,

Waterloo in 1976 and 1978, respectively. Since 1978 he has been with

no.3, pp. 945-952, July 1999.

Defence Research and Development Canada where he is currently a

[7] F. Patenaude, D. Boudreau and R. Inkol, “CFAR Detection Based on

senior scientist. He has made numerous contributions to the application

Windowed and Polyphase FFT Filter Banks for Channel Occupancy

of very large scale integrated circuit technology and digital signal

Measurements,” Proceedings of the 19th Biennial Symposium on

processing techniques to electronic warfare systems. He is an adjunct

Communications, Kingston, Ontario, Canada, pp. 339-343, May

professor with the Royal Military College and is a senior member of

1998.

the IEEE.

© 2007 ACADEMY PUBLISHER