Application of Biphase Complete Complementary Code for Ionospheric Sounding

This paper illustrates the processes carried out for the application of biphase complete complementary code (CCC) for ionospheric sounding to address the coherent interference problem in multi-station ionospheric sounding. An algorithm to generate the biphase CCC is described, and the detailed process of waveform construction and signal processing is presented. Characteristics of the autocorrelation and cross-correlation are analyzed through simulations, and the technical feasibility of the application of CCC is explored. Experiments of ionospheric sounding with the CCC are also implemented to verify performance. Results demonstrate that the CCC performs well in multi-station ionospheric sounding, and is capable of eliminating the coherent interference in the network of ionosondes, compared to the conventional complementary code.


Introduction
Activities of ionospheric sounding started from 1926 when remote sensing of the ionosphere by means of radio sounding was described [1]. With the development of digital technology, analog ionosondes evolved into digital instruments and custom computer circuits capable of digital integration and digital spectrum analysis. Reinisch et al. [2] developed the Digisonde Portable Sounder (DPS) with a digital signal processor in 1997. In order to establish a denser network of ionosondes, the Low-Cost Ionosonde was proposed by Stamper et al. [3]. The latest digital ionosonde Digisonde-4D [4] implements the advanced capabilities of new digital radio frequency (RF) circuitry and embedded computers while preserving the basic principles of Digisonde family. A network of ionosondes equipped with Digidonde-4D has been established in Europe to monitor travelling ionospheric disturbances (TIDs) by simultaneous vertical and oblique incidence sounding [5]. Nowadays, the universal software radio peripheral (USRP) is implemented in ionospheric sounding owing to its flexibility in data processing, less rigidity of hardware structures, and generally larger data transmission channel. Chirp sounding with software-defined radio (SDR) technology was proposed in 2003 [6]. Barona et al. [7] illustrated the construction of an ionosonde from a USRP with GUN-Radio and MATLAB. In addition, BPSK signal shaping and processing with a USRP was proposed by Ivanov et al. [8], which implies that SDR technology has received increased attention.
In 1985, Wuhan University (Wuhan, China) collaborated with the University of Paris-Sud (Paris, France) in the field of ionosphere research. Since then, Wuhan University has been devoted to ionospheric study and the development of ionospheric sounding systems. In 2003, the Wuhan Ionospheric Oblique Backscattering Sounding System (WIOBSS) was successfully developed by the Ionosphere Laboratory at Wuhan University. The system is able to operate at three different modes, including the vertical, oblique, and oblique backscatter sounding functions. In order to achieve frequency and clock synchronization, the Global Positioning System (GPS) was applied to the receiver. Recently, the frequency band of WIOBSS has been extended from high frequencies to very high frequencies to improve system performance [9]. Besides, Shi et al. [10] proposed an oblique-incidence ionosonde and applied it to the localization of ionospheric irregularities. Among the WIOBSS family, intrapulse phase-coded pulses are used as the emission signal so that the receiver can receive the echoes during the intrapulse periods. In order to achieve a high signal processing gain, long pseudorandom sequences such as the m-sequences and the almost perfect autocorrelation sequences are adopted [11]. Yao et al. [12] reported ionograms obtained by WIOBSS with different waveforms, and determined that the interphase-coded complementary code has the best performance. However, when two stations perform ionospheric sounding synchronously with the same complementary code, both vertical and oblique ionograms can be observed, thereby raising the difficulty of automatically extracting ionospheric parameters from the present ionograms using the current methods [13]. On the other hand, coherent interference can exist if the complementary code is non-orthogonal. Furthermore, if there are more than two stations operating, the coherent interference problem becomes worse. Research on coherent interference suppression has a long history; the methods are generally based on spatial averaging [14] for narrow band signals or spectral averaging [15] for broadband signals. Adaptive beamforming methods have been proposed to overcome the problem [16][17][18]. However, these methods aimed at the cancellation of the interference, are not appropriate for application in ionospheric sounding.
The complete complementary code (CCC) was first proposed by Suehiro et al. [19] to study the N-shift cross-orthogonal sequences. Additionally, Suehiro designed a multicarrier CDMA architecture based on the orthogonal complementary code [20], which achieves high spreading efficiency and offers multiple access interference (MAI) free operation. Recently, the CCC has also been applied in the MIMO system because of its ideal correlation property [21][22][23]. Quadri-phase complete complementary pairs were adopted in the topside ionosphere sounder system [24]. In this study, a biphase CCC is developed for application to ionospheric sounding in order to overcome the aforementioned weaknesses associated with the ionospheric sounding network. We firstly present a brief definition of CCC and its generation algorithm, and then describe the construction of the waveform and the details of signal processing. Characteristics of the autocorrelation and cross-correlation are subsequently analyzed through simulations, followed by discussions of the experimental results of ionospheric sounding with the application of the complete complementary code.

Complete Complementary Code
For N M-order sequence sets {(S 11 , S 12 , . . . , S 1M ), (S 21 , S 22 , . . . , S 2M ), . . . ,(S N1 , S N2 , . . . , S NM )}, each element S ij is a sequence of length L, and (S 11 , S 12 , . . . , S 1M ) is a subset. We can use a {M, N, L} sequence family to represent the sequence sets, where M is the number of sequences in a sequence set, N is the number of sequence sets and L is the length of the sequence S n,m . When the sum of autocorrelation functions of each sequence in the same subset is zero except for zero shift, and the sum of crosscorrelation of two sequences of the same index in each two subsets is zero at all shifts, the {M, N, L} sequence family is considered to be complete complementary code. In other words, the {M, N, L} sequence meets the following Equations (1) and (2).
where R S ij (τ) is the autocorrelation function of the sequence S ij , R S ij ,S kj (τ) is the crosscorrelation function of the sequence S ij and S kj , δ(τ) is a typical impulse function, and τ is the offset.

Waveform and Signal Processing
There are many kinds of generation algorithms for the complete complementary code [25][26][27][28]. Suchiro et al. [19] extended the idea of complementary sets to what is known as the complementary code, and proposed a method to construct complete complementary sequences of length N 2 . Han [28] expanded the method, and proposed a further method to generate a set of length MN (M, N ∈ Z + ; N ≤ M), and then an improved method to generate a set of length MN/P(M, N, P ∈ Z + ; N, P ≤ M) [27], which offers a practical way to shorten the length of sequences in the subsets.
In this study, the CCC is generated through the method proposed by Han and Suchiro [28], which uses the unitary matrix to derive the matrix composed of CCC. In order to be applied to the biphase modulation system, the Hadamard matrix [29] is used as an alternative of the unitary matrix. Figure 1 shows the flow chart of {N, N, N 2 } CCC's generation. Firstly, the number of sequence sets N should be determined (N = 2 k , k ∈ Z + ). Secondly, three N-order Hadamard matrices should be generated; they can be indexed by A, B, C. Then, a N 2 × N 2 matrix S can be generated by a recapitulate formula. Three different matrix operations, Kronecker product, vectorization, and the "Diag" operator, are used in the formula. At last, we should transform the matrix S into sequence sets. Each row of S is a sequence, and each row of N is a subset. Table 1 shows a set of {4, 4, 16} complete complementary sequences. In Table 1, + and −, respectively representing +1 and −1, the four subsets are indexed by letters A, B, C, D, and the sequences in a subset are indexed by numbers 0, 1, 2, 3.     An intra-pulse biphase coding scheme was used in the WISS, where +1 and −1 represent the phase 0 and π respectively. According to the periodic correlation characteristics of the complete complementary code, sequences should be modulated in a certain order. Figure 2 shows the single-cycle transmitting waveform, where A(t) is the control sequence used to open or close the transmit channel; s n (t) is a phase modulation signal based on the binary complete complementary code, T p represents the pulse width, and T r is the pulse repetition interval.
complementary code, sequences should be modulated in a certain order. Figure 2 shows the single-cycle transmitting waveform, where A(t) is the control sequence used to open or close the transmit channel; (t) is a phase modulation signal based on the binary complete complementary code, represents the pulse width, and is the pulse repetition interval. The waveform of a single cycle is represented by Equation (3) where , is the i-th symbol of the sequence , , is the carrier frequency, and u(t) is a rectangle pulse represented by Equation (4), Suppose that there are two sites respectively transmitting s (t) and s (t), which means that the signal will travel in different propagation paths. Therefore, the delay time will be different, and they can be represented by ∆ , ∆ . Correspondingly, the Doppler shifts of the echoes are ∆ , ∆ . The echo signal can be demonstrated by Equation (5).
where , denote the amplitude of different waveforms. By the matched filter of s (t), the final output s(t) can be represented by Equation (6), Due to the ideal auto-correlation property of CCC, the first part of Equation (6) is an impulse function, and the second part keeps zero in the entire time domain because of the orthogonality. Thus, we can extract the delay time ∆ and Doppler shift ∆ . Additionally, if the echo signal is compressed by the matched filter of s (t), ∆ and ∆ can be derived.

Feasibility Analysis
where S i n,m is the i-th symbol of the sequence S n,m , f is the carrier frequency, and u(t) is a rectangle pulse represented by Equation (4), Suppose that there are two sites respectively transmitting s 1 (t) and s 2 (t), which means that the signal will travel in different propagation paths. Therefore, the delay time will be different, and they can be represented by ∆τ 1 , ∆τ 2 . Correspondingly, the Doppler shifts of the echoes are ∆ f 1 , ∆ f 2 .
The echo signal can be demonstrated by Equation (5).
where A 1 , A 2 denote the amplitude of different waveforms. By the matched filter of s 1 (t), the final output s(t) can be represented by Equation (6), Due to the ideal auto-correlation property of CCC, the first part of Equation (6) is an impulse function, and the second part keeps zero in the entire time domain because of the orthogonality. Thus, we can extract the delay time ∆τ 1 and Doppler shift ∆ f 1 . Additionally, if the echo signal is compressed by the matched filter of s 2 (t), ∆τ 2 and ∆ f 2 can be derived.

Feasibility Analysis
In this section, characteristics of auto-and cross-correlation of the CCC and complementary codes are studied. Figure 3a shows the autocorrelation functions of two sets of non-orthogonal 16-bit complementary code and two subsets of a {4, 4, 16} complete complementary code. It can be noted from Figure 3a that the autocorrelation coefficients are zero except zero offset for the complementary code and CCC. Figure 3a indicates that both the complementary code and CCC have the ideal autocorrelation property. Since the coherent accumulation time of the {4, 4, 16} complete complementary code is twice of the complementary code, the autocorrelation gain is doubled in Figure 3a. Figure 3b shows the cross-correlation functions of the complementary code and CCC. The solid and dashed lines, respectively, represent the cross-correlation function between the two subsets of complete complementary code and the two non-orthogonal complementary codes. It can be seen from Figure 3b that the cross-correlation function of the CCC is better than that of the complementary code. Therefore, the ideal cross-correlation property of the complete complementary code makes it feasible for multiple stations to operate together without coherent interference.  Figure 3a. Figure 3b shows the cross-correlation functions of the complementary code and CCC. The solid and dashed lines, respectively, represent the cross-correlation function between the two subsets of complete complementary code and the two non-orthogonal complementary codes. It can be seen from Figure 3b that the cross-correlation function of the CCC is better than that of the complementary code. Therefore, the ideal cross-correlation property of the complete complementary code makes it feasible for multiple stations to operate together without coherent interference.  For pulse compression radars, an ambiguity function χ(∆τ, ∆f) is a two-dimensional function to measure the distortion of an echo signal, where ∆τ is the group delay of transmitted signal and ∆f is the Doppler shift. The narrowband ambiguity function is represented by Equation (7) in ionospheric sounding, where ( ) is the transmitted pulse waveform, and * ( − ) is the detected signal. The ambiguity function of the complete complementary code is the sum of ambiguity functions of each sequence in a subset [30], which can be represented by Equations (8) and (9), where and ∆τ = |∆τ − S |. Figure 4a shows the normalized ambiguity function of a {4, 4, 16} complete complementary code. In consideration of the pulse width of waveform in WISS, is set 25.6 μs. Figure 4b-g is For pulse compression radars, an ambiguity function χ(∆τ, ∆ f ) is a two-dimensional function to measure the distortion of an echo signal, where ∆τ is the group delay of transmitted signal and ∆ f is the Doppler shift. The narrowband ambiguity function is represented by Equation (7) in ionospheric sounding, where s(t) is the transmitted pulse waveform, and s * (t − τ)e j2π f t is the detected signal. The ambiguity function of the complete complementary code is the sum of ambiguity functions of each sequence in a subset [30], which can be represented by Equations (8) and (9), where and ∆τ 1 = ∆τ − ST p . Figure 4a shows the normalized ambiguity function of a {4, 4, 16} complete complementary code. In consideration of the pulse width of waveform in WISS, T p is set 25.6 µs. Figure 4b-g is different cross sections of the ambiguity function, where f is the Doppler shift. When f = 0, the complete complementary code keeps an ideal autocorrelation property. When f = 5, the sidelobe peak of the normalized ambiguity function is 0.00025, about −36 dB, indicating that the complete complement code still has excellent autocorrelation characteristics. When f = 200, the sidelobe peak rises to 0.01, approximately −20 dB. It is seen that when Doppler shift fluctuates at a large range, CCC maintains low sidelobe peak level. Table 2 shows the sidelobe peak of ambiguity function with different Doppler shifts. As the Doppler shift grows, the sidelobe peaks of the m-sequence and Barker code hardly change, while those of the complementary code and CCC decrease considerably. Since the Doppler shift is generally less than 5 Hz in ionospheric sounding [31], it infers that the Doppler shift has little effect on the application of the complete complementary code.  Table 2 shows the sidelobe peak of ambiguity function with different Doppler shifts. As the Doppler shift grows, the sidelobe peaks of the m-sequence and Barker code hardly change, while those of the complementary code and CCC decrease considerably. Since the Doppler shift is generally less than 5 Hz in ionospheric sounding [31], it infers that the Doppler shift has little effect on the application of the complete complementary code. Additionally, the number of subsets M and the length of sequence L cannot be too large. The longer the meta-code sequence, the larger the detection dead zone. Actually, the detection of the ionosphere is based on the smoothness of the ionosphere, and the state of the ionosphere is considered to be undisturbed within certain minutes. If the single cycle lasts too long, the detection validity can decrease significantly.  Additionally, the number of subsets M and the length of sequence L cannot be too large. The longer the meta-code sequence, the larger the detection dead zone. Actually, the detection of the ionosphere is based on the smoothness of the ionosphere, and the state of the ionosphere is considered to be undisturbed within certain minutes. If the single cycle lasts too long, the detection validity can decrease significantly.

Experiments and Discussions
In this study, the WISS, installed at Ningqiang (105.20 • E, 32.27 • N) and Wuhan (114.22 • E, 30.33 • N), are utilized to carry out multi-station ionospheric sounding experiments to investigate the performance of the CCC application. The system operated at the HF-band to accomplish vertical and oblique incidence sounding. The ionosondes operated in the swept-frequency mode. Three-wire dipole and log-periodic antennas were used to transmit and receive the signals, respectively. In the present experiment, the biphase intrapulse waveform was coded by a 16-bit complete complementary code or a {4, 4, 16} complete complementary code. All the stations synchronously transmitted and received the signals through the global positioning system (GPS). Table 3 shows the parameters of the emission waveform. Figure 5 illustrates the ionograms recorded at Ningqiang on 23 June 2018 when the ionosondes installed at Ningqiang and Wuhan synchronously transmitted the same 16-bit complementary code at 14:00 Beijing time (BJT). Parameters of the ionosphere at Ningqiang can be obtained from ionograms. For example, based on the delay time, we can derive the virtual height of the E and F2 layers (h'E = 140 km, h'F2 = 330 km). The red circles in Figure 5 show the reflection signal echoes of the F-layer at the middle of the propagation path from Wuhan to Ningqiang. It is seen from Figure 5 that the echoes of the vertical sounding at Ningqiang and oblique sounding from Wuhan overlapped on the ionograms. Obviously, these kinds of ionograms are difficult to be scaled automatically using the current software.

Experiments and Discussions
In this study, the WISS, installed at Ningqiang (105.20° E, 32.27° N) and Wuhan (114.22° E, 30.33° N), are utilized to carry out multi-station ionospheric sounding experiments to investigate the performance of the CCC application. The system operated at the HF-band to accomplish vertical and oblique incidence sounding. The ionosondes operated in the swept-frequency mode. Three-wire dipole and log-periodic antennas were used to transmit and receive the signals, respectively. In the present experiment, the biphase intrapulse waveform was coded by a 16-bit complete complementary code or a {4, 4, 16} complete complementary code. All the stations synchronously transmitted and received the signals through the global positioning system (GPS). Table 3 shows the parameters of the emission waveform.   Figure 5 show the reflection signal echoes of the F-layer at the middle of the propagation path from Wuhan to Ningqiang. It is seen from Figure  5 that the echoes of the vertical sounding at Ningqiang and oblique sounding from Wuhan overlapped on the ionograms. Obviously, these kinds of ionograms are difficult to be scaled automatically using the current software.  In the second experiment, the ionosondes installed at Ningqiang and Wuhan synchronously transmitted the different codes, non-orthogonal complementary codes, at 16:00 BJT. Figure 6 shows the ionograms recorded at Ningqiang on 14 July 2018. Similar to Figure 5, Figure 6 also displays the vertical sounding echoes from Ningqiang and oblique sounding echoes (in red circles) from Wuhan. The signals received at Ningqiang are processed by the local-field emission code. Due to the non-orthogonal complementary codes, the receiver at Ningqiang also recorded the oblique signal transmitted from Wuhan. Similar to the first experiment, the ionograms in this case are also difficult to scale automatically using the current tool. This therefore inspires us to further change the transmitted codes in the experiment to avoid the interference resulting from the oblique sounding signal from Wuhan.
In the third experiment, the complete complementary code was applied for ionospheric sounding. The transmitted code was a {4, 4, 16} CCC, as shown in Table 1. The ionosondes at Ningqiang and Wuhan, respectively, transmitted the subsets A and B to carry out ionospheric sounding. Figure 7 shows the ionograms recorded at Ningqiang on 23 June 2018. Figure 7a  In the second experiment, the ionosondes installed at Ningqiang and Wuhan synchronously transmitted the different codes, non-orthogonal complementary codes, at 16:00 BJT. Figure 6 shows the ionograms recorded at Ningqiang on 14 July 2018. Similar to Figure 5, Figure 6 also displays the vertical sounding echoes from Ningqiang and oblique sounding echoes (in red circles) from Wuhan. The signals received at Ningqiang are processed by the local-field emission code. Due to the non-orthogonal complementary codes, the receiver at Ningqiang also recorded the oblique signal transmitted from Wuhan. Similar to the first experiment, the ionograms in this case are also difficult to scale automatically using the current tool. This therefore inspires us to further change the transmitted codes in the experiment to avoid the interference resulting from the oblique sounding signal from Wuhan.
In the third experiment, the complete complementary code was applied for ionospheric sounding. The transmitted code was a {4, 4, 16} CCC, as shown in Table 1. The ionosondes at Ningqiang and Wuhan, respectively, transmitted the subsets A and B to carry out ionospheric sounding. Figure 7 shows the ionograms recorded at Ningqiang on 23 June 2018. Figure 7a,b correspond to the echoes processed by the subset A (vertical ionograms). Figure 7c,d correspond to the echoes processed by the subset B (oblique ionograms). Figure 7 clearly shows that the vertical and oblique echoes can be separated by use of different complete complementary code subsets. Based on the vertical ionograms, we can derive the virtual height of E and F2 layers in Ningqiang (h'E = 130 km, h'F2 = 290 km). Taking into account that the distance between Ningqiang and Wuhan is about 820 km, we can derive the ionospheric parameters in the midpoint of the propagation path according to the delay time in the oblique ionograms (h'E = 135 km, h'F2 = 310 km). Since the vertical and oblique echoes are distinguished by the matched filters of different subsets of CCC, the ionograms can be scaled automatically with the software. In the second experiment, the ionosondes installed at Ningqiang and Wuhan synchronously transmitted the different codes, non-orthogonal complementary codes, at 16:00 BJT. Figure 6 shows the ionograms recorded at Ningqiang on 14 July 2018. Similar to Figure 5, Figure 6 also displays the vertical sounding echoes from Ningqiang and oblique sounding echoes (in red circles) from Wuhan. The signals received at Ningqiang are processed by the local-field emission code. Due to the non-orthogonal complementary codes, the receiver at Ningqiang also recorded the oblique signal transmitted from Wuhan. Similar to the first experiment, the ionograms in this case are also difficult to scale automatically using the current tool. This therefore inspires us to further change the transmitted codes in the experiment to avoid the interference resulting from the oblique sounding signal from Wuhan.
In the third experiment, the complete complementary code was applied for ionospheric sounding. The transmitted code was a {4, 4, 16} CCC, as shown in Table 1. The ionosondes at Ningqiang and Wuhan, respectively, transmitted the subsets A and B to carry out ionospheric sounding. Figure 7 shows the ionograms recorded at Ningqiang on 23 June 2018. Figure 7a,b correspond to the echoes processed by the subset A (vertical ionograms). Figure 7c,d correspond to the echoes processed by the subset B (oblique ionograms). Figure 7 clearly shows that the vertical and oblique echoes can be separated by use of different complete complementary code subsets. Based on the vertical ionograms, we can derive the virtual height of E and F2 layers in Ningqiang (h'E = 130 km, h'F2 = 290 km). Taking into account that the distance between Ningqiang and Wuhan is about 820 km, we can derive the ionospheric parameters in the midpoint of the propagation path according to the delay time in the oblique ionograms (h'E = 135 km, h'F2 = 310 km). Since the vertical and oblique echoes are distinguished by the matched filters of different subsets of CCC, the ionograms can be scaled automatically with the software.

Conclusions
In this study, a biphase CCC has been proposed for application to multi-station ionospheric sounding system to overcome the coherent interference. We have developed an algorithm of CCC generation and constructed the corresponding emission waveform for realistic experiments. Our experimental results at the HF band show that the coherent interference can be effectively eliminated in multi-station ionospheric sounding by using CCC. Moreover, the separation of vertical and oblique echoes on the single ionogram make it feasible for automatically scaling of ionograms in multi-station ionospheric sounding system. Aiming to eliminate the coherent

Conclusions
In this study, a biphase CCC has been proposed for application to multi-station ionospheric sounding system to overcome the coherent interference. We have developed an algorithm of CCC generation and constructed the corresponding emission waveform for realistic experiments. Our experimental results at the HF band show that the coherent interference can be effectively eliminated in multi-station ionospheric sounding by using CCC. Moreover, the separation of vertical and oblique echoes on the single ionogram make it feasible for automatically scaling of ionograms in multi-station ionospheric sounding system. Aiming to eliminate the coherent

Conclusions
In this study, a biphase CCC has been proposed for application to multi-station ionospheric sounding system to overcome the coherent interference. We have developed an algorithm of CCC generation and constructed the corresponding emission waveform for realistic experiments. Our experimental results at the HF band show that the coherent interference can be effectively eliminated in multi-station ionospheric sounding by using CCC. Moreover, the separation of vertical and oblique echoes on the single ionogram make it feasible for automatically scaling of ionograms in multi-station ionospheric sounding system. Aiming to eliminate the coherent interference in the network of radio systems, this research has the important implication to locate ionospheric irregularities, monitor travelling ionospheric disturbances (TIDs), and even enhance the global navigation satellite systems (GNSS).
However, the number of subsets and the number of subset sequences are the primary limiting factors for the application of the complete complementary code. Since the number of subsets cannot exceed the number of meta-code sequences, increasing the number of subset sequences therefore means the exponential growth of a single probe cycle, which may subsequently enlarge the detection dead zone and degrade the applicability of the complete complementary code.