Interference Excision in Spread Spectrum Communications Using Adaptive Positive Time-Frequency Analysis

This paper introduces a novel algorithm to excise single and multicomponent chirp-like interferences in direct sequence spread spectrum (DSSS) communications. The excision algorithm consists of two stages: adaptive signal decomposition stage and directional element detection stage based on the Hough-Radon transform (HRT). Initially, the received spread spectrum signal is decomposed into its time-frequency (TF) functions using an adaptive signal decomposition algorithm, and the resulting TF functions are mapped onto the TF plane. We then use a line detection algorithm based on the HRT that operates on the image of the TF plane and detects energy varying directional elements that satisfy a parametric constraint. Interference is modeled by reconstructing the corresponding TF functions detected by the HRT, and subtracted from the received signal. The proposed technique has two main advantages: (i) it localizes the interferences on the TF plane with no cross-terms, thus facilitating simple ﬁltering techniques based on thresholding of the TF functions, and is an e ﬃ cient way to excise the interference; (ii) it can be used for the detection of any directional interferences that can be parameterized. Simulation results with synthetic models have shown successful performance with linear and quadratic chirp interferences for single and multicomponent interference cases. The proposed method excises the interference even under very low SNR conditions of − 10 dB, and the technique could be easily extended to any interferences that could be represented by a parametric equation in the TF plane.


INTRODUCTION
In spread spectrum (SS) communications, the message signal is modulated and spread over a wider bandwidth with a pseudonoise (PN) code also known at the receiver, and transmitted over the channel. The increase of the bandwidth yields a processing gain, defined as the ratio of the bandwidth of the transmitted signal to the bandwidth of the message signal, and it provides a high degree of interference suppression. However, there is a tradeoff between increasing the processing gain and the available frequency spectrum. In the case of a jammer with high power, the SS system may not be able to suppress the interference. Therefore, excising the interference prior to despreading the received signal is necessary to increase the performance of the system.
Most interference suppression techniques are designed to deal with narrowband interferences [1][2][3][4][5]. Among the time-domain approaches for narrowband interference excision, the most notable methods include adaptive notch filtering and decision-directed adaptive filtering techniques [6].
While SS systems can successfully reject narrowband interferences, their performance in rejecting wideband interferences is limited. In practical systems, it is not likely to transmit high-power wideband jamming signals due to the power limitations of the interference source. Additive white Gaussian noise can be considered as the only realizable wideband interference, which is very challenging to predict and excise. Therefore, substantial amount of research has been conducted on wideband interferences with narrowband instantaneous frequency elements such as FM signals. Most of these methods focus on suppressing the interference using TF distributions (TFDs) to localize the interference signals [7]. However, commonly used TFDs suffer from a tradeoff between the joint TF resolution and cross-term suppression. In [8], Amin proposed a method based on the Wigner-Ville distribution (WVD) of the signal, which represents the signal with precise TF localization; yet, the method is shown to suffer from cross-terms in the presence of multicomponent interferences. In the extension of this work [9], the authors use the Wigner-Hough transform ( Figure 1: Block diagram of a DSSS system. crossterms [9]; however, the system is shown to be sensitive to the signal model. Also, wavelet-based approaches provide a time-scale representation and perform poorly for narrowband activities in the high frequency range. In [10,11], different window length STFTs are used to localize the interference. In [12], the authors use a signal decomposition algorithm consisting of a chirp-based dictionary to represent linear chirp interferences on the TF plane. The chirp interferences can be modeled with few coefficients and the proposed method performs well with linear chirp interferences. However, the generalization of the system to include quadratic, hyperbolic, or sinusoidal FM interferences is not discussed. In [13], the instantaneous frequency of the interference is recursively estimated using the discrete evolutionary and Hough transforms, and the interference is subtracted from the signal by using the singular value decomposition of the de-chirped signal. In [14], the authors propose an adaptive TF exciser that decides the domain of the excision by evaluating both the time and frequency properties. This system performs well in the case of narrowband interferences. There are also the TF projection filtering techniques as proposed in [15,16]. A common characteristic of most of the interference excision algorithms is the continuing presence of bit errors even after the interference is suppressed. Earlier interference excision methods based on TFDs suffer from a tradeoff between the TF resolution and the TFD cross-terms [17][18][19]. They also perform the excision of limited type of interferences such as linear or sinusoidal interferences [9,12]. Considering these two disadvantages, we propose a new excision method based on constructing a positive TFD of the received SS signal using an adaptive signal decomposition technique, the matching pursuit (MP) algorithm [20]; followed by a line detection algorithm based on the HRT. By decomposing a signal into its components, the interaction between components can be kept under control and possibly eliminated. The decomposition will allow the construction of a cross-term free TFD by combining the TFDs of the individual components generated by the decomposition. Also, by using Gaussian functions as bases for decomposition, we can achieve a high TF resolution, since the Gaussian functions satisfy the equality in the uncertainty principle and provide optimal TF resolution [7]. We construct the TFD of the TF functions resulting from the MP, treat the TFD as an image, and detect the interfering signals using the HRT, which can detect any line satisfying a parametric equation. We then reconstruct a model of the interfering chirps using the TF functions and excise the reconstructed interference from the received signal.
The paper is organized as follows. In Section 2, the elements of a DSSS system are introduced. In Section 3, an adaptive signal decomposition technique based on MP and the construction of positive TFDs are explained. Then the HRT is introduced for the detection of chirps on the TF plane. This section ends with outlining the steps of the excision algorithm consisting of MP and HRT. In Section 4, the performance of the proposed system is evaluated in terms of jammer-to-signal-power ratio, bit error rate, and average chip error rate. The paper is concluded in Section 5.

DSSS SYSTEM
Let us consider a DSSS system as shown in Figure 1. In this system, the transmitter generates an SS signal which in turn is transmitted over a communications channel as a binary phase shift keying (BPSK) modulated signal. Additive channel noise as well as jamming signal act on the transmitted signal. At the receiver, the noise and interference corrupted signal is first demodulated. The "standard" SS receiver correlates the baseband SS signal with the synchronized PN sequence, and the resulting signal is processed and input into a threshold detector to estimate the transmitted binary data sequence.
Let b k = ±1 be the kth message symbol transmitted in a DSSS system such that where p k = [c 0 , . . . , c L−1 ] T for {k = 1, 2, . . .} is a PN sequence with a chip length L, c n = ±1 is the nth chip of the PN sequence, and w k is the SS signal. The received signal r k at the output of the BPSK demodulator will consist of the SS signal w k , the additive white Gaussian noise term n k , and the interference i k such that We use the notation r to refer to the received signal sequence: Similarly, we use the notations w, n, and i to refer respectively to the complete SS signal, noise, and interference sequences before they are separated into L-element vectors in the form w k , n k , and i k , for k = 1, 2, . . . . At the receiver, the received signal r is first synchronized and correlated with the same spreading signal p. To estimate b k , we use the PN sequence p k to despread r k , and integrate the result to generate the test statistic Λ:

S. Krishnan and S. Erküçük
Using the test statistic Λ, we estimate the message symbols as Let E w = w T w, E n = n T n, and E i = i T i. We define signalto-noise-power ratio (SNR) and jammer-to-signal-power ratio (JSR) as To illustrate the effects of increasing jammer power on the message symbol detection, we simulated the channel output with E n = 0, and a linear chirp as the jamming signal which sweeps the entire frequency spectrum of w. We changed the JSR values from 0 to 60 dB in 5 dB steps. To measure the performance of the SS signal, we despread the received signal r = w + i, with the PN sequence p, integrate the resulting sequence, and compare the result with a threshold to estimate the transmitted message symbols. The bit-errorrate (BER) results obtained from this simulation provide a measure of the built-in interference suppression, that is, selfexcision capability of the SS system. Figure 2 shows the BER values at different JSR levels. The results presented in Figure 2 show that the SS system was able to completely self-excise the interference for JSR < 10 dB, as manifested with BER = 0. The resistance of the system to interference decreased with increasing JSR. For JSR > 40 dB, we observed BER ≈ 50% indicating that the SS system cannot suppress any part of the interference.
From these observations we conclude that preprocessing of the SS signals is an essential step in expanding the operating range of SS systems in high-JSR environments. In particular, the preprocessing operations take the form of modeling the interference and excising from the SS signal before the despreading and detection steps as shown in Figure 1.

INTERFERENCE EXCISION ALGORITHM
The interference excision algorithm is a two-step process based on the matching pursuit (MP) algorithm and the Hough-Radon transform (HRT) [21].

Matching pursuit algorithm
The MP algorithm [20] is an adaptive signal decomposition technique that can decompose the signal into its TF functions. In MP, the signal x(n) of length N is decomposed into a linear combination of TF functions in {g γm (n)}, and can be represented as where The set {a m } are the expansion coefficients, and {g(n)} is the window function. K sm normalizes g(n). The scale factor s m and the temporal placement parameter p m control the width and the displacement of the window function, respectively. The parameters k m and φ m represent the frequency and the phase of the exponential function, respectively. k m allows the search for different frequencies at each scale. The discrete dictionary is limited with the set {γ m } = {(s m , p m , 2πk m /N), 1 < s m < N, 0 ≤ p m < N, and 0 ≤ k m < N}. One possible set of functions to be used in the dictionary is the set of Gaussian functions, where The equality in the uncertainty principle holds for Gaussian signals resulting in an optimal TF resolution [7]. Therefore, a dictionary consisting of Gaussian functions would result in fine TF resolution. In MP, the signal x(n) is projected onto the dictionary {g γm (n)} of TF functions with all possible window sizes, frequencies, and temporal placements. At each iteration, the best-correlated function g γm is selected from the dictionary and the remainder of the signal, which is called the residue, is further decomposed using the same iteration procedure. After M iterations, the signal x(n) can be represented as

EURASIP Journal on Wireless Communications and Networking
where R m x represents the residue of the signal x(n) after m iterations with R 0 x = x, such that The first term in (10) represents the first M Gaussian functions best matching the signal (we will refer to the first term as x (n)) and the second term (referred to as x (n)) represents the residue of the signal x(n). In order for the signal to be fully decomposed, the iteration process continues until all the energy in the residue signal is consumed. However, for some applications such as denoising, the signal does not need to be fully decomposed.
After the signal decomposition is achieved, the TFD W (n, w) may be constructed by taking the WVD [20] of the Gaussian functions represented in x (n): where W gγ m (n, w) is the WVD of the Gaussian function g γm (n) and W gγ m (n, w) takes discrete time and frequency values since {γ m } is a set of integers. The second term in (12) corresponds to the cross-terms of the WVD and should be rejected in order to obtain a cross-term free energy distribution of x (n) in the TF plane [20]. Therefore, the MP TFD which we denote with the symbol W (n, w) is given as The resulting MP TFD is a cross-term free distribution with high TF resolution.

The Hough-Radon transform
The directional interferences can be energy varying. Therefore, we require a directional element detector that can detect time-varying energy values. The line detector that can satisfy our needs is a detector that uses the combination of Hough and Radon transforms proposed in [22]. This detector has been mathematically proven to be an optimal detector as it provides the maximum likelihood identification of a chirp signal [17]. The combined Hough and Radon transform, the HRT, is an efficient tool to detect directional and time-varying energy components in the TF plane. We first discuss the Hough transform and the Radon transform, and then continue to discuss the advantages of using the combined HRT for TFDs.

The Hough transform
The Hough transform (HT) is a pattern recognition method for calculating the number of points that satisfy a parametric constraint. It is used in image processing applications such as object detection, texture analysis, character recognition, directional image analysis, and image compression. Although HT is mainly applied to straight-line detection, it can also be applied to other curves that can be described by equations [23]. Let the parametric constraint be represented as where U = (u 1 , u 2 , . . . , u K ) is a point in the space of possible features and Θ = (θ 1 , θ 2 , . . . , θ L ) is a point in the space of parameters. The parameter space is commonly referred to as Hough space. The constraint may represent a curve, a line or a surface depending on the interpretation of the feature point. Each point Θ 0 in the parameter space represents a constraint that is a particular instance of a curve, line or a surface. The constraint may be mapped into the Hough space by evaluating The parameter values consistent with the existence of a given feature point U 0 are curves that the particular point may lie on, and are given by Given a number of feature points that satisfy a constraint specified by the parameter Θ 0 , the sets generated by (17) for each feature will contain the point Θ 0 . Furthermore, (17) may be viewed as a hypersurface in a continuous space of parameters. The curves of features satisfying a particular constraint will intersect at the common point Θ 0 in the parameter space. The HT can be only applied to binary images.

The Radon transform
The Radon transform (RT) is a commonly used line detection technique in computer tomography [24]. The RT computes the projections of different angles of an image (TFD) or 2D data distribution i f (u, v) measured as line integrals along ray paths [24]: where θ is the angle of the ray path of integration, ρ is the distance of the ray path from the center of the image, and δ is the Dirac delta function. Equation (18) represents integration of i f (u, v) along the line ρ = u cos θ+v sin θ as illustrated in Figure 3. ρ denotes the distance of the perpendicular bisect from the origin, and θ denotes the angle spanned by the line. The RT adds up the pixel values in the given image along a straight line in a particular direction and at a specific displacement. The RT can be applied to both binary and graylevel images.

Combined Hough and Radon transform
The Hough and the Radon transforms are individually not adequate to detect directional elements with varying energy levels. The underlying principle of the HT is that it is a process for counting the number of pixels that satisfy parametric constraints in a binary image. This property may result in misdetection of some energy varying components. The RT may be seen as a special case of the HT for straight-line detection. While the RT can be applied to gray-level images, it does not encompass all possible variations of the HT. Considering the advantages and disadvantages of each transform, we use the combined HRT as proposed in [22]. Using the combined HRT, we can detect the pixels that form a parametric constraint in a gray-level image. These constraints can be straight lines or curves in the image of the TF plane. We will consider the TF plane as an image matrix which will replace the 2D data distribution i f (u, v) in the formulation of the RT given in (18). Let I be a K × N image matrix representation of the TF plane, where its elements I(k, n) represent the gray-level intensities. K is the number of rows corresponding to the number of frequency slots, and N is the number of columns that correspond to the number of time slots in the TFD. K and N vary according to the resolution of the TFD, the time duration, and the signal bandwidth. The formulation of the HRT for discrete data sets is given as follows: where and f (k, n, Θ) = 0 is the parametric constraint equation in the image plane. In general, the implementation of the HRT would require that we first determine a sufficiently fine search grid Θ s , for the parameter space, which will allow us to differentiate all parametric curves of the form given in (21) within the resolution limitations of the image matrix. This search grid functions as a quantized parameter space.
In the implementation of the HRT, the transform value R(Θ 0 ) at some Θ 0 ∈ Θ s contains the total energy in the pixels that satisfy the parametric constraint equation. Therefore, we can devise an HRT-based system to detect directional elements defined by parametric equations: the peak values of the HRT will yield the most likely parameter values. Figure 4 provides an overview of the proposed algorithm. After the interference excision, the "interference suppressed" SS signal is processed as before by first correlating with the synchronized PN sequence, integrating the resulting sequence, and estimating the transmitted data symbols using a threshold detector.

The excision algorithm
For the interference excision algorithm, we assume that the information on the number and type of interference signals is available. In particular, we assume that the interference signals are linear or quadratic FM signals which can be present simultaneously. Let τ ∈ {linear, quadratic} be the type of interference, and let M τ be the number of interference signals of type τ.
Step 1. The received signal r is modeled as a linear combination of Gaussian functions using the MP algorithm given in Section 3.1. Let r be the model generated by MP algorithm as in where g γm (n) are the Gaussian TF functions given in (8). The model order M is determined as the smallest positive integer which will make r = {r (0), r (1), . . .} satisfy the condition where N is the length of r.
Step 2. Formulate the parameter set G using the parameters of the Gaussian TF functions g γm , such that

EURASIP Journal on Wireless Communications and Networking
where k m and p m are the frequency and temporal placement parameters of the TF function g γm , respectively.
Step 3. Compute the cross-term free TFD of r using the WVD: where W gγ m (n, w) is the WVD of the Gaussian function g γm (n).
Step 4. Let I(k, n) be the K × N image matrix representation of W (n, w). For each interference type τ known to be present in the received signal, determine the corresponding quantized parameter space Θ τ and evaluate the HRT given in Section 3.2 and R(Θ τ ) using (19)-(21).
Step 5. For each interference type τ known to be present in the received signal, determine the M τ parameters τ } from the quantized parameter space Θ τ corresponding to the first M τ maxima of R(Θ τ ). Let Step 6. For each interference type τ known to be present in the received signal and for each Θ (m) where f τ (k, p, Θ) = 0 is the parametric constraint describing the interference of type τ and ΔΘ is the empirically determined confidence measure.
Step 7. For each interference type τ known to be present in the received signal and for each m ∈ {1, . . . , M τ }, construct the corresponding interference model as Step 8. Determine the interference excised SS signal by subtracting the interference models generated in Step 7 from the received signal: Before presenting the simulation results, the computational complexity of the proposed technique is briefly discussed. The search algorithm above uses both MP and HRT. Therefore, there are two factors that affect the computation time. For MP, the computation time will depend on how well the dictionary elements (Gaussian functions) given in (8) can match the interference signal. A conclusive comment on the effect of interference types (e.g., linear versus quadratic) on the computational complexity cannot be made as the representation of different parametric equations with Gaussian functions being not studied and compared in the literature. In [12], the authors used a chirp-based dictionary to reduce the computation time of linear chirps; however, this dictionary would fail to model other parametric equations. Also, it should be noted that the computation time will linearly increase with the number of interference signals. For HRT, the computation complexity of the search algorithm is O(N 2 ). To decrease the computational complexity of the HRT technique, variants of HRT such as the randomized Hough transform [25] can be considered. Also, a further study can be conducted on the optimization of HRT search algorithm, which is beyond the scope of this paper. In the following, the proposed excision algorithm is evaluated.

SIMULATION RESULTS AND DISCUSSION
The simulation results presented in this section are based on a DSSS system with L = 128 chips per message symbol b k . The transmitted message contained 100 message symbols. We assumed that the channel was nondispersive, and the received signal and the PN sequence were synchronized.

Performance measures
Bit error rate (BER). For the DSSS model used in this study, we process the received signal using the interference excision algorithm, and estimate the transmitted message symbols using the detector structure presented in Section 2. A comparison of the estimated message symbols { b k } with {b k }, and expressing the number of erroneous estimates as a percentage of the total number of message symbols yields the bit error rate. Chip error rate. We define the chip error as for n ∈ {0, . . . , L − 1} and k ∈ {1, 2, . . . }.

BER performance
To measure the performance of the DSSS system using the new interference excision algorithm developed in this paper, we evaluated the BER results for the following three interference scenarios, where we assumed the presence of (i) a single-component linear chirp, (ii) a single-component quadratic chirp, (iii) a multicomponent interference with linear and quadratic chirps.
The interferences were measured with JSR values in the range of 0 to 50 dB at 10 dB steps. We assumed the SNR to be 10 dB in each case. We suppressed the interference before despreading, using the proposed interference excision algorithm. We observed zero bit errors in all cases after the excision of single-component and multicomponent interferences. We repeated the same process for different SNR values in the range of −10 dB to 10 dB, and also recorded zero bit errors. One of the main reasons for having zero BER in these simulation runs is the accurate TF representation of interferences, and the successful detection by the HRT. A similar S. Krishnan and S. Erküçük  observation was made by Bultan et al. in [12], where they represent linear interferences with good TF localization using adaptive chirplet decomposition. However, they do not report any results on the excision of quadratic and/or multicomponent interferences. Other TFD-based methods reported bit errors for similar excision conditions [8,9,11].
As an example we plot in Figure 5 the MP TFD 1 of a multicomponent interference consisting of linear and quadratic chirps at JSR = 40 dB.

Chip error rate performance
We evaluated the DSSS system by calculating the percentage of chips received in error at various SNR values. Figures 6 and  7 show the simulation results for calculating the chip error rates for the JSR values 40 dB and 5 dB, respectively. Figure 6 shows the percentage of chips in error before and after the excision of single and multicomponent interferences. The JSR value of 40 dB is used because at this JSR level with no interference excision, the system BER is approximately 50 percent indicating that the system cannot suppress any part of the interference. It is observed that the excision of a single interference results in less chip error rate than the excision of a multicomponent interference. This is a result of the excision of a multicomponent interference introducing more noise than the excision of a single-component interference at the same power level. When the estimates of the interferences are excised from the SS signal, part of the SS signal in the vicinity of the interference localization are also suppressed. Therefore, multicomponent interferences are likely to introduce more residual noise. 1 Although the SS signal is also partly decomposed and represented on the TF plane, its lower energy compared to interferences makes it invisible on the MP TFD obtained by the WVD.    Figure 7 shows the results for the same experimental setup at JSR = 5 dB. JSR value of 5 dB is used since the system can suppress the interference partly without interference excision prior to despreading. In systems proposed by other researchers [8], the excision of the low-power interference degrades the performance of the system. In the case of the new interference excision algorithm we developed in this paper, the proposed system has substantially improved the chip  error rate. The results obtained show that the excision of the single interference results in less chip error rate, consistent with the results plotted in Figure 6.
For illustration purposes, in Figure 8 we provide the TFDs of the SS signal with a single interference (JSR = 5 dB) (left plot), the interference detected by HRT (middle plot), the interference excised SS signal (right plot).
The simulation results show that the proposed technique can be successfully used for excision of single-component and multicomponent chirp-like interferences using adaptive TFDs and the HRT. The new algorithm suppresses the interferences, while introducing an acceptable magnitude of noise, which can be overcome by the spreading gain. This results in zero bits in error after interference excision.

CONCLUSIONS
In this paper, we proposed a new interference excision algorithm and evaluated its performance in terms of the BER and chip error rates. The most striking observation resulting from the simulation studies is that there were no bit errors after the excision of single and multicomponent interferences at all JSR levels tested, that is, JSR ≤ 50 dB. Under similar test conditions, the algorithms developed in earlier studies reported bit errors with the notable exception of [12].
This highly desirable characteristic is the result of the following three factors.
(1) The model of the interference uses Gaussian functions, which provide optimal TF resolution within the limits of the uncertainty principle.
(2) The MP TFD uses WVD, which also localizes the components well and provides a high TF resolution. The modeling of the interferences as a linear combination of basis functions eliminated the cross-terms in the construction of the TFD for multicomponent interferences. Lack of cross-terms prevents undesired peaks in the HRT space, which may lead to incorrect parameter estimates.
(3) The HRT algorithm acts as an adaptive thresholding mechanism successfully determining the functions that model the interference.
Another important conclusion of this paper is that the proposed algorithm can excise any interference that can be modeled using a parametric equation. This is a consequence of the HRT being able to detect any directional element defined by parametric constraints. Other algorithms focus only on linear, or sinusoidal interferences, and break down, if there are nonlinear and/or multicomponent interferences.