IntroductionToSignalProcessing - DCC

Report
The Laser Interferometer
Gravitational-Wave Observatory
http://www.ligo.caltech.edu
Supported by the United States National Science Foundation
Introduction To Signal Processing & Data Analysis
Gregory Mendell
LIGO Hanford Observatory
LIGO-G1200759
Fourier Transforms
~
x ( ) 



1
x (t ) 
2
x (t )e



 i t
dt
i t
~
x ( )e d
  2f
e
it
 cost  i sin t
Digital Data
x
. . . . .
j=0
j=1 j=2
j=3
. .
…
j=4
j=N-1 T
t
Sample rate, fs:
f s  1 /(t )
Duration T, number of samples N, Nyquist frequency:
f sT  N ;
f Nyquist  f s / 2
Time index j:
t j  jt ;
Frequency resolution:
j  0,1,2,...N  1;
f  1 / T
Frequency index k:
fk  k / T ;
f kT  k ;
k  0,1,2,...N  1
Discrete Fourier Transforms
(DFT)
N 1
 2ijk / N
~
xk   x j e
j 0
1
xj 
N
N 1
2ijk / N
~
 xk e
k 0
The DFT is of order N2. The Fast Fourier Transform (FFT) is a
fast way of doing the DFT, of order Nlog2N.
DFT Aliasing
f f;
~
x k  ~
xk*
For integer m:
f  f  m fs ;
Useful Band :
~
xk  mN  ~
xk ;
~
xN k  ~
xk*
[0, f Nyquist ]  k  [0, N / 2]
Power outside this band is aliased into this band. Thus
need to filter data before digitizing, or when changing fs,
to prevent aliasing of unwanted power into this band.
DFT Aliasing
DFT Orthogonality
N 1
N 1
N 1
j 0
j 0
j 0

2ijk '/ N  2ijk / N
2ij ( k '  k ) / N
2i ( k '  k ) / N
e
e

e

e




j
N
1

r
S   r j ; rS   r j 1; rS   r j ' j 1   r j '  1  r N ; S 
1 r
j 0
j 0
j '1
j '0
N 1
 e
N 1
N 1
2i ( k '  k ) / N

j
j 0
N 1


2i ( k '  k ) / N N
1 e

1  e 2i ( k ' k ) / N
e
j 0
N 1
N
2ijk ' / N
e
1  e 2i ( k ' k )

 0, k  k ' ; N , k  k '
2i ( k '  k ) / N
1 e
 2ijk / N
 N kk '
Parseval’s Theorem
N 1
x
j 0
N 1
2
j
1

N
1
xj yj 

N
j 0
N 1
~
 xk
2
k 0
N 1
*
~
~
 xk y k
k 0
Correlation Theorem
N 1
c j'   x j y j j'
j 0
*
~
~
~
c k  xk y k
Convolution Theorem
N 1
C j '   x j y j ' j
j 0
~ ~~
C k  xk y k
Example: Calibration
•
•
•
•
•
•
Measure open loop gain G,
input unity gain to model
Extract gain of sensing
function C=G/AD from
model
Produce response function at
time of the calibration,
R=(1+G)/C
Now, to extrapolate for
future times, monitor
calibration lines in
DARM_ERR error signal,
plus any changes in gains.
Can then produce R at any
later time t.
A photon calibrator, using
radiation pressure, gives
results consistent with the
standard calibration.
DARM_ERR
G ( f )  C ( f ) D( f ) A( f )
1  G( f )
R( f ) 
C( f )
xext ( f )  R( f ) DARM _ ERR( f )
One-sided Power Spectral
Density (PSD) Estimation
Pk 
~
2 xk
2
t
2
T
Note that the absolute square of a Fourier Transform gives what
we call “power”. A one-sided PSD is defined for positive
frequencies (the factor of 2 counts the power from negative
frequencies). The angle brackets, < >, indicate “average value”.
Without the angle brackets, the above is called a periodogram.
Thus, the PSD estimate is found by averaging periodograms. The
other factors normalize the PSD so that the area under the PDS
curves gives the RMS2 of the time domain data.
Gaussian White Noise
n j  0;
~
nk
2
1

N
nj
N 1
2
1

N
N 1
N 1
n
j 0
2
j

2
~
 nk   n j N
k 0
2
2
j 0
Parseval’s Theorem is used above.
2
Power Spectral Density,
Gaussian White Noise
2

2
2 N 2t 2
2
 2 t 
Pk 
fs
T
P
2
fs
This ratio is constant,
independent of fs.
The square root of the PSD is the
Amplitude Spectral Density. It has
the units of sigma per root Hz.
2 2 1 N 2 2
2
2
RMS

Area





f

P


k
2 N
k 0 f s T
k 0
N /2
N /2
For gaussian white noise, the square root of the area
under the PSD gives the RMS of time domain data.
Amplitude and Phase of a
spectral line
x j  A cos(2ft j  0 )  A
fT  k
AN i0
~
xk 
e
2
e
2ift j  i0
e
2
2ift j i0
If the product of f and T is an
integer k, we call say this
frequency is “bin centered”.
For a sinusoidal signal with a “bin
centered” frequency, all the power
lies in one bin.
DFT Bin Centered Frequency
DFT Leakage
x j  A cos(2ft j  0 )  A
fT   ;
    k
e
2ift j  i0
e
2
2ift j i0
For kappa not an integer get sinc function response:
AN i0  sin(2 ) 1  cos(2 ) 
~
xk 
e 
i

2
2



2





For a signal (or spectral
2
2
2
A N sin ( ) disturbance) with a non2
~
xk 
bin-centered frequency,
2
2
4
 
power leaks out into the
neighboring bins.
DFT Leakage
Windowing reduces leakage
Hann Window:
1
 2j 
w j  1  cos

2
 N  1 
N
N
N
~
wk   k ,0   k ,1   k , 1 ; for
2
4
4
N
1
w
w
~
~
x j  w j x j  xk   ~
xk ' w
k k '
N k ' 0
1~ 1~
1~
w
~
xk  xk  xk 1  xk 1
2
4
4
N  N 1
Corrections to PSD
Am plitude Correction:
~x  2 ~x w
k
k
8 ~w 2
2
~
Energy Correction:
nk  nk
3
2
2
2
2
1AT
1 2 ( A / 2) T / 2 1 A (2T / 3)
2
SNR 


2
2 Sk
2 2 (3 / 8) S k
2
Sk
3
f  f
2
Example PSD Estimation
For iLIGO and aLIGO
PSD Statistics
n~  x  iy
2
~
n  x2  y2;
( x, y )dxdy 
if  x  y  0;  x 2  y 2  1
1  x2 / 2
e
2
1  y2 / 2
e
dxdy
2
r  x 2  y 2 ;   t an1 ( y / x); dxdy  rdrd
Rayleigh Distribution:
1 r 2 / 2
r 2 / 2
(r ,  )drd 
re
drd ;
(r )dr  re
dr
2
Chi-squared for 2 degrees of freedom:
1  / 2
2 1
  r ; d  rdr;
(  )d  e d
2
2
Histograms, Real and
Imaginary Part of DFT
For gaussian noise, the above are also gaussian distributions.
Rayleigh Distribution
Histogram of the square root of the DFT power.
Histogram DFT Power
Result is a chi-squared distribution for 2 degrees of freedom.
Chi-Squared Statistics
A chi-squared variable with ν degrees of freedom is
the sum of the squares of ν gaussian distributed
variables with zero mean and unit variance.
  x 2  y 2  z 2  ...;
r  x  y  z  ...;
2
2
2
2
if  x  y  z  ...  0
 x  y  z  ...  1
2
e.g., dxdydz r sin drdd ;
2
1/ 2   / 2
(  )d   e
2
2
  r ; d  2rdr
2
d

1
 / 2
2
1
(  )d  /2
 e
2 ( / 2)
d
Maximum Likelihood
Likelihood of getting data x for model h for Gaussian Noise:
1
P ( x | h) 
e
2  1
 ( x1  h1 ) 2
2 12
1
2  2
 ( x2  h2 ) 2
e
2 22
1
e
2  3
 x jhj 1 hjhj

 


2
  2
2
2
 j 
j
2 j j
j
j

e.g ., h j  f (h0 , cos ,  0 , )
2
( x j  h j )2
 2
 0,
h0
 2
 0,
 cos
 2
 0,
 0




 ( x3  h3 ) 2
2 32
...
Chi-squared
 2
0

Minimize Chi-squared = Maximize the Likelihood
Matched Filtering
 x j h j 1 h j h j   ~xk*h~k 1 h~k*h~k 

ln     2   2    
 

 j 

2

S
2
S
j
k
k
j
j 
k
k 


1
ln   ( x, h)  (h, h) This is called the log likelihood; x is
the data, h is a model of the signal.
2
A2
Maximizing over
e.g ., h  Ah;
ln   A( x, h)  (h, h)
2
a single amplitude
gives:
 ln  / A  ( x, h)  A(h, h)  0
~* ~
*~
2
~
xk hk
hk hk
( x, h )
1 ( x, h )
A
;
ln  max 
;
SNR  
/ 
( h, h )
2 ( h, h )
Sk
Sk
k
k
This can be generalized to maximize over more parameters, giving the
formulas for matched-filtering in terms of dimensionless templates.
Frequentist Confidence Region
Matlab simulation for signal:
•
Estimate
parameters by
maximizing
likelihood or
minimizing chisquared
•
Injected signal
with estimated
parameters into
many synthetic
sets of noise.
•
Re-estimate the
parameters from
each injection
s(t )  A cos2ft B sin 2ft
LIGO-G050492-00-W
Figure courtesy Anah Mourant, Caltech SURF student, 2003.
Bayesian Analysis
Bayes’
Theorem:
Likelihood:
Posterior
Probability:
Confidence
Interval:
P ( a ) P ( b | a )  P ( b) P ( a | b)
1
e
P( x | h ) 
2  1
( x1  h1 ) 2
2 12
1
2  2
P(h )
P( x | h )
P(h | x ) 
P( x )
hC
C   P( h | x )dx
0
Prior
e
( x2  h2 ) 2
2 22
...
Bayesian Confidence Region
s(t )  A cos2ft B sin 2ft
P( A, B; x)  exp(  2 / 2)
Matlab simulation for signal:
Uniform Prior:
x
LIGO-G050492-00-W
Figures courtesy Anah
Mourant, Caltech SURF student, 2003.
The End

similar documents