A&A 368, 369-376 (2001)
DOI: 10.1051/0004-6361:20000552

RFI excision using a higher order statistics analysis of the power spectrum

P. A. Fridman

ASTRON, Postbus 2, 7990AA Dwingeloo, The Netherlands

Received 7 September 2000 / Accepted 8 December 2000

Abstract
A method of radio frequency interference (RFI) suppression in radio astronomy spectral observations based on the analysis of the probability distribution of an instantaneous spectrum is described. This method allows the separation of the Gaussian component due to the natural radio source and the non-Gaussian RFI signal. Examples are presented in the form of computer simulations of this method of RFI suppression and of WSRT observations with this method applied. The application of real time digital signal processing for RFI suppression is found to be effective for radio astronomy telescopes operating in a worsening spectral environment.

Key words: power spectrum - probability distribution - RFI excision


1 Introduction

The increasing radio frequency interference (RFI) from transmissions and spurious signals of commercial telecommunication applications limit the real sensitivity of radio telescopes at certain frequencies. This radio-ecological environment at radio observatories appears to be worsening each year, while the observing systems get more sensitive. Further improvements in antenna quality and in receiver parameters do not give the expected results because of limitations due to RFI. Therefore it is very difficult to reach the design sensitivity, especially at low frequencies, even for radio telescopes situated in secluded areas (Spoelstra 1997; Kahlmann 1999; Cohen 1999).

Several methods of RFI excision have been proposed during recent years. The successful application of these methods depends on the type of radio telescope, type of observations, type of RFI, and there is no universal RFI mitigation algorithm.

1. Type of radio telescope - single dish, connected radio interferometer, VLBI. Single dish radio telescopes are especially vulnerable to RFI. Radio interferometric systems suffer from RFI to a lesser degree, (Thompson 1982), but auxiliary noise power measurements made at each antenna for calibration purposes are distorted by RFI;

2. Mode of observations - continuum or spectral. In continuum observations it is possible to sacrifice some part (5-10%) of the data stream polluted by a strong RFI and save the remained part without significant loss of sensitivity. In spectral observations RFI and signal-of-interest may occupy the same spectral region, so it is impossible to delete this particular part of the spectrum;

3. Type of RFI - impulse-like bursts, narrow-band or wide-band. Theoretical consideration and experimental data show that RFI may be basically represented as the superposition of two types of waveforms: impulse-like bursts and long narrowband and strongly correlated random oscillations (Middleton 1972,1977; Lemmon 1997). RFI of the first type can be effectively suppressed using fast real time thresholding before averaging at the total power detector or correlator (Fridman 1996; Weber et al. 1997). This is the processing in temporal domain.

RFI of the second type should be processed in the frequency domain, where there is a good contrast between a quasi-continuum spectrum of the system noise plus radio source signal and the narrow-band RFI signal. RFI in continuum observations may be excised using adaptive filtering methods (Fridman 1998), such that the parts of the spectrum affected by RFI are used with a smaller weighting than the other "unpolluted'' spectral regions during the averaging in the time-frequency plane. This sort of processing cannot be applied in spectral line observations, where RFI may contaminate the same spectral region that is of interest for the radio astronomer. In the case of the GLONASS navigational satellites (Ellington et al. 2000) the RFI was parametrically modelled and the model of RFI signal was subtracted from the radio telescope output so that the 1612 MHz OH spectral line was not distorted.

Adaptive interference cancellation (similar to adaptive noise cancellation, Widrow & Stearns 1985) can be useful for single dish observations with an auxuliary reference antenna, pointed "off-source''. The use of an auxiliary reference channel provides some benefits (Barnbaum & Bradley 1998), but not all radio telescopes are equipped or can be equipped with such a reference channel. Furthermore, the effectiveness of this adaptive noise cancelling strongly depends on the sensitivity of the reference channel. The reference channel should at least be as sensitive as the main channel for successful RFI suppression (Fridman 1997,2000).

A new method of post-correlation RFI removing was proposed (Briggs et al. 2000). Cross-correlation products of four signals (two main channels and two reference channels) are processed, using closure complex amplitude relations. RFI signal is estimated and subtracted from the outputs of main channels. So RFI suppression can be made off-line with the averaged data, see also references in (The Elizabeth and Frederick White Conference on Radio Frequency Interference Mitigation Strategies 1999).

Methods of spatial filtering in the temporal and frequency domains applicable for radio interferometers are discussed elsewhere (Leshem & van der Veen 1999). Spatial filtering methods use the difference in the direction-of-arrival of radio source and RFI. The RFI emission of well spatially localized sources could be suppressed in multielement radio interferometers using adaptive array philosophy, when zeros of a synthesized antenna pattern coincide with directions-of-arrival of undesirable signals (adaptive nulling techniques). This approach is popular while considering new projects of superlarge phased array radiotelescopes, (van Ardenne 1996; Bregman 2000). However, there are some limitations in using this techniques for the existing large radio interferometers like the WSRT, VLA, and GMRT: very sparse antenna arrays, phase-only computer control, complex impact on image postprocessing. Spatial filtering is effective when RFIs are strongly correlated at radio telescope antennas (like in a half-wavelength phased array). RFI at the adjacent antennas (hundreds meters separation) are highly correlated, but when the distance between antennas is more than 1 km (and especially for baselines 10-100 km) RFI at these points are less correlated, and effectiveness of the adaptive nulling is less, than for the continuous-like half-wavelength phased arrays. RFI at the sites of VLBI system, separated by hundreds and thousands of kilometers, are practically uncorrelated. So the impact of an RFI is an increased variance at correlator output, which is noticeable in the case of strong RFI (comparable with the system noise power).

In this paper, I consider a method of RFI excision based on the analysis of the probability distribution of the power spectrum at the radio telescope output. System and radio source noise generally has a Gaussian distribution with a zero mean. If a Fourier transformation is applied to such an "ideal'' signal, the real and imaginary components in every spectral bin are also Gaussian random values with zero mean. The instantaneous power spectrum, which is a square of magnitude of the complex spectrum, has an exponential distribution (or, in other terms, a chi-squared distribution with two degrees of freedom). Any changes of the input signal of this model, due to the presence of an RFI, yield a change of its statistics. The power spectrum in this case has a non-central chi-squared distribution with two degrees of freedom. Radio astronomy observations generally employ the averaging of the power spectrum over some integration period, which results in a conversion of the probability distribution of the averaged power spectrum into a Gaussian distribution.

Analysis of sample probability distribution before averaging allows an effective separation of these signals by means of real-time processing facilities employing digital signal processing. The method may be effective for single dish observations without a reference channel. The description of this separation by using the analysis of the power spectrum distribution will be presented below. I will also present some computer simulations, and some results from test observations at the Westerbork Synthesis Radio Telescope (WSRT).

2 Power spectrum probability distributions and moments

A radio telescope receiver output signal is described as follows:

\begin{displaymath}r(t)=x_{\rm sig}(t)+x_{\rm sys}(t)+x_{\rm RFI}(t),
\end{displaymath} (1)

where $x_{\rm sys}(t)$ is a system noise signal (receiver, feed, antenna), $x_{\rm sig}(t)$ is the radio source noise signal with some spectral features to be estimated, and $x_{\rm RFI}(t)$ is the RFI waveform. $x_{\rm sys}(t)$ and $x_{\rm sig}(t)$ are each supposed to be random processes with Gaussian probability distributions, zero mean and variances $\sigma _{\rm sys}^{2}$ and $\sigma _{\rm sig}^{2},$ respectively. However, $x_{\rm RFI}(t)$ is a narrow-band interference signal with variance $\sigma _{\rm RFI}^{2}$. After Fourier transformation on a time interval T, we get the instantaneous complex spectrum from Eq. (1):
$\displaystyle S(f)=\int_{0}^{T}r(t)\cos 2\pi ft{\rm d}t+ j\int_{0}^{T}r(t)\sin 2\pi ft{\rm d}t=$      
$\displaystyle =S\_re(f)+jS\_im(f).$     (2)

Probability distributions of the real part of the spectrum $S\_re(f)$ and the $\ $imaginary part $S\_im(f)$ are also Gaussian with a mean A(f) and a variance $\sigma _{n}^{2}(f)+{A}(f)_{n}^{2},$ where component A(f) is determined by the RFI and the $\sigma _{n}^{2}(f)$ is determined by system + signal spectral densities. A power spectrum (periodogram) is written as:

\begin{displaymath}P(f)_{T}=S(f)\overline{S(f)},
\end{displaymath} (3)

which corresponds to an envelope quadratic detector operation, the overbar denotes complex conjugation. Here P(f)T at each frequency bin f has the probability distribution of the square of the envelope of a randomly phased sine wave and a narrow-band Gaussian process (noncentral $\chi ^{2}$ distribution with two degrees of freedom):

\begin{displaymath}w(P)=\frac{1}{2\sigma _{n}^{2}}{\rm e}^{-\frac{P+{A}_{\rm RFI...
...eft(\frac{{A}_{\rm RFI}\sqrt{P}}{\sigma _{n}^{2}}\right),P>0,
\end{displaymath} (4)

which in the absence of RFI is an exponential distribution:

\begin{displaymath}w(P)=\frac{1}{2\sigma _{n}^{2}}{\rm e}^{-\frac{P}{2\sigma _{n}^{2}}}.
\end{displaymath} (5)

After averaging M periodograms (3), $P_{\rm aver}(f) = \sum_{m=1}^{M}P_{m}(f)$, we have for $P_{\rm aver}$ the non-central $\chi ^{2}$ distribution with 2M degrees of freedom (Whalen 1971):
$\displaystyle w(P_{\rm aver})$ = $\displaystyle \frac{1}{2\sigma _{n}^{2}}\left(\frac{P\sigma _{n}^{2}}{M{A}_{\rm...
...(-\frac{P}{2\sigma _{n}^{2}}-\frac{M{A}_{\rm RFI}^{2}}{2\sigma _{n}^{2}}\right)$  
    $\displaystyle \times I_{M-1}\left(\frac{{A}_{\rm RFI}\sqrt{PM}}{\sigma _{n}^{2}}\right)$ (6)

and in the absence of RFI:

\begin{displaymath}w_{0}(P_{\rm aver})=\frac{1}{(2\sigma _{n}^{2})^{M}\Gamma (M)}P^{M-1}{\rm e}^{-P^{2}/2\sigma _{n}^{2}}.
\end{displaymath} (7)

When the noncentral parameter ${A}_{\rm RFI}/\sigma_{n} >1$, both distributions (4) and (6) tend to acquire a Gaussian form. For big M (central limit theorem) the distribution (6) also is transformed into a Gaussian. This near-Gaussian distribution characterizes the power spectra at the outputs of common spectrum analyzers. Figure 1a shows how dramatically the distribution (4) is transformed, while ${A}_{\rm RFI}/\sigma _{n} $ grows from zero to 5.
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig1.eps}
\end{figure} Figure 1: The non-central $\chi ^{2}$ probability distribution as a function of a) ( ${A}_{\rm RFI}/\sigma _{n} $) and b) the order M
Open with DEXTER

Figure 1b represents this kind of transformation with periodogram averaging applied. Figure 2a shows several sections of the three-dimensional presentation of Fig. 1a.
  \begin{figure}
\par\includegraphics[width=8.7cm,clip]{fig2.eps}
\end{figure} Figure 2: a) The non-central $\chi ^{2}$ probability distribution as a function of a) the value of the noncentral parameter ( ${A}_{\rm RFI}/\sigma $) and b) the skewness and excess as a functions of ${A}_{\rm RFI}/\sigma _{n} $
Open with DEXTER

One observes that the presence of RFI changes the form of the distribution w(P) of the unaveraged periodogram. This property may be used for separating the (system + signal) noise component $\sigma _{n}^{2}(f)$ from the RFI component ${A}_{\rm RFI}(f)$.

The most promising method is to obtain an estimate of the sample probability distribution $\widehat{w}[P(f)]$ (histogram), and subsequently estimate the $\widehat{\sigma}_{n}^{2}(f)$ and $\widehat{A}_{\rm RFI}(f)$, which best fit this sample distribution (symbol $\widehat{~}$ denotes sample average). However, performing this processing in real-time with a sufficiently wide band (at the MegaHertz level) is a difficult task even for the most powerful modern Digital Signal Processors (DSPs).

A more attractive practical approach is to use the sample moments of the power spectrum. The first four statistical parameters, mean, variance, skewness and excess, for distribution (4) are:

mean = $\displaystyle 2\sigma _{n}^{2}+{A}_{\rm RFI}^{2},$ (8)
var = $\displaystyle 4\sigma _{n}^{4}+4\sigma _{n}^{2}{A}_{\rm RFI}^{2},$ (9)
skew = $\displaystyle \frac{2\sigma _{n}^{2}+3{A}_{\rm RFI}^{2}}{[\sigma _{n}^{2}+{A}_{\rm RFI}^{2}]^{3/2}}\sqrt{\sigma _{n}^{2}},$ (10)
exc = $\displaystyle 6(\sigma _{n}^{2}+2{A}_{\rm RFI}^{2})\frac{\sigma _{n}^{2}}{(\sigma _{n}^{2}+{A}_{\rm RFI}^{2})^{2}}\cdot$ (11)

In the absence of RFI, the $mean=2\sigma _{n}^{2}$, the variance = $4\sigma_{n}^{4}$, the skewness = 2, and the excess = 6. For strong RFI, the $skewness \rightarrow 0$ and the excess $\rightarrow 0$, which is a good indicator of the presence of RFI. Figure 2b shows this dependence of skewness and excess on the parameter ${A}_{\rm RFI}/\sigma _{n} $.

Having the sample values of $\widehat{mean},\widehat{var},\widehat{skew},\widehat{exc}$, this redundant system of nonlinear equations can be solved with respect to $\sigma _{n}^{2}$ and ${A}_{\rm RFI}^{2}$, with the obvios constraints $\sigma _{n}^{2}>0, {A}_{\rm RFI}^{2}>0$. Analysis of the "solution surface'' shows that the solution is unique. Even on the basis of only the first two equations, it is easy to derive that at each frequency

$\displaystyle \widehat{{A}_{\rm RFI}^{2}}=\sqrt{\widehat{mean}^{2}-\widehat{var}},$     (12)
$\displaystyle \widehat{\sigma }_{n}^{2}=(\widehat{mean}-\widehat{{A}_{\rm RFI}^{2}})/2.$     (13)

The value $2\widehat{\sigma }_{n}^{2}(f)$ as a function of frequency provides "clean'' power spectrum.

Sometimes radio astronomers use a cross-correlation spectrum instead of the auto-correlation spectrum, with the aim of suppressing large system noise spectral component. The probability distribution and the statistical parameters for this case are given in the Appendix.

3 Computer simulations

Some computer simulations have been done to investigate how effective the separation of the Gaussian noise component and the non-Gaussian RFI can be. Figure 3 represents a block diagram of the simulation program.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig3.eps}
\end{figure} Figure 3: A block diagram of the computer simulations
Open with DEXTER

The three signal components are as follows: (1) the system noise is a sequence of Gaussian random values passed through a passband Butterworth filter with cutoff frequencies [0.1, 0.9] and with order 7, (2) the signal noise imitating a spectral line is also represented as a sequence of Gaussian random values passed through a Butterworth filter with cutoff frequencies [0.3 0.35] and with order 5, and (3) the RFI waveforms consists of a superposition of two strong chirp (frequency modulated sinusoidal wave) signals. One of the RFI signals coincides in frequency with the position of the "spectral line''. These three signals are added and applied to a spectrum analyzer (FFT block), where an instantaneous complex spectrum is calculated. This operation imitates a discrete digital spectral analysis with L frequency channels. The generation of the mixture noise + RFI and of the Fourier transform was repeated M times. The upper part of the block diagram after the FFT block corresponds to the proposed RFI excision method, where statistical parameters (mean, variance, skewness, and excess) are calculated at each frequency bin for the M samples. Subsequently the estimates $\widehat{\sigma }_{n}^{2}(f_{l})$, for each frequency channel $l=1...\ L$ were calculated.

The lower part of the block diagram in Fig. 3 shows the common way of doing spectral analysis: calculating of the power spectrum $P(f)_{T}=S(f)\overline{S(f)}$ and subsequent averaging. At this point the two versions of the power spectrum $\widehat{\sigma }_{n}^{2}(f_{l})$ are available: one without RFI excision, and another with RFI excision. Figures 4-7 show the results of the some computer simulations.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig4.eps}
\end{figure} Figure 4: Computer simulations: a) the power spectrum with spectral line with no RFI; b) the input signal with no RFI; c) a two feature RFI power spectrum (logarithmic scale); d) the resulting input signal with RFI; e) the "cleaned'' power spectrum; g) the difference between a) and e)
Open with DEXTER

Figure 4a is an example of an input power spectrum without RFI with its input signal waveform in Fig. 4b. Figure 4c shows the "dirty'' power spectrum without any RFI excision on a logarithmic scale with the level of RFI about 30 dB above the system noise. Figure 4d is the waveform of the input signal with a strong RFI source on the same scale as the system noise. Finally, Fig. 4e shows the "clean'' power spectrum, where the "spectral line'' is clearly visible and the RFI is significantly suppressed, and Fig. 4g shows the "residuals'' or signal-of-interest distortions, that is the difference between "noRFI'' power spectrum a) and "clean'' power spectrum e). The rms level of these variations is approximately at the level of system noise ${\rm rms} \approx{\sigma^{2}_{\rm sys}}/\sqrt{M}$ in the absence of RFI, M=2.5 105.
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig5.eps}
\end{figure} Figure 5: Spectra of the statistical parameters for the model in Fig. 4c (with RFI): a) mean; b) variance; c) skewness; d) excess
Open with DEXTER

Figure 5 gives the spectra of four statistical parameters: mean and variance with a logarithmic scale, and skewness and excess with a linear scale. Figures 6 and 7 represent similar computer simulation, but with the RFI not coinciding with the "spectral line''.
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig6.eps}
\end{figure} Figure 6: Computer simulations: a) the power spectrum with spectral line with no RFI; b) the input signal with no RFI; c) a two feature RFI power spectrum (logarithmic scale); d) the resulting input signal with RFI; e) the "cleaned'' power spectrum; g) the difference between a) and e)
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig7.eps}
\end{figure} Figure 7: Spectra of the statistical parameters for the model in Fig. 6c (with RFI): a) mean; b) variance; c) skewness; d) excess
Open with DEXTER

Figures 7c and d show the remarkable characteristics of the skewness and the excess: they are "transparent'' for the Gaussian "spectral line'' component and react only to RFI presence, that is in the absence of RFI, the ${\rm skewness} = 2$, and the excess = 6, and for strong RFI both skewness and excess $\to0$. This is in harmony with (10), (11) and Figs. 1a, 2b, due to "Gaussianization'' of probability distribution function: skewness and excess are equal to zero for Gaussian distribution.

4 Test observations

Test observations were made with the WSRT to demonstrate the effectiveness of the proposed method with real radio astronomy signals. Digital spectral analysis was performed with the help of a SIGNATEC PMP-8 board with eight high performance Texas Instruments DSP TMS320C6201 processors (Fridman 2000). The baseband (0-1.25 MHz) signal at the analogue output of WSRT frequency conversion subsystem corresponding to one of the radio telescopes (RT6) is digitized in a 12-bit analogue-digital converter and then applied to the PMP-8, where the FFT and the averaged statistical parameters were calculated. The final presentation of the "dirty'' and the "clean'' spectra, and of their difference were done on a workstation.

Figures 89, and 10 show the results of three test observations of Galactic spectral lines: 21 cm HI lines in the direction of Crab and Orion, and an 18 cm OH line in the W3(OH) radio source.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig8_2.eps}
\end{figure} Figure 8: RFI excision using the WSRT with bandwidth = 1.25 MHz, halfsampling frequency = 1.563 MHz, and frequency resolution = 3.052 kHz: a) a neutral hydrogen spectral line (21 cm) in the Galaxy towards the Crab nebula with RFI superposed at 1420.397 MHz; b) the same spectrum after RFI excision simultaneously made with frame a; c) a magnification of the spectral line region (channels 30-100) with spectra of a) and b) superposed on a record without any RFI
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig9_1.eps}
\end{figure} Figure 9: RFI excision using the WSRT with bandwidth = 1.25 MHz, halfsampling frequency = 1.563 MHz, and frequency resolution = 3.052 kHz: a) a neutral hydrogen spectral line (21 cm) in the Galaxy towards the Orion nebula with RFI superposed at 1420.350 MHz; b) the same spectrum after RFI excision simultaneously made with frame a; c) a magnification of the spectral line region (channels 60-120) with spectra of a) and b) superposed on a record without any RFI
Open with DEXTER

An artificial continuous wave signal was generated with a Rohde & Schwarz SMR20 to serve as the RFI signal and was emitted in the direction of RT6 from the WSRT control building. The figures display the data with and without RFI excision applied. These three examples represent relatively strong interfering signals in astronomical data. The RFI suppression obtained in these observational spectra are respectively 17dB in Fig. 8, 13dB in Fig. 9, and 20dB in Fig. 10.

5 Conclusions

1.
Single dish spectral observations are most susceptible to RFI. Digital spectral analysis with real-time analysis of the power spectrum higher order statistics provides a promising way for narrow-band RFI excision taking into account the performance of modern DSP systems. The method can be also applied in continuum observations. There is no bandwidth limitations, in principle, as the algorithm works in frequency domain, after Fourier transformation. Several RFIs at different frequencies are processed in a parallel way, independently of each other, see Fig. 4;
2.
Computer simulations and test observations demonstrate the practical effectiveness of this method with respect to the narrow-band non-Gaussian RFI signals;
3.
The proposed method is a nonlinear procedure and the post-processing suppression gain (the RFI suppression $INR_{\rm input}/INR_{\rm output}$ ratio) depends on the input interference/noise ratio $INR_{\rm input} = {A}_{\rm RFI}^{2}/\sigma_{n}^{2}$ in a complex way. For large $INR_{\rm input}$ this gain can be estimated approximately as $10\log(INR_{\rm input})+5\log(M)$. The stronger RFI, the bigger the gain. The results from our tests with astronomical data show that the gain may range from 13 - 20 dB;
4.
The effectiveness of the method depends also on the temporal stability of INR: Eqs (8)-(11) must be solved under the assumption of constant A during the averaging. This limits the value for M and hence the averaging procedure should be divided on three stages: preliminary averaging - RFI subtraction - final averaging;
5.
The proposed method was tested for the autospectral, periodogram analysis, that is the power spectrum and its high-order statistics were calculated in real time after Fourier transform of the baseband signals. Most of the modern radio astronomy spectrum analyzers use auto- or cross-correlation techniques with 1-2 bit quantization and calculate only the first moment of power spectrum. I think that the direct periodogram spectral analysis with 12-16 bit sampling is more robust to strong RFI - less harmonic distortion products. This is also more suitable for implementation of other RFI mitigation methods using the same DSP hardware, as it is practiced in the RFI mitigation test system at WSRT.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig10_1.eps}
\end{figure} Figure 10: RFI excision using the WSRT with bandwidth = 2.5 MHz, halfsampling frequency = 3.125 MHz, and frequency resolution = 6.104 kHz: a) a hydroxyl (OH) emission line (18 cm) in the Galaxy towards W3(OH) at 1665.8 MHz with RFI superposed; b) the same spectrum after RFI excision simultaneously made with frame a; c) a magnification of the spectral line region (channels 250-300) with spectra of a) and b) superposed on a record without any RFI
Open with DEXTER

Appendix

For a product of two Gaussian values z=xy with variances $\sigma _{1}^{2},\sigma _{2}^{2}$, zero means and correlation coefficient $\rho $ the probability distribution function (PDF) is
$\displaystyle p(z)=\frac{1}{\pi \sigma _{1}\sigma _{2}\sqrt{(1-\rho ^{2})}}K_{0}\left[\frac{\vert z\vert}{\sigma _{1}\sigma _{2}(1-\rho ^{2})}\right]$      
$\displaystyle \times\exp \left[\frac{\rho z}{\sigma _{1}\sigma _{2}(1-\rho ^{2})}\right],$     (14)

and the corresponding characteristic function is

\begin{displaymath}f(\xi )=[1-2i\xi \sqrt{\sigma _{1}\sigma _{2}}+\xi ^{2}\sigma _{1}\sigma _{2}(1-\rho ^{2})]^{-1/2}.
\end{displaymath} (15)

Figure 11a shows a 3D presentation of how this PDF changes with rising $\rho $.
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig11.eps}
\end{figure} Figure 11: The product probability distribution with the correlation factor as a parameter, a) a 3D presentation; and b) sections for different correlation factors
Open with DEXTER

Figure 11b represents several sections of this dependence. It is clear that for $\rho =1$ this distribution is transformed into the distribution of a squared Gaussian value.

For a product of two Gaussian values z = xy with non-zero means the formula for the PDF is too complex, but the characteristic function is simple enough to calculate the necessary central moments and hence the skewness and the excess. The central moments are more convenient to calculate using cumulants: $\kappa_{r}=i^{-r}\frac{{\rm d}^{r}}{{\rm d}\xi^{r}}{\rm ln}f(\xi)\vert_{\xi=0}$. For variances equal to 1, the means a1, a2, and the correlation coefficient $\rho $, the characteristic function is (Deustch 1962):

$\displaystyle f(\xi )$ = $\displaystyle \frac{1}{[1-2i\xi \rho +\xi ^{2}(1-\rho ^{2})]^{1/2}}$  
    $\displaystyle \times\exp \left\{\frac{2i\xi a_{1}a_{2}-\xi ^{2}[a_{1}^{2}-2\rho a_{1}a_{2}+a_{2}^{2}]}{2[1-2i\xi \rho +\xi ^{2}(1-\rho ^{2})]}\right\}\cdot$ (16)

For the cross-spectral density at a particular frequency, we have to use $f_{\Sigma ^{2}}(\xi )=$ $f(\xi )^{2}$ due to two degrees of freedom, as for the auto-correlation spectrum. Thus the first two central moments corresponding to $f_{\Sigma ^{2}}(\xi)$ and to variances $\sigma _{1}^{2}$, and $\sigma _{2}^{2}$ (so that $\overline{x}=$ $a_{1}\sigma _{1}$, and $\overline{y}=$ $a_{2}\sigma _{2})$ are:

\begin{displaymath}m1_{\rm cor}=2\rho \sigma _{1}\sigma _{2}+2\overline{x}\overline{y},
\end{displaymath} (17)


$\displaystyle var_{\rm cor}=2\overline{x}^{2}\sigma _{1}\sigma _{2}+2\overline{y}^{2}\sigma _{1}\sigma _{2}+4\rho \overline{x}\overline{y}\sigma _{1}\sigma _{2}$      
$\displaystyle +2\sigma _{1}^{2}\sigma _{2}^{2}+2\rho ^{2}\sigma _{1}^{2}\sigma _{2}^{2}.$     (18)

For $\sigma _{1} = \sigma _{2} = \sigma$, and $\overline{x} = \overline{y} = A/\sqrt{2}$, we have the system of equations with respect to $\sigma$ and A:

\begin{displaymath}m1_{\rm cor}=2\rho \sigma ^{2}+A^{2},
\end{displaymath} (19)


\begin{displaymath}\mu _{2{\rm cor}}=var_{\rm cor}=2A^{2}\sigma ^{2}+2\rho A^{2}\sigma ^{2}+2\sigma ^{4}+2\rho ^{2}\sigma ^{4}.
\end{displaymath} (20)

For $\rho\to 1$, we find:
$\displaystyle m1_{\rm cor}$ = $\displaystyle A^{2}+2\sigma ^{2},$ (21)
$\displaystyle var_{\rm cor}$ = $\displaystyle 4A^{2}\sigma ^{2}+4\sigma ^{4},$ (22)

which is the same as for the auto-correlation spectrum. The following formulae are for the third and the fourth central moments of the sum of two products:
$\displaystyle \mu _{3{\rm cor}}$ = $\displaystyle 4\rho ^{3}\sigma ^{6}+12\rho \sigma ^{6}+6A^{2}\rho ^{2}\sigma ^{4}+12\rho A^{2}\sigma ^{4}+6A^{2}\sigma ^{4},$  
      (23)
$\displaystyle \mu _{4{\rm cor}}$ = $\displaystyle 24\rho ^{4}\sigma ^{8}+96\rho ^{2}\sigma ^{8}+24^{2}\sigma ^{8}+48A^{2}\rho ^{3}\sigma ^{6}+96A^{2}\rho ^{2}\sigma ^{6}$  
    $\displaystyle +96A^{2}\rho \sigma ^{6}+48A^{2}\sigma ^{6}+12A^{4}\rho^{2}\sigma ^{4}+24A^{4}\rho \sigma ^{4}+12A^{4}\sigma ^{4}.$  

Using these formulae for $\mu _{i}$, it is possible to calculate skewness
$\displaystyle skew_{\rm cor}$ = $\displaystyle \frac{\mu _{3{\rm cor}}}{\mu _{2{\rm cor}}^{3/2}}$  
  = $\displaystyle \frac{\sigma ^{4}\sqrt{2}(2\rho ^{3}\sigma ^{2}+6\rho \sigma ^{2}...
...\rho A^{2}+3A^{2})}{(\rho ^{2}\sigma ^{2}+\sigma ^{2}+\rho A^{2}+A^{2})^{3/2}},$  
      (24)

and excess
    $\displaystyle excess_{\rm cor}=\frac{\mu _{4{\rm cor}}}{\mu _{2{\rm cor}}^2}-3=$  
    $\displaystyle 3\sigma ^{2} \frac{(\rho ^{4}\sigma ^{2}+6\rho ^{2}\sigma ^{2}+\s...
...o A^{2}+2A^{2})}{(\rho ^{2}\sigma ^{2}+\sigma ^{2}+\rho A^{2}+A^{2})^{2}} \cdot$  
      (25)

When $\rho $ approaches 1, both Eqs. (25) and (26) reduce to Eqs. (10) and (11), respectively. Figures 12a,b show how the skewness and the excess depend on the RFI amplitude A with $\rho $ as a variable.
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{fig12.eps}
\end{figure} Figure 12: The skewness a) and the excess b) for a product with 2 degrees of freedom and $\rho $ as parameter
Open with DEXTER

When $\rho =1$ these curves are transformed into the curves in Figs. 2b,c corresponding to the sum of two squares of Gaussian values.

Acknowledgements
The author would like to thank W.A. Baan for helpful suggestions and valuable discussions.

References

 
Copyright ESO 2001