Transcript Slide 1

Statistics of the Spectral
Kurtosis Estimator
Gelu M. Nita and Dale E. Gary
New Jersey Institute of Technology
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
1
Population Spectral Kurtosis
The PSD estimate obtained by an N-points FFT
of a Gaussian time domain signal has an exponential distribution
 
p Pk 
1


e
Pk

; k 1
N
1
2
Population Spectral Kurtosis: statistical RFI discriminator
 k2
SK  2  1
k
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
2
Spectral Kurtosis Estimator
Power and squared power sums
M
S1   Pk
i 1
M
S2   Pk
2
i 1
Mean and variance unbiased estimators
2
MS

S
2
1
 k2 
M ( M  1)
1
k  S1
M
Spectral Kurtosis biased estimator
SK 
March 29, 2010
 k2
k
2
M  MS2 

 2  1
M  1  S1

RFI Mitigation Workshop, Groningen
The Netherlands
3
The SK Spectrometer
Key features:
Conceptual simplicity
Frequency channel independence
Straightforward FPGA implementation
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
4
Statistical thresholds for the rejection
of RFI outliers (M>>1)
(Nita et. al, 2007 PASP, 119, 805:827)
Assumed negligible bias
Hardware implementation of the SK
excision algorithm and initial tests revealed
that, although it performs generally well,
 
E SK  1
Variance (theoretically proven)
4
Var SK 
M
 
RFI symmetrical thresholds (1  3 )
(assumed SK normality)
6
1
M
March 29, 2010
•the lower threshold level is set too
low, failing to reject some RFI
contaminated channels
•the upper threshold level is set too
low, rejecting more RFI free channels
than statistically expected
CONCLUSION: More accurate RFI
rejection levels are needed for improved
and reliable performance.
RFI Mitigation Workshop, Groningen
The Netherlands
5
Characteristics of the SK Estimator
(Monte Carlo Simulations)
 
1
Bias: E SK  1 
M
Skewness:  1 
The probability distribution of the SK
estimator remains asymmetric even for a
fairly large accumulation number M,
approaching normality at a slower pace
than needed for practical applications.
10
M
Main goals of this study:
Redefine the SK estimator to
remove bias
Kurtosis Excess:
March 29, 2010
2 
250
M
Find an analytical expression for
probability distribution of the SK
estimator to allow accurate calculation
of its tail probabilities
RFI Mitigation Workshop, Groningen
The Netherlands
6
Starting Point
Spectral Kurtosis estimator
M  MS 2 
SK 
 2  1
M  1  S1

SK statistics is determined by the statistics of the
mean of squares to square of mean ratio
S2
2
P
MS2
M  k

2
2
S12
Pk
 S1 
 
M 
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
7
Joint Distribution (Monte Carlo):
Mean of Squares-Square of Mean



March 29, 2010
S2 and S12 are strongly
correlated random variables
Their (unknown) joint
distribution would be needed
to derive the distribution of
their ratio
Is there any work-around
approach?
RFI Mitigation Workshop, Groningen
The Netherlands
8
Joint Distribution (Monte Carlo):
Mean of Squares to Square of Mean Ratio - Square of Mean
Monte Carlo simulations suggest:


S2 /S12 and S12 are uncorrelated
random variables
Their individual distributions are
independent
This is the fundamental property that makes
the whole SK concept work by allowing
SK to have a unity value independently of
the power level
This property was analytically proven for
the exponential distribution based on
first principles (Nita & Gary, PASP, in
press)
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
9
Raw statistical moments of the
Mean of Squares to Square of Mean Ratio
Cov  f n , g n   n 2 f n 1 g n 1Cov  f , g 
Cov  x, y   E  xy   E  x  E  y 
 MS2  S1 2 
S
Cov  2 ,     E  2
M
 S1  M  
2
 MS2   S1  

  E  2  E     0

 S1   M  
 S2 
E
 MS 
 M 
 
2
E
 2 
 S  2 
S
 1  E 1 
 
 M  
  
March 29, 2010
n 1
2( n 1)
 MS n  S 2 n 
 MS2  S1 2 
 S1 
2  MS2 
2
1
Cov  2  ,     n  2   
Cov  2 ,     0
 S1   M  
 S1  M  
 S1   M 

n


S
 2 
E
 M  


n

 MS  
 




2

E 
 S 2  n 
 S 2  
 1   E   1   


  M   
RFI Mitigation Workshop, Groningen
The Netherlands
10
Analytical Results
(Nita & Gary, PASP, in press)
Raw moments of the square of mean
2 n
2n



M

1

2
n
!
Raw moments
of
the
mean
of
squares
to
square


S

  1 
  of mean ratio
E      
 
M
M

1
!


M 


   n   
 MS   M n  M  1!  n  n  2r !  M
r
2 
E 
t

   M  1  2n ! t n  r !
2
Raw
 r 0 s
 t 0
S moments
of the mean of square

 1  
M


 S 2  n     2 n  n  n  2 r  ! r 
E     
t 

n 
 M    M  t  r 0 r !  t 0
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
11
Redefinition of the SK Estimator
Bias of the originally proposed SK estimator
 M  MS2  
M   MS 2  
M
1
E SK  E 
 1
 E  2   1 
 2  1  
M
  M  1   S1   M  1
 M  1  S1
 
The unbiased SK estimator
M  1  MS 2 
SK 
 2  1
M  1  S1

 
E SK  1
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
12
Standard Moments of the
SK Estimator
 
1/  E SK  1
4M 2
4
 1 
2 

O 2 
 M  1 M  2  M  3 M  M 
32 4  M  2  M  3 5M  7 
1  3 
2
2
2
 M  1 M  4   M  5
2
3
2
4 3  M  2  M  3  M  98M  185M  78 
2  2 
2
2
 M  1 M  4   M  5 ( M  6)( M  7)
10
 1 
 O  3/2 
M
M 
256
 1 
 2  2  3 
O 2 
M
M 
 1  1 
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
13
Moment based approximation of the
SK estimator distribution using
Pearson Probability Curves
1   2  3

4  4 2  31  2 2  31  6 
2
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
14
Pearson Type IV PDF
 Pearson 1895; Nagahara 1999 
 


m  i m i 
2

1
x



2 
2

 
 x   

p( x) 
1  
  Exp   ArcTan 

1
a
a


a 







  m     m
2

 Heinrich 2004    
March 29, 2010
r (r  2) 1
16(r  1)  1 (r  2) 2
r
6(  2  1  1)
2 2  31  6
a
m
r2
2
  1/  (r  2) 2 1
1
2
2 6  r  1  1  r  2  


4
1
4
RFI Mitigation Workshop, Groningen
The Netherlands
15
Pearson IV CF and CCF
To compute the tail probabilities of the SK estimator, one needs to evaluate
the cumulative function (CF) and complementary cumulative function (CCF)
of the Pearson IV probability curve
x
P( x ) 
 p( x)dx;


1  P( x)   p( x)dx
x
Knowing the analytical expression for the Pearson IV PDF, the CF and CCF
can be computed analytically, by using the closed form expressions involving
Hypergeometric series provided Heinrich(2004) or Willink(2008).
Alternatively, CF and CCF can be computed by a simple numerical
integration.
The asymmetrical RFI thresholds are then chosen so as to provide
symmetric tails probabilities of rejecting true Gaussian signals of 0.13499%,
which are equivalent to the ±3 thresholds of a normal distribution.
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
16
RFI Threshold Computation Example
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
17
Pearson IV PDF vs. Monte Carlo Simulations
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
18
Time Domain Kurtosis Estimator
The DC (k=0) and Nyquist (k=N/2) frequency channels of a
FFT Spectrometer obey a different statistics.
The same statistics is obeyed by the power corresponding
to a narrow band time domain signal (FIR filter spectrometer)
 
p Pk 
1

P
Pk
1/2 

k
e
Following similar steps, an unbiased Time Domain Kurtosis (TDK) estimator
may be defined in this case to detect RFI contamination.
K
March 29, 2010
M  2  MS 2 
 2  1
M  1  S1

RFI Mitigation Workshop, Groningen
The Netherlands
19
Standard Moments of the
time domain Kurtosis Estimator
 
1/  E K  2
24 M 2
24
 1 
2 

O 2 
 M  1 M  4  M  6  M  M 
32 216  M  4  M  6  M  2 
1  3 
2
2
2
 M  1 M  8  M  10 
2
3
2
4 3  M  4  M  6   M  213M  474 M  368 
2  2 
2
 M  1 M  8 M  10  (M  12)(M  14)
6 6
 1 
 O  3/2 
M
M 
540
 1 
 2  2  3 
O 2 
M
M 
 1  1 
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
20
Kurtosis estimator Pearson IV PDF
and RFI thresholds (M>45)
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
21
Summary
The Spectral Kurtosis RFI excision algorithm recommends itself by
Conceptual simplicity
Frequency channel independence
Straightforward FPGA implementation (See Dale Gary’s following presentation)
Theoretically determined RFI thresholds for arbitrary integration time
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
22
Main References





Gelu M. Nita and dale E. Gary, “Statistics of The Spectral Kurtosis
Estimator”, 2010, PASP, in press
Gelu M. Nita, Dale E. Gary, Zhiwei Liu, Gordon J. Hurford, & Stephen M.
White, 2007, "Radio Frequency Interference Excision Using Spectral Domain
Statistics," Publications Of The Astronomical Society Of The Pacific, 119,
805.
Yuichi Nagahara, 1999, "The PDF and CF of Pearson type IV distributions
and the ML estimation of the parameters," Statistics & Probability Letters, 43,
(1999), page 251.
Joel Heinrich, 2004, "A Guide to the Pearson Type IV Distribution," Collider
Detector at Fermilab internal note 6820, 2004, http://wwwcdf.fnal.gov/publications/cdf6820 pearson4.pdf
Willink, R., 2008, "A Closed-Form Expression For The Pearson Type IV
Distribution Function,“ Aust. N. Z. J. Stat. 50(2), 199, 205
March 29, 2010
RFI Mitigation Workshop, Groningen
The Netherlands
23