- Email: [email protected]
Polynomial-phase signal analysis using stationary moments
SIGNAL PROCESSING Signal Processing
54 (1996) 239-248
Polynomial-phase signal analysis using stationary moments A. Ferrari*, C. Theys, G. Alengrin Lahoratoiw 13s. CNRS, UniversitP de Nice Sophia-Antipolis, 41 Bd NupolPon III, 06041 Nice Cedes. Frunw Received
22 December
1995: revised
13 June 1996
Abstract This article proposes a method for polynomial-phase signals analysis relying on their time-invariant high-order moments. The set of these moments is first characterized in the noiseless and noisy case. It is demonstrated that the only identifiable phase parameter is the highest-order coefficient, the estimation requiring moments of order at least the double of the phase degree. Next, the consistency of the estimated moments calculated by sample averages is proved. Finally, an estimation algorithm based on the moments having for lags the multiples of a given ‘root’ is developed. Performances, evaluated by computer simulations, show the influence of the different parameters and the efficiency of the method. Zusammenfassung Der Beitrag schllgt eine Methode vor, um Signale polynomialer Phase unter Beriicksichtigung ihrer zeitinvarianten Momente hoherer Ordnung zu analysieren. Die Menge dieser Momente wird zunachst Rir den rauschfreien und den verrauschten Fall charakterisiert. Es wird gezeigt, dab der einzige identifizierbare Phasen-Parameter der Koeffizient hochster Ordnung ist, und die Schlzung Momente erfordert, deren Ordnung mindestens das doppelte des Phasengrades ist. Danach wird die ijbereinstimmung der geschatzten Momente durch Mittelung bewiesen. Abschliegend wird ein Schatzalgorithmus auf der Basis von Momenten, ausgehend von einem Vielfachen einer gegebenen ‘root’, entwickelt. Das Verhalten, durch RechnerSimulationen bewertet, zeigt den EinfluS der verschiedenen Parameter und die Leistungsfahigkeit der Methode. RCsumC Cet article propose une methode d’analyse de signaux a phase polynomiale utilisant leurs moments d’ordre superieurs qui sont invariants dans le temps. L’ensemble de ces moments est tout d’abord caracterise dans le cas non bruit6 et bruit& On demontre en particulier que seul le coefficient du terme de plus haut degre est identifiable et ce uniquement a partir de moments dont l’ordre est superieur ou tgal au double du degre de la phase. La consistance des estimees de ces moments, calculees par moyenne temporelle est ensuite proude. Enfin, un algorithme d’estimation utilisant les moments dont les retards sont les multiples d’une ‘racine’ donnee est developpe. Les performances Cvaluees par simulations montrent l’influence des differents paramttres ainsi que I’efficacitt de cette mtthode.
K~Jww~‘s: Polynomial-phase
* Corresponding
signals; High-order
author. Tel.: +33-93217953;
moments;
fax: +33-93212054;
Parameter estimation
e-mail: [email protected].
0165-1684/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved PZI: SO165-1684(96)001 IO-7
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
240
1. Introduction This article addresses the estimation of the parameters of a polynomial-phase signal. This problem, encountered in different radar systems, is usually solved using a time-frequency analysis or phase-only algorithms. For a detailed introduction on this topic see, for example, [5]. This paper investigates a parametric solution based on statistical properties of this signal: although it is clearly non-stationary, some of its high-order moments are shift-invariant. This property is exploited for the retrieval of the phase parameters. The paper is organized as follows. In Section 1, the conditions verified by the delays of these ‘stationary’ moments are derived in the noiseless and noisy case. Section 2 is devoted to the statistical study of the estimated moments. In particular, consistency of the sample estimates is proved. Finally, in Section 3 an estimation algorithm is derived and its performances analyzed by computer simulations.
Proposition 1. The signal x, in (1) is (II,. . . , ZQ)stationary if and only if the following condition, referred as H( p, M), holds: Wp,M):
Vq E
(11 - Ii,,)
f:
{l,...,
M - 1},
= 0.
(3)
k=l
In this case the stationary moments equal M2p,xU1,.
. . , 12p)
= A2P exp
j aM f:
(1;1”
(4)
k=l
Proof. The 2pth-order moments equal A2P exp {j$} with
4=
2 9,aq((n -k k=l
lk)’
M&n;
11,.. . , Z2p)
- (n + lk+p)‘)
q=l
2. High-order moments of polynomial-phase signals 2.1.
The noiseless case
We define a noisy polynomial-phase
signal y, as
yn =G+G,
(1)
where a0 is a random variable uniformly distributed in the interval [0,2rc) and e, is a complex-independent identically distributed (iid) noise. The 2pth-order moment of the signal z, is defined as
M2,,&;
11,. . . ,
(2) The signal z, is said (11,. . . , 12P)-stationary if Mz,,,(n; 11,.. . ,1+) is shift-invariant, i.e., M,,,(n; 11, . . . ,lzp) is not a function of it. In this case, the index n will be omitted in M,,(n; 11,.. . ,lzp). It is worthy to note that the definitions obtained with a different number of conjugated and unconjugated terms lead to a 2pth-order moment of x, that equals zero.
This polynomial of the variable n is not a function of n if and only if ‘dm = 1,. . . , M, b, = 0. bl = 0 implies CE, aq(‘f)ckp,l (Ii-’ - I:;;) = 0. This equation has to be verified Va,. This implies H(p,M). Conversely, it can be easily checked that H(p,M) implies Vm = 1,. . . ,M, b, = 0. The coefficient bo equals Ego a4 Ckp,l (1: - lz,,). Application of H(p,M) to this equality terminates the proof of Eq. (4). 0 This result suggests the following remarks: - The only identifiable parameter of a polynomialphase signal from its stationary moments is the coefficient of highest degree, aM. The estimation of all the ak can be performed using the iterative
241
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
algorithm consisting in: estimate UM, multiply the signal by exp {-j aM@‘}, estimate CIM-_I, . . . . ~ H(p,M) is a set of M - 1 non-linear equations for 2p integer variables. A trivial solution consists in choosing for (l,+t , . . . , l+,) a permutation of (11,. . . , lP). Unfortunately, this solution, according to (4), leads to a moment independent of a,+,,. An important point is then to find non-trivial solutions of H(p,M). An analytic solution of the diophantine system H( p,M) is well beyond the scope of this paper. The solution consisting in the use of nested loops together with constraints on the delays to reduce the computation will be discussed further. To illustrate these results, we apply Proposition 1 to the cases M = 2,M = 3, assuming without loss of generality that 11 = 0. Example 1 (x, is u chirp: M = 2). - lf p = 1,12 must verify: 12 = 0. The signal is only (0,0)-stationary and Mz,(O, 0) = A*. _ If p = 2, (12,1x, 14) must verify 12- 13- 14 = 0. The signal is (0,13+14,13, /4)-stationary and M+JO, 13+ 14,139 14)
= A4
exP{j
M2,wxQA
12,.
> 12~4)
(-l)“jMQM
fi
lM+k
(5)
k=l
Proof. Consider the two p x p diagonal matrices C = Diag{lr, . . , Ip}, D = Diag{ I,+, , . , lzp}. The condition H( p, M) can be rewritten in terms of these matrices as H(p,M):
Vq E {I,...,
M-
l},
tr{Cq} = tr{Dq},
(6)
where tr{M} is the trace of the matrix M. If we apply the Hamilton-Cayley theorem, [3], to these matrices, we obtain cp - c,cp-I
- c*cp-2
- .” - cp-,c
- cpI = 0, (7)
DP - dlDP-’
- d2DP-2 -
. - d,_,D
- dPI = 0,
2Q21314)).
(8)
Example 2 (x, is an hyperbolic phase signal: M=3). - If p = 1, 12 must verify: 12 = 0 and 1: = 0. The signal is (O,O)-stationary and Mz,JO, 0) = A’. - Ifp=2, (12,13,14)mustverify12-ls-&=Oand 1: - 13- 1: = 0. These two equations imply 1314 = 0. The signal is (0,14,0,14) and (0,13,ls,O)-stationary and M+,(O, LO, 14) = &JO, k,k 0) = A4. - If p = 3, (12,13,14,15,1~) must verify 12+ 13 - 1415 - 16 = 0 and 1; + 1: - 1: - 1: - 1: = 0. An analytical solution to H(3,3) is not available; however nontrivial solutions exist, see Table 1. Using the two previous equalities, it is possible to show that the six order stationary moments of x, equals A6 exp{-j 3u3141s16). These results proposition.
of order 2M of x, equal
are generalized
by the following
Proposition 2. The coeficient aM of the signal in ( 1) is only identifiable from the stationary moments of order higher or equal to 2M. The stationary moments
where ck and dk are the coefficients of the characteristic polynomial of C and D. The Faddeev method [3] allows the recursive computation of these coefficients using the relations kCk = tr{Ck} - cl tr{c”-‘}
- . . - ck-1 tr{C},
k = l,...,p, kdk = tr{Dk} - dl tr{Dk-‘}
(9) - . ‘. - dk_1 tr{D),
k = l,...,p,
(10)
with CO= 1 and do = 1. The recursive H( p, M) to these equations gives VqE{l,...,
min(M-l,p)},
cq=dq.
application
of
(11)
Consider first the case p < M. The trace of the difference of the two equalities (7) and (8) together with (11) and N(p,M) gives tr{Cp} =tr{DP}. If we now multiply Eq. (7) by C and Eq. (8) by D, the trace of the difference of these two equalities will give tr{C P+* } - tr{DJ’+‘}. Iterating, we then obtain
242
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
‘dq, tr{Cq} = tr{Dq}. The application of this result with q=M to (4) terminates the proof of the first assertion of the proposition. Consider now the case p = M. In this case, the trace of the difference of the two equalities (7) and (8) together with ( 11) and H(M, M) gives tr{C”}
- tr{oM}
(A)~-’
=
(12)
fi lk - (-l)“-ifi k=l
lp+k
0,13)M2,e(l3
+M,,x(n;
l3 +
+Mz,x(~
0,14)M2,,(13
+Mz,x(n;
13 +
+
l4914)
14, ~4)M2,e(o, +
13) 14,l3)
14,~3)M2,e(o,
(15)
14).
x, being only (O,O)-stationary, stationarity of this moment, independent of e,, requires 13 = 14 = 0. This solution being useless, stationarity will be obtained using the independence of the noise samples. In fact, the constraint 13 # 0 and 14 # 0 will imply the nullity of all the noise second-order moments:
k=l
(13) Application the proof.
+M,,x(n;
v’I3
of this result with Ii = 0 to (4) terminates 0
#
0914
M4,Jo,
#
0,
13 +
14,13,14)
=
M4,x(o,l3
+
14,139 14).
(16)
The noisy case
2.2.
In the noisy case, we will consider the moments of the signal yn at lags for which x, is (11,. . . , IQ,)stationary. For brevity, we introduce the ordered k-uple 2, constituted by a subset of k elements of (I 1,. . . , lp). 9: is a k-uple of (I,+,, . . ., lQ),?& the complementary of yk in (Ii, . . . , lp) and 3: the complementary of 9: in (Z,+i, . . . , 1~). Consequently, we will use the condensed notation M.&n; &, 2:) for the 2kth-order moment of z,,. Using the circularity of e, and the hypothesis on a0 distribution, we can expand the 2pth-order moment of yn as M2p,y(%
2;)
Unfortunately, the (Ii,. . . , IQ,)-stationarity of xB does not imply the (_&, _%‘f )-stationarity of x,, and then the (Ii , . . . , 12P)-stationarity of yn. To illustrate this difficulty we consider the following example.
Example 3 (x, is a chirp: M = 2, and p = 2). If we expand M+y at the stationary Example 1, we obtain M4,@, =
13 + M4,x(0,/3
14,139 +
lags of M+,
+
M4,e(0,13
+
to the general
Proposition 3. For all solutions of H( p, M) verifying (1 I,..., (II,...
lp} n {4,+1,.. , 12p)=
MzM,~(O,
M2,,x(h
., hp} = 0 , . . . , 12p).
we have MQ,,~ Moreover, if p =M,
12,. . . ,12~)
= A2M eXp
1)“j
(-
hfLZ,+f fi
l,,,f+k
.
(17)
k=l
It is interesting by now to investigate the possibility of using cumulants instead of moments in order to solve at the same time the additive noise problem and the multicomponent signal case. The properties of cumulants, [l], imply that for xn defined in (l),
A = cum(x,+h,.
=
~,.M,x(n;
. . ,JW,,~,*+~,+,,
yp,
2;
),
. . . ,~n*,~~,)
(18)
(19)
see
l4) 14,137 l4)
This result can be easily extended case using (14).
I49 13914)
with y,, = cum( ej@, . . . , ejQ, e-jao, . . . , e--jao). Assume that y,, is formed by m noisy polynomialphase signals of order M, xi, with independent initial phases. The additive property of the cumulants of the
A. Ferrari et al. /Signal
sum of independent
processes implies
If each component is (11, . . . , IQ))-stationary, (19) implies that C&,,_,,(n; LZ”,2; ) is not a function of n and
Meanwhile,
this solution is not tractable because is defined from the M2k,Jn; P”,,9:) mentioned above, are not necessary
C,.,(P& 2:) which, as stationary. 2.3.
The set of stationary moment delays
=
0.
In the case of a parabolic-phase signal an explicit solution exists, see Eq. (16), but, for higher degrees,
243
analytical solutions are not available. We suggest to find them using 2M - 1 loops to test condition H(A4,M). In this context, it is important in order to reduce the computational cost of the scan to define a set of minimal size containing the lk. This is done in a way similar to [7]. The moment is simultaneously symmetric with respect to the first and the last A4 parameters lk. The smallest lk can then be moved to position 1 or position M+l;infactifmin{lr,...,l~}=min{l~+~,...,/~~} the two sets are not disjoints. In this last case, the obvious relation M~MJJ&, 2;) = MzM,J_Y~, Y,)* allows to move the smallest lk to position 1. From the definition of stationarity, if (11,12,. . . ,L~M) verifies H(M,M),Vr,(r + ll,r + 12,. . . , r + 12~) also verifies H(M,M). According to this, the first and smallest delay can be subtracted to the others leading to delays that are all nonnegative for the first M and positive for the last M. Both these sets can then be reordered in increasing order giving
1, = 0,
We discuss in this section the computation of the delays lk leading to a stationary moment of y,. We will limit our study to the minimal order moment p = M. According to the previous sections, we need to find the 2M-uple verifying H(M,M), { 11,. . . , 1~) fl {1A4+1,...,12M}
Processing 54 (1996) 239-248
0<1/2< ... 6lM,
161&f+, < ...
dlZM. (22)
An important remark is that if (Ii, 12,. . ,Z~M) verifies H(M,M) and {It,., ., 1~) n {l~+l,..., 12~) = 0,Vr # O,{rll,rl2 ,..., rZZ*} also verifies these two conditions. A direct consequence is that the set of interest can be reduced to only its coprime members. An element of this set will be called a root of order M. Table 1 gives the 11 first roots of order 3 and 4.
Table I Eleven first roots of order 3 and 4 12
13
14
3 4 5 5 6 6 7 7 7 9 9
3 5 7 8 9 11 7 8 II 9 10
I I I 2 I 2 I 2 I 3 I
15
16
-3141510
12
13
14
1s
16
17
lx
4w6h18
I
4 6 8 9 10 I2 9 IO 12 12 12
-12 -36 -7 -108 -120 -216 -108 -180 -180 -324 -216
3 4 5 5 5 6 6 7 7 8 8
4 7 5 10 12 7 I3 9 11 9 II
7 II 10 I5 17 I3 19 16 18 17 I9
1 1 I I 2 1 1 I 2 2 I
1 2 2 3 2 3 4 4 3 3 5
6 9 8 12 15 10 I.5 I2 I5 14 I4
6 10 9 I4 I5 12 18 15 I6 I5 I8
144 720 516 2016 3600 1440 4320 2880 5760 5040 5040
2 3 2 4 3 4 3 5 3 6
A. Ferrari et al. /Signal Processing 54 (1996) 239-248
244
3. Statistical properties of the sample estimates The stationary moments of the signal under scope may be estimated on the basis of a single realization of the time series y,, n = 1,. . . , N. The sample moment function can be estimated as
where ZM = max{_YP U 9;}, I, = min{ZP U 3.) and N, = N - 1~ + 1,. We address in this section the statistical properties of these estimates. For this sake, let
(26) The term of this sum corresponding to the empty sets _5&=0and A?:=0 equals M2,&Zpp,2;). Hence, dM$,,( TP, 9’) equals the sample mean of (26) where the sum ranges only over the non-empty subsets LYEu 9: # 0. Consequently, the expansion of the squared modulus of this sum provides a first expression of ~dM$JYP,Ip,*)~2. Using the notation (2), the expectation of this expansion reduces to the sum over the subsets _%,,49,* and Y,,_??z of -
1
with _CZkUL?z # 0 and LZlUYz
M$,,J-%2; I- ~wG% 3; ).
M,,e(u;2k
u {v
-
ZJ +
y;},
(24)
It is straightforward that the sample moment function (23) is an unbiased estimator of MQ,,J 5$,, 9; ), hence
E{dM&,(~p,2; )I = 0.
8
N-1~
c N,’ u,v=l-I”,
=
#
=5gLf; u {u-
u+
up/%}>.
(27)
(25)
The next question is the convergence in probability of M$J_!Z& 2;) to M+JYp, 9;). The answer is given by the following proposition.
Proposition 4. The sample estimate M$,(5$,, 2’;) tends to MzP,,(YP, L$*) in the mean square sense as N goes to infinity. Proof. In [2] the authors derived the asymptotic properties of the joint moment sample averages for a set of not necessarily stationary complex processes. These results are obtained under certain mixing conditions. However, these hypotheses imply the absolute summability of the signal cumulants that is evidently not verified in our case, see (19). To derive the consistency, the inner term of the sample estimate is first expanded as
In this expression the order of the moments has been omitted for notational simplification and b + d denotes the set obtained adding b to all the elements of d. Next, this sum is upper bounded by the sum of the absolute value of each term. The absolute value of the moment of x, being equal to a power of A, the consistency of the sample moment is determined by the convergence to zero of each term:
$ s
N2 IM,,,(L&U{U-u+y;}, u,v=l-I”,
2; u {u- u + %})\. Substituting <
$
(28)
r 4 u - u and t 5 u in (28) gives
N--lu+l,,-l
c
IM.,e(=% u {z+ T;},
’ r=-Nfh-I,,+1 9;
u (7 +
=%)>I.
(29)
A. Ferrari et al. /Signal
The noise e, being assumed iid, for IrI > 1~ - I,, the high-order moment of e, can be factored as
where c(M = (-l)“hriaM
M.,,(Yk u {r + T:}> 2;
u {r + 4pl))
= M.,e(-%, ~;)M.,&; = M.,e(% The stationarity
y:,
=$W..eW;,
(30)
21).
(31)
n $pp=
8 and hence _Ykn li”y* = 0 and -4”,* n 21 = 8. This implies that Eq. (3 1) equals 0 and that the sum in (29) 0 tends to 0 when N tends to infinity. In order to analyze the performances of an algorithm using the stationary moments, a derivation of the asymptotic covariances of the sample estimates jn+riX NE{ldM$,,(~~,
yz)12}
(32)
would be useful. For this purpose, a development analogous to the one above can be driven. However, after the substitution of the indices r and t in (27), the calculus terminates with
(33) for which a closed expression
derivation
This section is devoted to the derivation of an algorithm for the estimation of the coefficient aM of the signal (1). This estimation will rely on a matching between estimated stationary moments and their theoretical expression. Consider a root of order M: (0,12,. . . , 22~). The estimation procedure will exploit the moments associated to the first L multiples of this root. For purpose of notational simplification we use the following notation: m(r) a MzM,~JO, r12,. . , rl2M), Proposition
r = 1,. . . ,L.
(34)
3 implies
m(r) = A2M exp{jPa,+f},
r = 1,. . .,L,
lM+k.
(36)
To avoid a function minimization, the algorithm will consist of a least-squares fit between the angles of the estimated moment k(y), Y = 1,. . . ,L, and their theoretical value. According to (35) and (36) this solution is
The previous expression assumes phase unwrapping during the computation of the angles. To avoid this problem we will require that /rM XM/ < n, Y = 1,. . . , L. Thus, the condition required to avoid any ambiguity about aM is
For notational simplification, this algorithm will be referred in the sequel by SMF (stationary moments fitting). 4.2. Computer simulations
is not available.
4. Estimation algorithm
4. I. Mathematical
fi k=l
y/)
of the moment requires Yp
245
Processing 54 (1996) 239-248
(35)
In order to evaluate the performances of SMF, computer simulations using an hyperbolic-phase signal (M = 3) of N = 20 samples have been drawn. The phase coefficients are: a3 = 0.0005, a2 = 0.007, al = 0.003 and a0 = 1. The mean square error (MSE) of a^3 has been estimated using 500 noise realization for each signal to noise ratio. The first question is the determination of the parameter L. The first root of order 3, (3,3,1,1,4), has been selected and MSE curves associated to various L have been computed. If we arbitrarily assume that an average on at least N/3 samples is required in (23) the maximum value of L is given by the integer part of 2N/(31~ - 31,). Fig. 1 gives the results of these simulations: the MSE decreases as L increases. If we establish a parallel between L and the maximum correlation lag used in traditional spectral estimation algorithms [4], this result leads to the same conclusion. Consequently, the maximum value
A. Ferrari et al. 1 Signal Processing 54 (1996) 239-248
246
80-
I
I
3oO
5
10
15
I
25
30
Fig. 1. MSE for 63 using SMF. The root is (3,3,1,1,4)
of L will be retained in the following simulations. The next problem is the selection of the root. Fig. 2 represents the MSE curves for the first three roots (3,3,1,1,4), (4,5,1,2,6), (5,7,1,3,8) with, respectively, L = 4,3 and 2. The root giving the lower MSE is (4,5,1,2,6). This result illustrates the dual influence of the parameter L and the deviation 1~ - 1,. The MSE decreases when both these quantities will increase. In this sense, the pair (4,5,1,2,6) wit1 lo - I, = 5 and L = 3 represents a compromise with respect to the other pairs for which lil,~- 1, = 3, L = 4 and 1~ - 1, = 7, L = 2. Finally, the performances of SMF have been compared to the Cramer-Rao lower bound derived in [6] and to a parametric method for analysis of constant-amplitude polynomial-phase signals: the discrete polynomial-phase transform (DPPT) [5]. The first step of this method is a transformation of the signal into a single tone signal, DYM(x,,r). This is performed iteratively by A4 - 1 phase differentiations. At each step, the phase differentiation of the current signal is obtained by multiplying the sample at instant n by
35
1
40
and L = 2,3 and 4.
the conjugated sample at instant n - z. With the notation xsq = x if q is even and xsq = x* if q is odd, the authors of [5] have demonstrated that
M-l Li JJ
[_&(n
-
qT)]CyL’)
(38)
q=o =
A2M-’ expQ$} expCjM!?-‘aMn},
(39)
with $ = (M - I)!tM-‘a~_, - 0.5(M - l)M!z”‘aM. In the noisy case, the coefficient UM is then estimated from the global maximizer of the periodogram of D%4(Y,, 7). The framework derived in this paper allows a novel interpretation of the DPPT. Using a classical manipulation, the squared modulus of the discrete Fourier transform of DYM(yn, r), n = (M - l)r + 1,. . . ,N, can be rewritten as the discrete Fourier transform of its estimated autocorrelations. This shows that the DPPT retrieves aM from the estimated values of the 2”-order
A. Ferrari et al. / Signal Processing
247
54 (1996) 239-248
1OOr 95-
90-
- - (4,5,1,2,6), L=3 - -
Y
5o0
Fig. 2. MSE for 4
I
5
using SMF
According
15
these
(40) moments
1 25
30
35
40
with (3,3,1,1,4) and L = 4, (4,5,1,2,6) and L = 3, (5,‘?‘,1,3,8)and L = 2
N - (A4 - 1)~. to (39)
L=2
I
10
moments:
/sJ <
(5,7,1,3,8).
equal,
in the
noiseless case, A2” exp{ -jM!?-‘uMS} and are obviously stationary. An immediate consequence of this result is that ‘v’A4, H(2”-‘,M) always has a solution. An essential point for the effectiveness of this algorithm is the stationarity of the moments (40) in the noisy case. For a reason of simplicity we will consider only the case A4 even. From (38), the different delays occurring in these moments are {s - (2k + l)z, -2kz; k = 0,. ..,M/2 - 1) for the unconjugated samples and {s - 2kz, -(2k + 1)~; k = 0,. . , M/2 - 1} for the conjugated samples. Applying Proposition 3, stationarity will then require s # kz with lkl
- The moment corresponding to s = 0 is not stationary. However, this lag does not affect the maximiser of the periodogram and consequently it does not affect the estimation of aM. - For the other lags, stationarity requires 7 > N - (M 1)~ or equivalently 7 b N/M. Following the same reasoning that was made above for the ‘optimal’ root in the SMF, we can conclude that z must be as small as possible and consequently we can consider that its best value is N/M. This result is in perfect agreement with the optimal value obtained in [5]. Fig. 3 gives the Monte-Carlo simulations results and the Cramer-Rao lower bound. The DPPT has been implemented using an optimization algorithm initialized by the estimate given by an FFT of the transformed signal zero-padded to 29. The parameter 7 equals N/IV. For SNR > 10 dB, performances of the proposed algorithm are very close to the CRLB and performances of DPPT correspond to the values predicted by the theoretical analysis derived in [5]. For SNR < 10 dB, DPPT exhibits the traditional breakdown, whereas the efficiency of SMF remains almost constant. This result could be explained by the fact
A. Ferrari et al. 1 Signal Processing 54 (1996) 239-248
248
90-
80P E
/
,.I’
50I’ /’ /’ 40
0
I
I
I
5
10
15
I
I
SN&3,
25
Fig. 3. MSE for a^3 using SMF, DPPT and Cramer-Rao
that the moments order involved are, respectively, 2M for DPPT and 2A4 for SMF, these last ones being estimated with better precision.
[3] [4]
5. Conclusion A new framework for analysis of polynomial-phase signals is derived. It is based on shift-invariance properties of certain signal high-order moments. Simulations have been presented to demonstrate the effectiveness of the approach.
References [l] D.R. Brillinger, Time Series: Data Analysis and Theory, Holden Day, San Francisco, CA, 198 1. [2] A.V. Dandawate and G.B. Giannakis, “Asymptotic theory of mixed time averages and &h-order cyclic-moment and
[5]
[6]
[7]
30
35
40
lower bound.
cumulant statistics”, IEEE Trans. Inform. Theory, Vol. 41, No. 1, January 1995, pp. 216-232. F.R. Gantmacher, Theory of Matrices, Vol. 1, Chelsea Publishing Co., New York, 1959. S.M. Kay, Modern Spectral Estimation, Signal Processing Series, Prentice-Hall, Englewood Cliffs, NJ, 1987. S. Peleg and B. Friedlander, “The discrete polynomial-phase transform”, IEEE Trans. Signal Process., Vol. 43, No. 8, August 1995, pp. 1901-1912. S. Peleg and B. Porat, “The Cramer-Rao lower bound for signals with constant amplitude and polynomial phase”, IEEE Trans. Signal Process., Vol. 39, No. 3, March 1991, pp. 749-752. P. Tichavsky and A. Swamy, “Statistical characterization of sample fourth-order cumulants of a noisy complex sinusoidal process”, IEEE Trans. Signal Process., Vol. 43, No. 7, July 1995, pp. 1620-1630.
54 (1996) 239-248
Polynomial-phase signal analysis using stationary moments A. Ferrari*, C. Theys, G. Alengrin Lahoratoiw 13s. CNRS, UniversitP de Nice Sophia-Antipolis, 41 Bd NupolPon III, 06041 Nice Cedes. Frunw Received
22 December
1995: revised
13 June 1996
Abstract This article proposes a method for polynomial-phase signals analysis relying on their time-invariant high-order moments. The set of these moments is first characterized in the noiseless and noisy case. It is demonstrated that the only identifiable phase parameter is the highest-order coefficient, the estimation requiring moments of order at least the double of the phase degree. Next, the consistency of the estimated moments calculated by sample averages is proved. Finally, an estimation algorithm based on the moments having for lags the multiples of a given ‘root’ is developed. Performances, evaluated by computer simulations, show the influence of the different parameters and the efficiency of the method. Zusammenfassung Der Beitrag schllgt eine Methode vor, um Signale polynomialer Phase unter Beriicksichtigung ihrer zeitinvarianten Momente hoherer Ordnung zu analysieren. Die Menge dieser Momente wird zunachst Rir den rauschfreien und den verrauschten Fall charakterisiert. Es wird gezeigt, dab der einzige identifizierbare Phasen-Parameter der Koeffizient hochster Ordnung ist, und die Schlzung Momente erfordert, deren Ordnung mindestens das doppelte des Phasengrades ist. Danach wird die ijbereinstimmung der geschatzten Momente durch Mittelung bewiesen. Abschliegend wird ein Schatzalgorithmus auf der Basis von Momenten, ausgehend von einem Vielfachen einer gegebenen ‘root’, entwickelt. Das Verhalten, durch RechnerSimulationen bewertet, zeigt den EinfluS der verschiedenen Parameter und die Leistungsfahigkeit der Methode. RCsumC Cet article propose une methode d’analyse de signaux a phase polynomiale utilisant leurs moments d’ordre superieurs qui sont invariants dans le temps. L’ensemble de ces moments est tout d’abord caracterise dans le cas non bruit6 et bruit& On demontre en particulier que seul le coefficient du terme de plus haut degre est identifiable et ce uniquement a partir de moments dont l’ordre est superieur ou tgal au double du degre de la phase. La consistance des estimees de ces moments, calculees par moyenne temporelle est ensuite proude. Enfin, un algorithme d’estimation utilisant les moments dont les retards sont les multiples d’une ‘racine’ donnee est developpe. Les performances Cvaluees par simulations montrent l’influence des differents paramttres ainsi que I’efficacitt de cette mtthode.
K~Jww~‘s: Polynomial-phase
* Corresponding
signals; High-order
author. Tel.: +33-93217953;
moments;
fax: +33-93212054;
Parameter estimation
e-mail: [email protected].
0165-1684/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved PZI: SO165-1684(96)001 IO-7
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
240
1. Introduction This article addresses the estimation of the parameters of a polynomial-phase signal. This problem, encountered in different radar systems, is usually solved using a time-frequency analysis or phase-only algorithms. For a detailed introduction on this topic see, for example, [5]. This paper investigates a parametric solution based on statistical properties of this signal: although it is clearly non-stationary, some of its high-order moments are shift-invariant. This property is exploited for the retrieval of the phase parameters. The paper is organized as follows. In Section 1, the conditions verified by the delays of these ‘stationary’ moments are derived in the noiseless and noisy case. Section 2 is devoted to the statistical study of the estimated moments. In particular, consistency of the sample estimates is proved. Finally, in Section 3 an estimation algorithm is derived and its performances analyzed by computer simulations.
Proposition 1. The signal x, in (1) is (II,. . . , ZQ)stationary if and only if the following condition, referred as H( p, M), holds: Wp,M):
Vq E
(11 - Ii,,)
f:
{l,...,
M - 1},
= 0.
(3)
k=l
In this case the stationary moments equal M2p,xU1,.
. . , 12p)
= A2P exp
j aM f:
(1;1”
(4)
k=l
Proof. The 2pth-order moments equal A2P exp {j$} with
4=
2 9,aq((n -k k=l
lk)’
M&n;
11,.. . , Z2p)
- (n + lk+p)‘)
q=l
2. High-order moments of polynomial-phase signals 2.1.
The noiseless case
We define a noisy polynomial-phase
signal y, as
yn =G+G,
(1)
where a0 is a random variable uniformly distributed in the interval [0,2rc) and e, is a complex-independent identically distributed (iid) noise. The 2pth-order moment of the signal z, is defined as
M2,,&;
11,. . . ,
(2) The signal z, is said (11,. . . , 12P)-stationary if Mz,,,(n; 11,.. . ,1+) is shift-invariant, i.e., M,,,(n; 11, . . . ,lzp) is not a function of it. In this case, the index n will be omitted in M,,(n; 11,.. . ,lzp). It is worthy to note that the definitions obtained with a different number of conjugated and unconjugated terms lead to a 2pth-order moment of x, that equals zero.
This polynomial of the variable n is not a function of n if and only if ‘dm = 1,. . . , M, b, = 0. bl = 0 implies CE, aq(‘f)ckp,l (Ii-’ - I:;;) = 0. This equation has to be verified Va,. This implies H(p,M). Conversely, it can be easily checked that H(p,M) implies Vm = 1,. . . ,M, b, = 0. The coefficient bo equals Ego a4 Ckp,l (1: - lz,,). Application of H(p,M) to this equality terminates the proof of Eq. (4). 0 This result suggests the following remarks: - The only identifiable parameter of a polynomialphase signal from its stationary moments is the coefficient of highest degree, aM. The estimation of all the ak can be performed using the iterative
241
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
algorithm consisting in: estimate UM, multiply the signal by exp {-j aM@‘}, estimate CIM-_I, . . . . ~ H(p,M) is a set of M - 1 non-linear equations for 2p integer variables. A trivial solution consists in choosing for (l,+t , . . . , l+,) a permutation of (11,. . . , lP). Unfortunately, this solution, according to (4), leads to a moment independent of a,+,,. An important point is then to find non-trivial solutions of H(p,M). An analytic solution of the diophantine system H( p,M) is well beyond the scope of this paper. The solution consisting in the use of nested loops together with constraints on the delays to reduce the computation will be discussed further. To illustrate these results, we apply Proposition 1 to the cases M = 2,M = 3, assuming without loss of generality that 11 = 0. Example 1 (x, is u chirp: M = 2). - lf p = 1,12 must verify: 12 = 0. The signal is only (0,0)-stationary and Mz,(O, 0) = A*. _ If p = 2, (12,1x, 14) must verify 12- 13- 14 = 0. The signal is (0,13+14,13, /4)-stationary and M+JO, 13+ 14,139 14)
= A4
exP{j
M2,wxQA
12,.
> 12~4)
(-l)“jMQM
fi
lM+k
(5)
k=l
Proof. Consider the two p x p diagonal matrices C = Diag{lr, . . , Ip}, D = Diag{ I,+, , . , lzp}. The condition H( p, M) can be rewritten in terms of these matrices as H(p,M):
Vq E {I,...,
M-
l},
tr{Cq} = tr{Dq},
(6)
where tr{M} is the trace of the matrix M. If we apply the Hamilton-Cayley theorem, [3], to these matrices, we obtain cp - c,cp-I
- c*cp-2
- .” - cp-,c
- cpI = 0, (7)
DP - dlDP-’
- d2DP-2 -
. - d,_,D
- dPI = 0,
2Q21314)).
(8)
Example 2 (x, is an hyperbolic phase signal: M=3). - If p = 1, 12 must verify: 12 = 0 and 1: = 0. The signal is (O,O)-stationary and Mz,JO, 0) = A’. - Ifp=2, (12,13,14)mustverify12-ls-&=Oand 1: - 13- 1: = 0. These two equations imply 1314 = 0. The signal is (0,14,0,14) and (0,13,ls,O)-stationary and M+,(O, LO, 14) = &JO, k,k 0) = A4. - If p = 3, (12,13,14,15,1~) must verify 12+ 13 - 1415 - 16 = 0 and 1; + 1: - 1: - 1: - 1: = 0. An analytical solution to H(3,3) is not available; however nontrivial solutions exist, see Table 1. Using the two previous equalities, it is possible to show that the six order stationary moments of x, equals A6 exp{-j 3u3141s16). These results proposition.
of order 2M of x, equal
are generalized
by the following
Proposition 2. The coeficient aM of the signal in ( 1) is only identifiable from the stationary moments of order higher or equal to 2M. The stationary moments
where ck and dk are the coefficients of the characteristic polynomial of C and D. The Faddeev method [3] allows the recursive computation of these coefficients using the relations kCk = tr{Ck} - cl tr{c”-‘}
- . . - ck-1 tr{C},
k = l,...,p, kdk = tr{Dk} - dl tr{Dk-‘}
(9) - . ‘. - dk_1 tr{D),
k = l,...,p,
(10)
with CO= 1 and do = 1. The recursive H( p, M) to these equations gives VqE{l,...,
min(M-l,p)},
cq=dq.
application
of
(11)
Consider first the case p < M. The trace of the difference of the two equalities (7) and (8) together with (11) and N(p,M) gives tr{Cp} =tr{DP}. If we now multiply Eq. (7) by C and Eq. (8) by D, the trace of the difference of these two equalities will give tr{C P+* } - tr{DJ’+‘}. Iterating, we then obtain
242
A. Ferrari et al. / Signal Processing 54 (1996) 239-248
‘dq, tr{Cq} = tr{Dq}. The application of this result with q=M to (4) terminates the proof of the first assertion of the proposition. Consider now the case p = M. In this case, the trace of the difference of the two equalities (7) and (8) together with ( 11) and H(M, M) gives tr{C”}
- tr{oM}
(A)~-’
=
(12)
fi lk - (-l)“-ifi k=l
lp+k
0,13)M2,e(l3
+M,,x(n;
l3 +
+Mz,x(~
0,14)M2,,(13
+Mz,x(n;
13 +
+
l4914)
14, ~4)M2,e(o, +
13) 14,l3)
14,~3)M2,e(o,
(15)
14).
x, being only (O,O)-stationary, stationarity of this moment, independent of e,, requires 13 = 14 = 0. This solution being useless, stationarity will be obtained using the independence of the noise samples. In fact, the constraint 13 # 0 and 14 # 0 will imply the nullity of all the noise second-order moments:
k=l
(13) Application the proof.
+M,,x(n;
v’I3
of this result with Ii = 0 to (4) terminates 0
#
0914
M4,Jo,
#
0,
13 +
14,13,14)
=
M4,x(o,l3
+
14,139 14).
(16)
The noisy case
2.2.
In the noisy case, we will consider the moments of the signal yn at lags for which x, is (11,. . . , IQ,)stationary. For brevity, we introduce the ordered k-uple 2, constituted by a subset of k elements of (I 1,. . . , lp). 9: is a k-uple of (I,+,, . . ., lQ),?& the complementary of yk in (Ii, . . . , lp) and 3: the complementary of 9: in (Z,+i, . . . , 1~). Consequently, we will use the condensed notation M.&n; &, 2:) for the 2kth-order moment of z,,. Using the circularity of e, and the hypothesis on a0 distribution, we can expand the 2pth-order moment of yn as M2p,y(%
2;)
Unfortunately, the (Ii,. . . , IQ,)-stationarity of xB does not imply the (_&, _%‘f )-stationarity of x,, and then the (Ii , . . . , 12P)-stationarity of yn. To illustrate this difficulty we consider the following example.
Example 3 (x, is a chirp: M = 2, and p = 2). If we expand M+y at the stationary Example 1, we obtain M4,@, =
13 + M4,x(0,/3
14,139 +
lags of M+,
+
M4,e(0,13
+
to the general
Proposition 3. For all solutions of H( p, M) verifying (1 I,..., (II,...
lp} n {4,+1,.. , 12p)=
MzM,~(O,
M2,,x(h
., hp} = 0 , . . . , 12p).
we have MQ,,~ Moreover, if p =M,
12,. . . ,12~)
= A2M eXp
1)“j
(-
hfLZ,+f fi
l,,,f+k
.
(17)
k=l
It is interesting by now to investigate the possibility of using cumulants instead of moments in order to solve at the same time the additive noise problem and the multicomponent signal case. The properties of cumulants, [l], imply that for xn defined in (l),
A = cum(x,+h,.
=
~,.M,x(n;
. . ,JW,,~,*+~,+,,
yp,
2;
),
. . . ,~n*,~~,)
(18)
(19)
see
l4) 14,137 l4)
This result can be easily extended case using (14).
I49 13914)
with y,, = cum( ej@, . . . , ejQ, e-jao, . . . , e--jao). Assume that y,, is formed by m noisy polynomialphase signals of order M, xi, with independent initial phases. The additive property of the cumulants of the
A. Ferrari et al. /Signal
sum of independent
processes implies
If each component is (11, . . . , IQ))-stationary, (19) implies that C&,,_,,(n; LZ”,2; ) is not a function of n and
Meanwhile,
this solution is not tractable because is defined from the M2k,Jn; P”,,9:) mentioned above, are not necessary
C,.,(P& 2:) which, as stationary. 2.3.
The set of stationary moment delays
=
0.
In the case of a parabolic-phase signal an explicit solution exists, see Eq. (16), but, for higher degrees,
243
analytical solutions are not available. We suggest to find them using 2M - 1 loops to test condition H(A4,M). In this context, it is important in order to reduce the computational cost of the scan to define a set of minimal size containing the lk. This is done in a way similar to [7]. The moment is simultaneously symmetric with respect to the first and the last A4 parameters lk. The smallest lk can then be moved to position 1 or position M+l;infactifmin{lr,...,l~}=min{l~+~,...,/~~} the two sets are not disjoints. In this last case, the obvious relation M~MJJ&, 2;) = MzM,J_Y~, Y,)* allows to move the smallest lk to position 1. From the definition of stationarity, if (11,12,. . . ,L~M) verifies H(M,M),Vr,(r + ll,r + 12,. . . , r + 12~) also verifies H(M,M). According to this, the first and smallest delay can be subtracted to the others leading to delays that are all nonnegative for the first M and positive for the last M. Both these sets can then be reordered in increasing order giving
1, = 0,
We discuss in this section the computation of the delays lk leading to a stationary moment of y,. We will limit our study to the minimal order moment p = M. According to the previous sections, we need to find the 2M-uple verifying H(M,M), { 11,. . . , 1~) fl {1A4+1,...,12M}
Processing 54 (1996) 239-248
0<1/2< ... 6lM,
161&f+, < ...
dlZM. (22)
An important remark is that if (Ii, 12,. . ,Z~M) verifies H(M,M) and {It,., ., 1~) n {l~+l,..., 12~) = 0,Vr # O,{rll,rl2 ,..., rZZ*} also verifies these two conditions. A direct consequence is that the set of interest can be reduced to only its coprime members. An element of this set will be called a root of order M. Table 1 gives the 11 first roots of order 3 and 4.
Table I Eleven first roots of order 3 and 4 12
13
14
3 4 5 5 6 6 7 7 7 9 9
3 5 7 8 9 11 7 8 II 9 10
I I I 2 I 2 I 2 I 3 I
15
16
-3141510
12
13
14
1s
16
17
lx
4w6h18
I
4 6 8 9 10 I2 9 IO 12 12 12
-12 -36 -7 -108 -120 -216 -108 -180 -180 -324 -216
3 4 5 5 5 6 6 7 7 8 8
4 7 5 10 12 7 I3 9 11 9 II
7 II 10 I5 17 I3 19 16 18 17 I9
1 1 I I 2 1 1 I 2 2 I
1 2 2 3 2 3 4 4 3 3 5
6 9 8 12 15 10 I.5 I2 I5 14 I4
6 10 9 I4 I5 12 18 15 I6 I5 I8
144 720 516 2016 3600 1440 4320 2880 5760 5040 5040
2 3 2 4 3 4 3 5 3 6
A. Ferrari et al. /Signal Processing 54 (1996) 239-248
244
3. Statistical properties of the sample estimates The stationary moments of the signal under scope may be estimated on the basis of a single realization of the time series y,, n = 1,. . . , N. The sample moment function can be estimated as
where ZM = max{_YP U 9;}, I, = min{ZP U 3.) and N, = N - 1~ + 1,. We address in this section the statistical properties of these estimates. For this sake, let
(26) The term of this sum corresponding to the empty sets _5&=0and A?:=0 equals M2,&Zpp,2;). Hence, dM$,,( TP, 9’) equals the sample mean of (26) where the sum ranges only over the non-empty subsets LYEu 9: # 0. Consequently, the expansion of the squared modulus of this sum provides a first expression of ~dM$JYP,Ip,*)~2. Using the notation (2), the expectation of this expansion reduces to the sum over the subsets _%,,49,* and Y,,_??z of -
1
with _CZkUL?z # 0 and LZlUYz
M$,,J-%2; I- ~wG% 3; ).
M,,e(u;2k
u {v
-
ZJ +
y;},
(24)
It is straightforward that the sample moment function (23) is an unbiased estimator of MQ,,J 5$,, 9; ), hence
E{dM&,(~p,2; )I = 0.
8
N-1~
c N,’ u,v=l-I”,
=
#
=5gLf; u {u-
u+
up/%}>.
(27)
(25)
The next question is the convergence in probability of M$J_!Z& 2;) to M+JYp, 9;). The answer is given by the following proposition.
Proposition 4. The sample estimate M$,(5$,, 2’;) tends to MzP,,(YP, L$*) in the mean square sense as N goes to infinity. Proof. In [2] the authors derived the asymptotic properties of the joint moment sample averages for a set of not necessarily stationary complex processes. These results are obtained under certain mixing conditions. However, these hypotheses imply the absolute summability of the signal cumulants that is evidently not verified in our case, see (19). To derive the consistency, the inner term of the sample estimate is first expanded as
In this expression the order of the moments has been omitted for notational simplification and b + d denotes the set obtained adding b to all the elements of d. Next, this sum is upper bounded by the sum of the absolute value of each term. The absolute value of the moment of x, being equal to a power of A, the consistency of the sample moment is determined by the convergence to zero of each term:
$ s
N2 IM,,,(L&U{U-u+y;}, u,v=l-I”,
2; u {u- u + %})\. Substituting <
$
(28)
r 4 u - u and t 5 u in (28) gives
N--lu+l,,-l
c
IM.,e(=% u {z+ T;},
’ r=-Nfh-I,,+1 9;
u (7 +
=%)>I.
(29)
A. Ferrari et al. /Signal
The noise e, being assumed iid, for IrI > 1~ - I,, the high-order moment of e, can be factored as
where c(M = (-l)“hriaM
M.,,(Yk u {r + T:}> 2;
u {r + 4pl))
= M.,e(-%, ~;)M.,&; = M.,e(% The stationarity
y:,
=$W..eW;,
(30)
21).
(31)
n $pp=
8 and hence _Ykn li”y* = 0 and -4”,* n 21 = 8. This implies that Eq. (3 1) equals 0 and that the sum in (29) 0 tends to 0 when N tends to infinity. In order to analyze the performances of an algorithm using the stationary moments, a derivation of the asymptotic covariances of the sample estimates jn+riX NE{ldM$,,(~~,
yz)12}
(32)
would be useful. For this purpose, a development analogous to the one above can be driven. However, after the substitution of the indices r and t in (27), the calculus terminates with
(33) for which a closed expression
derivation
This section is devoted to the derivation of an algorithm for the estimation of the coefficient aM of the signal (1). This estimation will rely on a matching between estimated stationary moments and their theoretical expression. Consider a root of order M: (0,12,. . . , 22~). The estimation procedure will exploit the moments associated to the first L multiples of this root. For purpose of notational simplification we use the following notation: m(r) a MzM,~JO, r12,. . , rl2M), Proposition
r = 1,. . . ,L.
(34)
3 implies
m(r) = A2M exp{jPa,+f},
r = 1,. . .,L,
lM+k.
(36)
To avoid a function minimization, the algorithm will consist of a least-squares fit between the angles of the estimated moment k(y), Y = 1,. . . ,L, and their theoretical value. According to (35) and (36) this solution is
The previous expression assumes phase unwrapping during the computation of the angles. To avoid this problem we will require that /rM XM/ < n, Y = 1,. . . , L. Thus, the condition required to avoid any ambiguity about aM is
For notational simplification, this algorithm will be referred in the sequel by SMF (stationary moments fitting). 4.2. Computer simulations
is not available.
4. Estimation algorithm
4. I. Mathematical
fi k=l
y/)
of the moment requires Yp
245
Processing 54 (1996) 239-248
(35)
In order to evaluate the performances of SMF, computer simulations using an hyperbolic-phase signal (M = 3) of N = 20 samples have been drawn. The phase coefficients are: a3 = 0.0005, a2 = 0.007, al = 0.003 and a0 = 1. The mean square error (MSE) of a^3 has been estimated using 500 noise realization for each signal to noise ratio. The first question is the determination of the parameter L. The first root of order 3, (3,3,1,1,4), has been selected and MSE curves associated to various L have been computed. If we arbitrarily assume that an average on at least N/3 samples is required in (23) the maximum value of L is given by the integer part of 2N/(31~ - 31,). Fig. 1 gives the results of these simulations: the MSE decreases as L increases. If we establish a parallel between L and the maximum correlation lag used in traditional spectral estimation algorithms [4], this result leads to the same conclusion. Consequently, the maximum value
A. Ferrari et al. 1 Signal Processing 54 (1996) 239-248
246
80-
I
I
3oO
5
10
15
I
25
30
Fig. 1. MSE for 63 using SMF. The root is (3,3,1,1,4)
of L will be retained in the following simulations. The next problem is the selection of the root. Fig. 2 represents the MSE curves for the first three roots (3,3,1,1,4), (4,5,1,2,6), (5,7,1,3,8) with, respectively, L = 4,3 and 2. The root giving the lower MSE is (4,5,1,2,6). This result illustrates the dual influence of the parameter L and the deviation 1~ - 1,. The MSE decreases when both these quantities will increase. In this sense, the pair (4,5,1,2,6) wit1 lo - I, = 5 and L = 3 represents a compromise with respect to the other pairs for which lil,~- 1, = 3, L = 4 and 1~ - 1, = 7, L = 2. Finally, the performances of SMF have been compared to the Cramer-Rao lower bound derived in [6] and to a parametric method for analysis of constant-amplitude polynomial-phase signals: the discrete polynomial-phase transform (DPPT) [5]. The first step of this method is a transformation of the signal into a single tone signal, DYM(x,,r). This is performed iteratively by A4 - 1 phase differentiations. At each step, the phase differentiation of the current signal is obtained by multiplying the sample at instant n by
35
1
40
and L = 2,3 and 4.
the conjugated sample at instant n - z. With the notation xsq = x if q is even and xsq = x* if q is odd, the authors of [5] have demonstrated that
M-l Li JJ
[_&(n
-
qT)]CyL’)
(38)
q=o =
A2M-’ expQ$} expCjM!?-‘aMn},
(39)
with $ = (M - I)!tM-‘a~_, - 0.5(M - l)M!z”‘aM. In the noisy case, the coefficient UM is then estimated from the global maximizer of the periodogram of D%4(Y,, 7). The framework derived in this paper allows a novel interpretation of the DPPT. Using a classical manipulation, the squared modulus of the discrete Fourier transform of DYM(yn, r), n = (M - l)r + 1,. . . ,N, can be rewritten as the discrete Fourier transform of its estimated autocorrelations. This shows that the DPPT retrieves aM from the estimated values of the 2”-order
A. Ferrari et al. / Signal Processing
247
54 (1996) 239-248
1OOr 95-
90-
- - (4,5,1,2,6), L=3 - -
Y
5o0
Fig. 2. MSE for 4
I
5
using SMF
According
15
these
(40) moments
1 25
30
35
40
with (3,3,1,1,4) and L = 4, (4,5,1,2,6) and L = 3, (5,‘?‘,1,3,8)and L = 2
N - (A4 - 1)~. to (39)
L=2
I
10
moments:
/sJ <
(5,7,1,3,8).
equal,
in the
noiseless case, A2” exp{ -jM!?-‘uMS} and are obviously stationary. An immediate consequence of this result is that ‘v’A4, H(2”-‘,M) always has a solution. An essential point for the effectiveness of this algorithm is the stationarity of the moments (40) in the noisy case. For a reason of simplicity we will consider only the case A4 even. From (38), the different delays occurring in these moments are {s - (2k + l)z, -2kz; k = 0,. ..,M/2 - 1) for the unconjugated samples and {s - 2kz, -(2k + 1)~; k = 0,. . , M/2 - 1} for the conjugated samples. Applying Proposition 3, stationarity will then require s # kz with lkl
- The moment corresponding to s = 0 is not stationary. However, this lag does not affect the maximiser of the periodogram and consequently it does not affect the estimation of aM. - For the other lags, stationarity requires 7 > N - (M 1)~ or equivalently 7 b N/M. Following the same reasoning that was made above for the ‘optimal’ root in the SMF, we can conclude that z must be as small as possible and consequently we can consider that its best value is N/M. This result is in perfect agreement with the optimal value obtained in [5]. Fig. 3 gives the Monte-Carlo simulations results and the Cramer-Rao lower bound. The DPPT has been implemented using an optimization algorithm initialized by the estimate given by an FFT of the transformed signal zero-padded to 29. The parameter 7 equals N/IV. For SNR > 10 dB, performances of the proposed algorithm are very close to the CRLB and performances of DPPT correspond to the values predicted by the theoretical analysis derived in [5]. For SNR < 10 dB, DPPT exhibits the traditional breakdown, whereas the efficiency of SMF remains almost constant. This result could be explained by the fact
A. Ferrari et al. 1 Signal Processing 54 (1996) 239-248
248
90-
80P E
/
,.I’
50I’ /’ /’ 40
0
I
I
I
5
10
15
I
I
SN&3,
25
Fig. 3. MSE for a^3 using SMF, DPPT and Cramer-Rao
that the moments order involved are, respectively, 2M for DPPT and 2A4 for SMF, these last ones being estimated with better precision.
[3] [4]
5. Conclusion A new framework for analysis of polynomial-phase signals is derived. It is based on shift-invariance properties of certain signal high-order moments. Simulations have been presented to demonstrate the effectiveness of the approach.
References [l] D.R. Brillinger, Time Series: Data Analysis and Theory, Holden Day, San Francisco, CA, 198 1. [2] A.V. Dandawate and G.B. Giannakis, “Asymptotic theory of mixed time averages and &h-order cyclic-moment and
[5]
[6]
[7]
30
35
40
lower bound.
cumulant statistics”, IEEE Trans. Inform. Theory, Vol. 41, No. 1, January 1995, pp. 216-232. F.R. Gantmacher, Theory of Matrices, Vol. 1, Chelsea Publishing Co., New York, 1959. S.M. Kay, Modern Spectral Estimation, Signal Processing Series, Prentice-Hall, Englewood Cliffs, NJ, 1987. S. Peleg and B. Friedlander, “The discrete polynomial-phase transform”, IEEE Trans. Signal Process., Vol. 43, No. 8, August 1995, pp. 1901-1912. S. Peleg and B. Porat, “The Cramer-Rao lower bound for signals with constant amplitude and polynomial phase”, IEEE Trans. Signal Process., Vol. 39, No. 3, March 1991, pp. 749-752. P. Tichavsky and A. Swamy, “Statistical characterization of sample fourth-order cumulants of a noisy complex sinusoidal process”, IEEE Trans. Signal Process., Vol. 43, No. 7, July 1995, pp. 1620-1630.