Je tenais tout d'abord à remercier Dr Rodolphe Weber pour m'avoir permise d'effectuer ce travail sous sa direction, pour ses conseils, pour sa patience et pour son temps, et Dr Albert Jan Boonstra, pour sa disponibilité et son aide tout au long du parcours, et surtout pour son excellente thèse, qui fut ma bible pendant 3 ans. Sans vous deux jamais ce travail n'aurait été possible. Et même si aujourd'hui je suis passée dans le « camp ennemi », je ne me fais aucun soucis pour la radioastronomie, car tant qu'il y aura des personnes comme vous avec vos inestimables qualités techniques et surtout humaines, il sera toujours possible d'observer les astres, même si le spectre deviendra de plus en plus pollué. En espérant que vous pourriez faire profiter plus de personnes de votre encadrement de qualité, de votre dévouement au travail et de votre grandeur d'âme.

Ensuite, j'adresse mes plus vifs remerciements aux rapporteurs, Pr Van Ardenne et Pr Gosselin, d'abord pour votre temps, ensuite pour la qualité de vos rapports, car vos remarques et suggestions ont été extrêmement précieuses. Je remercie également Pr Dudock De Wit pour avoir accepté de présider le jury de ma soutenance.

Mes sincères remerciements vont à TOUS les membres d'ASTRON, qui ont fait de mon séjour une expérience inoubliable. Je remercie particulièrement Mme et Mr Noordam, car il est rare de trouver des personnes avec de si grandes qualités humaines, prêts à ouvrir leur maison pour accueillir les passants du moment, qui viennent de différents horizons, autour d'un bon bridge, et de grandes discussions. Sans oublier Mme Hugen, ma propriétaire. Même si la plupart de nos conversations se faisaient par des signes, elle m'a montré par sa gentillesse et par sa chaleur la vraie nature des Hollandais, un pays et un peuple désormais très chers à mon coeur.

Je remercie Mr Jennane d'avoir dirigé ma thèse, Mr Harba de m'avoir donné de m'avoir offert l'opportunité de travailler sur ce sujet au LESI, ainsi que Mme Mounaim Rousselle et Mr Leger. J'adresse ma reconnaissance à TOUS mes collègues d'ISS Prisme, sans oublier les anciens du LESI dont Mme Feng, pour l'ambiance au labo, pour les fous rires et pour surtout pour leur aide tout au long de ces trois années. Je tenais particulièrement à remercier Mme Elizabeth Rowley Jolivet, pour la correction de ma thèse. Merci pour vos précieuses remarques, vos explications était si claires que ça devenait agréable et ludique de reprendre vos corrections. Ma reconnaissance et mon amitié vont spécialement à l'équipe signal, à Mlle Ait Allal, pour sa bonne humeur et son soutien, en lui souhaitant beaucoup de succès pour la suite de ses travaux, à Mr Ravier pour ses conseils et son suivi, et à Mme Jabloun pour son temps, même si je n'ai pas eu le plaisir de rester longtemps après sa venue.

Enfin, je remercie ma famille qui a cru en moi, et qui m'a apporté toute l'aide dont j'avais besoin pour en arriver là, je ne serai jamais arrivée à ce diplôme sans vos encouragement, votre confiance et votre soutien mes chers parents, j'espère avoir été à la hauteur de vos attente. Mes remerciements vont ensuite à ma deuxième famille, mes chers amis, qui m'ont aidée à surmonter les moments difficiles, et sans qui, je ne sais pas si ce travail aurait aboutit un jour : vous n'hésitiez jamais à être là pour moi que ce soit pour le soutien morale, l'accueil le week end, l'aide sur LaTex, et les articles IEEE, etc...

Introduction

Les astres ont toujours été une source de crainte, de passion, d'inspiration et de quête pour l'Homme à travers les âges. Cette quête de connaissance de l'inconnu avait, à l'antiquité, pris forme de science de prédictions écliptiques, telle que la prédiction des éclipses de Vénus par les civilisations Maya, ou la prédiction des éclipses solaires et lunaires à l'époque grecque alexandrine. Les premières théories sur l'Univers débutèrent par la théorie géocentrique, défendues par Ptolémée, Aristote et Hipparque. Cette théorie établit un modèle physique de l'Univers centré sur la Terre, qui est immobile, tandis que les autres astres (dont le Soleil et la Lune) ont des trajectoires circulaires autour de ce centre. La théorie géocentrique fut, par la suite, contredite par la théorie héliocentrique, évoquée au 13 e siècle par le savant Al-Toosi de l'observatoire de Maragheh (Perse). Après des mesures des trajectoires des astres à l'aide de l'Astrolabe (un des premiers instruments de mesures astronomiques), il déduit que la Terre tourne autour du Soleil. Cette théorie héliocentrique fut reprise par Copernic au 16 e siècle, puis complétée par Kepler et Galilée au début du 17 e siècle. D'ailleurs, c'était l'amélioration apportée par Galilée aux lunettes astronomiques (qui se comportent comme une longue vue marine, mais avec un facteur de grossissement plus important)en 1609 qui a permis d'obtenir les premières précisions sur la théorie héliocentrique. Entre 1684 et 1687, Newton a reussi à établir, d'une part, les lois gravitationnelles qui appuyaient cette théorie, et d'autre part, il a inventé le télescope réflecteur avec des miroirs réfléchissants qui a permis d'améliorer davantage la qualité des observations, et d'ouvrir une fenêtre optique pour l'étude des astres.

Longtemps après, Karl Jansky découvrit en 1932 l'existence d'émissions radios provenant de la Voie Lactée qui brouillaient les communications transatlantiques. Ce fut ainsi que commença une nouvelle ère d'observations astronomiques via les radiofréquences (figure 1). En effet, le premier radiotéléscope, proprement dédié à la réception de signaux cosmiques, a été construit en 1938 par Grote Reber. Ce dernier a, par ailleurs, entrepris une série de travaux de recherche qui ont permis :

• la confirmation de l'observation faite par Jansky,

• l'établissement d'une carte radiofréquence du ciel,

• la découverte d'un nombre important de signaux cosmiques à basse énergie, ce qui a permis en 1950 de différencier entre les émissions thermiques et les synchrotrons 1 .

Cette dernière découverte a été le pas précurseur dans l'expansion et la multiplication des constructions des radiotéléscopes à travers le monde à l'issue de la 2 nde guerre mondiale. A Introduction titre d'information, le plus grand radiotéléscope de nos jours, nommé l'Arecibo radio telescope, se trouve au Porto Rico. Mais pour répondre aux besoins scientifiques de plus en plus grandissants, des instruments plus sensibles et plus performants ont dû être envisagés, notamment l'utilisation des réseaux phasés d'antennes, essentiellement en interférométrie.

Figure 1: Fenêtres d'observation astronomique dans le spectre du rayonnement électromagnétique. Les télescopes optiques permettent d'observer dans la fenêtre du visible et de l'infrarouge. Les radiotélescopes permettent d'étudier les astres dans le domaine radio.

En 1993 [START_REF]Square kilometer array[END_REF], lors d'une réunion de l'Union des Radio Sciences Internationale (URSI), un groupe de travail international fut créé pour le développement des grands télescopes. Le 10 août 2000, il fut décidé la construction du plus grand radiotélescope du monde de par sa surface collective (de l'ordre du kilomètre carré), et le plus performant de par sa fenêtre d'observation (de 70 MHz à 10 GHz), connu sous le nom de Square Kilometer Array (SKA) (figure 2). Les tâches scientifiques qu'aura à traiter SKA seront [START_REF] Schilizzi | Preliminary specifications for the square kilometer array[END_REF]:

• L'époque de réonisation,

• l'évolution des galaxies, la Cosmologie et l'énergie noire,

• l'origine et l'évolution du magnétisme cosmique,

• tests des champs puissants de gravité à l'aide des Pulsars et des trous noirs,

• étude de l'origine de la vie, Figure 2

: Vue d'artiste du radiotélescope Square Kilometer Array (SKA). Cette image représente la zone centrale du radiotélescope (5 à 10 km de diamètre). Au premier plan, le réseau d'antennes phasé, observant dans la bande 70MHZ-1200 MHz. Il y aura 250 stations de ce type dans un rayon de 200km. Arriere plan : antennes paraboliques de 12m-15m de diamètre, fonctionnant dans la bande de 700MHz-10GHz. Le foyer pourra être un récepteur simple ou un réseau de récepteurs positionnés sur le plan focal. Il y aura environ 1200 paraboles dans la zone centrale et 1200 autres réparties sur des bras spiralés jusqu'à une distance 3000 km. Les deux sites potentiels sont l'Australie et l'Afrique du Sud. Avec l'autorisation de SKA telescope organisation.

Introduction

• et l'exploration de l'inconnu ... Les tâches d'études et de développement ont été réparties sur plusieurs entités : Les Etats Unis, L'Europe, et le reste du Monde (l'Australie, le Canada, l'Afrique du Sud, la Chine et l'Inde). Les partenaires européens se sont regroupés dans le projet Square Kilometer Array Design Studies (SKADS) [START_REF]Square kilometer array desing studies[END_REF]. Ce projet regroupe les Pays Bas (représentés principalement par ASTRON et l'Université de Groningen), la Grande Bretagne (les Universités de Cambridge, de Manchester et d'Oxford), la Suède (l'Université de Chalmers), l'Allemagne (l'Institut Marx Planck pour la Radio Astronomie), l'Italie (l'Observatoire de Bologne), l'Espagne (l'Observatoire National Astronomique d'Espagne) et la France (représentée par l'Observatoire de Paris et l'Université d'Orléans).

L'Université d'Orléans, de par cette thèse, participe à l'amélioration de l'élimination de brouilleurs. En effet, plus un instrument est sensible, et plus il sera vulnérable face aux interférences. De plus, les dernières décennies connaissent une explosion sans précédent des moyens de télécommunications. En outre, des développements en télécommunications seront mis en oeuvre dans un futur proche tels que l'avènement des télévisions et radios numériques, ou la possibilité de téléphoner et même d'avoir des connexions à internet depuis les transports aériens. Il faudrait alors des méthodes de suppressions de brouilleurs plus robustes, et une différentiation plus précise entre les signaux utiles (les signaux cosmiques) et les signaux d'origine entropique (figure 3). Le projet SKADS se base sur l'étude de propositions possibles à apporter au projet SKA. Dans ce but, cette thèse proposera des méthodes exploitant la propriété de cyclostationnarité présente sur la plupart des brouilleurs. En effet, selon [START_REF] Rohlfs | Tools of radio astronomy[END_REF], les signaux cosmiques sont une somme d'émissions indépendantes, et en utilisant le théorème de centrale limite, ces signaux sont gaussiens, blancs et stationnaires. La plupart des signaux de télécommunications, sont, quant à eux, cyclostationnaires, offrant une possibilité de différenciation entre les signaux de télécommunications et les signaux utiles (signaux cosmiques), lors du traitement spatial des interférences (pour des radiotélescopes à réseaux phasés, tels que SKA)).

Cette thèse est une contribution à l'amélioration des observations pour les radiotélescopes à réseaux phasés en présence d'interférences. L'originalité de cette thèse repose sur l'utilisation de la séparation spatiale entre les sources cosmiques et les brouilleurs issus des télécommunications en se basant sur la cyclostationnarité de ces derniers. Cette thèse s'inscrit dans le cadre du projet SKADS décrit ci-dessus, elle s'appuie sur des techniques de traitement spatial existantes, et qui seront adaptées au contexte radioastronomique. Les performances des techniques proposées seront évaluées à travers des simulations sur des données synthétiques et/ou réelles.

Le premier chapitre abordera les principes généraux de la cyclostationnarité. Des exemples seront donnés sur les brouilleurs observés classiquement en radioastronomie.

Dans le second chapitre, un modèle général de signaux reçus dans un radiotélescope à réseau phasé sera proposé. Les matrices de corrélations dans l'approche classique et cyclostationnaire seront définies. La décomposition en sous-espaces signal et bruit des matrices de corrélation étant une étape clé dans le traitement multidimensionnel des signaux, elle sera expliquée pour l'approche classique et cyclostationnaire.

Le troisième chapitre étudiera les techniques principales existantes en télécommunications pour le traitement d'antennes. Plus spécialement les méthodes basées sur la cyclostationnatrité Figure 3: Spectrogramme de la bande [START_REF] Dandawate | Nonparametric polyspectral estimators for kth-order (almost) cyclostationary processes[END_REF]

-300 MHz de l'environnement du radiotélescope LOFAR (annexe A.1). Certaines interférences ont disparu du spectre (quelques chaînes TV). Néanmoins, les bandes vides actuelles (170-230 MHz) seront progressivement occupées par des transmissions issues des radio et TV numériques (DAB et DVB). Avec l'autorisation de Albert-Jan BOONSTRA et ASTRON.

Introduction seront expliquées afin de définir des pistes de développement de méthodes de suppression d'interférences cyclostationnaires pour notre application.

La détection d'interférences est une technique permettant d'identifier et de masquer les observations brouillées. Trois méthodes de détection d'interférences cyclostationnaires seront présentées dans le quatrième chapitre, et seront comparées dans leurs performances aux méthodes de détection classiques existantes en radioastronomie.

Dans le cinquième chapitre, une méthode d'estimation et de soustraction d'interférences cyclostationnaires sera présentée. Cette technique permet d'éliminer le brouilleur, tout en limitant les distorsions introduites sur l'observation. Des simulations seront proposées sur des signaux synthétiques et réels.

Enfin, le filtrage spatial est une méthode générale pouvant être appliquée dans quasiment tous les cas d'observation. Son algorithme sera présenté dans le sixième chapitre avec des exemples de traitement d'interférences synthétiques et réelles en utilisant des observations issues du radiotélescope LOFAR (cet instrument précurseur de SKA est décrit en annexe A.1).

• Introduction

The present thesis is a contribution to observation improvements for phased array radio telescopes, in presence of radio frequency interferers (RFIs). The originality of the study is the use of the cyclostationarity property, in order to improve the spatial separation between cosmic sources and telecommunication signals. This thesis is part of the square kilometre array design study (SKADS) project [START_REF]Square kilometer array desing studies[END_REF]. The proposed approaches are based on existing spatial processing techniques, which will be adapted to the radio astronomy context. The performance of the techniques presented will be evaluated through simulations on synthetic and/or real data A distinction between cosmic signals and man-made signals can be achieved using cyclostationarity. This property is due to the redundancy of the periodic characteristics involved in the modulation schemes of the telecommunication signals (carrier frequency, baud rate and coding scheme). The first chapter will present the general principles of cyclostationarity. Examples on common interferers encountered in radio astronomy will be given.

In the second chapter, a general model of received signals in phased array radio telescopes will be proposed. Correlation matrices will be defined within respectively the classic and the cyclostationary approaches. Decomposition of the correlation matrices into signal and noise subspaces, which is a key step in multidimensional signal processing, will be explained for both classic and cyclostationary approaches.

The third chapter is an overview on the main existing array processing techniques. More specifically, cyclostationarity-based methods will be explained, in order to define RFI mitigation improvements for our application.

Interference detection techniques allow the identification and the blanking of corrupted observations. Three RFI detection algorithms based on cyclostationarity will be presented in the fourth chapter, and will be compared to classic detection techniques in radio astronomy.

In the fifth chapter, an estimation and subtraction technique of cyclostationary RFIs will be presented. This technique allows RFI suppression, with limited distortions on the observation. Simulations will be proposed on synthetic and real data. Finally, spatial filtering, which is a general method that can be applied on almost all observation cases, will be presented in the sixth chapter. Examples of synthetic and real RFI processing will be presented.

Chapter 1

General definitions of the cyclostationarity property

Most telecommunication signals present a hidden periodicity due to the periodic characteristics involved in the signals construction related to e.g.carrier frequency, baud rate, coding scheme...). These parameters are usually scrambled and hidden by the randomness of the message to be transmitted. However, by using a cyclostationary approach, this hidden periodicity can be recovered, thus making identification of the telecommunications signal possible.

Cosmic signals, in contrast, according to [START_REF] Rohlfs | Tools of radio astronomy[END_REF], are considered to be Gaussian centered, white and stationary. Thus, the separation criterion will be the cyclostationarity assumption of the radio frequency interferences (RFIs) against the stationary assumption of the cosmic sources, called the signals of interest (SOI). This chapter will present the cyclostationarity concept.

The earliest explanation of cyclostationarity in signal processing appeared in Gardner's book [START_REF] Gardner | Introduction to random processes with applications to signals and systems[END_REF]. In [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF], Gardner presented his nonprobabilistic statistical theory of cyclostationarity, while Giannakis [START_REF] Dandawate | Statistical tests for presence of cyclostationarity[END_REF] [START_REF] Dandawate | Nonparametric polyspectral estimators for kth-order (almost) cyclostationary processes[END_REF] [START_REF] Dandawate | Asymtotic theory of mixed time averages and kth-order cyclic-moment and cumulant statistics[END_REF], defined cyclostationarity within the framework of stochastic processes. In other words, Gardner adopted the infinite time averaging notation < . > ∞ instead of the classic expectation operation E {.}. In accordance with Gardner, we will proceed in this thesis with the time-space theory using the notation < . > ∞ instead of the probability-space theory. At certain points in the study, however, we will employ E {.} when necessary. In this section we will also consider the asymptotical case, represented by infinite time averaging. An exhaustive overview on cyclostationarity theory and applications can found in [START_REF] Serpedin | Bibliography on cyclostationarity[END_REF] [START_REF] Gardner | Cyclostationarity: Half a century of research[END_REF].

In the first section, we will describe the cyclostationarity principle and will derive the different mathematical tools used to express this property. In the second section, the cyclostationarity will be illustrated, for common telecommunication interferences. In the third part of this chapter, we will present a practical technique to retrieve the hidden periodicities mentioned above.

Cyclostationarity principle

This section will define cyclostationarity for readers unfamiliar with the concept, since all the document is based on this principle.

CYCLOSTATIONARITY PRINCIPLE

Consider the autocorrelation function of a signal r(t) defined by:

R r,r * (t, τ ) ≜ r(t + τ 2 )r * (t - τ 2 )
∞ .

(1.1)

• If R r,r * (t, τ ) = Rr, r * (τ )
, the signal is said to be wide-sense stationary,

• If R r,r * (t + T 0 , τ ) = R r,r * (t, τ ), then R r,r * (t, τ ) exhibits a second order periodicity, and the process is said to be wide-sense cyclostationary. We define the corresponding fundamental frequency and its harmonics as cyclic frequencies. These cyclic frequencies are noted α = k/T 0 , with k ∈ ℤ.

• Moreover, if there is more than one fundamental cyclic frequency, the signal is called polycyclostationary.

Figure 1.1 illustrates the cyclostationarity principle for a binary signal. In particular, it shows that telecommunications signals can be modeled as stationary or cyclostationary by respectively randomizing or not randomizing the phase.

(b) Autocorrelation representation when the temporal origin is seen as a random parameter: the process is stationary. (c) Autocorrelation representation when the temporal origin is constant: the process becomes cyclostationary.

Consider a discrete-time signal r(t). If r(t) is a cyclostationary signal, then R r,r * (t, τ ) is periodic with period T 0 and can be decomposed into its Fourier series:

R r,r * (t, τ ) = +∞ α R α r,r * (τ )e j2παt , with α = k T 0 and k ∈ ℤ, (1.2)
where R α r,r * (τ ) are the Fourier coefficients defined by:

R α r,r * (τ ) ≜ R r,r * (t, τ )e -j2παt T 0 . (1.3)
As R α r,r * is a cyclic (sinusoidal) weighting of the autocorrelation function, with the factor e -j2παt , we call R α r,r * the cyclic autocorrelation function. Moreover, the cyclic autocorrelation function R α r,r * can be interpreted [START_REF] Gardner | Cyclostationarity in communications and signal processing[END_REF] as the crosscorrelation between two signals v(t) and w(t), representing the signal r(t) with respectively frequency-shifts of α 2 and -α 2 :

R α r,r * (τ ) = R v,w * (t, τ ) with v(t) = r(t)e -jπαt
and w(t) = r(t)e jπαt .

(1.4)

Thus, a more generalized definition of the cyclic correlation function can be expressed with:

R α r,r * (τ ) = r(t + τ 2 )r * (t - τ 2 )e -j2παt ∞ . (1.5)
Note that if α = 0 or widesense stationary, we obtain R α r,r * (τ ) = R r,r * (τ ). In some cases (see section 1.2), for instance amplitude modulation (AM) signals, the cyclic conjugated correlation function provides another cyclostationary representation of the signal r(t). It is defined by:

R α r,r (τ ) = r(t + τ 2 )r(t - τ 2 )e -j2παt ∞ . (1.6)
From these definitions, we can see that R α r,r * (τ ) and R α r,r (τ ) are two interesting tools to discriminate between cyclostationary and stationary signals. Indeed:

• If R α
r,r * (τ ) and R α r,r (τ ) are zero for α = 0 then r(t) is stationary. In other words, for cosmic signals R α r,r * (τ ) and R α r,r (τ ) = 0.

• If ∃α = 0, such as R α r,r * (τ ) = 0 or R α r,r (τ ) = 0, then r(t) is cyclostationary.

Note that in some papers [START_REF] Gardner | Spectral correlation of modulated signals: Part I -analog modulation[END_REF] [40], the authors prefer the use of the cyclic or the cyclic conjugated spectrum, which are just the Fourier transform of the cyclic and the cyclic conjugated autocorrelation functions (figure 1.1).

The cyclostationary properties of some modulations , considered as the most common interferences, will be derived in the next section.

Cyclostationary properties of some common RFIs

In this section, the statistics of a non-exhaustive list of telecommunications signals will be studied, in order to illustrate the cyclostationarity property defined in the previous section. In subsections 1.2.1, 1.2.2 and 1.2.3, we will first define some calculation tools, then will apply the results on the modulation schemes.

Amplitude modulation

Consider a carrier signal r(t) = e j(2πf 0 t+ϕ 0 ) where f 0 and φ 0 are respectively the carrier frequency and phase.

The cyclic correlation function is:

R α r,r * (τ ) = e j(2πf 0 (t+ τ 2 )+ϕ 0 ) e -j(2πf 0 (t-τ 2 )+ϕ 0 ) e -j(2παt) ∞ = e j2πf 0 τ e -j(2παt) ∞ =      0 if α = 0, e j2πf 0 τ if α = 0.
(1.7)

The cyclic conjugated autocorrelation function of r(t) is:

R α r,r (τ ) = e j(2πf 0 (t+ τ 2 )+ϕ 0 ) e j(2πf 0 (t-τ 2 )+ϕ 0 ) e -j2παt ∞ = e j2ϕ 0 e j2π(2f 0 -α)t ∞ =      0 if α = 2f 0 , e j2ϕ 0 if α = 2f 0 .
(1.8)

Let us now consider r(t) to be the product of two independent signals w(t) and x(t): r(t) = w(t)x(t). The cyclic autocorrelation function of this product is given by ( [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF] pages 400 and 401):

R α r,r * (τ ) = β R α-β x,x * (τ )R β w,w * (τ ).
(1.9)

A similar formula can be derived for the cyclic conjugated autocorrelation function:

R α r,r (τ ) = β R α-β x,x (τ )R β w,w (τ ).
(1.10)

• Application to amplitude modulation

If a(t) ∈ ℝ is a stationary signal, and f 0 and φ 0 are respectively the carrier frequency and phase, the expression of an amplitude modulated signal is:

r(t) = a(t) e j(2πf 0 t+ϕ 0 ) . x(t) (1.11)
According to equation (1.9):

R α r,r * (τ ) = β R α-β x,x * (τ )R β a,a (τ ). (1.12)
Since a(t) is stationary, the cyclic autocorrelation function of a(t) is:

R β a,a (τ ) =      0 if β = 0, R a,a (τ ) if β = 0, (1.13)
where R a,a (τ ) is the autocorrelation function of a(t).

Considering equations (1.13) and (1.7), equation (1.12) becomes:

R α r,r * (τ ) = R α x,x * R a,a (τ ) =      0 if α = 0, e j2πf 0 τ R a,a (τ ) if α = 0. (1.14)
Similarly, the cyclic conjugated correlation function of an AM modulated signal will be:

R α r,r (τ ) = R α x,x R a,a (τ ) =      0 if α = 2f 0 , e j2ϕ 0 R a,a (τ ) if α = 2f 0 .
(1.15)

Since R α r,r * (τ ) is zero for all values of α except when α = 0, no cyclostationary property can be stated in this way. However, the cyclic conjugated autocorrelation function reveals that R α r,r (τ ) is non zero for the cyclic frequency α = 2f 0 , and thus, r(t) is cyclostationary. Figure 1.2 shows the cyclic conjugated autocorrelation function and its spectrum, for an amplitude modulated signal.

Binary phase shift keying

Consider the baseband digital signal q(t) given by:

q(t) = k a k h(t -kT sym -t 0 ) with k ∈ ℤ, (1.16)
where a k ∈ ℝ are independent values chosen among a set of centered fixed levels with corresponding power equal to σ 2 a , h(t) is the impulse response of the emission filter, T sym is the symbol duration and t 0 is an unknown time delay. The autocorrelation function of q(t) is given by:

R q,q * (t, τ ) = q(t + τ 2 )q * (t -τ 2 ) ∞ = σ 2 a k h(t + τ 2 -kT sym -t 0 )h * (t - τ 2 -kT sym -t 0 ).
(1.17)

From (1.17), it is straightforward to see that R q,q * (t + T sym , τ ) = R q,q * (t, τ ). Thus, q(t) is cyclostationary with cyclic frequency 1/T sym . As R q,q * (t, τ ) is cyclic, it can be decomposed into its Fourier series:

R q,q * (t, τ ) = l∈ℤ R l Tsym q,q * (τ )e j2π l Tsym t , (1.18)
with:

R 1 Tsym q,q * (τ ) = l T sym Tsym R q,q * (t, τ )e -j2π l Tsym t dt = σ 2 a T sym Tsym k h(t + τ 2 -kT sym -t 0 )h * (t - τ 2 -kT sym -t 0 )e -j2π l Tsym t dt = σ 2 a T sym +∞ -∞ h(t + τ 2 -t 0 )h * (t - τ 2 -t 0 )e -j2π l T t dt = σ 2 a T sym e -j2π l Tsym t 0 +∞ -∞ h(t + τ 2 )h * (t - τ 2 )e -j2π l Tsym t dt r l Tsym h,h * (τ ) (1.19)
Consequently, the cyclic autocorrelation function will be:

R α q,q * (τ ) =        σ 2 a Tsym e -j2παt 0 r α h,h * (τ ) if α = l Tsym with l ∈ ℤ, 0 if α = l Tsym .
(1.20)

Another demonstration can be found in [START_REF] Gardner | Exploitation of spectral redundancy in cyclostationary signals[END_REF], where the author has employed a time averaging approach in order to prove the result given by equation (1.20).

• Application to BPSK modulation

Binary phase shift keying is a modulation that uses two phases which are separated by 180 • and so can also be termed 2-PSK (i.e. a k in equ. (1.16) takes the values +a or -a with a ∈ ℝ). The expression of a BPSK modulated signal is:

r(t) = q(t) e j2πf 0 t+jϕ 0 . x(t) (1.21)
Using equation (1.9), the cyclic autocorrelation function will be as follows:

R α r,r * (τ ) = β R α-β x,x * (τ )R β q,q * (τ ). (1.22)
As shown with equations (1.20) and (1.7), the cyclic autocorrelation function of q(t) and x(t) will be:

R β q,q * (τ ) =        σ 2 a Tsym e -j2πβt 0 r β h,h * (τ ) if β = l Tsym with l ∈ ℤ, 0 if β = l Tsym .
(1.23)

R γ x,x * (τ ) =      0 if γ = 0, e j2πf 0 τ if γ = 0 ⇒ α = β.
(1.24)

From equations (1.23) and (1.24), we deduce that the cyclic autocorrelation function of a BPSK modulated signal will be:

R α r,r * (τ ) =        σ 2 a Tsym e -j2παt 0 e j2πf 0 τ r α h,h * (τ ) if α = l Tsym , 0 if α = l Tsym .
(1.25)

Figure 1.3 shows the spectrum of the cyclic autocorrelation function for a BPSK modulated signal.

In the same way, we evaluate both R β q,q (τ ) and R γ x,x (τ ) in order to deduce the cyclic conjugated autocorrelation function of BPSK modulated signals:

R β q,q (τ ) =        σ 2 a Tsym e -j2πβt 0 r β h,h (τ ) if β = l Tsym with l ∈ ℤ, 0 otherwise. (1.26) R γ x,x (τ ) =      0 if γ = 2f 0 , e j2ϕ 0 if γ = 2f 0 ⇒ α = 2f 0 + β. (1.27)
Consequently, the cyclic conjugated autocorrelation function of a BPSK modulated signal will be:

R α r,r (τ ) =        σ 2 a Tsym e j2ϕ 0 e -j2π(α-2f 0 )t 0 r l Tsym h,h (τ ) if α = 2f 0 + l Tsym , 0 otherwise .
(1.28)

Figure 1.4 shows the spectrum of the cyclic conjugated autocorrelation function of a BPSK modulated signal.

M-phase shift keying and quadrature amplitude modulations

Let q(t) be the signal with the following expression:

q(t) = k c k h(t -kT sym -t 0 ), (1.29)
where c k = a k + jb k , a k and b k are independent values chosen among the same set of centered fixed levels, with corresponding power of

σ 2 a = σ 2 b = σ 2 c /2, h(t)
is the impulse response of the emission filter, T sym is the symbol duration and t 0 is an unknown time delay.

The same approach as detailed in section 1.2.2 will be used, however:

c k c * l ∞ = a k a l ∞ + b k b l ∞ + j [ b k a l ∞ -a k b l ∞ ] =      0 if k = l, σ 2 c if k = l.
(1.30)

Therefore, the autocorrelation function of q(t) will be:

R q,q * (t, τ ) = σ 2 c T sym k h(t + τ 2 -kT sym -t 0 )h * (t - τ 2 -kT sym -t 0 ) (1.31)
Thus, the cyclic autocorrelation function will be: .32) where:

R α q,q * (τ ) =        σ 2 c Tsym e -j2παt 0 r α h,h * (τ ) if α = l Tsym with l ∈ ℤ, 0 otherwise . ( 1
r l Tsym h,h * (τ ) = +∞ -∞ h(t + τ 2 )h * (t - τ 2 )e -j2π l
Tsym t dt.

To calculate the cyclic conjugated autocorrelation function, we have:

c k c l ∞ = a k a l ∞ -b k b l ∞ + j [ b k a l ∞ + a k b l ∞ ] = 0 ∀ k, l, (1.33)
therefore: R q,q (t, τ ) = 0 ⇒ R α q,q (τ ) = 0 ∀ α.

(1.34)

• Application to the m-PSK and QAM modulations

The statistics of the quadrature amplitude modulations (QAM) are similar to the phase shift keying modulations (m-PSK). A m-PSK or QAM modulated signal expression is:

r(t) = q(t) e j2πf 0 t+jϕ 0 . x(t) (1.35)
Using equation (1.9), the cyclic autocorrelation function will be given by:

R α r,r * (τ ) = β R α-β x,x * (τ )R β q,q * (τ ).
(1.36)

CHAPTER 1. GENERAL DEFINITIONS OF THE CYCLOSTATIONARITY PROPERTY Similarly to section (1.2.2), we obtain:

R β q,q * (τ ) =        σ 2 c Tsym e -j2πβt 0 r β h,h * (τ ) if β = l Tsym with l ∈ ℤ, 0 otherwise. (1.37) R γ x,x * (τ ) =      0 if γ = 0, e -j(2πγτ ) if γ = 0 ⇒ α = β.
(1.38)

From equations (1.37) and (1.38), the cyclic autocorrelation function will be:

R α r,r * (τ ) =        σ 2 c Tsym e -j2παt 0 e -j(2πf 0 τ ) r α h,h * (τ ) if α = l Tsym , 0 otherwise.
(1.39)

R β q,q (τ ) = 0 ∀ β ⇒ R α r,r (τ ) = 0 ∀ α (1.40)

Continuous phase modulations

CPM modulations are generally used in modern satellite and mobile communications due to their the ability to control the power and bandwidth efficiencies and their constant-modulus properties [START_REF] Napolitano | Cyclic spectral analysis of continuous-phase modulated signals[END_REF]. CPM schemes include Gaussian minimum shift keying (GMSK) signals used in GSM communication systems, and continuous phase frequency shift keying (CPFSK) signals which are commonly used in wireless modems.

Continuous phase frequency shift keying

Frequency shift keying (FSK) signals are nonlinearly modulated digital signals [START_REF] Proakis | Digital communications[END_REF] [START_REF] Ferréol | Second-order blind separation of first and second-order cyclostationary sources -application to AM, FSK, CPSK and deterministic sources[END_REF]. CPFSK is the most commonly used variation of FSK modulations. The expression of such modulated m-FSK signals is given by:

r(t) = σ F SK k e j(θ k +2πf d a k (t-kTsym)) h(t -kT sym )e j(2πf 0 t+ϕ 0 ) , ( 1.41)
where: σ 2 F SK is the power of the FSK signal r(t), θ k is the phase of the symbol k, f d is the peak frequency deviation, a k are the transmitted M-ary symbols taking their values in the alphabet ±1, ±3, ..., ±(m -1), T sym is the symbol duration, h(t) is a rectangular pulse of duration T sym and amplitude 1, and f 0 and φ 0 are respectively the carrier frequency and phase. The FSK modulation expression is equivalent to expression (1.41). However, for CPFSK signals, we consider that θ k represents the accumulation of all symbols up to (k -1)T sym and thus:

θ k = 2πf d T sym k-1 q=-∞ a q .
1.2. CYCLOSTATIONARY PROPERTIES OF SOME COMMON RFIS

Figure 1.5: Spectrum of a QPSK (4-PSK) modulated signal cyclic autocorrelation function with a rectangular window. The spectrum is non-zero for cyclic frequencies α = l

Tsym , l ∈ ℤ.

CHAPTER 1. GENERAL DEFINITIONS OF THE CYCLOSTATIONARITY PROPERTY

The CPFSK signals are cyclostationary [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF] [2] [START_REF] Ferréol | Second-order blind separation of first and second-order cyclostationary sources -application to AM, FSK, CPSK and deterministic sources[END_REF] and both the cyclic and cyclic conjugated autocorrelations functions are non zero for:

R α r,r * (τ ) = 0 ⇒ α = l Tsym with l ∈ ℤ R α r,r (τ ) = 0 ⇒ α = 2f 0 + (2l + 1)f d with 0 ≤ l ≤ m-2 2 .
(1.42)

Gaussian minimum shift keying

GMSK is a continuous-phase frequency-shift keying modulation scheme. It is similar to standard minimum-shift keying (MSK); however the digital data stream is first shaped with a Gaussian filter before being applied to a frequency modulator. This has the advantage of reducing sideband power. However, the Gaussian filter increases the modulation memory in the system and causes intersymbol interference. GMSK schemes correspond to a CPM signal with infinite-width pulse, the expression of GMSK transmitted signal is given by [START_REF] Napolitano | Cyclic spectral analysis of continuous-phase modulated signals[END_REF]:

r(t) = 1 √ 2πδT sym e - t 2 2δ 2 T 2 sym ⊗ π 2T sym h( t T sym
)e j(2πf 0 t+ϕ 0 ) , (1.43) where δ = √ ln 2/2πB (B is generally equal to 0.3), f 0 and φ 0 are respectively the carrier frequency and phase and h(t/T sym ) is defined by:

h( t T sym ) ≜      1 if |t| ≤ Tsym 2 , 0 otherwise .
GMSK signals are cyclostationary [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF] [2] [START_REF] Napolitano | Cyclic spectral analysis of continuous-phase modulated signals[END_REF] and both the cyclic and cyclic conjugated autocorrelations functions are non zero for:

R α r,r * (τ ) = 0 ⇒ α = l Tsym with l ∈ ℤ R α r,r (τ ) = 0 ⇒ α = 2f 0 ± 1 2Tsym .
(1.44)

Orthogonal frequency division multiplexing modulations

OFDM modulations belong to the family of multi-carrier transmission schemes. OFDMlike techniques, and especially OFDM/OQAM modulations, are employed in the European digital radio broadcasting (DAB) and digital terrestrial TV broadcasting (DVBT), which are problematic interferers in radio astronomy. The general OFDM/OQAM transmitted signal is expressed by [START_REF] Ciblat | A blind frequency offset estimator for OFDM/OQAM systems[END_REF]:

r(t) =   Q-1 q=0 +∞ l=-∞
x q,l g x (t -lT sym ) + jy q,l g x (t -lT sym + T sym 2 ) e j2π q Tsym (t-lTsym)

  e j2π∆f 0 t , (1.45) which can be interpreted as QAM-type symbols, of duration T sym /Q (T sym is the duration of the total OFDM symbol, and Q is assumed to be known), shaped by the pulses g a (t) and g b (t) (usually square-root raised cosine), translated to the sub carrier q/T sym . ∆f 0 is the carrier frequency offset.

PRACTICAL RETRIEVAL OF CYCLIC FREQUENCIES

For OFDM/OQAM signals, the cyclic frequencies are [START_REF] Bölcskei | Blind estimation of symbol timing and carrier frequency offset in wireless OFDM systems[END_REF] α = {q/Q} q=0,Q-1 and the cyclic conjugated frequencies [START_REF] Ciblat | A blind frequency offset estimator for OFDM/OQAM systems[END_REF] α = {α 0 + q/Q} q=0,Q-1 where α 0 = ∆f 0 mod 1.

Conclusion

Other demonstrations for the following modulations are illustrated in [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF]. Spectral correlations have been detailed for analog modulated signals in [START_REF] Gardner | Spectral correlation of modulated signals: Part I -analog modulation[END_REF], and for digital modulations in [START_REF] Gardner | Spectral correlation of modulated signals: Part II -digital modulation[END_REF]. Detailed developments on the spectrum of the exemplified autocorrelation functions dedicated to radio astronomy can be found in [START_REF] Bretteil | Estimation spectrale en temps réel de sources radioastronomiques en présence de brouilleurs cyclostationnaires: Application au radiotélescope décimétrique de Nançay[END_REF].

It has been shown that most modulation schemes present some cyclostationary properties. The corresponding cyclic frequencies can be obtained through the specifications of the transmitter or can be measured from the data itself. In the next section we will present a practical technique to retrieve these cyclic frequencies, for any kind of modulation scheme.

Practical retrieval of cyclic frequencies

When doing RFI site measurement, it could be interesting to retrieve the modulation characteristics of the RFI from the measured data. The expressions of the cyclic and cyclic conjugated autocorrelation function of a cyclostationary signal, r(t), are given by equations (1.5) and (1.6):

R α r,r * (τ ) = r(t + τ 2 )r * (t -τ 2 )e -j2παt ∞ R α r,r (τ ) = r(t + τ 2 )r(t -τ 2 )e -j2παt ∞ , (1.46)
We derived a finite sample length of these formulas, assuming discrete time signals:

R α r,r * (τ ) = r(t + τ 2 )r * (t -τ 2 )e -j2παt L ≜ 1 L L-1 k=0 r(k + τ 2 )r * (k - τ 2 )e -j2παt
(1.47)

R α r,r (τ ) = r(t + τ 2 )r(t -τ 2 )e -j2παt L ≜ 1 L L-1 k=0 r(k + τ 2 )r(k - τ 2 )e -j2παt
(1.48)

For τ fixed, these equations correspond to the Fourier transform of respectively r(k + τ /2)r * (k -τ /2) and r(k + τ /2)r(k -τ /2) over L samples.

In practice, by choosing τ = 0 and α = k/L, k = 0, .., L -1, we obtain a simple analysis tool based on the Fast Fourier Transform (FFT). In this case, equations (1.47) and (1.48) can be seen as a spectral analysis of the |r(k)| 2 and r 2 (k).

If r(t) is cyclostationary, |r(k)| 2 and r 2 (k) will contain periodicities buried in the noise. Hence, spectral lines related to r(t) cyclic frequencies will pop up in the spectrum obtained through equations (1.47) and (1.48). The accuracy of the cyclic frequency estimation will depend on the time series length, L. In [START_REF] Weber | Data preprocessing for decametre wavelength exoplanet detection: an example of cyclostationary rfi detector[END_REF], we derived the detection performances of this approach.

This level of analysis is sometimes not adequate to determine exactly the kind of modulation behind these cyclic frequencies. For example, m-PSK or QAM modulations with same baud rate but different bit rates cannot be differentiated. Then, more complex methods must be applied.

In [START_REF] Dandawate | Statistical tests for presence of cyclostationarity[END_REF], Dandawate and Giannakis proposed a more sophisticated approach based on detection theory where spectral lines can be stated as cyclic frequencies with a given false alarm probability. In fact, after Fourier transform, any kind of nonlinear transformations may provide some spectral lines corresponding to the modulation characteristics. In particular higher order cyclic approaches may improve modulation identification [START_REF] Marchand | Multiple hypothesis modulation classification based on cyclic cumulants of different orders[END_REF].

We applied our technique on the Westerbork synthesis radio telescope (WSRT, see A.2) data. Figure 1.6 represents the spectrum of the received signal for one antenna. The received RFI is a GPS (global positioning system) satellite signal, BPSK modulated. Thus, its theoretical cyclic and cyclic conjugated frequencies are respectively α = l/T sym , l ∈ ℤ and α = 2f 0 + l/T sym (f 0 and T sym are respectively the carrier frequency and the BPSK symbol duration). Figure 1.7 shows the results obtained with this technique. The theoretical cyclic frequencies can be retrieved on the measured spectra obtained from equations (1.47) and (1.48).

Conclusion

For single dish radio telescopes, in the cyclostationary case, the autocorrelation of the received signal can be written as:

R(t, τ ) = R s (τ ) + R RF I (t, τ ), ( 1
R s (τ ) = t/R RF I =0 R(t, τ )dt.
(1.50)

3. RFI estimation: under certain conditions, it is possible to estimate the RFI power spectrum from R(t, τ ). This estimate can then be subtracted from the power spectrum of the measured signal to obtain an estimate of the SOI power spectrum [START_REF] Bretteil | Estimation spectrale en temps réel de sources radioastronomiques en présence de brouilleurs cyclostationnaires: Application au radiotélescope décimétrique de Nançay[END_REF]. The last two techniques can be used even if the RFI and the SOI share the same timefrequency area.

The one dimensional cyclostationary tools and approaches presented in this chapter will be extended next to the multidimensional sensor case. First we will propose a multidimensional sensor model corresponding to our radio astronomical context, then we will extend the previous cyclostationary concept to the proposed model.

Chapter 2

Phased array model description

A phased array consists of an array of similar elements that are excited in or out of phase with each other to control beam steering or scanning and radiation patterns without moving the array [START_REF] Hansen | Phased array antennas[END_REF]. The array is assumed to be in free space. The purpose of this chapter is to present a multidimensional data model for phased array radio telescopes, which will be used in the manuscript. Some related tools will also be described in order to simplify the analysis of the proposed interferences mitigation techniques discussed in the following chapters.

The first section will present the assumptions made on the array and the received signals.

A general model of the telescope output will be therefore proposed. Based on this model, the correlations between the sensors output will be discussed in the second section. Then, the correlation matrix model will be extended to the cyclic approach, using the cyclostationary properties of the radio frequency interferers (RFIs). The correlation matrices contain useful information on the inpinging signals and the array response, that can be extracted using subspace decomposition. This process will be explained in the third section for both the classic and the cyclic cases. The chapter will end with some highlights resulting from the comparison between the classic and the cyclic approaches.

Data model

Single source model

We assume the sky as a collection of discrete point sources, which can be of an astronomical or a telecommunications nature. The astronomical sources are complex white Gaussian Noise with statistics CN (0, σ 2 s ), and the RFIs are defined by their modulation properties. We consider one discrete source in the sky u(t), received with a delay on the i th sensor of an array. The delay τ θu,i is measured relative to a reference sensor. The output of each sensor, illustrated in figure 2.1, is expressed by:

z i=1...M (t) = γ i,θu (t)u(t -τ θu,i ) + n i (t), (2.1)
where:

• γ i,θu (t) is the sensor beam response including the electronic chain response. γ i,θu (t) is direction dependent and time varying due to motion of the source. For cosmic sources, this response can be different from one sensor to another. Nevertheless, it can be known through calibration (which involves estimating γ i,θu (t) by pointing the telescope at a CHAPTER 2. PHASED ARRAY MODEL DESCRIPTION known single point source). For RFI sources, this response is usually unknown since the signal enters through the side lobes, which are not well characterized (single beam feed array case). For aperture arrays, this gain can be known if the RFI enters the main lobe.

• n i (t) is the independent identically distributed (i.i.d) complex Gaussian Noise with statistics CN (0, σ 2 i ). It is assumed to be uncorrelated between sensors. Certain assumptions on the received signal must be taken into consideration:

• Broad-band / narrow-band assumptions: A signal is considered as narrow-band if the signal is band-limited (by its own nature or by sub-band filtering [START_REF] Boonsta | Radio frequency interferences mitigation in radio astronomy[END_REF]) to a frequency band ∆f 0 centered at f 0 , with ∆f 0 /f 0 << 1 [START_REF] Proakis | Digital signal processing: Principles, algorithms and applications[END_REF], or ∆f 0 << (2πτ max ) -1 (where τ max is the maximum propagation delay between the array elements) [94] [44].

Let (t) be the complex envelope of the source u(t). The complex envelope τ θu,i (t) of a τ θu,i delayed signal will be:

Under the narrow-band assumption:

τ θu,i (t) ∼ = (t)e -j2πf 0 τ θu,i .
Thus, if we consider the received signal u(t -τ θu,i ) expression:

u(t -τ θu,i ) = (t -τ θu,i )e -j2πf 0 τ θu,i e j2πf 0 t , τ θu,i (t) (2.2)
under the narrow-band assumption, we obtain:

u(t -τ θu,i ) ∼ = (t)e j2πf 0 t e -j2πf 0 τ θu,i . u(t) (2.3)
Hence, if the narrow-band condition holds, equation (2.1) becomes:

z i (t) = γ i,θu (t)e -j2πf 0 τ θu,i u(t) + n i (t). a i,u (t) (2.4)
In this thesis, except when specified otherwise, we will consider block processing and will assume that the parameters of our model are constant over the block duration. In particular, a i,u (t) will be constant over the block duration. For example [START_REF] Boonsta | Radio frequency interferences mitigation in radio astronomy[END_REF], the block duration in the Westerbork synthesis radio telescope (WSRT -the Netherlands) varies from 10 ms up to 100ms, while the stationarity time of an astronomical source is 10s, due to the Earth's rotation. If a i,u (t) is not constant over the block duration, we can use a multipath model (see the multipath effects below). Thus, the M sensors array model will be:

z(t) = a u u(t) + n(t), (2.5)
where: z(t) = [z 1 (t)...z M (t)] T is the array output vector, a u = [a 1,u ...a M,u ] T is the spatial signature of the source u and n(t) = [n 1 ...n M ] T is the system noise vector with independent Gaussian complex entries.

Except in section 6.5, the narrow-band hypothesis on the received signals will be assumed herein.

• Near-field / far-field considerations: τ θu,i can be expressed in a closed form from the sensor array geometry. If we consider the far field hypothesis (i.e the emitted spherical wavefront can be considered flat at the sensor array), and we assume that the propagation delays are determined by a source in direction l u and that each sensor has a position in the array represented by p i , then cτ θu,i = p T i l u , where c is the propagation speed (figure 2.2).

Following [START_REF] Rohlfs | Tools of radio astronomy[END_REF], the far-field assumption can be expressed by the satisfied condition :

d source > 2D 2 array λ wavelength , ( 2.6)
where d source is the distance between the array and the source, D array the array aperture diameter and λ wavelength the signal wavelength.

This signal can be an interference (r(t)) or a cosmic source (s(t)). The position of each sensor is given by p i . n i is the additive noise. The array output vector z(t) is constituted by all the sensors'

outputs z i=1,...,M (t).

DATA MODEL

The far-field assumption generally holds for cosmic sources [START_REF] Rohlfs | Tools of radio astronomy[END_REF]. Hence, once the sensor array has been calibrated (i.e γ i,θu is known), a u is completely linked with the source direction and the sensor array geometry. Then, beamforming or sky mapping can be performed. This is not always the case for RFI. When the far-field hypothesis is not valid, the nearfield is assumed, and we consider that the different delays are unknown as the direction of the RFI may not be known. Thus, a u can again be considered as a random complex vector without any predictable structure. With this model, the sensor array geometry does not have to be known.

• Multi path effects: Except for sciltillation (multi paths caused by ionospheric effects), the cosmic sources multi paths can be neglected, due to the fact that the ground plane is neither flat nor perfectly conductive (which may have led to a reflexion of the signals) [START_REF] Boonsta | Radio frequency interferences mitigation in radio astronomy[END_REF]. However, the multipath effects cannot always be neglected for RFIs. We therefore obtain a multiple sources configuration (see next section).

In conclusion, under the narrow-band assumption:

z(t) = a u u(t) + n(t), (2.7)
where u(t) represents either a cosmic source (i.i.d complex Gaussian) or a RFI (modulated signal), n(t) is the system noise vector, considered as independent complex Gaussian entries and a u is the spatial signature of u, and can be considered as a random complex vector without any predictable structure because of one or more of the following reasons:

• the sensor array is uncalibrated (valid for cosmic sources or RFIs),

• the signal enters through the side lobes (RFI case only),

• the far field hypothesis is not valid (RFI case only).

Since calibration and imaging are beyond the scope of this thesis, we do not need a more precise model for a u . The structured model will only be used when simulating sky maps or beamforming. In fact, our model is quite general. Thus, having a calibrated array and/or the far field hypothesis are just particular cases of our model.

Multiple sources model

If we consider receiving K inpinging sources u k=1...K (t), under the narrow-band assumption, equation (2.4) becomes:

z i (t) = K k=1 γ i,θu k (t)e -j2πf 0 τ θu k ,i u k (t) + n i (t), = K k=1 a i,u k u k (t) + n i (t).
(2.8)

If A u = [a u 1 , ...a u K ] is a M × K matrix where each a u k = [a 1,u k ...a M,u k ] T
is the spatial signature of the corresponding k th signal, and u(t) = [u 1 (t)...u K (t)] is the K × 1 vector of the K received signals, the output of the telescope in this case will be:

CHAPTER 2. PHASED ARRAY MODEL DESCRIPTION z(t) = A u u(t) + n(t),
(2.9)

where n(t) is the M × 1 system noise vector with independent Gaussian complex entries. Among the K sources, we consider receiving K 1 RFIs that are cyclostationary, with the same cyclic frequency α and K 2 other sources, including RFIs that are not cyclostationary at α and cosmic sources which are white Gaussian, as represented by figure 2.3.1

Figure 2.3: Separation between received signal types

For simplicity's sake, we will consider next that we do not receive cyclostationary interferences with cyclic frequencies different than α. From all the above-mentioned considerations, the telescope output vector z(t) can be modeled in complex baseband form as:

z(t) = A r r(t) + A s s(t) + n(t) , RFI cosmic sources system noise (2.

10)

where A r = a r 1 , ..., a r k1 , ..., a r K 1 is a M × K 1 matrix where each a r k1 = [a 1,r k1 ...a M,r k1 ] T is the spatial signature of the corresponding k1 th RFI, and

r(t) = [r 1 (t)...r K 1 (t)] T is the K 1 × 1 vector of the K 1 cyclostationary interferers, A s = a s 1 , ..., a s k2 , ..., a s K 2 is a M × K 2 matrix where each a s k2 = [a 1,s k2 ...a M,s k2 ]
T is the spatial signature of the corresponding k2 th cosmic source, and s(t) = [s 1 (t)...s K 2 (t)] T is the K 2 × 1 vector of the K 2 white, Gaussian, independent cosmic source signals. n(t) is the system noise vector, considered as independent complex Gaussian entries.

This model is a multidimensional extension of the single source model (see the conclusion to section 2.1.1). We recall that A r and A s contain non structured spatial signatures. We also assume that:

• Each incoming cosmic sources signal has a power of σ 2 s k2 , with k2 = 1...K 2 . However, the K 1 signals are not independent of each other, with respect to the multipath issue or coupling issues.

• The system noise on each sensor has a power of σ 2 n i , with i = 1...M . • The K 1 signals are independent of the K 2 signals and the system noise (the K 2 signals are also independent of the system noise).

The following section will study the correlations between the antennas using the telescope output expressed by equation (2.10).

Correlation Matrices

In radio astronomy, Sky images are obtained through the correlation matrices measured by phased array radio telescopes. Indeed, all the spatial and statistical information of the cosmic sources are contained in these matrices.

From our previous model, in the next subsection, we will derive the expression of the correlation matrix, the cyclic correlation matrix and the cyclic conjugated matrix. In radio astronomy, the correlation matrices are calculated with τ = 0. However, in this section, we will also give the expressions of the different correlation matrices for τ = 0, since for the cyclostationary case, the maximum of information may occur for τ = 0. Moreover, it will be necessary to employ expressions with τ = 0 when studying the broad band case (section 6.5).

Classic correlation matrix

The classic correlation matrix is expressed by:

R(τ ) = z(t + τ 2 ).z H (t - τ 2 ) ∞ , ( 2.11)
with:

z(t) = A r r(t) + A s s(t) + n(t).
(2.12)

According to the independence assumption of r(t), s(t) and n(t), we obtain that:

r k s * l ∞ = 0 r k r * l ∞ = 0 s k n * l ∞ = 0 however s k s * l ∞ = 0 r k n * l ∞ = 0 n k n * l ∞ = 0 (2.13)
(Note that we do not assume for now the independence of cosmic sources and system noise in order to have a general correlation matrix model, i.e for cases where the system noise is correlated XXX) With this assumption, the classic correlation matrix expression will be:

R(τ ) = A r R r (τ )A H r + A s R s (τ )A H s + R n (τ ). = R RF I (τ ) + R Cosmic (τ ) + R N oise (τ ) (2.14)
defining :

• R r (τ ) ≜ r(t + τ 2 ).r H (t -τ 2 ) ∞ , • R s (τ ) ≜ s(t + τ 2 ).s H (t -τ 2 ) ∞ , • R n (τ ) ≜ n(t + τ 2 ).n H (t -τ 2 )
∞ .

An illustration of the skymaps using R(0) is given by figure 2.4. Maps have been obtained with a whole sky scanning.

Cyclic correlation matrix

The cyclic correlation matrix is expressed by:

R α (τ ) = z(t + τ 2 ).z H (t - τ 2 )e -j2παt ∞ . ( 2.15)
All the contributions from the non-cyclostationary signals, and the cyclostationary signals that have a cyclic frequency α ′ = α tend to zero. From these considerations, equation (2.15) reduces to:

R α (τ ) = A r R α r (τ )A H r = R α RF I (τ ), (2.16)
where

R α r (τ ) ≜ r(t + τ 2 ).r H (t -τ 2 )e -j2παt ∞ .

Cyclic conjugated correlation matrix

In some cases (see section 1.1), when the cyclic correlation matrix is zero for all α = 0 (example: AM modulated signals), we will have to use the conjugated correlation matrix, defined by:

R α (τ ) = z(t + τ 2 ).z T (t - τ 2 )e -j2παt ∞ .
(2.17)

Similarly to the cyclic correlation function, because of the non cyclostationarity of cosmic sources and system noise, (2.15) reduces to:

R α (τ ) = A r R α r (τ )A T r = R α RF I (τ ), (2.18) with R α r (τ ) ≜ r(t + τ 2 ).r T (t -τ 2 )e -j2παt
∞ .

CORRELATION MATRICES

4: Skymaps with real and synthetic data. The upper figures represent real observations using the LOFAR initial test stations (ITS) (see appendix A.1, by courtesy of Albert-Jan Boonstra and ASTRON). The lower figures represent the corresponding simulations that we have obtained using the proposed model. The skymaps on the left represent an observation with a strong received RFI (cosmic sources are not visible in the map). Figures on the right represent the same observation without RFI. 3 sources appear in the map. For the lower figures, the correlation matrices were obtained firstly with

Independent signals case

If we consider now that all the signals are independent (independence of all the sources), R r , R s , R n , R α r (τ ) and R α r (τ ) reduce to diagonal matrices. Each element of the diagonal is, depending on the specific case, either the autocorrelation, cyclic autocorrelation or cyclic conjugated autocorrelation function of each source (see definitions in section 1.1). Expressions (2.14), (2.16) and (2.18) become respectively:

R(τ ) = A r diag R r 1 ,r * 1 (τ ), ..., R r K 1 ,r * K 1 (τ ) A H r +A s diag R s 1 ,s * 1 (τ ), ..., R s K 2 ,s * K 2 (τ ) A H s +diag (σ 2 1 , ..., σ 2 M ) R α (τ ) = A r diag R α r 1 ,r * 1 (τ ), ..., R α r K 1 ,r * K 1 (τ ) A H r R α (τ ) = A r diag R α r 1 ,r 1 (τ ), ..., R α r K 1 ,r K 1 (τ ) A H r (2.19)
where diag(.) is the diagonal operator, that converts a vector into a diagonal matrix, where the off-diagonal elements are set to zero.

Conclusion

The classic and cyclic correlation matrices described in this section are structured according to the characteristics of the RFIs, the astronomical signals and the system noise. These different structures lead to a possible separation between the signals. Indeed, from the correlation matrices, we can estimate the spatial signatures of the incoming signals. These spatial signatures will be used later in all the RFI mitigation algorithms that will be presented. Their estimation is based on a subspace decomposition of the corresponding correlation matrices. This is the topic of the next section.

Subspace decomposition

In order to remove the RFIs from the received data, we will have to estimate the spatial signature of the undesired signals. Each received signal is identified by its spatial signature in the received data model (section 2.1). When spatial signatures of the interfering sources are unavailable, we can estimate them using a subspace decomposition of the corresponding correlation matrix. This process will be detailed in this section, for R(τ ), R α (τ ) and R α (τ )

(respectively the classic, cyclic and cyclic conjugated correlation matrices).

In the following chapters, we will use the notation "subspace decomposition" to refer to an eigenvalue decomposition (EVD) for full-rank matrices and a singular value decomposition (SVD) for cyclic correlation matrices. In this section, for explanatory purposes, we will distinguish the two sorts of decomposition (EVD and SVD), but in the rest of the document, we will use the general term of eigenvalues for both SVD and EVD (subspaces decomposition) to avoid too many notations for the reader. Further explanation can be found in [START_REF] Friedberg | Linear Algebra[END_REF] and [START_REF] Krim | Two decades of array signal processing research[END_REF] for the classic EVD, and in [START_REF] Gardner | Cyclostationarity in communications and signal processing[END_REF] for the cyclic SVD.

SUBSPACE DECOMPOSITION

Classic eigenvalues decomposition

Consider the correlation matrix expressed in equation (2.14). As R N oise is assumed to be always a full rank matrix, R is full rank as well, and can be decomposed into eigenvalues. If U contains the eigenvectors u i : U = [u 1 , ..., u M ], and Λ contains the corresponding eigenvalues

λ i : Λ = diag(λ 1 , ..., λ M ) (sorted in decreasing order λ 1 ≥ λ 2 ≥ ...λ M ≥ 0): R(τ ) = A r R r (τ )A H r + A s R s (τ )A H s + R n (τ ) = UΛU H , (2.20)
where U is a M ×M unitary matrix(UU H = U H U = I). A physical interpretation of eigenvalue decomposition can be that the eigenvectors give an orthogonal set of directions present in the correlation matrix sorted in decreasing order of dominance. The eigenvalues give the power of a beamformer matched in that direction. Without further assumptions, there is no link between these directions and the spatial signatures contained in A r and A s . Now, we consider that the noise is i.d.d (independent and identically distributed) and spatially white 2 with the same power on each antenna equal to σ 2 N oise . If we assume a large number of samples:

λ 1 ≥ λ 2 ≥ ...λ K ≥ λ K+1 ≃ λ K+2 ≃ ... ≃ λ M ≃ σ 2 N oise . (2.21)
where K is the number of inpinging signals

(K = K 1 + K 2 consisting of K 1 RFIs and K 2 cosmic sources).
From the EVD, it is possible to find an orthonormal basis which forms two subspaces [START_REF] Friedberg | Linear Algebra[END_REF] [48]:

• Signal subspace generated by the K first eigenvectors

U s = [u 1 , ..., u K ],
• Noise subspace generated by the M -K smallest eigenvalues:

U N = [u K+1 , ..., u M ].
Noise subspace eigenvectors are orthogonal to the K vectors of the signal subspace.

However, the RFI spatial signatures matrix A r cannot be estimated from U s . The only way for A r and U s to span the same subspace is to neglect the cosmic contributions in equation (2.20) (i.e: R(τ

) = A r R r (τ )A H r + σ 2 N oise δ(τ )I)
, where δ(.) is the Dirac delta function on τ .

Cyclic and cyclic conjugated singular values decomposition

As R α and R α are non hermitian, we perform a singular values decomposition (SVD) instead of an EVD. Moreover, R α and R α have a rank less than or equal to K 1 [START_REF] Xu | Direction of arrival estmation via exploitation of cyclostationarity -a combination of temporal and spatial processing[END_REF].

If Λ Cyclic is the singular values matrix, U Cyclic the input basis vectors matrix and V Cyclic the output basis vectors matrix, the correlation matrix (eq.(2.16)) will be decomposed into:

R α (τ ) = A r R α r (τ )A H r = R α RF I (τ ) = U Cyclic Λ Cyclic V H Cyclic (2.22)
2 If the noise matrix contains unequal elements on the diagonal (i.e not identical for all sensors), a prewhitening can be done by calibration procedure. We will call N the estimated noise matrix [START_REF] Boonsta | Radio frequency interferences mitigation in radio astronomy[END_REF]. The corrected model R w to which the EVD will be applied is :

R w = N -1 2 R N -1 2 .
Here, the estimated singular values will be:

λ 1 ≥ λ 2 ≥ ...λ K 1 ≥ λ K 1 +1 = λ K 1 +2 = ... = λ M = 0. (2.23)
The K 1 dominant singular values are related only to the RFIs directions. Thus, span(U Cyclic ) = span(A r ). No assumptions have been made on the calibration of the noise, as it does not affect the estimation of the K 1 strongest singular values.

However, since the spatial signatures of the interferers are not orthogonal to each other, the columns of the signal subspace formed by the K 1 input basis vectors

U cyclic = [u 1 , ..., u K 1 ]
do not correspond exactly to the columns of A r . In that case, u 1 will point to a direction common to all the interferers, u 2 will point to the remaining direction orthogonal to u 1 , u 3 will point to the remaining direction orthogonal to u 1 and u 2 etc. However, in the case of a single received interferer K 1 = 1 (which is a common case, see footnote in section 2.1.2), A r is equal to a r 1 . Thus, we obtain for the cyclic case, using equation (2.19), that:

R α (τ ) = R α r 1 ,r * 1 (τ )a r 1 a H r 1 = λ 1 u 1 u H 1 (2.24)
It is straightforward to see that:

λ 1 = R α r 1 ,r * 1 (τ ) a r 1 2 u 1 = a r 1 a r 1 -1 (2.25)
Similar results are obtained for Rα .

Conclusion

With negligible cosmic sources, and with a calibrated system (i.e: R(τ ) ≈ R RF I (τ ) + σ 2 N oise δ(τ )I), the classic and cyclic approaches are equivalent, since the signal subspace derived from the subspace decomposition spans the same subspace as the RFI spatial signatures. However, as we make the model more complex with uncalibrated noise, or add cosmic sources, only the cyclic approach can estimate the RFI subspace. Consequently, the cyclic approach will be more robust than the classic approach. Figure 2.5 summarizes the previous results.

Conclusion

In this chapter, a global model with related tools has been presented, for narrow-band received signals. The study of the correlations between the sensors' outputs has shown that, the cyclic correlation matrices asymptotically contain only information on the RFIs. Indeed, by using the subspace decomposition of the cyclic correlation matrices, we obtained an interesting result: the estimated signal subspace spans the same subspace as that defined by the RFI spatial signature vectors, regardless of the presence of cosmic sources and/or uncalibrated system noise.

For our general model, the signal subspace in the classic case will depend on the interferers and the cosmic sources. In the cyclic case, it will depend only on the RFIs.

Thus, more robust methods can be derived using the cyclostationarity principle, compared to the classic approaches. From the cyclic correlation matrices (R α and Rα ) subspace decomposition, eigenvalues and subspaces related to the RFI will only be estimated and used to design three different RFI mitigation approaches. Chapter 4 will describe the cyclo-detection approach for phased array radio telescopes. In some cases, the particular structure of the cyclic correlation matrices (see equations (2.24) and (2.25)) allows an estimation of the RFI correlation matrix (R RF I ), in order to remove this RFI contribution from the global correlation matrix (R). This method will be explained in chapter 5. Chapter 6 will describe cyclic spatial filtering. It is a general method that can be extended to broad-band received signals, also discussed in chapter 6.

Since the proposed techniques are based on existing array processing methods, an overview of main phased array techniques will be presented in the next chapter. Figure 2.6 summarizes the cyclic approaches, compared to the classic approaches.

Figure 2.6: Classic (using R) and cyclic (using R α ) approaches for RFI mitigation in phased arrays. In the presence of one cosmic source and one RFI source for instance, two dominant eigenvalues result from the subspace decomposition of R. From the cyclic subspace decomposition results only one dominant eigenvalue. The presence of the cosmic source in the estimation of the spatial signature may affect the classic RFI mitigation techniques (detection and spatial filtering). In addition to the classic approaches, that may be improved using cyclostationarity, the estimation and subtraction technique may be proposed to process a certain class of interferers.

Chapter 3

Overview of multidimensional signal processing

This chapter presents an overview of existing array processing techniques. The aim of this chapter is to extract useful approaches that can be applied to our RFI mitigation issue.

We will first start with a description of existing multidimensional RFI mitigation techniques in radio astronomy, pointing out the limitations of each method and the interest of employing the cyclostationarity concept.

The second part of the section will present some well-known algorithms that are employed in telecommunications applications such as source localization, signal extraction and nulling. Depending on the "philosophy" behind each algorithm, these techniques can be divided into two kinds of methods :

• Beamforming methods which are based on a weighted filtering of the signals, according to an a-priori knowledge of the array characteristics and the properties of the signals of interest.

• High-resolution approaches, based on signal and noise subspaces determination.

For telecommunications algorithms, we will modify the general data model given in section 2.1 as follows, by neglecting the cosmic sources part, in a M sensors network receiving K sources:

z(t) = Ar(t) + n(t), (3.1)
where z(t) is the array output vector, A = [a 1 , ..., a K ] is a M ×K matrix where each a k=1,...,K = [a 1 k ...a M k ] T is the steering vector of the corresponding source, and r(t) = [r 1 (t)...r K (t)] and n(t) = [n 1 (t)...n M (t)] are respectively the sources signals and the spatially white Gaussian system noise. Note that the term "steering vectors" is used instead of "spatial signature", as the algorithms presented are based on the assumption of a structured vector. For example, an equi-spaced linear array will have the following steering vector:

a k=1...K = 1, e -j(2π f c ∆)sin(θ k ) ...e -j(2π f c ∆)(M -1)sin(θ k ) T , (3.2)
where f is the center frequency of the signal, c the propagation speed and ∆ the distance between the antennas (these parameters are assumed to be known). Note that if K = 1, A = a.

Existing methods in radio astronomy

Since the introduction of the phased array concept in radio astronomy, several multidimensional RFI mitigation approaches had been proposed [START_REF] Boonstra | Rfi mitigation methods strategy[END_REF]. These methods can be regrouped into RFI detection methods, subspace approaches and adaptive beamforming. The latter will be studied in 3.3.1.

In the following approaches, two important assumptions are made:

• the astronomical source signal power must be much less than the interferer powers.

• the system noise must be whitened.

Thus, the received data model will be:

z(t) = A r r(t) + n(t) , RFI system noise (3.3)
where

A r = a r 1 , ...a r K 1 is a M × K 1 matrix of the K 1 spatial signature of each corresponding k1 RFI signal, and r(t) = [r 1 (t)...r K 1 (t)] is the K 1 × 1 vector of the K 1 interferers. n(t)
is the M × 1 system noise vector with independent identically distributed Gaussian complex entries.

Classic detection

Since the detection is based on a binary problem, the observations can be expressed by two testing hypotheses H 0 and H 1 :

• H 0 : there is no interference,

• H 1 : there is at least one interferer detected.

Let T be the detection criterion, γ the detection threshold and d the detection decision. The detection hypotheses yield:

T ≤ γ → d = 0 → H 0 , T ≥ γ → d = 1 → H 1 .
(3.4)

In the case of an array radio telescope, the basic idea beyond the detection criterion is to maximize the power coming from the RFI by focusing the array in that direction. Then, it is compared with the noise power, σ 2

n . An important consideration is that the system noise is assumed calibrated (i.e. R noise = σ 2 n I). Assuming that the spatial vector is known, we can beamform the array toward the RFI. We define y(t) to be the output of the matched beam former in the direction of z(t), y(t) = |a r | -1 a H r z(t). The corresponding power can be expressed by

P y = |a r | -1 a H r Ra H r .

EXISTING METHODS IN RADIO ASTRONOMY

The Neyman-Pearson detector based on the data matrix Z = [z(1), ...z(L)] (with L the length of the time series) considers the estimated data correlation matrix and is given by: [47]

T (Z) = Lσ 2 N oise a H r Ra r a H r a r > H 1 < H 0 γ. (3.5)
where γ is a power threshold, which is related to a given false alarm probability. The Neyman-Pearson detector is generally used in radio astronomy multichannel detection because of its independence of the telescope characteristics and its independence of the required a-priori knowledge of the RFI presence probabilities or the cost of correlation.

Depending on the process applied on the correlation matrix R = z(t)z H (t)

L
, where L is the length of the time series, several multi-channel detectors have been proposed. Among these methods, the following can be mentioned [START_REF] Boonstra | Rfi mitigation methods strategy[END_REF]:

The eigenvalue test [7] [54] [55] [51] [86]

If a telescope array is calibrated, and if the location or direction of an interferer is known, then the spatial signature of the interferer is also known. In that case, T (Z) (eq. (3.5)) will compare the received power in the direction a r to the noise power σ 2 N oise , which is known as the matched spatial filter detector.

If the interferer spatial signature is unavailable, a r can be estimated by the eigenvector u 1 corresponding to the largest eigenvalue of R.

Thus, as a r R a H r = λ1 for a r = u 1 , the test statistic in (3.5) takes the form of a matched spatial filter:

T (Z) = Lσ 2 N oise λ1 > H 1 < H 0 γ (3.6)
A threshold γ can be set by simulations or using the maximum eigenvalue estimate as defined in [START_REF] Edelman | Eigenvalues and condition numbers of random matrices[END_REF]

γ = σ 2   1 + M L   2 (3.7)
The norm-based detector [START_REF] Leshem | Multichannel detection of gaussian signals with uncalibrated receivers[END_REF] The selected norm for interferences detection is the Frobenius norm. This detector requires a pre-whitening of the correlation matrix. The test statistic applied on the pre-whitened correlation matrix R w is:

T (Z) = || R w || 2 F > H 1 < H 0 γ (3.8) Note that || R w || 2 F = λ 2 i
, where λ i are the eigenvalues of R w .

The generalized likelihood ratio test detector [54] [52] [53]

The eigenvalue analysis showed that in the absence of interferences all eigenvalues are asymptotically equal to σ 2 N oise . If the noise power σ 2 N oise is known, we can apply the (generalized) likelihood ratio test (GLRT), which leads to a test statistic given by:

T (Z) = -LM log M i=1 λi σ 2 N oise = -LM log | R| σ 2 N oise (3.9)
where λi is the i th eigenvalue of R. An interference is detected if T (Z) > γ.

If the noise power is unknown as well, the Minimum Description Length (MDL) detector [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF] can be used, even if, at present, some improvements need to be made to the classic approach, such as determination of the false alarm rate.

Spatial filtering

To overcome the continuous interferences issue (i.e observations where the RFIs cannot be masked on the time-frequency plane), filtering methods are applied. One filtering method proposed for radio astronomy is spatial filtering [START_REF] Leshem | Multichannel interference mitigation techniques in radio astronomy[END_REF] [7] [START_REF] Raza | Spatial filtering of RF interference in radio astronomy[END_REF] The spatial filtering is a subspace approach based on nulling the power received from undesired directions (RFIs incoming directions). This nulling is achieved by projecting out the undesired direction, which is such that PA r = 0. The projector P is defined by:

P ≜= I -A r (A H r A r ) -1 A H r (3.10)
where A r is the M × K 1 matrix containing the K 1 spatial signatures of the RFIs. When available, these spatial signatures are removed from the observations by applying the projector on the correlation matrix R:

R cleaned = PRP (3.11)
When the spatial signatures of the interferers are unknown, they can be estimated by a subspace analysis of the correlation matrix, which can be written in terms of eigenvalues and eigenvectors as

R = UΛU H (3.12)
where U is a M × M matrix containing the eigenvectors, and Λ is a M diagonal matrix containing the eigenvalues. Assuming that the astronomical contribution is small, the eigenvalue decomposition yields M -K 1 eigenvalues that are equal to the system noise power σ 2 N oise (see section 2.3.1). Thus, span(A r ) = span(U K 1 ), where U K 1 is formed by K 1 dominant eigenvectors. As U K 1 and the noise subspace U N (formed by the M -K 1 eigenvectors corresponding to the noise eigenvalues) are orthogonal, the projector expression can be written as:

P = U N U H N , (3.13)
and thus, the projector can be built using only the noise subspace, without requiring spatial signature estimation.

Conclusion

Classic detection and spatial filtering techniques have been developed by assuming that cosmic sources are neglected and that the telescope is calibrated. If one of these conditions is not satisfied, the false alarm rate will increase (for detection), and the interferer subspace estimation will be distorted (for spatial filtering).

To overcome these issues and obtain better performances (a more accurate detection of the interferers and a robust estimation of the RFI subspace), we can use the cyclostationary properties of the interferer. Cyclostationarity is a widely used concept in telecommunications approaches to improve existing techniques. The next sections will present an overview of the 3.3. BEAMFORMING METHODS general array processing algorithms, in order to explore some resulting principles for multichannel RFI mitigation.

Beamforming methods

Multidimensional processing techniques were first proposed for application in radio communications in [START_REF] Frost | An algorithm for linearly constrained adaptive array processing[END_REF], [START_REF] Winters | Optimum combining in digital mobile radio withco-channel interferernce[END_REF] and [START_REF] Swales | The performance enhancement of multibeam adaptive base station antennas for cellular land mobile radio systems[END_REF]. The aim of these algorithms is the extraction of the signal of interest, and the determination of the direction of arrival of the undesired signal (for radio astronomy purposes). Beamforming algorithms can be divided into two categories: adaptive beamforming techniques, which employ minimization of the output power subject to linear constraints, in order to protect the desired signal from being canceled with the minimization process, and property restoral techniques (SCORE, CAB), which adapt the array to restore the properties of known signals of interest to the array output.

Classical and adaptive beamforming

The generalized method consists in a weighted summation of each channel. The weights represent time-lags in the time domain and a phase difference in the frequency domain. The output of each sensor z(t) is multiplied with the conjugate of the weights:

y(t) = w H z(t). (3.14)
In order to estimate the weights vector w, several cost functions were defined, allowing the estimation of the direction of arrival of the sources. The most commonly employed cost function is the power function of the output of the beamformer :

P = |y(t)| 2 ∞ = w H Rw, (3.15)
where R is the correlation matrix of the array output.

Fixing w = a(θ) (which is the basic method), the DOA of the sources is given by the value of θ which maximizes the output power :

P (θ) = a H (θ)Ra(θ).
(3.16)

The adaptive beamforming process [START_REF] Capon | High resolution frequency-wavenumber spectrum analysis[END_REF] (also known as the Capon method or the minimum variance distortionless response (MVDR) method) uses the constraint w H a(θ) = 1, which leads to a minimization of the received power in all directions except the direction of observation. The filter is the weights vector min w w H Rw with the constraint w H a(θ) = 1. The solution is given by :

w M V DR = R -1 a(θ) a H (θ)R -1 .a(θ) (3.17)
The Capon filtering power will be:

P (θ) = 1 a H (θ)R -1 a(θ) . (3.18)
Several other adaptive beamforming methods have been derived from the capon principle, among which, [START_REF] Van Trees | Optimum array processing: Part IV of detection, estimation, and modulation theory[END_REF] [80] and [START_REF] Li | Robust adaptive beamforming[END_REF], where adaptive beamforming approaches robust to small sample size problems and steering vector estimation errors have been presented. Adaptive beamforming techniques have also been employed in radio astronomy [START_REF] Fridman | DSP experimental system for radio frequency interference mitigation at radio telescope[END_REF] [25] [START_REF] Ellingson | Beamforming and interference cancelling with very large wideband arrays[END_REF] [3].

Cyclic adaptive beamforming technique

The adaptive beamforming algorithm proposed in [START_REF] Castedo | An adaptive beamforming technique based on cyclostationary signal properties[END_REF] [START_REF] Castedo | Linear constrained adaptive beamforming using cyclostationaritysignal properties[END_REF]is based on the property of cyclostationarity, using a particular class of nonlinear transformations that generates spectral lines. A zero-mean complex signal r(t) generates a spectral line at a frequency α after passing through the nonlinearity (.) p if it has a non zero p th order cyclic moment defined by:

m α pr = r p (t)e -j2παt ∞ = 1 T lim T →∞ T /2 -T /2 r p (t)e -j2παt dt, (3.19)
which is a generalization of the second order definitions given in section 1.1. The cyclic moment m α pr represents the spectral line at α contained in the PSD of z p (t). As shown in [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF] and section (1.3), the cyclic moments correspond, under some assumptions, to the coefficients of the Fourier series expansion of the statistical moments.

The algorithm proposed in [START_REF] Castedo | Linear constrained adaptive beamforming using cyclostationaritysignal properties[END_REF] exploits these spectral lines by choosing the coefficients w of a narrow band antenna array, where y(t) = w H z(t) is the beamforming output, which minimizes:

J = |e j2παt -y p (t)| 2 .
(3.20)

The proposed cost function is therefore the mean square error between a complex exponential and the beamforming output after the (.) p nonlinearity, where the p and α are selected according to the order and the frequency of the spectral lines generated by the signal to be extracted (i.e the existence of non-zero p th order cyclic moment at α).

In order to minimize the cost function J, the computed weights will be:

w(t + 1) = w(t) + µǫ * (t)y p-1 (t)z(t), (3.21)
where µ is the algorithm step size and ǫ(t) = e j2παt -y p (t) is the error signal whose variance is minimized via the gradient (eq. (3.20).

A blind adaptive beamforming approach, that does not require knowledge of the cyclic frequencies, based on the gradient-based technique, was presented in [START_REF] Hung | A new adaptive beamformer for underwater communications[END_REF]).

Self-coherence restoral technique

Cyclic adaptive beamforming was first presented in [START_REF] Agee | Spectral self-coherence restoral: A new approach to blind adaptive signal extraction using antenna arrays[END_REF], where three self-coherence restoral (SCORE) algorithms were introduced: the Least Square SCORE, the Cross SCORE and the Auto-Score. The spectral self-coherence property of a signal is defined as the existing (non-zero) correlation with a frequency-shifted and probably correlated version of the signal for certain discrete values of frequency shift. The spectral self coherence functions were developed in [START_REF] Gardner | Statistical spectral analysis: A nonprobabilistic theory[END_REF] and [START_REF] Gardner | Introduction to random processes with applications to signals and systems[END_REF], where it was shown that complex wide-sense cyclostationary and wide-sense almostcyclostationary waveforms exhibit spectral self-coherence or conjugate spectral self-coherence at discrete multiples of the time periodicities of the waveform statistics. Indeed, a scalar waveform r(t) is said to be spectrally self-coherent (eq. (3.22)) and respectively spectrally conjugated self-coherent (eq. (3.23)) at frequency α if the correlation between r(t) and r(t) (respectively r * (t)) frequency shifted by α is non zero for certain lag τ :

ρ α rr * (τ ) ≜ r(t+ τ 2 ).r * (t-τ 2 )e -j2παt ∞ |r(t+ τ 2 )| 2 ∞ |r * (t-τ 2 )e j2παt | 2 ∞ ≜ R α rr * (τ ) Rrr(0) = 0, (3.22
)

ρ α rr (τ ) ≜ r(t+ τ 2 ).r(t-τ 2 )e -j2παt ∞ |r(t+ τ 2 )| 2 ∞ |r(t-τ 2 )e j2παt | 2 ∞ ≜ R α rr (τ )
Rrr(0) = 0, (3.23)
Note that, for some authors and for our document, the definition of the conjugated correlation is the inverse of the definition given in [START_REF] Agee | Spectral self-coherence restoral: A new approach to blind adaptive signal extraction using antenna arrays[END_REF].

The spectral self coherence restoral concept consists in finding the weights of the beam former output that suppress all undesired signals (including the noise), by minimizing a least square function, using the self-coherence property. The resulting weights are:

w SCORE = g SCORE R -1 aR RF I ,
where R r is the signals of interest correlation matrix (note that knowledge of this matrix is required), g SCORE = [a ( * ) ] H Cρ α ss ( * ) (τ )e -j2πτ is a gain constant, and C is the control vector. This algorithm is presented as the least square (LS) SCORE method. Two other variants of the SCORE algorithm, the cross SCORE and the auto SCORE, were developed in [START_REF] Agee | Spectral self-coherence restoral: A new approach to blind adaptive signal extraction using antenna arrays[END_REF]. Some improvements have been made to the initial SCORE algorithms: To avoid mismatch in the cyclic frequency resulting from the carrier frequency offset or the Doppler shift, some algorithms were proposed in [START_REF] Lee | Robust adaptive array beamformig for cyclostationary signals under cycle frequency error[END_REF] and [START_REF] Lee | Efficient robust adaptive array beamformig for cyclostationary signals[END_REF] to estimate the true cyclic frequency for the SCORE methods. In [START_REF] Du | Simple and pratical cyclostationary beamforming algorithms[END_REF], a gradient based algorithm was proposed by maximizing the self-coherent coefficients constrained with a normalized beamforming weight, thus increasing the convergence speed.

Cyclic adaptive beamforming algorithm

The main attractive quality of the CAB algorithm [START_REF] Wu | A fast algorithm for adaptive beamforming of cyclic signals[END_REF], [START_REF] Wu | Blind adaptive beamforming for cyclostationary signals[END_REF] compared to the SCORE method, is its faster implementation. When the array receives only one desired signal, the array output data vector z(t) and its frequency shifted version z(t)e j2παt contain the signal components r(t) and r(t)e j2παt that have a high correlation value at the computation cyclic frequency α. Thus, the weight vector of the beam former w CAB is obtained from the singular vector corresponding to the largest singular value cyclic correlation matrix R α .

Variants of the CAB algorithm, the Constrained CAB and the Robust CAB methods were also proposed in [START_REF] Wu | A fast algorithm for adaptive beamforming of cyclic signals[END_REF] to, respectively, suppress interference and increase robustness.

Conclusion

Adaptive beamforming methods suffer from:

• The cost of the matrix inversion,

• the risk of several minimals (or maximals, depending the defined cost function),

• sensitivity to pointing errors and uncertainties in aperture characteristics

• the required knowledge of the properties of the signal of interest.

Moreover, cyclic adaptive beamforming methods seems to be quite complex for our aim to find a general RFI mitigation approach. In addition, the determination of the step size is a an additional problem when using this gradient-based algorithm.

The SCORE methods are based on adapting the array in order to maximize the cyclic fluctuations of the correlation function. The major advantage of the SCORE algorithms over conventional methods is that the only necessary parameters are the cyclic frequencies of the desired signals. Nevertheless, the methods suffer from slow convergence and require the adaptation of two beam formers simultaneously. In addition, the methods fail completely when operating under environments which contain more than one spectrally self-coherent signal, at the computation frequency. Despite these drawbacks, the spectral lines concept is attractive, and a detector based on this principle, was proposed in [START_REF] Bretteil | Estimation spectrale en temps réel de sources radioastronomiques en présence de brouilleurs cyclostationnaires: Application au radiotélescope décimétrique de Nançay[END_REF], for a single dish radio telescope.

The CAB algorithm is easy to implement. It only requires the extraction of the main eigenvector obtained from the subspace decomposition of the cyclic correlation matrix. This eigenvector will be used as the beamformer weights vector. However, while these algorithms are efficient for one signal, in the presence of more than one signal with the same cyclic frequency, certain conditions have to be verified. Besides, the accurate knowledge of the array, and the calibration assumption do not fit our purpose of proposing a general method regardless to the array characteristics. Since this observation is valid for all the beamforming techniques presented, we will not use these principles in the present thesis.

High resolution methods

This section is an overview of the main high resolution methods. These methods are based on a subspace decomposition of the correlation matrix R in order to define the signal subspace vectors and the noise subspace vectors. The estimation of these subspaces is therefore used to minimize or maximize a cost function. We distinguish between two categories of high resolution approaches: linear methods, and subspace decomposition approaches.

Linear methods

Linear methods are based on a determination of the signal subspace and the noise subspace without decomposition of the correlation matrix into values. The assumption common to all these methods is that the steering vectors matrix A must be full rank, and contain at least K linearly independent lines that form a full K rank matrix. This submatrix is noted A linear and the corresponding output data vector z linear to the K used lines of A linear . A compl is the complementary submatrix containing the lines of A that have not been involved in the construction of A linear . z compl is the corresponding output data vector.

The determination of the signal subspace and the noise subspace, according to each method, makes it possible to build an orthogonal projector . This projector is subsequently used to compute a cost function (3.24). Minimizing this cost function gives the direction of arrivals (DOA) of the desired signals. (3.24)

Fuller details about the linear methods can be found in [START_REF] Marcos | Méthodes linéaires haute résolution pour l'estimation de directions d'arrivée de sources. performances asymptotiques et complexité[END_REF].

Bearing estimation without eigen decomposition method

The BEWE technique was first presented in [START_REF] Yeh | Simple computation of projection matrix for bearing estimations[END_REF] and was then extended for low signal to noise ratios [START_REF] Yeh | Projection approach for bearing estimations[END_REF] with a statistical study in [START_REF] Stoica | Statistical analysis of a subspace method for bearing estimation without eigendecomposition[END_REF].

In the presence of spatially white Gaussian noise, a (M -K) × K matrix is defined by :

U s = z compl (t)z H linear (t) ∞ = A compl rA H linear , ( 3.25)
Assuming that A compl is a full rank matrix, which implies that M > 2K, U s and A compl span the same subspace. The orthogonal projector of the BEWE method is therefore defined by:

BEW E = I -U s U H s . (3.26)

The propagator

The propagator method was published in [START_REF] Munier | Spatial analysis using new properties of the crossspectral matrix[END_REF] and then developed in [START_REF] Marcos | On a high resolution array processing method nonbased on the eigenanalysis approach[END_REF], [START_REF] Marcos | Source bearing estimation and sensor positioning with the propagator method[END_REF] and [START_REF] Marcos | The propagator method for source bearing estimation[END_REF].

Using the definition given for all linear methods, the propagator is the unique K × (M -K) operator defined by:

P H A linear = A compl or A H Q = 0 with Q ≜ (= def ) [P -I M -K ] T .
(3.27)

From these relations, it appears that Q span an orthogonal subspace to A. Therefore, if Q s = I K P H T then: span {Q s } = span {A}. The propagator P consequently can define the noise subspace represented by Q and the signal subspace represented by Q s .

The determination of Q and Q s to build the projector is given by estimating the propagator P.

In order to estimate P, we define two submatrices of the correlation matrix, G = z(t)z linear ∞ and H = z(t)z compl ∞ , such that R = [G, H]. The projector is therefore defined by (3.28)

P = (G H G) -1 G H H. (3.28)
The orthogonal projector of the Propagator method is defined by 3.29

propagator = I -Q s Q H s . (3.29)

Subspace method without eigen decomposition

The SWEDE method was presented in [START_REF] Eriksson | On-line subspace algorithms for tracking moving sources[END_REF]. The steering vectors matrix is reorganized such that A = [A 1 A 2 A 3 ] T , where A 1 and A 2 are K × K dimension matrices and A 3 M -2K × K matrix. Strong assumptions were made : M > 3K, A is a full rank matrix and A 1 and A 2 are non singular matrices.

Two unique K × (M -2K) operators, P 1 and P 2 were defined:

P 1 = A -H 1 A H 3 and P 2 = -A -H 2 A H 3 . (3.30)
Consequently, the matrix W, defined by 3.31, spans an orthogonal subspace to A, and therefore can define the noise subspace.

W =    P 1 0 0 P 2 I I    . (3.31)
The noise subspace estimation is given by the matrix V SW EDE = WD, where D is determined such that span {V SW EDE } ⊥ span {A}, which corresponds to the noise subspace. The signal subspace V s-SW EDE will be therefore an orthogonal subspace to V SW EDE . The choice of D gives a set of subspace methods, corresponding to different versions of the SWEDE technique, that do not require an eigenvalue decomposition of the correlation matrix to determine the signal and noise subspaces.

The SWEDE orthogonal projector of the method is given by (3.32):

SW EDE = I -V s V H s . (3.32)
Some other versions of the SWEDE algorithm are also presented in [START_REF] Eriksson | On-line subspace algorithms for tracking moving sources[END_REF], using the sensors output vector z(t) or the correlation matrix R.

Conclusion

To avoid subspace decomposition, linear methods make certain strong assumptions, such as A compl and A linear must be full rank. The performance of the method resides in the choice of z compl and z linear . Indeed if the elements of z linear correspond to nearby sensors, A linear can be singular. If the sources are close, the elements of z linear must be chosen from spaced sensors.

However, the proposed projector seems to be interesting for our radio astronomy context, as we project out the signal subspace (instead of the noise subspace). This principle is similar to the spatial filter (section 3.1.2), for which we will propose some improvements using cyclostationarity in chapter 6.

Subspace decomposition methods

Subspace decomposition methods yield an estimation of the desired signals DOA, with a theoretically infinite resolution, independently of the signal to noise ratio. This attractive property has given rise to a number of publications. A non-exhaustive list of the main subspace decomposition approaches is presented below.

Multiple signal classification

The MUSIC method was proposed by [START_REF] Bienvenu | Optimality of high resolution array processing using eigenstructure approach[END_REF] and [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] (with a more detailed study in [START_REF] Schmidt | A signal subspace approach to multiple emitter location and spectral estimation[END_REF]). The algorithm uses the orthogonality property between the sources steering vectors and the noise subspace eigenvectors.

Using the structured formula of the steering vector a(θ) and the eigenvectors resulting from the subspace decomposition of the correlation matrix(R) (see 2.3.1):

θ k = θ ⇒ a H (θ)u i=1,M = 0, (3.33)
where θ k is the direction of arrival of the desired (or undesired) source. Consequently, the power norm expressed by eq.(3.34) has to be minimized, according to a variation of θ in order to find the direction of the source.

p(θ) = |U H N a(θ)| 2 = a H (θ)U N U H N a(θ), (3.34)
where U N is the matrix containing the noise subspace eigenvectors. Some variants of the MUSIC algorithm, known as the Root MUSIC algorithms, can be found in [START_REF] Barabell | Improving the resolution performance of eigenstructure based direction finding algorithms[END_REF], [START_REF] Swindlehurst | An analysis of the subspace fitting algorithm in presence of sensor errors[END_REF] and [START_REF] Wong | Root MUSIC based azimuth elevation angle of arrival estimation with uniformly spaced but arbitrarily oriented velocity hydrophones[END_REF].

Estimation of signal parameters via the rotational invariance technique

The ESPRIT method was proposed by [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF]. The ESPRIT method is a narrow-band source localization algorithm, requiring the assumption of a two similar subnetworks sensors network that can be deduced one from the other by translation.

Consequently, the ESPRIT algorithm allows a direct estimation of the sources direction of arrival, using the signal subspace properties. When applicable, this method avoids the costly numerical search of minimals, or maximals of a function (such as for the classical beamforming methods and the MUSIC algorithm).

Consider the Vandermonde matrix V (which is a matrix with the terms of a geometric progression in each row), whose dimensions correspond to the noise subspace dimensions. V up and V down are respectively the matrices that contain the M -1 first lines and the M -1 last lines of V . In the same way, U s-up and U s-down are matrices that contains respectively the M -1 first lines and the M -1 last lines of the signal subspace vectors matrix U s . We can verify that :

V up = V down J, (3.35)
with J = diag e -j(2π f c ∆)sin(θ 1 ) ...e -j(2π f c ∆)(M -1)sin(θ K ) . Note that the columns of V and U s form two bases of the same subspace. Consequently, an invertible matrix G, which is the translation matrix from one base to the other, verifies:

V up = U s G. ( 3
U s-up = Φ ESP RIT U s-down , ( 3.37)
where Φ ESP RIT is a spectral matrix, whose poles are θ k=1...K , defined by:

Φ ESP RIT = GJG -1 . (3.38)
Some other methods use this principle of two (or more) similar subnetworks that are related by a translation or a rotation relations, such as covariance difference methods [START_REF] Paulraj | Eigenstructure methods for direction of arrival estimation in the presence of unknown noise fields[END_REF], [START_REF] Tuteur | A new method for signal detection and estimation using the eigenstructure of the covariance difference[END_REF] and [START_REF] Prasad | A transform-based covariance differencing approach for some classesof parameter estimation problems[END_REF], and the spatial smoothing methods [START_REF] Reddi | On a spatial smoothing technique for multiple source location[END_REF] and [START_REF] Shan | On spatial smoothing for direction-of-arrival estimation of coherent signals[END_REF].

Cyclic subspace decomposition methods

The MUSIC algorithm was adapted to cyclostationary sources in [START_REF] Gardner | Simplification of MUSIC and ESPRIT by exploitation of cyclostationarity[END_REF] and [START_REF] Schell | Cyclic MUSIC algorithms for signal-selective direction estimation[END_REF], and named "cyclic MUSIC". The correlation matrix that has to be decomposed in this algorithm, instead of R, is R α , or R α when available (see 2.2.2 for definitions ). As explained in 2.3.2, if the desired or undesired sources are cyclostationary, the result of the subspace decomposition of the cyclic correlation matrix is a signal subspace depending on only the cyclostationary sources (at the computation cyclic frequency α), thus allowing an estimation of the steering vectors of these sources. Some improvements to the cyclic-MUSIC algorithm are given in [START_REF] Izzo | An interference tolerant algorithm for localization of cyclostationary signal sources[END_REF] and [START_REF] Liu | Unitary cyclic DOA algorithm for coherent cyclostationary signals[END_REF].

Other improvements using the extended cyclic correlation matrix (which is a combination of both cyclic and cyclic conjugated correlation matrix) can be found in [START_REF] Chargé | Traitement d'antennes pour les télécommunications, localisation des sources et autocalibration[END_REF].

The Cyclic-ESPRIT algorithm, presented in [START_REF] Gardner | Simplification of MUSIC and ESPRIT by exploitation of cyclostationarity[END_REF], also employs the cyclic correlation matrix instead of the classic one. It is based on a generalized eigenvalue (when receiving one cyclostationary signal) estimated from the cyclic correlation matrix R α , and the cyclic output matched correlation matrix (corresponding to a correlation between each sensor and a matched sensor in the same plane displaced by a distance ∆):

R α M atched = R α Φ with Φ = e j2πf ∆ sin(θ)/c , (3.39)
where f is the center frequency of the signal, c the speed of propagation and ∆ the distance (these parameters must be known). The generalized eigenvalue will then be λ = 1/Φ * . The direction of arrival θ can be therefore estimated.

Conclusion

For subspace methods, the direction finding procedure requires structured steering vectors, which is not always the case in practice. These methods are sensitive to the errors committed in the modeling of the signals. The theoretically infinite resolution of these algorithms requires exact knowledge of the correlation matrix, which is not realistic, as only an estimate of the correlation matrix is available. The obtained zeros are in fact minimals. For non-cyclostationary methods, an additional issue is the determination of the number of eigenvalues involved in the extraction of respectively the system and noise subspaces. Some rank tracking methods have been proposed to avoid this issue in [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF], [START_REF] Yang | Projection approximation subspace tracking[END_REF] and [START_REF] Abed-Meraim | Fast orthonormal past algorithm[END_REF]. Despite the cost of a subspace decomposition of the correlation matrix, this limitation can be compensated by the interest of the information extracted from the eigen structures.

Conclusion

Several improvements can be made to the classic RFI mitigation methods in radio astronomy using cyclostationarity. For example, the detection methods fail in the presence of a significant cosmic source, due to a non-differentiation between the sources of interest and the interferers, leading to a wrong RFI detection. For subspace approaches, calibration and noise whitening are important and costly steps, which can be avoided by using subspace decomposition of the cyclic correlation matrix instead of the classic correlation matrix. The proposed improvements can be achieved using recent results in telecommunication array processing. Multidimensional telecommunication techniques fall into two categories : beamforming techniques and high resolution methods.

Beamforming methods suffer from several limitations. For classical and adaptive beamforming techniques, the most relevant problem is the computational cost and the extensive a priori knowledge required of the system parameters . Indeed, the correlation matrix can be inaccurately estimated due to limited number of data samples, and the knowledge of the steering vector can be imprecise due to constrained direction errors or imperfect array calibration. The Score methods convergence is very slow, and the algorithms fail in the presence of several self-coherent signals at the computational frequency α. In the rest of this study, among the proposed RFI mitigation approaches, we will not employ beamforming techniques.

For RFI mitigation in the radio astronomy context, high resolution methods seem to perform better than beamforming techniques. Despite the strong assumptions required for linear methods, the orthogonal projector can be interesting after some modifications, to avoid the high risk of errors committed by choosing very close lines, or the presence of nearby sensors or sources. Moreover, it has been shown in subspace decomposition approaches that the spatial signature of the received cyclostationary signal can be estimated via the dominant eigenvector, regardless of knowledge or the architecture of the array. Additionally, the projector is quite similar to the spatial filter in radio astronomy. Thus, improving the projector using the cyclic subspace decomposition results amounts to improving the spatial filter using cyclostationarity. The proposed approach will be described in chapter 6.

In short, the main result obtained from the study of the telecommunications methods described is the improvement provided by the use of the cyclostationarity principle. In the following chapters, our proposed array processing techniques, based on cyclostationarity, will be compared to existing radio astronomy multichannel RFI mitigation methods. In the next chapter, we will present three detection algorithms based on cyclostationarity, which are improvements on the classic detectors presented in this chapter.

Chapter 4

Cyclic detection for phased array radio telescopes

For some applications in radio astronomy, such as time-frequency representation of the observations, the polluted time-frequency slots can be flagged (with a flag (1 or 0) which marks a time-frequency slot as corrupted with radio frequency interferers (RFI)), blanked (timefrequency slots are replaced by fake data, usually zero) or excised (time-frequency slots are removed) from the data without completely losing astronomical information. These techniques are applied when the RFI do not share continuously the same time frequency slots with the signal of interest (astronomical sources). Prior to the flagging/blanking/excision process, a detection procedure must indicate if the time frequency slot is polluted or not.

In detection theory, the problem is rephrased as a hypothesis testing issue:

• H 0 : there is no interference, so the reception model will be: z(t) = A s s(t) + n(t),

• H 1 : at least one interferer is present in the band z(t) = A r r(t) + A s s(t) + n(t), where A r and A s are the spatial signature matrices of respectively the interferers r(t) and the cosmic sources s(t), and n(t) is the system noise vector (see equation (2.10)).

The detection techniques can be conceptually implemented at any point of the processing line. However, such techniques are traditionally applied in radio astronomy to post correlation data (i.e. R).

For example, if we consider the correlation matrix R as the output product of the radio telescope array, the basic approach is to split the whole observation into K smaller independent observations R k , where k is the observation index1 . Consequently, a clean version of R would be:

R clean = k/R k ∈H 0 R k
Now, depending on the degree of knowledge on the involved signals, more or less optimal approaches can be derived. The objective is to define a detection criterion T (z) and a threshold γ.

In the first section, we will define our performance analysis tool. In the second section, from the classical detector described in section 3.1.1, we will derive three cyclic methods based on the cyclic (conjugated) correlation matrix (in this section, the delay τ will be systematically equal to zero). The first one requires knowledge of the cyclic frequency, the second one is a blind approach based on the Frobenius norm. Simulation results assessing the performances of both approaches will also be provided. An application of the blind detector on the WSRT A.2 will be shown as well.

Fisher criterion definitions

The aim of the detection algorithms is the rejection of the interferences so as to obtain clean observations. However, this rejection has to be accurate enough to avoid either rejecting the information of interest (cosmic sources), known as false alarm detection, or not missing the RFIs, known as miss detection.

In order to describe the ability of the detectors to separate between H 0 and H 1 , the Fisher criterion [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF] [29] has been used. It is defined by:

F = (E H 1 [T (z)] -E H 0 [T (z)]) 2 V ar H 1 [T (z)] + V ar H 0 [T (z)] , ( 4.1)
where E[.] and V ar [.] are respectively the mathematical expectation and the variance of the detection criterion T (z). In fact, the Fisher criterion describes the distance between H 0 and H 1 obtained from a given detection criterion T (z). The interesting point is that it does not require the complete statistics knowledge of T (z) under the hypothesis H 0 and H 1 . Mean and variance are sufficient to assess the performance of the criterion T (z) (see figure 4.1)

This approach is well adapted to our RFI scenario where the temporal characteristics of the RFI modulation scheme must be taken into account and cannot be simplified by a Gaussian model.

In the next section, we will present detectors based on cyclostationarity.

Cyclic detection approaches

As shown in section 2.2, the cyclic and cyclic conjugated correlation matrices R α and R α depend asymptotically only on the information related to the cyclostationary interferer. Thus, we can expect better estimations of the RFI indicators, which will reduce the false alarms and will increase the detection rate, compared to classic detectors.

In the classic approach (see section 3.1.1), different criteria based on the eigenvalues of the classic correlation matrix were proposed (assuming a specific model R = R RF I + σ 2 n I). From this principle, the first approach is to apply the same procedure to R α or R α , but in our case, with no hypothesis on the cosmic sources or the system noise. However, these cyclic techniques require prior knowledge of the cyclic frequency α. The cyclic frequency can be known, for telecommunication signals, through previous observations, or via the technical documents. However, we can blindly search for cyclic frequencies (i.e we do not know α), by extending our single dish detector described in section 1.3 to the multidimensional case. All these proposed algorithms will be compared with one another, and with the classic approach.

Description of the eigenvalue approach

In [START_REF] Gardner | Simplification of MUSIC and ESPRIT by exploitation of cyclostationarity[END_REF], it was shown that the rank of the cyclic correlation matrix R α is equal to the number of cyclostationary signals received .Thus, the rank can be estimated from the number of the non-zero eigenvalues. Indeed, if λ i (where i = 1...M ) are the eigenvalues of R α , then λ i=K 1 ...M = 0 (see section 2.3.2). A simple cyclo detector algorithm will be to compare the maximum eigenvalue with a given threshold.

The subspace decomposition detection method is summarized by figure 4.2. From an implementation point of view, subspace decomposition has a heavy computational load. The next subsection presents another approach which is computationally less costly. In the classic approach, and for the corresponding model (R = R RF I + σ 2 n I), frob 2 (R) is similar to the GLRT detector (see section 3.1.1). We transpose this approach to R α with our general model (R = R RF I + R Cosmic + R N oise ).

Description of the norm-based detector

2: The subspace decomposition detection algorithm: After a subspace decomposition of the cyclic correlation matrix R α , we compare the eigenvalues to the threshold in order to obtain the number of K 1 dominant eigenvalues corresponding to the number of the K 1 received cyclostationary RFIs with cyclic frequency α.

Cyclic frequency known

Under H 0 (no interferers: r(t) = 0), the Frobenius norm of the the cyclic correlation matrix frob 2 (R α ) = 0, as R α = 0. Under the H 1 hypothesis,

frob 2 (R α ) = K 1 k=1 λ 2 k ,
where λ k are the eigenvalues resulting from the subspace decomposition of R α (see section 2.3.2), since M -K 1 other eigenvalues are zero. Thus, the cyclic norm-based detector is asymptotically not sensitive to general detection issues (calibration problems and presence of cosmic sources), as only the presence of the RFI is the determining criterion for distinguishing between H 0 and H 1 .

Blind detector for cyclic frequencies retrieval

Knowledge of the cyclic frequencies is a key element for all the RFI mitigation techniques described in this document. Now, we will define a blind detector which can be used when the modulation scheme is unknown. We consider the M × M instantaneous correlation matrix defined by: R(t) = z(t) * z H (t).

(4.

3

)
The algorithm is defined by the following steps:

MHz, which correspond to the GPS cyclic frequencies (multiples of the BPSK symbol duration T sym ). Right hand figure: retrieval of the cyclic conjugated cyclic frequencies: Spectral lines appear for corresponding BPSK cyclic conjugated frequencies, which correspond to twice the carrier frequency 2f 0 and additional multiples of the symbol rate 1/T sym . WSRT data by courtesy of Albert-Jan Boonstra and ASTRON.

1. We consider a set of L successive instantaneous matrices and we stack them in a data cube. The cube size is M × M × L. Each M × M slice represents an instantaneous correlation matrix. The third dimension represents time.

2. Then, we apply a L bin FFT to each vector [z k (t)z * l (t)] t=0,..L-1 (i.e. one row in the third dimension of the data cube). We obtain another data cube where each slice (M × M ) is related to one cyclic frequency α = l/L, l = 0, .., L -1.

3. Finally, we compute the square of the Frobenius norm (frob 2 ) for each slice. We obtain a monodimensional signal which will provide spectral lines at frequencies corresponding to RFI cyclic frequencies:

blind detector (α) = M k=1 M l=1 z k (t)z * l (t)e -j2παt L 2 .
This approach can be seen as a generalization of the retrieval method proposed in section 1.3 for one-dimensional signal. It can also be used as a technique to analyze cyclostationary signals. Obviously, the same approach can be applied to the cyclic conjugated case.

We have tested this blind detector on the Westerbork telescope (A.2). The method has been applied on the same GPS interferer as that used in section 1.3. The GPS are BPSK modulated signals, and thus, the theoretical cyclic and cyclic conjugated frequencies are respectively α = l/T sym , l ∈ ℤ and α = 2f 0 + l/T sym (f 0 and T sym are respectively the carrier frequency and the BPSK symbol duration). Figure 4.3 confirms the result obtained previously in figure 1.7. However, due to the multidimensional approach, the spectral line levels are higher.

CYCLIC DETECTION APPROACHES

The next section is a performance study of the described methods. They will be compared to classic detection methods to show the improvements afforded by the use of the cyclostationarity hypothesis. In particular, the effect of finite length time series will be simulated.

Cyclic detector performances

In this part, the performances of the different cyclic detection algorithms proposed will be studied, and compared to the classic detectors for phased array radio telescopes. Monte Carlo simulations will be performed on different reception scenarios, including a random sorting of the inpinging directions. The Fisher criterion will be used to evaluate the performances of the detectors.

For the simulations, the following signal model will be considered:

H 0 → z(t) = √ 1 -ρa s s(t) + √ ρn(t) H 1 → z(t) = σ r a r r(t) + √ 1 -ρa s s(t) + √ ρn(t) (4.4)
where:

• The spatial signatures a r and a s are chosen randomly : a r or a s = e jφ with φ random value ∈ [0, 2π].

• The RFI is an AM signal : r(t) = b(t)e j2πf 0 t+jφ 0 . b(t) is Gaussian white noise N (0, 1), φ 0 is a random phase and f 0 = 0.3 in normalized frequency.

• The system noise n(t) is assumed independent and identically distributed (i.i.d) N (0, 1) for the moment

• The parameter ρ defines the contribution ratio between the cosmic source and the system noise in the generated signal. If ρ = 1, the scenario involves no cosmic source contribution. The parameter ρ is also related to the signal power to noise power ratio (SNR) through the following equation: SNR= 10 log 10 ( 1 ρ -1). For example, If ρ = 0.9 (10 % contribution of the cosmic source in the total noise power), then SNR = -9.54 dB.

• We define the interference power to cosmic source plus noise ratio (INR) as: INR= σ 2 r (in our simulation, the cosmic source plus noise power is always equal to 1).

The number of sensors in the array is M = 8. The interference to noise ratio (INR) varies from -30 dB to 5 dB. The different correlation matrices will be estimated over a finite set of samples (L = 1024, 8192 and 65536 samples). For each set of parameters (INR, ρ, L), 500 trials are averaged to obtain the final Fisher value.

As the classic method is affected by the presence of cosmic sources and uncalibrated telescope, the three main reception scenarios are based on these issues. To study the effects of each problem separately, we assume first a calibrated system (i.i.d noise) receiving one interferer, but no cosmic source (i.e. ρ = 1). Secondly, a cosmic source will be included in the generated signal with ρ = 0.9, still with a calibrated system noise. Finally, the uncalibrated noise case will be studied in a scenario with no cosmic source.

No cosmic sources, calibrated noise

We compared first the dominant eigenvalue extracted from the subspace decomposition of the classic and the cyclic correlation matrices. Figure 4.4 shows the simulation results. Similar performances were obtained for both the classic and the cyclic approaches, as described in section 2.3.3.

: Fisher criterion comparison (classic vs cyclic approach) for subspace based detectors in the absence of cosmic sources (ρ = 1). The numbers on the right are the sample number, L, used to estimate the classic and cyclic correlation matrices. For example, both detectors perform well down to INR= -15 dB for L = 1024 samples.

The Frobenius norm-based detector, for both the classic and the cyclic case, was compared for the two approaches in the same reception scheme (one interferer and calibrated system noise). We also added a detector based on the Frobenius norm of the cyclic correlation matrix diagonal (diag(R α )) with the cyclic and classic Frobenius norm. The Frobenius norm of diag(R α ) limits the computation complexity.

The performances of the three detectors are presented in figure 4.5. Once again, we obtained similar performances for both the classic and the cyclic full approaches, as described in section 2.3. Even though the diag(R α ) detector performs less well than the classic and cyclic approaches, its simple implementation offers an attractive compromise between performance and required resources. Moreover, by comparing figure 4.4 with figure 4.5, it can be seen that norm based detectors perform a few dB lower that subspace detectors whatever the approach is (cyclic or classic, see figure 4.6). Indeed since the Frobenius norm is the sum of all the eigenvalues, all estimation errors will be added for this detector.

: Fisher criterion comparison (classic vs cyclic approach) for the Frobenius norm based detectors, in the absence of cosmic source (ρ = 1). The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

Cosmic source, calibrated noise

In this part, the same simulations are performed in order to compare the classic and the cyclic detectors behavior in the presence of a cosmic source (ρ = 0.9). In other words, the source power represents a 10% contribution to the total noise power. In figure 4.7, the dominant eigenvalue is extracted, from both classic and cyclic correlation matrices by subspace decomposition. The sensitivity of the classic approach to the presence of a cosmic source is clearly shown, whereas the cyclic approach remains robust even if the scenario departs from the "no cosmic source" scenario. Thus, for the cyclic detector, the separation between H 0 and H 1 can be achieved, for 65536 samples, down to an INR of -24 dB.

: Fisher criterion comparison between the Frobenius norm based detector and the subspace based detector, in the classic approach, in the absence of cosmic sources (ρ = 1). The system noise is i.i.d. The correlation matrices have been estimated over

: Fisher criterion comparison (classic vs cyclic approach) for the subspace based detectors, in the presence of a cosmic source with 10% contribution. The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

In figure 4.8, similar simulations have been set up for the Frobenius norm detectors, and confirm the robustness of the cyclic approach. At low INR, an interesting point is that the classic Frobenius norm detector improves its detection performance compared to the classic subspace approach (see figure 4.9). The reason is that the Frobenius detector behaved more as a power detector than as a detector based on spatial discrimination. It will therefore be less sensitive to any spatial mismatch. In other words, the reason why the classic Frobenius norm detector is less efficient is the "no cosmic scenario". It is also the reason why it performs better in the "presence of cosmic source scenario".

: Fisher criterion comparison (classic vs cyclic approach) for the Frobenius norm based detectors, in the presence of a cosmic source with 10% contribution (ρ = 0.9). The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

In figure 4.10, we pushed the cosmic contribution up to 50%. This is not a realistic scenario, but it provides some information on the detectors limitations. The bad effect on the classic detectors is confirmed. However the interesting point is that the cyclic subspace detector is also affected by the presence of a strong source. Once again, the norm based detectors (right hand figure in 4.10) appear to be more robust for both the classic and the cyclic case, compared with the subspace detectors.

Uncalibrated noise

The third part of the simulations addresses the calibration problem for classic and cyclic detectors. It is assumed that the array receives only one interferer while the system noise is uncalibrated. To generate an uncalibrated system, the following steps has been computed:

• A power variation parameter ρ was set to 0.2.

• For each antenna, a signal was computed, whose values are taken randomly between ± √ ρ, say b i , i = 1, ..., M .

• The uncalibrated noise was generated on each antenna (d i (t) = b i . * n i (t), where . * is the element-wise multiplication).

• The total noise power σ 2 N = M i=1 σ 2 i was computed, where the σ 2 i are the power of each corresponding d i (t).

• The corresponding power of the RFI σ 2 r was evaluated for each INR:

σ 2 r = (IN R * σ 2 N )/M .
• The RFI (the AM signal) was generated with power of σ 2 r .

• The same simulation steps as in the previous section were then followed (classic and cyclic correlation matrices subspace decomposition, Fisher criterion computation, etc...).

To evaluate the detectors performances, the Fisher criterion was computed for both classic and cyclic approaches (figure 4.11). It shows that calibrated noise is an important assumption for the classic approach. Even for large INR (0 to 5 dB), the extracted main eigenvalue (the subspace detection criterion) is corrupted in the classic case by the noise uncalibration. Once again, the robustness of the cyclic approach is demonstrated. Since array calibration needs clean observations, the cyclostationary approach could be an interesting way to provide such information to the array calibration process.

Conclusion

Different detectors based on cyclostationarity have been proposed, and their performances analyzed. It has been shown that asymptotically they are not sensitive to general detection issues (calibration problems and presence of cosmic sources), as only the presence of the RFI is the determining criterion for distinguishing between H 0 and H 1 .

Conclusion

In this chapter, three new detectors based on cyclostationarity have been proposed:

• The first algorithm is based on the subspace decomposition of the cyclic correlation matrix.

• The second algorithm is based on the Frobenius norm of the cyclic correlation matrix.

• The third algorithm is a blind detector for cyclic frequencies retrieval.

These algorithms have demonstrated by simulation their robustness against the presence of cosmic sources or uncalibrated noise. These results are based on the fact that the cyclic correlation or cyclic conjugated correlation matrices depend asymptotically only on the cyclostationary RFI.

In practice, cyclostationary detectors could be useful when looking for clean frequency bands prior to the array calibration. Another interest would be radio spectrum analysis and RFI modulation characteristics retrieval, during RFI site monitoring and survey. Moreover, the ability of such detectors to discriminate between stationary and cyclostationary signals could also be used for Pulsar2 observations. In particular, giant pulsar pulses are interesting phenomena for radio astronomers but are difficult to observe since most of the time they are completely random. By combining a classical detector together with a cyclic one, the observations could be flagged only in the presence of giant pulses even in the presence of RFI. In this case, the cyclic detector is used to identify false alarm detections triggered by RFI signals

The ability of the cyclostationary approach to perform efficient subspace decomposition will be used in the following sections to derive first an "estimation and subtraction" method (chapter 5), then a cyclic spatial filtering method (chapter 6).

Chapter 5 Estimation and subtraction method

The raw output product of a phased array radio telescope is the classic correlation matrix, R(τ = 0). With this classic correlation matrix, images or skymaps can be derived1 .

As explained in the introduction of the previous chapter, in the presence of radio frequency interferers (RFI), a clean version of R can be obtained by discarding or blanking intermediate correlation matrices R k detected as being polluted. However, such an approach will also remove part of the information and hence, it will decrease the sensitivity of the observations. Moreover, this detection and blanking process does not take full advantage of the spatial processing possibilities of phased array radio telescopes.

In this chapter, we propose an estimation and subtraction (E & S) approach. It is based on knowledge of the RFI modulation model and on the ability of the cyclostationary approach to extract the unknown parameters of this model from observations. Each kind of modulation requires a specific study. In this thesis, only the AM and the BPSK modulations have been evaluated. Thus, in the next section, we will present the algorithm for these two modulations. As the method is especially dedicated to interferometry, skymaps simulations will be given to illustrate the method. Finally, an example of spectrum estimation using data obtained with the low frequency array (LOFAR) will be presented.

The E & S algorithm

Figure 5.1 illustrates the principle of the E & S method. The aim of the method is the estimation of the interferer correlation matrix, R RF I , using the cyclic correlation matrix of the data, R α , which contains interferers only (R α RF I ). R RF I will be then subtracted from the global correlation matrix R, to arrive at R cleaned .

From equation (2.25) in section 2.3.2, we have shown that the subspace decomposition of the cyclic correlation matrix R α or the cyclic conjugated correlation matrix R α provides straightforward information on the spatial signature, but is also linked with the monodimensional cyclic correlation of the RFI R α RF I , assuming only 1 interferer:

λ 1 = R α r 1 ,r * 1 (τ ) a r 1 2 u 1 = a r 1 a r 1 -1 (5.1)
In the following, we will investigate how to use these formulae in the case of an AM interferer and a BPSK interferer. • Amplitude modulated interferers:

Assume that we receive 1 AM modulated interferer with a carrier frequency f 0 and carrier phase ϕ 0 . In this case, the global correlation matrix will be:

R = σ 2 r a r a H r +R Cosmic + R N oise , R AM (5.2)
where σ 2 r is the power of the AM interferer, a r is the spatial signature of the RFI, R Cosmic and R N oise are respectively the cosmic sources and the system noise correlation matrices.

From section 1.2.1, we know that AM signals express their cyclostationarity properties only through the cyclic conjugated correlation. Thus, by combining equations (1.15) and (2.24), the asymptotic expression of the cyclic conjugated correlation matrix will be (for τ = 0):

R α = σ 2 r e j2φ 0 a r a T r = R α AM , (5.3)
From equation (5.2), we have: R AM = σ 2 r a r a H r . In order to reconstruct R AM , we need to extract σ 2 r , φ 0 and a H r from (5.2) and (5.3) :

• We saw in section 2.3 that, using eigenvalues analysis, we can estimate a r . This estimate is noted a r .

• We can estimate the RFI power σ 2 r with the following operation:

σ 2 r = 1 M 2 frob(R α ./ a r a T r ),
where ./ is the element-wise division.

THE E & S ALGORITHM

• We can estimate the carrier phase

φ 0 = 1 2M 2 arg   M i=1 M j=1 R α ./ a r a T r i,j  
Thus, in the case of AM modulation, we can estimate the correlation matrix of the interferer using one of the following two expressions:

R AM = σ 2
r a r a H r .

(5.4)

R AM = e -j2 φ 0 R α a * r a T r a * r -1 a H r .
(5.5)

• BPSK modulated interferers:

As R α is non zero for BPSK modulated signals, the same processing can be done for 1 received BPSK interferer with the following emission characteristics: carrier frequency f 0 , carrier phase ϕ 0 , emission filter with impulse response h(t) ∈ ℝ and impulse duration T sym .

Using the expressions of a BPSK signal statistics given in section 1.2.2, we obtain:

R(τ ) = σ 2 r T sym r h,h * a r a H r +R Cosmic + R N oise , R BP SK (5.6) with: r h,h * = +∞ -∞ h(t)h * (t)dt.
The cyclic conjugated correlation matrix will be (asymptotically, for α = 2f 0 + l/T sym , l ∈ ℤ):

R α = σ 2 r T sym e j2φ 0 e -j2π l Tsym t 0 r l Tsym h,h a r a T r = R α BP SK . ( 5.7)
with t 0 is the reception time delay, and

r l Tsym h,h = +∞ -∞ h(t)h(t)e -j2π l
Tsym t dt and l ∈ ℤ.

Since h(t) ∈ ℝ, r h,h * = r h,h . For l = 0, equation (5.7) becomes:

R α=2f 0 = σ 2 r T sym e j2πφ 0 r h,h a r a T r = R α BP SK . ( 5.8)
Thus:

σ 2 r T sym r h,h = G = 1 M 2 frob(R 2f 0 ./ a r a T r ).
(5.9)

Applying the same method as for the AM signals, we can estimate R BP SK with either:

R BP SK = G a r a H r (5.10)
or:

R BP SK = e -j2 φ 0 R α a * r a T r a * r -1 a H r .
(5.11) Algorithms 1 and 2 summarize the two E & S algorithms for AM and BPSK signals with cyclic frequency α = 2f 0 . The aim of the two algorithms is the estimation of the RFI correlation matrix from the cyclic correlation matrix, in order to subtract it from the global correlation matrix. The two approaches differ in the required parameters for the estimation of the RFI correlation matrix.

In the BPSK case, if we assume knowledge of the emission filter h(t), it is possible to use other cyclic frequencies α = 2f 0 + l/T sym . With this assumption, it is also possible to base the E & S method on the cyclic correlation matrix at α = l/T sym for BPSK, m-PSK and QAM signals. In the next section, we will apply the E & S technique on sky imaging.

Algorithm 1 First algorithm of the estimation and subtraction method

Decompose R 2f 0 into subspaces, Estimate the RFI spatial signature ( a r = u r 1 ) Estimate amplitude of the RFI σ 2 r or G:

σ 2 r (or G)= 1 M 2 frob(R α .
/ a r a T r ) Estimate the correlation matrix relative to the RFI:

R RF I = G a r a H r Subtract R RF I from R: R cleaned = R -R RF I
Algorithm 2 Second algorithm of the estimation and subtraction method Decompose R α into subspaces, Estimate the RFI spatial signature ( a r = U r 1 ) Estimate the carrier phase φ 0 :

φ 0 = 1 2M 2 arg   M i=1 M j=1 R α ./ a r a T r i,j  
Estimate the correlation matrix relative to the RFI:

R RF I = e -j2 φ 0 R α a * r a T r a * r -1 a H r Subtract R RF I from R: R cleaned = R -R RF I

Simulation results on the E & S method

In this section, we will illustrate the E & S technique through simulations on skymaps. For the phased array, we will use the ITS configuration with M = 60 antennas (see appendix A.1).

Figure 5.2 shows a corrupted observation compared to an interference free observation. We will compare the cleaned results to this reference figure. As there is no such equivalent technique to the E & S method using the classic correlation matrix R, we cannot compare the performance of this approach with an existing one. However, error maps (cleaned mapreference map without RFI) will be computed in order to evaluate the E & S method.

SIMULATION RESULTS ON THE E & S METHOD

One interferer

We performed the E & S method on a simulated observation, containing 1 white Gaussian cosmic source with steering vector a s 1 pointing towards the zenith (i.e. here the spatial signatures are structured relative to the ITS configuration), and 1 BPSK modulated RFI whose carrier frequency f 0 = 0.3 (normalized frequency), its baud rate is 1/T sym = 1/8, the emission filter is rectangular, its power is σ 2 r 1 and its steering vector is a r 1 . The cyclic correlation matrix that has been computed is R α=2f 0 . Three cases of INR were selected: 10 dB (strong interferer), 0 dB (weak interferer) and -10 dB (very weak interferer). The signal subspace was estimated using the strongest eigenvalue. Figure 5.3 shows the results of the simulation. For all INR values, the directions of the RFIs have been correctly estimated, and the corrupting signal has been subtracted from the correlation matrices.

Presence of several interferers

We performed the E & S method on a simulated observation containing the same Gaussian cosmic source, and 3 BPSK modulated RFIs with different emission characteristics (different carrier frequencies and baud rates). The corresponding steering vectors of each one of the 3 BPSK signals are a r 1 , a r 2 and a r 3 .

The first BPSK carrier frequency is f 1 0 = 0.3 and the baud rate is 1/T sym = 1/8. The second BPSK carrier frequency is f 2 0 = 0.48 and the baud rate is 1/T sym = 1/16. The third BPSK carrier frequency is f 3 0 = 0.25 and the baud rate is 1/T sym = 1/32. The carrier frequencies are normalized. The emission filters for the 3 BPSK are rectangular. The interferers have the same power σ 2

r 1 = σ 2 r 2 = σ 2 r 3 .
Three cyclic correlation matrices were computed: R α=2f 1 0 , R α=2f 2 0 and R α=2f 3 0 . For each

3: Estimation and subtraction method with 1 BPSK interferer. Left to right: The upper figures represent initial skymaps for respectively a strong (INR= 10 dB), a weak (INR= 0 dB) and a very weak interferer (INR= -10 dB). The correlation matrix has been computed over 8192 time samples. The midline figures represent the observations after applying the E & S method: the RFI has been correctly subtracted and is no longer visible on the skymaps. The lower figures represent the error maps.

REAL DATA RESULTS USING THE E & S METHOD

(b) The observation after applying the E & S method using 3 different cyclic frequencies: the RFIs are almost completely removed, and the shape of the map is fairly well preserved. (c) The error map.

cyclic correlation matrix, the eigenvector corresponding to the strongest eigenvalue was extracted. The INR is 0 dB (all the signals have the same level of power).

Figure 5.4 shows the results of the simulation. Since each column of the estimated 3 dimensions signal subspace can be identified as a r 1 , a r 2 and a r 3 , the estimated directions correspond quite closely to the directions of the RFIs.

Real data results using the E & S method

In this section, we present the test of our E & S method applied on LOFAR (see appendix A.1) data. Since only the output signal from 8 stations (i.e. M = 8) was available, we will provide a frequency spectrum rather than a skymap, as the number of sensors is too small, leading to low resolution of the maps. The estimation results are based on a real AM signal received in the 0 -20 MHz band of LOFAR. The different steps of our test are:

• First, the AM signal was filtered from the output data vector of one antenna, at the central frequency 16.7 MHz (figure 5.5) using a passband filter of 200 KHz bandwidth. The same process was carried out on the 8 antennas.

• After filtering the signal from each antenna, we built a new output data vector z AM (containing the 8 obtained signals), and then computed the classic instantaneous correlation matrix of z AM .

• Using the blind detector based on the Frobenius norm (described in section 4.2.2), we estimated the cyclic frequency α = 0.15 (normalized frequency).

• Then, we computed the cyclic correlation matrix of z AM , to extract the eigenvector corresponding to the dominant eigenvalue.

• The carrier phase is therefore estimated (according to the second E & S algorithm). We obtained φ AM = 0.8147.

• We applied the RFI estimation method described in [START_REF] Bretteil | Estimation spectrale en temps réel de sources radioastronomiques en présence de brouilleurs cyclostationnaires: Application au radiotélescope décimétrique de Nançay[END_REF] to estimate the power spectrum of z AM from the instantaneous cross-correlation function R zx,z * y (t, τ ) (where x and y refer to two different antennas).

• The AM signal spectrum is then estimated from the Fourier transform of the cyclic cross-correlation function at α (figure 5.6). To confirm the results obtained, we added a fake source (which is not cyclostationary at α), and repeated the process. Figure 5.7 shows that the AM spectrum was correctly estimated, without being affected by the presence of the fake source.

Conclusion

In this chapter, we have proposed an estimation and subtraction (E & S) approach. It is based on knowledge of the RFI modulation signature, and on the ability of the cyclostationary approach to extract the unknown parameters of this signature from observations. However, this technique is limited by its applicability. While it can be used on AM or BPSK, for m-PSK or QAM signals, some knowledge of the signals characteristics is required such as the emission filter characteristics, which can be known through literature. Further research is necessary to conclude about the applicability of this technique to other modulation schemes.

Another point for further study would be to evaluate the exact distortion of the cosmic source signal induced by such an approach. In particular, it could be interesting to compare it with distortions induced by spatial filtering methods, since spatial filtering is a more general technique that can be applied on almost all interferers. In the next chapter, we will see how classic spatial filtering can be extended to cyclic spatial filtering.

Chapter 6

Cyclic spatial filtering

Cyclic spatial filtering is an improvement on the classic spatial filtering detailed in 3.1.2. It is a general method that can be applied either on narrow-band or broad-band signals, and on signals for which either the cyclic correlation matrix R α or the cyclic conjugated correlation matrix R α is available.

We will start the chapter with a presentation of the general cyclic spatial filtering algorithm. The error on the estimation of the spatial signatures will be computed, for both classic and cyclic spatial filtering. The performances of the two approaches will be compared via the study of the remaining INR after filtering, in comparison with the initial INR, for small arrays and for a very large number of sensors. Skymap simulation results will be presented, to illustrate the impact of spatial filtering. Results using LOFAR real data will be presented to test the practical effectiveness of cyclic spatial filtering compared to the classical approach. The chapter concludes with an extension of the algorithm to the broad-band case, with an illustration using a real broad-band interferer received by LOFAR.

General algorithm of cyclic spatial filtering

Consider receiving K 1 interferers (cyclostationary with cyclic frequency α) and K 2 cosmic point sources (white Gaussian) with a M sensors array. The system noise is independent Gaussian complex entries. In accordance with section 2.3.2, we can estimate the spatial signatures using subspace decomposition of the estimated cyclic correlation matrix over L time samples, or the cyclic conjugated correlation matrix computed over L time samples, and τ = 0:

R α = z(t)z H (t)e -j2παt L = U cyclic Λ cyclic V H cyclic , R α = z(t)z T (t)e -j2παt L = U cyclic Λ cyclic V H cyclic , ( 6.1)
where z(t) = A r r(t) + A s s(t) + n(t) is the telescope output vector (see section 2.1 for definitions), and α is the cyclostationary parameter which characterizes the RFI signals. Since the cosmic sources signal s(t) and the system noise signal n(t) are not cyclostationary with the cyclic frequency α, the cyclic correlation matrix will depend asymptotically on RFIs only:

R α = R α RF I + R α Cosmic + R α N oise . 0 0 R α = R α RF I + R α Cosmic + R α N oise . 0 0 (6.2)
Consequently, the signal subspace formed by the K 1 largest singular vectors

U r K 1 = [u 1 ...u K 1 ] (with U r K 1 = U cyclic or U cyclic depending on the decomposition of R α or R α
) will span the same subspace as the RFI spatial signatures A r .

Once these spatial signatures have been estimated, the RFI signals can be filtered out by applying a projector on the telescope outputs:

z cleaned (t) = P cyclic z(t), (6.3)
where P cyclic is the cyclic spatial projector defined by

P cyclic = I -U r K 1 U H r K 1
and I is the M × M identity matrix Since A r spans the same subspace as U r K 1 , A r can be written as a linear combination of U r K 1 =⇒ A r = U r K 1 ϑ (ϑ is an arbitrary normalized vector). This leads to the following result when applying the projector on A r :

P cyclic A r = (I -U r K 1 U H r K 1 )U r K 1 ϑ = 0. (6.4)
It follows that equation (6.3) becomes, using equation (6.4):

z cleaned (t) = P cyclic (A s s(t) + n(t)) (6.5) This result represents a cleaned output signals vector. The bias induced by P cyclic in Equ(6.5) can be removed as explained in [START_REF] Van Der Tol | Performance analysis of spatial filtering of rf interference in radio astronomy[END_REF].

Alternatively, the same projector can be applied directly on the correlation matrix, similarly to the classic approach described in section 3.1.2. Thus, the cleaned correlation matrix will be obtained by: R cleaned = P cyclic RP cyclic , (

where:

R = A r R r A H r + A s R s A H s + R n . = R RF I + R Cosmic + R noise (6.7)
The next section will demonstrate by simulations the effectiveness of this approach through a comparison with the classical approach.

Performances of the spatial filtering method

In this section we will compare the performances of the classic and the cyclic approaches. We consider the case where we have K 1 = 1 interferer and K 2 = 1 cosmic source, using an array of M = 8 (subsections 6.2.2 and 6.2.2) or M = 1000 (subsection 6.2.3) antennas. The RFI is a BPSK modulated signal, with the following characteristics: the carrier frequency is f 0 = 0.3 (normalized frequency), the baud rate is 1/T symbol = 1/8, the emission filter is rectangular, its power is σ 2 r 1 and its spatial signature is a r 1 . The cosmic source power is σ 2 s 1 . Without loss of

Algorithm 3 General Algorithm of cyclic spatial filtering

Decompose R α (or R α ) into subspaces Estimate the RFI spatial signature ( A r = U r K 1 ) Calculate the cyclic spatial projection operator

P cyclic = I -U r K 1 U H r K 1
Apply this operator on the sensors output z(t) to filter out the RFI: z cleaned (t) = P cyclic z(t) The obtained cleaned output signal will be:

z cleaned (t) = P cyclic (A s s(t) + n(t))
generality, the system noise power in the simulations for each antenna has been set to σ 2 n . We consider the correlation matrix defined by:

R = σ 2 r 1 a r 1 a H r 1 + σ 2 s 1 a s 1 a H s 1 + σ 2 n I . R r 1 D (6.8)
The first subsection will study the estimation errors of the spatial signature, in the presence of a cosmic source. In the second subsection, the effects of the projector on the filtered correlation matrix are analyzed. Finally, the third subsection examines the behavior of the classic and cyclic approaches in the case of a very large number of sensors, which will be achieved for the SKA telescope.

The inner product as performance indicator

Performances are analyzed through the inner product, which is the product of the true RFI spatial signature a r 1 with the estimated one (the strongest eigenvector u 1 ), in accordance with the variations in the INR. This inner product, a H r 1 u 1 , is equal to 1 if the spatial signature is accurate.

In figure 6.1, the performance and the limitations of both classic and cyclic spatial filtering for the 50% source case are shown. A 50% source contribution indicates that the source to noise ratio (SNR) is 0 dB, which means that the system noise and the cosmic source have the same power.

For IN R > 0 dB, the main eigenvector extracted using the two methods is very close to the true spatial signature of the interferer, since the total noise power contribution ( the cosmic source + system noise) is (much) less than the power of the RFI. Thus, the power in the first eigenvalue is dominated by the RFI, and the corresponding eigenvector relates to a direction vector pointing in the direction of the RFI.

As the INR decreases, the performance of the classic approach drops off, since the array signature vector is altered by the cosmic source contribution. Thus, below an INR of -5 dB, the spatial signature of the interferer is badly estimated. This decreasing performance is consistent with the results shown in figure 6.2.

The cyclic approach is more robust to this cosmic source contribution.

The remaining INR

We define the Interference to Noise Ratio before any spatial filtering by:

IN R bef ore = tr(R r 1 ) tr(D) = σ 2 r 1 σ 2 s 1 + σ 2 n . (6.9)
We also define the Interference to Noise Ratio after spatial filtering by:

IN R af ter = tr( PR r 1 P) tr(D) = σ 2 r 1 tr( Pa r 1 a H r 1 P) σ 2 s 1 + σ 2 n , ( 6.10)
where P is the estimated projector obtained through either the classic approach P classic or the cyclic approach P cyclic . Figure 6.2 shows the computed IN R bef ore as a function of the IN R af ter . The projector P is estimated from time series with length L = 8192. We have also used different ratios between the cosmic source power and the system noise power in the total noise power contribution. For example, a 50% contribution indicates that the system noise and the cosmic source will have the same power.

When the projector is badly estimated, the RFI is not filtered out, then, apart from a small bias, IN R bef ore = IN R af ter . This bias (≃ 0.6 dB) can be explained by the fact that even though the interferer spatial signature is wrongly estimated, the projector will remove a subspace dimension and, hence, a small fraction of the RFI power. When the number of antennas increases, this offset decreases.

When we have only system noise (i.e. 0% source), the two methods yield similar simulation results. This is due to the fact that in the classic (asymptotic) approach, the estimated eigenvector does not depend on the (whitened) noise power, assuming equal noise powers for each of the antennas.

However as we add a cosmic source, differences between the two approaches become apparent. This result can be explained by the fact that in the classic spatial filtering case, unlike the cyclic case, the estimated largest eigenvector is influenced by the source. Similar results are observed when the previous equal noise powers assumption is not satisfied (see the last 2 curves on figure 6.2).

For large IN R bef ore , we retrieve an extreme case [79] [7]:

IN R af ter = 1 L 1 + 1 M IN R bef ore . ( 6.11)
For L = 8192, the theoretical limit (Equ.6.11) is -39.09 dB, which is close to that reached in the simulation:-39.83 dB.

Very large number of sensors

The simulation presented in figure 6.3 was obtained with a 1000 sensors array. We considered the sources as a very weak contribution in the total noise power (1%). We kept the same method for estimating the INR before and after applying the projector, and the length for the time series was similar to that in the previous subsection (L = 8192). In the case of a very large number of sensors (M = 1000), the array will be more sensitive to the received signals. We obtain better performances for both classic and cyclic approaches compared to the 8 sensors configuration. However, the cyclic approach in this simulation clearly outperforms the classic one. This quite satisfying result indicates that, in the case of an array with a very large number of sensors (SKA for instance), the greater sensitivity to RFI when using the cyclic technique may lead to a better filtering, and thus, an enhanced robustness of the method.

For L = 8192, the limit reached in the simulation is -39.17 dB, which is very close to the theoretical limit (Equ.6.11): -39.1339 dB. The simulations above show that, for low INR, and in the presence of a cosmic source, the cyclic method outperforms the classic one. In the next section, we will demonstrate with simulations (skymaps) and experimentally (using real LOFAR telescope data), that for large INR the cyclic method also has advantages over the classic spatial filtering approach.

Simulation results on the SF method

In this section we will use the ITS configuration to build skymaps in order to illustrate the impact of the cyclic spatial filtering method (using the estimated cyclic correlation matrix R α ). We will compare it to classic spatial filtering (using the classic correlation matrix R).

Figure 6.4 provides the initial skymaps with and without RFI. The latter will be used as a reference skymap.

One interferer

We performed classic spatial filtering (SF) and cyclic spatial filtering (CSF) on an observation containing:

• 1 white Gaussian cosmic source with spatial signature a s 1 ,

• 1 strong BPSK modulated RFI. The carrier frequency is f 0 = 0.3 (normalized frequency), baud rate is 1/T symbol = 1/8, the emission filter is rectangular, the power of the RFI power is σ 2 r 1 and its spatial signature is a r 1 .

The cyclic correlation matrix that was computed is R α=2f 0 .

Three cases of INR have been selected: 10 dB (strong interferer), 0 dB (weak interferer) and -10 dB (very weak interferer). The signal subspace was estimated using the strongest eigenvalue. Figure 6.5 shows the results of the simulation. For strong interferers, the maximum eigenvalue of R and R α is related to the RFI only. The estimated signal subspace is a good estimate of the interferer spatial signature. However, as the INR decreases, the performance of classic spatial filtering drops. Indeed, as the power of the interferer declines, the power of the cosmic source becomes relatively strong. The extracted eigenvector, which is mostly related to the strongest signal (here the cosmic source) gives a wrong estimation of the RFI direction. For very weak interferers, skymaps show that the cosmic source is completely removed. For CSF, as the strongest eigenvector is always related to the RFI only, the method performs well even in the presence of very weak interferers. The results obtained in the previous section are confirmed with the skymaps.

Presence of several interferers

We performed the SF and the CSF methods on an observation containing 1 white Gaussian cosmic source with spatial signature a s 1 and 3 BPSK modulated RFIs with different emission characteristics; their respective spatial signatures are a r 1 , a r 2 and a r 3 .

The first BPSK carrier frequency is f 1 0 = 0.3 and the baud rate is 1/T sym = 1/8. The second BPSK carrier frequency is f 2 0 = 0.48 and the baud rate is 1/T sym = 1/16. The third BPSK carrier frequency is f 3 0 = 0.25 and the baud rate is 1/T sym = 1/32. The carrier frequencies are normalized. The emission filters for the 3 BPSK are rectangular. The interferers have the same power σ 2 r 1 = σ 2 r 2 = σ 2 r 3 . The INR is 0 dB (all the signals have the same level of power). The three cyclic correlation matrix that were computed are: R α=2f 1 0 , R α=2f 2 0 and R α=2f 3 0 . For each cyclic correlation matrix, the eigenvector corresponding to the strongest eigenvalue was extracted, and merged to form 3-dimension subspace. For R, the signal subspace was formed using the 3 main eigenvectors resulting from the subspace decomposition of R. Figure 6.6 shows the results of the simulation. For the classic spatial filtering method, the 3 main eigenvectors extracted from the subspace decomposition of R contain a mix of all the signals since no specific directions pop up. Thus, all the signals are partially filtered out. For cyclic spatial filtering, as each column of the 3-dimensions estimated signal subspace can be identified as a r 1 , a r 2 and a r 3 , the estimated directions correspond very closely to the directions of the RFIs. We compared the case of different cyclic frequencies using CSF with the Estimation and Subtraction approach (chapter 5). The E & S method generates less distortion than cyclic spatial filtering but it needs particular configurations as explained in chapter 5. The advantage of the CSF method is that it can be applied on practically all kinds of RFIs using either R α or R α .

Real data results using the SF method

We have applied classic and cyclic spatial filtering to real observations using the LOFAR radio telescope (see appendix A.1).

In figure 6.7, we have observed in the 160 -240 MHz LOFAR band, which contains a very strong transmitter (a pager) at 170 MHz with an INR of 47 dB. The array configuration consisted of M = 8 LOFAR antennas. The cyclic frequency, α, of the pager was first estimated from the data, using the blind detector described in section 4.2.2. We obtained α = 0.1221 in normalized frequency (figure 6.7). Figure 6.8 represents the eigenvalues obtained from the classic and the cyclic correlation matrices (resp. R and R α ) which were derived from baseband data of the 8 antennas. This figure shows that the interferer signal subspace can be fairly well estimated using one dimension in the cyclic decomposition, whereas it needs two dimensions in the classic one. The more dimensions that are used to remove the interferer, the more information about the cosmic sources

(b) Results of the application of classic spatial filtering: the 3 RFIs have been filtered but part of the cosmic source is also removed. (c) Results of CSF with 3 different interferers (3 different cyclic frequencies): the interferers have been filtered out. (d) Comparison with the E & S method: in addition to removing the RFI, the shape of the map has been preserved.

BROAD-BAND ISSUE

that is thrown away in the same way. We therefore used only the eigenvector corresponding to the strongest eigenvalue to build the projector for both methods. Figure 6.9 demonstrates the effect of the projector on the pager. For both approaches (clasic and cyclic) only one dimension has been removed. Using the cyclic method, the pager has been removed more effectively than by using the classic approach. The main reason is that the spatial filter is applied to uncalibrated data; unlike the classic spatial filter, the cyclic method in not dependent on calibration.

Broad-band issue

In all the previous sections and chapters, we assumed that the narrow-band assumption holds. In this section, it will be shown that the spatial filtering method can be generalized to the broad-band case. The modifications made to the algorithm will be first explained, then an illustration using LOFAR real data will be presented.

Broad-band algorithm

In the case where the narrow band assumption is no longer verified, the algorithm described in 6.1 has to be modified. The proposed approach is based on an algorithm from Xu and Figure 6.9: Spectrum of one antenna output after applying cyclic and classic spatial filtering (M = 8, L = 65536).

BROAD-BAND ISSUE

Kailath [START_REF] Xu | Direction of arrival estmation via exploitation of cyclostationarity -a combination of temporal and spatial processing[END_REF]. For purposes explanation, we assume that there is no cosmic source and no system noise.

If the received interferer is a broad-band signal, then: .12) where γ i,θr is the sensor beam response and the electronic chain response (see section 2.1.1 and figure 2.1).

z i (t) i=1...M = γ i,θr r(t -τ θr,i ). ( 6
According to equation (2.2), we can express r(t -τ θr,i ) with: r(t -τ θr,i ) = (t -τ θr,i )e -j2πf 0 τ θr ,i e j2πf 0 τ θr ,i t , τ θr ,i (t) (6.13) where τ θr ,i (t) is the complex envelope expression of the τ θr,i delayed signal r(t -τ θr,i ) Its cyclic correlation with cyclic frequency α is [START_REF] Xu | Direction of arrival estmation via exploitation of cyclostationarity -a combination of temporal and spatial processing[END_REF]:

R α z i (τ ) = R α r (τ )e -j2πατ θr ,i . (6.14)
This equation is just a straightforward application of the time delay and phase shift duality transposed to cyclic frequency. Let us now define the global output cyclic correlation z α (τ ):

z α (τ ) = R α z 1 (τ ), R α z 2 (τ ), ..., R α z M (τ ) T = a(α)R α r (τ ), (6.15)
where a(α) = γ 1,θr e -j2πατ θr ,1 , ..., γ M,θr e -j2πατ θr ,M T is called the pseudo-spatial signature 1 .

A similar expression was obtained when expressing the sensor output with the narrow band condition (eq.(2.4)): z(t) = a r r(t), (6.16) with a r = γ 1,θr e -j2πf 0 τ θr ,1 , ..., γ M,θr e -j2πf 0 τ θr ,M T .

Thus, the RFI mitigation techniques applied to z(t) can be transposed to z α (τ ). However, the key point is that in the general case (cosmic sources and/or uncalibrated noise), z α (τ ) content, asymptotically, only information on the RFI.

The sample covariance matrix will be:

R sample = z α (τ )(z α (τ )) H N ,
where N is the number of delays τ (in practice τ = kT S , k = 0, ..N -1. T S is the sampling period).

By applying a subspace decomposition on this sample covariance matrix, we can estimate the pseudo spatial signature a(α) using the eigenvector u r corresponding to the dominant eigenvalue. Thus, if the array is calibrated (i.e. γ i,θr = γ θr , for all i = 1, ..M ), then: u r = 1 γ θr a(α) = e -j2πf 0 τ θr ,1 , ..., e -j2πf 0 τ θr ,M T

u(f k ) = [u r ] f k α
Calculate the cyclic spatial projection operators:

P cyclic (f k ) = I -u(f k )u H (f k )
Apply this operator either on the sensors output z(t) to filter out the RFI: z cleaned (t) = P cyclic (f k )z(t), or on the correlation matrix of the output vector R = z(t)z H (t) L : R cleaned = P cyclic (f k )RP cyclic (f k ).

Figure 6.11 shows that the interferer signal subspace can be fairly well estimated using the broad band approach, whereas it needs at least 6 dimensions in the classic one ! If we use all these dimensions to remove the interferer in the classic case, almost all the information about the cosmic sources will be thrown away. However, from the subspace decomposition of the classic correlation matrix, it can be seen that, due to the spread of the eigenvalues, and to the power variation of the DAB from one antenna to another, a pre-calibration of the data must be achieved before spatial filtering, in order to remove one dimension from the classic correlation matrix using the broad-band spatial filter. Figure 6.12 shows the effect of the projector on the DAB using the cyclic method. In fact, the projector was computed for only one frequency f k chosen as the center frequency of the 4 MHz band, and was applied as it is on the 4 MHz band. The DAB was removed efficiently even though the narrow band condition is not completely fulfilled in this example. Note that the input data was calibrated before applying the algorithm since the power of the DAB varies from one antenna to another.

Conclusion

Cyclic spatial filtering is a general method that can be applied to almost all kinds of interferers whether they are narrow-band or broad-band. Except for the broad-band case, the technique is independent of calibration issues and the presence of cosmic sources. The use of cyclostationarity in the spatial filter described in this chapter leads to better performances compared to existing methods. Moreover, it has been shown that in practice, the multipath issue can be overcome when using the cyclic correlation matrix instead of the classic one. Indeed, the more dimensions that are removed, the more useful data that is lost. The number of corrupted dimensions in the classic correlation matrix is in practice superior to the dimensions that have to be removed from the cyclic correlation matrix, even for very strong interferers, while performance simulations have shown that the two techniques are equivalent for high INR. For broad-band interferences, the structure of the cyclic correlation matrix makes it possible to build a pseudo-data correlation matrix, which takes into account the broad-band characteristics of the interferer. Thus, only one dimension (if we receive one broad-band interferer) is removed. These kinds of approaches seem to be very promising and should have priority for further testing and research in the framework of SKA.

However, these good performances are obtained if and only if the cyclic frequency is accurately estimated, using the cyclic frequency finder algorithm described in 4.2.2, or retrieved from the interferer description (which can be found in the literature).

Furthermore, the spatial filter adds deep zeros in the skymaps, while the E & S method does not change the maps, but can be applied only under certain conditions described in the previous chapter. The appropriate method has to be chosen according to the application.

Conclusion

La présente étude a été menée dans le cadre de SKADS [START_REF]Square kilometer array desing studies[END_REF] pour élaborer des méthodes de suppression d'interférences. L'originalité du travail repose sur l'utilisation de la propriété de cyclostationnarité des signaux de télécommunication dans le cas du traitement spatial de ces interférences. Cette étude s'appuie sur le calcul de la matrice de corrélation cyclique ou cyclique conjuguée qui correspond à une mesure de périodicités dans les matrices de corrélations instantanées obtenues par corrélation successive du réseau d'antennes. Dans ce document, plusieurs approches ont été proposées, cependant elles reposent toutes sur le fait que ces matrices de corrélation cyclique ou cyclique conjugué ne dépendent asymptotiquement que de l'information spatiale relative aux brouilleurs. En effet, l'étude des sous espaces signal et bruits générés par la décomposition de la matrice de corrélation cyclique, a montré que le sous espace signal couvre le même sous-espace que celui défini par les signatures spatiales des interférences, même en présence de sources cosmiques, et dans le cas où le bruit système n'est pas calibré.

Le problème de détection a été abordé, des améliorations par rapport aux techniques actuelles dites classiques (n'utilisant pas la cyclostationnarité) ont été proposées, et comparées dans leurs performances. Outre la capacité à identifier des observations polluées ou non, les techniques de détection proposées peuvent déterminer le nombre de ces brouilleurs, paramètre important dans le traitement spatial des brouilleurs. De plus, lorsque la fréquence cyclique n'est pas déterminée, nous pouvons la retrouver en aveugle. La détection peut également être une étape précédant les deux autres méthodes décrites dans ce document, et qui nécessitent la connaissance préalable de la fréquence cyclique : la méthode d'estimation et de soustraction et le filtrage spatial.

Dans la méthode d'estimation et de soustraction, le spectre du brouilleur est estimé et soustrait directement des observations. Cette technique a l'avantage, en interférométrie, de ne pas modifier les images du ciel. Néanmoins, le problème de cette méthode est qu'elle n'est applicable que sur certains types de brouilleurs (AM -BPSK voire QPSK avec quelques informations supplémentaires), où la matrice de corrélation des brouilleurs peut être estimée à partir de la matrice de corrélation cyclique des observations. Pour les autres types de modulations, une étude de la cyclostationnarité à des ordres supérieurs serait nécessaire.

Pour les autres types de brouilleurs, notamment les brouilleurs large bande, la méthode de filtrage spatial cyclique a été proposée. C'est une technique générale de suppression des brouilleurs, sans contraintes d'application.

De façon générale, les performances des méthodes proposées sont meilleures que celles des méthodes classiques, sous les hypothèses où ces dernières ne sont pas applicables : la noncalibration du système (bruit système non blanc) et la présence de sources cosmiques fortes. Les futurs instruments seront plus sensibles, et de ce fait, la contribution de ces sources ne sera plus négligeable. De plus, la calibration des réseaux phasés étant problématique, l'application des méthodes cycliques permettrait de s'affranchir de ces contraintes.

Cette thèse a montré l'intérêt d'une approche cyclostationnaire dans le cas d'un radiotélescope à réseau phasé. Cependant, un certain nombre d'approfondissements seraient nécessaires. Ainsi, une étude approfondie de l'impact de ces méthodes sur les images du ciel serait nécessaire, notamment en lien avec les problèmes de calibration. D'autre part, des tests exhaustifs sur les données réelles seraient également très intéressants pour définir plus précisément le domaine d'applicabilité des techniques proposées. En effet, lors de cette étude, nous n'avions pas d'éléments suffisants pour étudier ces images en pratique, et nous nous sommes donc basés sur des simulations. De futures recherches sur les effets du filtrage spatial sont en cours, et des conclusions plus précises pourront alors être données.

Le radiotélescope LOFAR rentrant en phase opérationnelle, il sera un très bon précurseur de ce que sera le radiotélescope SKA. En effet, de par sa sensibilité, sa surface effective et son architecture, LOFAR doit permettre de tester et de prolonger concrètement les approches proposées dans ce mémoire. Une thèse en continuité de ce travail est d'ailleurs d'ores et déjà prévue.

• Conclusion

The present study has been achieved in the framework of the SKADS project [START_REF]Square kilometer array desing studies[END_REF], in order to define some RFI mitigation tools. The originality of our work is based on the use of the cyclostationarity property of telecommunications signals when applying spatial interferences mitigation processing. This study is based on the cyclic and cyclic conjugated correlation matrices, which correspond to a periodicity measurement in the instantaneous correlation matrices, resulting from successive correlation of the array. In this document, several approaches have been proposed. However, they all depend on the fact that the cyclic and cyclic conjugated correlation matrices contain asymptotically only spatial information related to the interferers. From the subspace decomposition of the cyclic correlation matrices, we obtained an interesting result: the estimated signal subspace spans the same subspace as that defined by the RFI spatial signature vectors, regardless to the presence of cosmic sources and/or uncalibrated system noise.

The detection issue has been studied, and some improvements on the existing "classic" approaches (that do not employ cyclostationarity) have been proposed, and their performance compared. In addition to the ability to identify the polluted or unpolluted observations, the detection techniques proposed can determine the number of interferers, which is an important parameter in RFI mitigation. Moreover, when the cyclic frequency cannot be determined, it can be retrieved blindly. Detection can also be a previous step to the other two techniques described in this document, and which require prior knowledge of the cyclic frequency: the estimation and subtraction technique and spatial filtering.

In the estimation and subtraction approach, the interferer spectrum can be estimated and subtracted directly from the observations. The main advantage of such a technique in Sky imaging is that it does not alter the skymaps. However, the main drawback of this approach is that it can be applied only on certain kinds of interferers (AM, BPSK and QPSK with some additional information), for which the RFI correlation matrix can be estimated from the cyclic correlation matrix of the observations. For other modulation schemes, a high order cyclostationarity study is required.

For other kinds of interferers, and more specifically for broad-band interferers, the spatial filtering technique has been proposed. It is a general RFI mitigation method, without application constraints.

Overall, the proposed approaches perform better than the classic approaches, when working under the hypotheses where the classic approach is not applicable, namely an uncalibrated system (the system noise is not white), and in the presence of strong cosmic sources. As future instruments will be more sensitive, it will no longer be possible to disregard the contribution of these sources. Moreover, the problematic system calibration may be avoided when applying the cyclic techniques. This thesis has shown the interest of cyclostationarity for phased array radio telescopes. However, further investigation is still required on the impact on skymaps, and especially the impact on calibration issues. Moreover, exhaustive tests on real data would be interesting so as to define more precisely the applicability context of the proposed techniques. In fact, in the present study, insufficient elements were available to evaluate these skymaps in practice, and the techniques were therefore tested only on simulated data. Further research on spatial filtering will be done, and more precise conclusions may therefore be provided.

Since the LOFAR radio telescope is in the roll-out phase, it will be a good demonstrator for the square kilometer array. Thanks to its sensitivity, effective collecting area and architecture, the LOFAR will allow the approach proposed in this thesis to be tested and improved. Another PhD topic is already underway with this aim in view. Indeed, one of the major improvements that SKA will bring to radio astronomy is that the field of view, i.e. the region of the sky that can be imaged in a single observation, will be much larger than what is currently possible. This increase will be achieved by placing a receiver array in the focus of each parabolic dish of the WSRT, instead of the single receiver element that the current system employs.

Symbols and operators

Spatial processing of cyclostationary interferers for phased array radio telescopes

This thesis is a contribution to observation improvements for phased array radio telescopes, in the presence of radio frequency interferers (RFIs). The originality of the study is the use of the cyclostationarity property, in order to improve the spatial separation between cosmic sources and telecommunication signals. This thesis is part of the European SKADS project, which aims to improve RFI mitigation techniques for future instruments in radio astronomy.

We have proposed three spatial processing techniques: detection, estimation and subtraction and spatial filtering. The performance of the techniques presented have been evaluated through simulations on synthetic and/or real data, and compared to existing approaches.

3. 3 . 3

 33 Self-coherence restoral technique . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Cyclic adaptive beamforming algorithm . . . . . . . . . . . . . . . . . . 3.3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 High resolution methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Linear methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Subspace decomposition methods . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Random binary signal. (a) Temporal representation.(b) Autocorrelation representation when the temporal origin is seen as a random parameter: the process is stationary. (c) Autocorrelation representation when the temporal origin is constant: the process becomes cyclostationary.

Figure 1 . 2 : 20 1. 2 .

 12202 Figure 1.2: (a)Cyclic conjugated autocorrelation function of an amplitude modulated signal (b) Its Fourier transform (F T ): named the cyclic spectrum. Both representations are non-zero only for the cyclic frequency α = 2f 0 . The modulation term a(t) is a rectangular filtered white noise.

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Spectrum of a BPSK modulated signal cyclic autocorrelation function, considering a rectangular window h(t). The spectrum is non-zero for cyclic frequencies α = lTsym , l ∈ ℤ.

Figure 1 .

 1 Figure 1.5 represents the spectrum of the cyclic autocorrelation function of a quaternary phase shift keying (QPSK) modulated signal. The cyclic conjugated autocorrelation function of mPSK and QAM is zero for all values of α, since, as shown by equation (1.34):

Figure 1 . 6 :Figure 1 . 7 :

 1617 Figure 1.6: Spectrum of received signal by one WSRT antenna. The spectrum shows, in addition to the system noise, the presence of an RFI which is a GPS signal. GPS signals are BPSK modulated, and characterized with a carrier frequency f 0 and a symbol duration T sym . WSRT data by courtesy of Albert Jan Boonstra and ASTRON.

Figure 2 . 1 :

 21 Figure 2.1: A sensor array radio telescope configuration. The sensor receives a signal u(t -τ θu,i ) from the direction θ u . This signal is affected by both the sensor beam response and the electronic chain response ( γ i,θu (t)). The received signal plus the additive system noise constitute the sensor output. The array output vector z(t) is constituted by all the sensors' outputs z i=1,...,M (t).

Figure 2 . 2 :

 22 Figure 2.2: A phased array radio telescope configuration in the far-field hypothesis. M sensors receive the signal u(t) from l u . This signal can be an interference (r(t)) or a cosmic source (s(t)). The position of each sensor is given by p i . n i is the additive noise. The array output vector z(t) is constituted by all the sensors' outputs z i=1,...,M (t).

  the spatial signature y coordinate of the spatial signature

Figure 2 .

 2 Figure 2.4: Skymaps with real and synthetic data. The upper figures represent real observations using the LOFAR initial test stations (ITS) (see appendix A.1, by courtesy of Albert-Jan Boonstra and ASTRON). The lower figures represent the corresponding simulations that we have obtained using the proposed model. The skymaps on the left represent an observation with a strong received RFI (cosmic sources are not visible in the map). Figures on the right represent the same observation without RFI. 3 sources appear in the map. For the lower figures, the correlation matrices were obtained firstly withK 1 = 1 and K 2 = 3, for the observation with RFI (left figure), then with K 1 = 0 and K 2 = 3 for the observation without RFI (right figure). These skymaps are based on the ITS antennas positions (i.e for each source we have defined a structured spatial signature).

  Figure 2.4: Skymaps with real and synthetic data. The upper figures represent real observations using the LOFAR initial test stations (ITS) (see appendix A.1, by courtesy of Albert-Jan Boonstra and ASTRON). The lower figures represent the corresponding simulations that we have obtained using the proposed model. The skymaps on the left represent an observation with a strong received RFI (cosmic sources are not visible in the map). Figures on the right represent the same observation without RFI. 3 sources appear in the map. For the lower figures, the correlation matrices were obtained firstly withK 1 = 1 and K 2 = 3, for the observation with RFI (left figure), then with K 1 = 0 and K 2 = 3 for the observation without RFI (right figure). These skymaps are based on the ITS antennas positions (i.e for each source we have defined a structured spatial signature).

Figure 2 . 5 :

 25 Figure 2.5: Comparison between the estimated subspaces resulting from the EVD of R and the SVD of R α . For our general model, the signal subspace in the classic case will depend on the interferers and the cosmic sources. In the cyclic case, it will depend only on the RFIs.

3. 4 .

 4 HIGH RESOLUTION METHODS P (θ) = a H (θ) a(θ).

4. 1 .

 1 FISHER CRITERION DEFINITIONS

CHAPTER 4 .

 4 CYCLIC DETECTION FOR PHASED ARRAY RADIO TELESCOPES

Figure 4 . 1 :

 41 Figure 4.1: The Fisher criterion principle. Left hand figure: The T (z) criterion is not able to completely separate the two hypotheses H 0 and H 1 . Miss detections and false alarms will occur. The corresponding Fisher criterion is low (here, F = 4.5). Right hand figure: The T (z) criterion completely separates the two hypotheses H 0 and H 1 . The detection is perfect. The corresponding Fisher criterion is higher (here, F = 32).

  Consider a M ×M matrix A = [A kl ] k,l=1...M .The Frobenius norm is defined by:

Figure 4 .

 4 Figure 4.2: The subspace decomposition detection algorithm: After a subspace decomposition of the cyclic correlation matrix R α , we compare the eigenvalues to the threshold in order to obtain the number of K 1 dominant eigenvalues corresponding to the number of the K 1 received cyclostationary RFIs with cyclic frequency α.

  frobenius norm of the cyclic conjugated correlation matrix

Figure 4 . 3 :

 43 Figure 4.3: The norm-based detector. Left hand figure: retrieval of the cyclic frequencies: Spectral lines appear for frequencies that are multiples of 1.023 MHz, which correspond to the GPS cyclic frequencies (multiples of the BPSK symbol duration T sym ). Right hand figure: retrieval of the cyclic conjugated cyclic frequencies: Spectral lines appear for corresponding BPSK cyclic conjugated frequencies, which correspond to twice the carrier frequency 2f 0 and additional multiples of the symbol rate 1/T sym . WSRT data by courtesy of Albert-Jan Boonstra and ASTRON.

Figure 4 . 4

 44 Figure 4.4: Fisher criterion comparison (classic vs cyclic approach) for subspace based detectors in the absence of cosmic sources (ρ = 1). The numbers on the right are the sample number, L, used to estimate the classic and cyclic correlation matrices. For example, both detectors perform well down to INR=-15 dB for L = 1024 samples.

Figure 4 . 5

 45 Figure 4.5: Fisher criterion comparison (classic vs cyclic approach) for the Frobenius norm based detectors, in the absence of cosmic source (ρ = 1). The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

Figure 4 . 6

 46 Figure 4.6: Fisher criterion comparison between the Frobenius norm based detector and the subspace based detector, in the classic approach, in the absence of cosmic sources (ρ = 1). The system noise is i.i.d. The correlation matrices have been estimated over L = 1024, 8192 and 65536 samples.

  Figure 4.6: Fisher criterion comparison between the Frobenius norm based detector and the subspace based detector, in the classic approach, in the absence of cosmic sources (ρ = 1). The system noise is i.i.d. The correlation matrices have been estimated over L = 1024, 8192 and 65536 samples.

Figure 4 . 7

 47 Figure 4.7: Fisher criterion comparison (classic vs cyclic approach) for the subspace based detectors, in the presence of a cosmic source with 10% contribution. The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

Figure 4 . 8

 48 Figure 4.8: Fisher criterion comparison (classic vs cyclic approach) for the Frobenius norm based detectors, in the presence of a cosmic source with 10% contribution (ρ = 0.9). The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices.

Figure 4 . 9 :

 49 Figure 4.9: Fisher criterion comparison between subspace detector and Frobenius norm detector in the classic approach. The cosmic source contribution is 10% (ρ = 0.9). The numbers on the right are the number of samples, L, used to estimate the classic correlation matrices.

Figure 4 . 10 :

 410 Figure 4.10: Classic and cyclic approaches in the presence of a significant cosmic source. A 10% source contribution corresponds to ρ = 0.9 and to a SNR of -9.54 dB. A 50% contribution corresponds to ρ = 0.5 and to a SNR of 0 dB. Left hand figure: subspace detection using the dominant eigenvalue extracted from the subspace decomposition of the classic (R) and cyclic (R α ) correlation matrices. Right hand figure: Frobenius norm detectors.

Figure 4 . 11 :

 411 Figure 4.11: Classic and cyclic approaches in the presence of uncalibrated noise. The correlation matrices were estimated over 1024, 8192 and 65536 samples. Left hand figure: subspace detection using the dominant eigenvalue extracted from the subspace decomposition of R (classic) and R α (cyclic). Right hand figure: Frobenius norm of R and R α detector comparison.

Figure 5 . 1 :

 51 Figure 5.1: The Estimation and Subtraction principle: We estimate the RFI correlation matrix R RF I using the cyclic correlation matrix R α , as R α ∼ = R α RF I . Assuming then we can extract R RF I from R α , we subtract the resulting estimate R RF I from the global correlation matrix R.

Figure 5 . 2 :

 52 Figure 5.2: Initial skymaps simulated from the ITS configuration with M = 60 antennas (see appendix A.1). Left hand figure: the reference skymap simulating an ideal sky with no interferer. Only the cosmic source is visible at the zenith on the skymap. Right hand figure: the cosmic source is completely drowned in the presence of a strong interferer.

Figure 5 .

 5 Figure 5.3: Estimation and subtraction method with 1 BPSK interferer. Left to right: The upper figures represent initial skymaps for respectively a strong (INR= 10 dB), a weak (INR= 0 dB) and a very weak interferer (INR= -10 dB). The correlation matrix has been computed over 8192 time samples. The midline figures represent the observations after applying the E & S method: the RFI has been correctly subtracted and is no longer visible on the skymaps. The lower figures represent the error maps.

Figure 5 . 4 :

 54 Figure 5.4: Skymaps for 3 interferers with different emission characteristics, for INR= 0 dB and 8192 time samples: Left to right: (a) Observation with the presence of 3 RFIs. (b) The observation after applying the E & S method using 3 different cyclic frequencies: the RFIs are almost completely removed, and the shape of the map is fairly well preserved. (c) The error map.

Figure 5 . 5 :

 55 Figure 5.5: Spectrum of the AM signal. Left hand figure: TF plane of the observed 0 -100 MHz band with LOFAR. Right hand figure: from the observed spectrum, we filtered the AM signal using a passband filter of 200 KHz bandwidth.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Spectrum of the AM signal, compared to the estimated spectrum using the RFI estimation process based on cyclostationarity.

  Figure 6.1 clearly shows that the RFI spatial signature is still well estimated at IN R = -15 dB.

Figure 6 . 1 :

 61 Figure 6.1: Simulation results of the inner product between the estimated spatial signature and the true one, for 8192 samples and 8 antennas.

Figure 6 . 2 :

 62 Figure 6.2: Simulation results of the INR before and after applying spatial filtering, for 8192 samples and 8 antennas. The first 8 curves satisfy the equal noise powers assumption. For the last two curves, the noise powers fluctuate within 20% between antennas.

Figure 6 . 3 :

 63 Figure 6.3: Simulation results of the INR before and after applying spatial filtering, for 8192 samples and 1000 antennas.

6. 3 .Figure 6 . 4 :

 364 Figure 6.4: Initial skymaps simulated from the ITS configuration with M = 60 antennas (see appendix A.1). Left hand figure: Simulated single source skymap of ideal sky without interferer, the SNR = 0 dB. Right hand figure: the cosmic source is completely drowned in the presence of a strong interferer. Here a BPSK with f 0 = 03 (normalized frequency) and INR = 10 dB.

6. 4 .Figure 6 . 5 : 93 CHAPTER 6 Figure 6 . 6 :

 46593666 Figure 6.5: Classic and cyclic spatial filtering with 1 BPSK. Upper line: from left to right: skymaps for respectively a strong (INR= 10 dB)a weak (INR= 0 dB) and a very weak interferer (INR= -10 dB). The correlation matrix was computed over 8192 time samples. Midline: the classic spatial filtering results on the skymaps. The performance declines as the INR decreases. The lower figures represent the observation after applying the cyclic spatial filtering method: the RFI has been correctly subtracted and is no longer visible on the skymaps

Figure 6 . 7 :

 67 Figure 6.7: Strong RFI (a pager) received with the LOFAR. Left hand figure: TF plane of the observed 160 -240 MHz band with LOFAR. Right hand figure: cyclic frequencies retrieval using the blind detector based on the Frobenius norm of the instantaneous correlation matrix, described in section 4.2.2.

Figure 6 . 8 :

 68 Figure 6.8: The eigenvalue decomposition of the classic and the cyclic correlation matrices estimated from real data from the LOFAR telescope. A strong transmitter is present in the dataset (see figure 6.7). M = 8 antennas were used and the correlation matrices were estimated over L = 65536 samples.

Figure 6 . 11 :Algorithm 4

 6114 Figure 6.11: The eigenvalue decomposition of the classic correlation matrix and the sample cyclic correlation matrix estimated from LOFAR DAB data. A broadband transmitter is present in the dataset (see figure 6.9). M = 8 antennas were used and the correlation matrices were estimated over N = 8192 delays.

Figure 6 . 12 :

 612 Figure 6.12: Broad-band approach. Left hand figure: TF plane of the observed 200 -300 MHz band with LOFAR. Right hand figure: spectrum of one antenna output after applying cyclic spatial filtering (M = 8, N = 8192).

A. 2 .

 2 Figure A.2: ITS parameters for interferometry: performing a beam scan using the antennas location allows LOFAR to observe the whole sky. By courtesy of Albert-Jan Boonstra and ASTRON.

  Cette thèse est une contribution à l'amélioration des observations pour les radiotélescopes à réseaux phasés en présence d'interférences. L'originalité de cette thèse repose sur l'utilisation de la séparation spatiale entre les sources cosmiques et les brouilleurs issus des télécommunications en se basant sur la cyclostationnarité de ces derniers. Cette thèse s'inscrit dans le cadre du projet européen SKADS pour l'amélioration des techniques de suppression d'interférences en radioastronomie pour les futurs instruments d'observations. Nous avons proposé trois techniques de traitement d'interférences : la détection, l'estimation et la soustraction, et le filtrage spatial. Les performances des techniques proposées ont été évaluées à travers des simulations sur des données synthétiques et/ou réelles, et comparées aux techniques existantes.Mots clés : cyclostationnarité, traitement d'antennes, réseau phasé, traitement d'interférences, radioastronomie.

  Keywords: cyclostationarity, array processing, phased array, RFI mitigation, radio astronomy.Institut PRISME -Equipe ISS, Polytech'Orléans 12 Rue de Blois BP 674445067 Orleans Cedex 2

  Denotes vector, a = [a kl ] k,l=1...M A Denotes matrix, A = [A kl ] k,l=1...M E {.} Expected value (.) * Conjugate < . > ∞ Infinite time averaging < . > ∞ = 1 T lim.dt for nonprobabilistic approach, or < . > ∞ = E {.} for stochastic approach < . > L Finite time averaging for continuous signals < . > L = 1 Frobenius norm of the K × L matrix A, frob(A) =

	Acronyms
	A, a RFI	Denote scalars Radio frequency interferers
	a SOI	Signal of interest
	INR	Interference to noise ratio
	SNR	Signal to noise ratio
	SKA	Square kilometer array
	T /2 SKADS Square kilometer array design and study
	LOFAR	Low frequency array		T →∞	-T /2
	WSRT	Westerbork synthesis radio telescope
	ITS GPS DAB	Initial test station Global positioning system for discrete signals, < . > L = 1 L Digital audio broadcasting	L/2	.	L	L/2 -L/2 .dt,
	-L/2 Digital video broadcasting -terrestrial Convolution operator Amplitude modulation Transpose operator Binary phase shift keying Conjugate transpose operator Qua Matrix inverse Quadrature amplitude modulation Element-wise multiplication Offset quadrature amplitude modulation Element-wise division Element-wise multiplication Complex conjugate operator Frequency shift keying The M × M identity matrix Continuous phase frequency shift keying diag(.) Diagonal operator, converts a vector into a diagonal matrix ⊗ DVBT (.) T AM (.) H BPSK (.) -1 QPSK . * QAM ./ OQAM (.) CPM I M ×M FSK CPFSK tr(.) GMSK Gaussian minimum shift keying Trace operator det(.) OFDM Orthogonal frequency division multiplexing Determinant a √ a H a FFT Fast Fourier transform EVD Eigenvalue decomposition Vector norm, |a| = frob(A) K k=1 L SVD Singular value decomposition DOA Direction of arrival l=1 e . SCORE Self-coherent restoral Exponential log LS Least square e-based logarithm log 10 CAB Cyclic adaptive beamforming 10-based logarithm δ(.) BEWE Bearing estimation without eigen decomposition Dirac delta function. SWEDE Subspace method without eigen decomposition	A 2 kl
	MUSIC	Multiple signal classification
	ESPRIT Estimation of signal parameters via the rotational invariance technique
	E&S	Estimation and subtraction method
	SF	Spatial filtering
	CSF	Cyclic spatial filtering

Du fait du mouvement des électrons, tous les corps émettent un rayonnement thermique caractéristique de leur température. Cependant, il existe également des émissions non thermiques, beaucoup plus intenses, provenant de particules chargées se déplaçant dans les champs magnétiques galactiques et intergalactiques. Lorsque l'énergie d'une particule est si élevée que sa vitesse tend vers celle de la lumière, l'émission radioélectrique de ces particules est appelée radiation synchrotron.

τ θu,i (t) = (t -τ θu,i )e -j2πf 0 τ θu,i .

In practice, only one RFI is present within a given bandwidth. Thus, K 1 = 1 in general, and K

only refers to the number of received cosmic sources. However, for some interferers, such as GPS signals, or in the presence of multipath, K 1 > 1.

The dump time for R will depend on the stationarity conditions of the observation, typically 10s, and the dump time for R k will depend on the stationarity conditions of the RFI, typically 1ms[START_REF] Clerc | High performance receiver for rfi mitigation in radio astronomy : Application at decameter wavelengths[END_REF] [START_REF] Baan | Radio frequency interference mitigation at the westerbork synthesis radio telescope: Algorithms, test observations, and system implementation[END_REF]

Pulsars are rapidly rotating highly magnetized neutron stars which produce radio beams that sweep the sky like a lighthouse. If the beam is oriented towards the Earth, it will produce periodic pulses which can be measured with radio telescopes and dedicated backends.

Roughly, the Fourier transform of R will transform correlation information (i.e. visibilities in interferometry) into sky image intensities[START_REF] Leshem | Multichannel interference mitigation techniques in radio astronomy[END_REF].

In fact, in the presence of cosmic sources or system noise, the same expression will be derived since the cyclic approach removes all the components which are not cyclostationary at cyclic frequency α.

Remerciements

Assuming knowledge of α and f 0 , We can derive from u r , another phase vector : u(f 0 ) = [a(α)] f 0 α = e -j2πf 0 τ θr ,1 , ..., e -j2πf 0 τ θr ,M T We propose to use the broadband approach to retrieve u r and then to construct a projector P cyclic (f 0 ), as defined in section 6. 1. In other words, the observed band will be split into K smaller sub-bands centered at frequency f k , k = 0, ..., K -1. The sub-band bandwidth will fit the narrow band assumption. From the broad-band analysis, we will extract the u r vector. Then, we will derive K cyclic projectors P cyclic (k), k = 0, .., K -1 which will be applied individually on each sub-band (see figure 6.10). Figure 6.10: Broad-band spatial filtering. The observed band will be split into K smaller sub-bands centered at frequency f k , k = 0, ..., K -1. The sub-band bandwidth will fit the narrow band assumption. From the broad-band analysis, K cyclic projectors P cyclic (k) will be derived, and which will be applied individually on each sub-band, or on each correlation matrix R(f k ).

Real data results

We have observed in the 222 -226 MHz LOFAR band, which contains a broad-band interferer (Digital Audio Broadcasting (DAB)) at 223 -225 MHz with an INR of 5 dB. The array configuration consisted of M = 8 LOFAR antennas. Figure 6.11 shows the eigenvalues obtained from the classic correlation matrix R = z(t)z H (t) N , and the sample cyclic correlation matrix R sample = z α (τ )(z α ) H (τ )

N

, with N = 8192. The cyclic frequency, α, of the observed German DAB has been retrieved from the literature [START_REF] Hoeher | Performance of an RCPCcoded OFDM-based digital audio broadcasting (DAB) system[END_REF] Appendix A

Existing phased array radio telescope A.1 The low frequency array

The LOFAR [START_REF] Website | Low frequency array for radio astronomy[END_REF] (figure A.1) is an instrument built by ASTRON, dedicated to interferometry. The core is situated in Exloo (the Netherlands) with stations located across the Netherlands, and extensions to other countries (Germany, France, Great Britain, Sweden, and possibly Poland and Ukraine). In LOFAR, antennas are grouped in so-called stations in which the signals from over a hundred antennas are combined using phased-array beamforming. The beamformed signals of many stations are combined centrally by correlating them. It is currently in the roll-out phase with a compact core area comprising 32 stations and 45 remote stations, and operates in the band 30-240 MHz. The LOFAR radio telescope can represent a demonstrator for SKA, since the predicted total effective collecting area will be up to 1 square kilometer. Additionally, LOFAR is currently the most sensitive instrument in radio astronomy.

• Initial test station

The Initial Test Station (ITS) is the first configuration of the LOFAR. It has been operational since December 2003. The ITS is composed of 60 inverse V-shaped dipoles. The position of the dipoles allows the whole sky to be mapped, as shown by figure A.2.

A.2 The Westerbork synthesis radio telescope

The WSRT (figure A.3) [START_REF] Website | The netherlands institute for radio astronomy[END_REF], located in Westerbork, Midden-Drenthe, the Netherlands, is a synthesis radio telescope built by ASTRON. It is an aperture synthesis interferometer that consists of a linear array of 14 antennas arranged on a 2.7 km East-West line. Ten of the telescopes are on fixed mountings while the remaining four dishes are movable along two rail tracks. The telescope has been operational since 1970.

The telescopes in the array can operate at several frequencies between 120 MHz and 8.3 GHz. The WSRT is often combined with other telescopes around the world to perform Very Long Baseline Interferometry (VLBI) observations, being part of the European VLBI Network.

The WSRT is a demonstrator of the square kilometer array (SKA) within the project APERTIF. This project aims to increase the field of view of the WSRT with a factor 25.

Figure A.1: LOFAR station at the Effelsberg site (Germany). Bottom left : the low frequency array (30 MHz-80MHz). Top the high frequency array (110 -240 MHz). Each black box is itself an array of high frequency antennas. Copyright of the Max Planck Institute for radio astronomy, ASTRON and LOFAR.

Figure A.3: The Westerbork synthesis radio telescope (WSRT) -the Netherlands. The grey building in the bottom right corner of the picture is the EMBRACE (Electronic MultiBeam Radio Astronomy ConcEpt) dome. EMBRACE is an aperture-plane phased-array prototype using tiles of Vivaldi feed antennas. It is part of the FP6 European Project SKADS [93] and it has installations at Westerbork and Nancay.

List of Figures

Spectrogramme de la bande 20-300 MHz de l'environnement du radiotélescope LOFAR (annexe A.1). Certaines interférences ont disparu du spectre (quelques chaînes TV). Néanmoins, les bandes vides actuelles (170 -230 MHz) seront progressivement occupées par des transmissions issues des radio et TV numériques (DAB et DVB).

Avec l'autorisation de Albert-Jan BOONSTRA et ASTRON. . . . .

Random binary signal. (a) Temporal representation. (b) Autocorrelation representation when the temporal origin is seen as a random parameter: the process is stationary. (c) Autocorrelation representation when the temporal origin is constant

Spectrum of a BPSK modulated signal cyclic autocorrelation function, considering a rectangular window h(t). The spectrum is nonzero for cyclic frequencies α = l

Tsym , l ∈ ℤ. . . . . . . . . . . . . . . . . .

LIST OF FIGURES

Spectrum of the cyclic conjugated autocorrelation function for a BPSK modulated signal, with h(t) a rectangular window. The spectrum is non-zero for cyclic frequencies

Skymaps with real and synthetic data. The upper figures represent real observations using the LOFAR initial test stations (ITS) (see appendix A.1, by courtesy of Albert-Jan Boonstra and ASTRON). The lower figures represent the corresponding simulations that we have obtained using the proposed model. The skymaps on the left represent an observation with a strong received RFI (cosmic sources are not visible in the map). Figures on the right

The norm-based detector. Left hand figure: retrieval of the cyclic frequencies: Spectral lines appear for frequencies that are multiples of 1.023 MHz, which correspond to the GPS cyclic frequencies (multiples of the BPSK symbol duration T sym ). Right hand figure: retrieval of the cyclic conjugated cyclic frequencies: Spectral lines appear for corresponding BPSK cyclic conjugated frequencies, which correspond to twice the carrier frequency 2f 0 and additional multiples of the symbol rate 1/T sym . WSRT data by courtesy of Albert-Jan

Fisher criterion comparison (classic vs cyclic approach) for subspace based detectors in the absence of cosmic sources (ρ = 1). The numbers on the right are the sample number, L, used to estimate the classic and cyclic correlation matrices. For example, both detectors perform well down to INR=

Fisher criterion comparison (classic vs cyclic approach) for the Frobenius norm based detectors, in the absence of cosmic source (ρ = 1)

. The numbers on the right are the number of samples, L, used to estimate the classic and cyclic correlation matrices. . . . . .

LIST OF FIGURES

Fisher criterion comparison between the Frobenius norm based detector and the subspace based detector, in the classic approach, in the absence of cosmic sources (ρ = 1). The system noise is i.i.d. The correlation matrices have been estimated over

Estimation and subtraction method with 1 BPSK interferer. Left to right: The upper figures represent initial skymaps for respectively a strong (INR= 10 dB), a weak (INR= 0 dB) and a very weak interferer (INR= -10 dB). The correlation matrix has been computed over 8192 time samples. The midline figures represent the observations after applying the E & S method

LIST OF FIGURES

A.3 The Westerbork synthesis radio telescope (WSRT) -the Netherlands. The grey building in the bottom right corner of the picture is the EMBRACE (Electronic MultiBeam Radio Astronomy ConcEpt) dome. EMBRACE is an aperture-plane phased-array prototype using tiles of Vivaldi feed antennas. It is part of the FP6 European

Project SKADS [START_REF]Square kilometer array desing studies[END_REF]