Cyclic Spectral Features Extracting of Complex Modulation Signal Based on ACP Method

Published Online June 2011 in MECS

Complex modulation technique is widely used in low probability of intercept wideband signal [1], which makes Non-cooperative intercept receivers difficulty in looking for these signals. The conventional signal processing methods such as energy detection and power spectral density analysis cannot provide details about the complex modulation parameters. Recent research [2] have shown that intercept receiver can increase its processing gain by implementing cyclo-stationary signal processing method so as to raise the effect of features extraction.

This paper investigates an improved cyclic spectral feature extracting method based on averaged cyclic periodogram The key parameters (bandwidth B , code length Nc , code rate Rc , codes period tm , modulation period T and carrier frequency fc etc.) of complex modulated signals are estimated on the condition of unknown prior knowledge. The detail procedures of the method are presented and the extracting effects on the frank coded and costas coded signals are evaluated, since the two typical complex modulated signals have the firstand second-order statistic periodicity [9].

• II.  M ODELING OF S IGNAL P ARAMETERS

The complex modulated signal intercepted by noncooperate receiver can be written as

A t) = s( t) + < t)                             (1)

where n ( t ) is addictive white Gaussian noise (AWGN), s ( t ) is complex modulated signal. Two typical signals are modeling as below.

• A. Frank Coded Signal

Frank coded signal is one type of polyphase-coded signal which can be written as[4]

s ( t ) = A exp{ / 2 p f c t + f t ) + q_c\}             (2)

where A is the margin of s ( t ) , f is carrier frequency, q _c is initial phase. We start from a pulse of duration T .The pulse is divided into N bits of identical duration t_b = T / N ,the code rate is R ^ = 1 / T ,and each bit is assigned with a different phase value a_t , k = 1,2, L , N .The phase envelope is

N ft) = ^ a/ect(/- ktb) ,which can be presented in k=1

Table.1

a Frank code is used with M sub-codes, a carrier frequency of f = 1000 Hz and sampling frequency of f = 7 kHz with SNR=0dB [10]. The Frank code is a polyphase code with each sub-code phase defined as fu=2P (i -1,(/'1f                       (3)

where i = 1,2, L M , j = 1,2, L M , N_c = M . The number of codes is M , total code length is N_c and code period is t_m [11]. The ambiguity function of Frank coded signal can be written as

m = -•

| g„ (t) | Wh (t, n +

B. Costas Coded Signal

Costas coded signal is divided into N bits of identical duration b = T / N , the frequency code sequence is a = { a 1 , a , , L , а м } ,the code length is M . The complex envelop of Costas coded signal can be described as[6]

d t) =

where

1 M , ^ a um[^-

( m - 1)/ b ]

g„ (t) = 7 oT%(t)^tt — t) ejpnt/T»(5)

•

Wh(t, n) =o   wh(t)w*(t- t)e j2pntdt(6)

- • and

X t) = Wh (t) ◊ %(t)

% where %(t) is periodic extend of x(t) , wh (t) is the

where

um (t)

and

= I exP(T2.pvmt), = t0

0 £ t £ tb elsewhere

a

window function with length T .Fig.1 shows the ambiguity function of 16 elements Frank coded signal.

Figure 1. Ambiguity function of 16 elements Frank coded signal

The AF of Frank coded signal can be written as

1 M

c^ t, n) =    a exp[j2,p(m - 1)mb{F mm(t, m)

^M m = 1

M

+ a  Fmn[ t - (m - n)b, mi}

where

^F mn ⁽ t , m = ^{(1 -} ^ )^^^ ^ ^exP(-. j b ^- j 2 P f n t ) >| t | ^£ t b t_b     a

Figure 2. Zero-Doppler cut of 16 elements Frank coded signal AF

Fig.3 shows the ambiguity function of 7 elements Costas coded signal.

Figure 3. Ambiguity function of 7 elements Costas coded signal

The zero-doppler cut AF of Frank coded signal is shown in Fig.2, which is also the autocorrelation function of Frank coded signal. The correlation function is zero for displacements of multiples of M since the rows of the Frank matrix are orthogonal. The autocorrelation has a magnitude of unity for displacement of one more or one less than multiples of M . The Frank code exhibits good properties of relatively low sidelobes.

Figure 4. Zero-Doppler cut of 7 elements Costas coded signal AF

The zero-doppler cut AF of Costas coded signal is shown in Fig.4, which is also the autocorrelation function of Costas coded signal. We can conclude the first null is located at T / M = t_fc / M = t_fc / 7 .This implies that the pulse compression of a Costas signal is M . To reduce the sidelobe pedestal, M will have to be increased. Besides, the value M / t_b is approximately the spectral width of the Costas signal[6].

We can summarize both Frank and Costas coded signals’ parameters and the interrelationship among them like:

b = 1/ b

T = N ◊ b                           (15)

Rc = 1 / T in Frank coded signal and

B = M / t b, T = N ◊ b R c = 1 / t

• III.    C YCLIC S PECTRAL F EATURES BASE ON ACP M ETHOD

• A.    A General Quadratic Form of ACP Method

Given an finite length signal sequence { Xn ]} L = 0 , any non-parametric estimator $ 2 x ( f , a, L ) of the cyclic spectrum S 2 j X f , a can be deduced from the general quadratic form[7].

$2x(f, a,L) = Da L- p = 0

L - 1

1a ql [p qXp q = 0

X

X[qi * e" j2P f+ ba)pD e^P f - ba) qD where QL is a suitably chosen positive semi-definite kernel such as to preserve the interpretation of $2 x( f, a, L) as a power density and in particular $2x( f a, L) 3 0 for a = 0 .The double Discrete Time

Fourier Transform (DTFT) of kernel Q_L which can be used as y(l , h)

L - 1 L - 1

y L (l , h) = D ² a a Q l [ P q i e jpl ⁽ p ^- q ^{) D} e j ² ph q D (18) p= 0 q = 0

Different kernel Q_L result in different spectral estimators. Let { w [ n |} N w 0¹ be a positive and smooth N -long data-windows and let w^n | = w [ n - kR be its shifted version by R samples. Then the double DTFT of kernel Q_L is

L - 1 L - 1

y(l , h) = D ² a a Q_l [ P q i e j ^{2 pl ( p- q ) D} e -j ² phqD p= 0 q = 0

= W( l) W(l - h) * DR (h) e-jphRD1^ K-1’ with    W( f)   the DTFT of   w[ni   and

D R D у) = Islnl p/ R D L ) .

K K sinpl R D)

The averaged cyclic periodogram is obtained by introducing (19) into (17)

K

$ 2 x ( f, a , L) = -L- a X^ k ( f + ba ) X ( k ( f - ba )* (20) K D _t = 1 N w         N w

R + N „ - 1

with   X(k (f) =D a  wk[nШnie j2pfnD   the w        n = R short-time DTFT of the kth weighted sequence { w^[n^ni} R+ Nw-1                             and

K = |( L - N w ) / Ri + 1 .Formula(20) can be very efficiently implemented since the parallel two liner FFT processing faster than double DTFT in (17).Imposing N_w to be a power of 2, the ACP method is more efficient. In practice application, we can use the modulo[8] cyclic spectrum density function (or the modulus square) after the method of integration along the frequency axis, that is defined by spectral correlation model function of the marginal distribution cycle frequency.

P s ( a, L ) =q⁺ ^ I $ 2 x ( f , a, L ) | df              (21)

This reduces the frequency of the same cycle phase between different frequency spectral of random phase cancellation, concentration energy, highlighting the relative intensity of spectral lines, improving the detection probability and estimation accuracy.

• B.    Cyclic Spectral Features of Frank Coded Signal

The ACP method estimation results of a Frank coded signal are shown in Fig.5 and Fig.6, which represent signal parameters in bi-frequency plane respectively. Fig.1 illustrates the measurement of fc andB ; Fig.2 illustrates the measurement of fc and Rc [14]. In Fig.5, we can get B by measuring the distance of coordinates of peak value on frequency axis, namely B = 1000 Hz, and get carrier frequency by projecting the peak value of support area on cyclic frequency axis, namely fc = Fc / 2 ; In Fig.6, we can get Rc by measuring the distance between main and sub- support areas on cyclic frequency axis. We can calculated more parameters in the basis of B , Rc and fc ,

Cycle frequency (Hz)

Figure 5. Cyclic spectral features of Frank coded signal in main support area

Figure 6. Cyclic spectral features of Frank coded signal in sub-support area

• C.    Cyclic Spectral Features of Costas Coded Signal

The ACP method estimation results of a Costas coded signal are shown in Fig.7 and Fig.8, which represent signal parameters in bi-frequency plane respectively. Fig.7 illustrates the measurement of B ; Fig.8 illustrates the measurement of f_m and R_c [14]. In Fig.7, we can get B by measuring the distance of coordinates of peak value on frequency axis, In Fig.8 we can get R_c by measuring the distance between main and sub- support areas on cyclic frequency axis. and get frequency modulation sequence by projecting the peak value of support area on cyclic frequency axis; We can see the cross terms[9] exists between frequency 4.7KHz and 5.6KHz, which mean the two frequency in the frequency modulation sequence are adjacent, but we can’t get the order of the frequency sequence. We can calculated more parameters in the basis of B , R_c and f_c , which will discuss in next section.

Figure 7. Cyclic spectral features of Costas signal in main support area

Figure 8. Cyclic spectral features of Costas signal in sub-support area

With -10dB noise added, measurement ability in Fig.5, 6, 7 and 8 remains fairly robust [15] since AWGN is not correlated and is suppressed in this cyclic spectral analysis.

• IV.    C YCLIC F EATURES E XTRACTION OF C OMPLEX M ODULATED S IGNALS BASED ON ACP M ETHOD

The block diagram of the complex modulation extraction algorithm that uses ACP method results is shown in Fig. 9.

• A.    ACP Method Processing on Frank Coded Signal

The first step is to crop main and sub-support area[11]. After cropping the area of support, an adaptive threshold operation is performed to reduce the amount of noise present on bifrequency plane[10]. This is followed by scanning axis for detection of maximum magnitude on bifrequency plane. The scanning results are used to measure the bandwidth B , carrier frequency f_c , subsequently calculate the code length N modulation period t_m and modulation period T .

An approximation to the PDF is constructed by scanning both the i -and j -axis and finding the magnitude above threshold corresponding to each index in main and sub-support areas simultaneously. The horizontal scan gives the i -index values as

Figure 9. Block diagram of complex modulation signal feature extraction processing using ACP method

N x,   a S(,j)                             (22)

j = 1

and

PDF. = . x ^/=(1: N                           (23)

a x j = (1: N )

This energy distribution enables the location of the support areas in the bifrequency plane to be determined. In main support area, the scan on both i - and j - axis from the lower index values to the higher index values to a level above -6dB gives i ₁₁ and j ₁₁ ; from the higher index values to the lower index values to a level above -6dB gives i ₁₂ and j ₁₂ . In sub-support area, the scan remains the same as in main support area and give i ₂₁ , i ₂₂ , j ₂₁ , j ₂₂ . Thus

2 f = ^{i 21} 2 ⁱ 2 ²                                  (24)

and

B = j ¹¹ 2 j iz                               (25)

We introduce(24) and (25) into (15),the modulation time t_b ,code length N and modulation period T can be obtained.

• B.    ACP Method Processing on Costas Coded Signal

For Costas coded signal, the ACP algorithm processing steps remains the same on the whole, but the key parameter fc in Frank coded signal becomes frequency modulation sequence fm [11].In main support area of Costas coded signal, the scan on both i - and j -axis from the lower index values to the higher index values to a level above -6dB gives i11 and j11 ; from the higher index values to the lower index values to a level above -6dB gives i1M and j1M . In sub-support area, the scan remains the same as in main support area and givei21 i22 , j2M , j2M . Thus

2 f = i ! i ⁺ i !( i + 1) c ₂

and

B = j i1 ^{+ j}M

1 ^M

R„              (^ - i )

^c    m - 1 i " ₂ ² i   ²⁽ i ^{- 1)}

We introduce(26), (27)and (28)into (16), the modulation time t_b ,code length N and modulation period T can be obtained[12].

In terms of a noisy background, the scan method described above will degrade the signal’s PDF. The procedure can be modified to use an adaptive noise filter[13].

• V.  T ESTING THE A LGORITHM

To evaluate the algorithm’s accuracy, the parameter extraction algorithm is tested in this section.

• A.    ACP Method Processing on Frank Coded Signal

The algorithm checks the signal for the existence of the noise and uses cyclic spectral analysis method to extract the parameters. The signals with no noise can analyzed with 64Hz of frequency resolution but the signals with -6dB noise cannot analyzed with it due to the computational difficulties. The signal with -6dB noise do not have a code rate solution, since the resolution of the frequency should be smaller than the largest parameter being measured or extracted [14].When code length N_c is

9 and 16, the corresponding code rate R_c are

R = — = ---¹---= 111.1 Hz > 64 Hz   (29)

Nt b,    9 * 10 ^{- 3} s

R = — = ---¹---= 62.5 Hz < 64 Hz   (30)

c   Nt6   16 * 10- 3s from (29)and (30)We see the frequency resolution D f can achieve 64Hz when Nc £ 9 ,while N3 16 , D f can’t achieve 64Hz,,neither the code Rc ,the effectiveness of ACP is decided by code length Nc .

Three Frank coded signals are analyzed in the test, the key parameters are f_c , f_s , B , N , and the range of SNR is -10~0dB.The normalized RMS error (NRMSE)[15] is used to verify the effectiveness of parameter extracting. If $ k ( k = 1,2, L , N ) is the estimate value of a quantity whose exact value is a , then NRMSE is defined by:

NRMES = ![1 ^ k - a ^J ¹/² a ) Nl = 1

and gives an indication of how good a measurement is relative to the size of the thing being measured.

200 times Monte Carlo simulation tests are carried out in the condition of different SNR. The parameters used to generate the Frank coded modulations are f_s =7000Hz, f_c =1000Hz, B =1000Hz, and N_c =4, 9, 16. The results of parameters extraction are showed in Table 1 The NRMSE of tests are showed in Table 2, The NRMSE curves are plotted in Fig. 10 and Fig.11.

Figure 10. f_c NRMSE curves of Frank coded signal

TABLE I.     T URE AND M EASURED P ARAMETERS OF F RANK

S IGNAL

Signal	True Parameters
	fc	B	Rc	tb	T	Nc
1	1K	1k	250	1	4	4
2	1K	1k	111	1	9	9
3	1K	1k	62.5	1	16	16
Signal	Measured Parameters
1	1001	1012	246	0.99	4.07	4.11
2	998	1231	116	0.81	8.62	10.61
3	1073	957	61.5	1.04	16.3	15.56

Figure 11. Bandwidth B NRMSE curves of Frank coded signal

TABLE II.    NRMSE OF K EY P AMAMETERS OF F RANK S IGNAL

Signal	fc	B	Rc	t µ b	T µ	N µ
1	0.081	0.055	0.013	0.036	0.013	0.087
2	0.073	0.049	0.010	0.032	0.012	0.127
3	0.072	0.031	0.009	0.029	0.010	0.053

For the test results, The NRMSE curve of carrier frequency f_c and bandwidth B are shown in Fig.6, The curves of code rate R_c and modulation period T are just like the curve of f_c [17],which are influenced by resolution of cyclic frequency D a ;the curves of modulation time t_b is just like the curve of B ,which is influenced by the resolution of frequency D f .

The test results shows that the NRMSE curves of B and t_b are more remarkable than that of the other parameters. The NRMSE of code length N_c is decided by both D a and D f , the results is more complicated. The NRMSE of other parameters decline with the increase of N , since the long code length can obtain good processing gain.

• B.    ACP Method Processing on Costas Coded Signal

First we evaluate the computation complexity of ACP method on feature extraction of Costas signals. The code length is N , ACP method compute two parallel FFT, and each computation follow the linear logarithm complexity[18], means O ( 2 N l og₂ N ) .The traditional cyclic spectral estimator compute double FFT, the computation follow the quadratic order complexity, mean O ( N ) . O ( 2 N log₂ N ) is lower order than O ( N ) in general[19], and the computation complexity of ACP method declines significantly as code length N increases. Meanwhile, the resolution restriction is just the same as Frank coded signal showed in (29) and (30), hence the effectiveness of ACP is decided by code length N .

200 times Monte Carlo simulation tests are carried out in the condition of different SNR. The parameters used to generate the Costas coded modulations are f_s =7000Hz, f_c =1000Hz, B =1000Hz, and f m = {4,7,1,6,5,2,3} Hz . The NRMSE curves are plotted in Fig. 12 .

Figure 12. True and Measured Parameters of sig5 and sig6

in Fig.12,The curves of the key parameters preserve the same declining type. The NRMSE of modulation time t_b , code length N and modulation period T is larger than the other parameters, which means the former three parameters are easily influenced by AWGN[20]. In addition, the ACP method processing results of Costas coded signal has the more limitations than Frank coded signal. The effectiveness of extracting limited by the resolution of frequency and cyclic frequency on both situations, while ACP method can’t determine the order of Costas coded signal frequency modulation sequence.