Transcript 2-Danaher

Signal Processing For
Acoustic Neutrino Detection
(A Tutorial)
Sean Danaher
University of Northumbria,
Newcastle UK
Sean Danaher
What is Signal Processing?
• Used to model signals and systems where there is
correlation between past and current inputs/outputs
(in space or time)
• Two broad categories: Continuous and Sampled
processes
• A host of techniques Fourier, Laplace, State space,
Z-Transform, SVD, wavelets……
• Fast, accurate, robust and easy to implement
• Sadly also a big area. Engineers spend years
learning the subject (Not long Enough!)
• Smart as Physicists are I can’t cover a 3+ year
syllabus in 40 minutes!
• Notation (not just i and j)!
Sean Danaher
How Can Signal Processing be used
in Acoustic Neutrino Detection?
At Virtually all stages of the process
1. Accurate Parameterisation of distributions with no amenable
analytical form e.g. Shower Energy distribution profiles
2. Speed up computation e.g. Acoustic integrals by many orders of
magnitude
3. Design, Understand, simulate Hydrophones, Microphones and
other acoustic transducers (and amplifiers)
4. Design, Understand, simulate filters both analogue and digital
(tailored amplitude, phase response, minimise computing time)
5. Estimate and recreate noise and background spectra
6. Design Optimal Filters e.g. Matched Filters
7. Design Optimal Classification Algorithms (Separating Neutrino
Pulses from Background)
8. Improve Reconstruction accuracy……………
Will cover 1-5 in this talk……
Sean Danaher
Singular Value Decomposition
( )( )( )
W
L
VT
We can get an approximation
of the original data by setting
the L values to zero below a
certain threshold
Similar techniques are used in
statistics CVA, PCA and Factor
Analysis
Based on Eigenvector Techniques
Good SVD algorithms exist in ROOT and MATLAB
Sean Danaher
Decompose the
Data matrix into 3
matrices.
When multiplied
we get back the
original data.
W and V are
unitary. L contains
the contribution
from each of the
eigenvectors in
descending order
along main
diagonal.
SVD II
Somewhere between vector
five and seven we get the
“Best” representation of the
data
•Better than the observation as
noise filtered out
•Highly compressed (only 1-2%
of original size e.g. 6/500x1.5)
SVD done on
Noisy data
Sean Danaher
•Have basis vectors for data so
can produce “similar” data+++
Building the Radial Distribution
4 examples chosen have
maximum variation in shape
Hadronic Component
using Geant IV
SVD works well
with the radial
distribution
•Three vectors sufficient
to fit the data
•Can use a linear
mixture of these three
vectors as input to the
Acoustic MC
We Reproduce both the shape and variation
Sean Danaher
Acoustic Integrals
P (t ) 
Slow decay
fast thermal
energy deposition
d2
dt 2
Sean Danaher

E0  d
 (r / c  t ) dv
4 Cp r dt
Band limited by
•Water Properties
•Attenuation
Pressure (Pa)
P (arbitrary)
Method I
0.1
0.5
0.05
0
0

-0.2
0.20
2
-0.05
3
4
0
Time (ps)
5
-0.1-6
-4
-2
0
2
Time (s)
4
6 -4
x 10
5
6
•Throw a number if MC points ~107
•Assume the energy at each point in
the cascade deposited as a Gaussian
Distribution.
•Calculate distance to observer from
each point
•Propagate the Gaussian to the
observer (an analytical solution is
known for this)
•Sum over all the points
•Provided the width of the Gaussian is
small in comparison to the shower
radius we will get the pressure pulse
Sean Danaher
7
8
Expensive
5 flops per time point
If using 1024 points
on time axis
5*1024*1e7=5.1e10 flops
9
10
11
3 x 10
4
2.5
Frequency
-0.5-5
2
1.5
1
0.5
00
2
4
6
8
Longitudinal Distance (m)
10
12
Radial Distribution
Need to do this thousands of
times for different distances,
angles, distributions, etc.
12
The Convolution Integral
Given a signal s(t) and a system with an
impulse response h(t) (response to (t)) then
y(t) to an arbitrary input is given by

y (t )   s(t )h(t   )d

s(t) Bipolar
Acoustic Pulse
h(t) Impulse response
high pass filter (SAUND)
y(t) Electrical
Pulse
We can use this to process all the MC points in the cascade simultaneously
Sean Danaher
A Few examples
RC circuit
Mixed signal
This situation is similar
to our acoustic integral
*
Sean Danaher
=
2d
Convolution Theorem
Convolution in the time domain is
Multiplication in the frequency domain
(and visa versa)
Will speed the process still further
• Convolution very expensive computationally
(o=n2)
– Convert the signal and impulse response to the
frequency domain (Fourier Transform)
– Provided the number of points n=2m very efficient
FFTs o=n log n
– Multiply and take Inverse Transform
Sean Danaher
Acoustic Integrals II
•Assume each point produces a Delta Function
•The sampled equivalent of a Delta function is 1
•The integral can be done using a histogram function
(modern histogram algorithms are very efficient)
•The derivative can be done in the frequency domain
(d/dt  i No overhead)
•We then convolve the response with the entire signal
P (t ) 
E0  d

4 Cp dt
Sean Danaher


E0  d
 (r / c  t ) dv
4 Cp r dt
E0  d
 (r / c  t )dv 
E xyz (t )
r
4 Cp dt

Only Histogram depends on the
number of MC points
P [t ] 
n 511

n 512
E xyz ( )w ( )D( )Aw ( )e inT
FT of Exyz
Derivative j
Blackman window Water Attenuation
Histogram: 103 flops
Scale histogram by 1/d: 103
flops (only needed in near field)
FFT nLogn c 104 flops
Multiplication c 4x103 flops
IFFT nLogn c 104 flops
Sean Danaher
6 orders of
magnitude less
than approach 1
Limitation: Will not work if volume of the
integral is so great that attenuation
varies significantly (Large 100+m
source in near field)
State Space Analysis
Mode 1
MIMO systems
Mixed Mechanical/electrical models etc
Matrix Based Method
Quantum Mechanics
Mode 14
Mode 100
X = AX + BU
Y = CX + DU
All modern control
algorithms use SS
methods.
 1 
 
d
 2 
i
 H
 3 
dt
 
   4 
 . 
i   H
 
 . 
 
 n
Pictured 20000 state
simulation of a square acoustic membrane:
(100x100 masses 2 states x and v)
A Matrix 400x106 elements
Sean Danaher
SS Implementation
States
Inputs
States
A
B
Outputs
C
D
A is of size States  States
B is of size States  Inputs
C is of size Outputs  States
D is of size Outputs  Inputs
Sean Danaher
States are degrees of
freedom of the system.
Things that store energy
•Capacitor Voltages
•Inductor currents
•Positions and Velocities
of masses
dV
 Cx
dt
di
if x  IL  VL  L  Lx
dt
if x  position  x  velocity
if x  Vc  IC  C
if x  velocity  x  acceleration
Simple Example LC
Ic  C
dVc
dt
VL  ()L
VL  Vc
1

0

C
A

 1 0 


 L

IL  IC
x1  VC Cx1  x2
x 2  IL
dIL
dt
y1  x1
Lx2   x1 y 2  x2
Response to Initial Conditions
1.5
Voltage
1
0.5
0
Amplitude
-0.5
-1
-1.5
1.5
Current
1
0.5
0
-0.5
-1
-1.5
0
5
Sean Danaher
10
15
20
Time (sec)
25
30
35
40
 1 0
C

0
1


Maths simple to
implement
Example shows
the behaviour
with an initial
1V on the
Capacitor (1F)
and 1A flowing
through the
inductor (1H)
Speaker/Microphone
Frame
Cone
Magnet
assembly
N
S
Magnet
assembly
Coil
d 2x
dx x

Mechanically Force produced by current Fe  Bi  m 2  r
x =x, (position)
dt
dt s
1
x2=dx/dt, (velocity)
x3=i, (current)
Loss to air term
Electrically coil moving in a magnetic field
Creates an EMF
V B
dx
di
 L  Ri
dt
dt
Bode Diagram
-60
x
2.5 10
Magnitude (dB)
-80
-120
Amplitude
-140
-160
90
Phase (deg)
Step Response




 0
 0 0
1
0 
2




1 r
B 
1
X
X  
0  U 1.5
 sm m



m
m
1




1

B

R
 0

 0

0.5
L

L
L 

0
Y  ( 0 1 0) X

(0 0)U
-100
0
-90
-1802
10
-4
3
10
4
10
Frequency (rad/sec)
Sean Danaher
5
10
6
10
-0.5
0
0.5
1
1.5
Time (ms)
2
2.5
3
Simple Hydrophone Model
x1  VC , x2  v, x3  x
F=KVc
Mass
m
Heart of Hydrophone
Piezo electric crystal
Piezo electric effect
K

 1
0 

 1
  RC  C
0

F=x/S


RC


1 
r
K
B   0 1


A
 m
Sm 
m
 0 0


0 
1
0



F=ma







1

 1
0
0 0



D R
C R
F=rv




 0 0
 0 1 0
Omni works
in breathing
mode
40cm
Vin  Vc
R
30cm
I  I0e
Ic
Vin

x

Kv
Vc
“Erlangener” water tank (Niess)
Simulated Hydrophone signal
3rd
order simulation.
Do we need higher?
Amplitude (Arbitrary)
25
=5e-3
20
15
10
5
0
Gaussian
cross-section
-5
-10
-1525
Sean Danaher
= 0.09
30
35
40
Time (eq. cm)
45
50
Acoustic Simulation
•Bipolar Pulse
•1000 element simulation
 0

 -8
 0

4
A = 
 0

 .
 .

 0
1
0
0 0 0 0 0 0
4 0 0 0 0 0
.
.
. .
. .
.
.
0 0 1 0 0 0 0
0 -8 0 4 0 0 0
.
.
. .
. .
.
.
•Two types of medium; Heavy and light
0 0 0 0 1 0 0 . . . .
•pulse reflects against walls
.
. . . . . . . . . .
•mismatch between sections
.
. . . . . . . . . .
•velocity on each section
0 0 0 0 0 0 0 0 0 1 0
•amplitude change
Enhance to 2D?
•dispersion
Embed Hydrophone models etc.?
Sean Danaher
. 0

. 0
. 0

. 0
. 0

. .
. .

-2 0 
Butterworth Filter
1930s
1
| H ( ) |
1   2n
Stable
Good Frequency
Response
Buterworth Impulse Response
1.2
1st order
2nd order
1
4th order
Note Impulse response
•Output Delay
•Oscillation
Output
0.8
16th order
32nd order
0.6
0.4
0.2
0
-0.2
0
Sean Danaher
8th order
10
20
30
Time(s)
40
50
Implementation
Passive Filters
1st Order Familiar RC
2nd Order LRC
3rd Order “”
At acoustic frequencies easiest
to use op-amps. Expensive to
make high precision lossless
inductors
1
s4 +2.6131s3 +3.4142s 2 +2.6131s+1
1
1
 2
 2
s +0.7654s+1 s +1.8478s+1
=
ADC
R
| h(s ) | 2
, ADC  1  a
s  (3  ADC )s  1
Rb
+
Rb
Sean Danaher
+
Ra
=
-
Cascade 2nd Order
sections to make
Butterworth high
order filters
Recovering Phase information
e.g. SAUND Data
It is trivial to design
digital filters which
have a constant
group delay d/d
and hence no phase
distortion
0.15
Amplitude (Arbitrary)
0.1
S+N
causal
non causal
0.05
0
-0.05
-0.1
0
0.2
0.4
0.6
time (s)
0.8
1
1.2
If however we know the
filter response e.g.
Butterworth we can run the
data through the filter
backwards.
-3
x 10
This increases the order of the filter by a factor of 2
Sean Danaher
The Digital Filter
The digital filter is simple! Based upon
sampled sequences
a0 y [n]  a1y [n  1]  a2 y [n  2]  ....  b0u[n]  b1u[n  1]  b2u[n  2]  ....
Negative coefficients =>Causal
10
2
6
1
y [n ]   x[n ]  x[n  1]  5
2
9
Crude low-pass 8
Called FIR or
5
0
MA
8
4
Simple moving
average Filter
Sean Danaher
5
6
4
5.5
7
8.5
6.5
2.5
4
6
Perfect Integrator 6
8
y[n]  y[n  1]  x[n] 9
7
Called IIR
2
or
4
AR
9
9
4
9
6
14
23
30
32
36
45
54
58
67
The Z Transform and
Sampled Signals
X (z) 
n 

x [ n ]z  n
n 
if x  1, 3, 4  X ( z )  1  3z 1  4z 2
if x  1, a, a ,...  X ( z ) 
2
n 

n 0
az 1

n

1
1  az 1
Geometric sequences are of prime interest because
u(t )  eit  u[t ]  ei nt is a geometric sequence
Sean Danaher
A few z transforms
   /10,  1
1
0
-1
-1
0
1
3
2
1
0
-1
-2
-3
-4
0
h( z ) 
5
10
15
20
25
 
y [n]  x[n]  2 cos   y [n  1]  y [n  2]
 10 
30
   /14,  1.05
1
0
-1
10
5
0
-5
-10
-15
0
Poles must be
inside the unit
circle for stability
5
10
15
20
25
30
   / 5,  0.95
1
0
-1
2
1.5
1
0.5
0
-0.5
-1
-1.5
0
Sean Danaher
1
 
z 2  2 cos   z 1  1
 10 
Frequency Response
20
15
10
5
5
10
15
20
25
30
0
0
0.2
0.4
0.6
0.8
Fraction of Nyquist frequency
1
Simple Transfer Function
z  0.5
z  0.5
y [n ]  x [n ]  0.5 x[n  1]
 0.5 y [n  1]
h( z ) 
3
2.5
|H()|
2
1.5
2
0.5
00
0.2
0.4
0.6
0.8
frequency (Nyquist=1)
1
60
(degrees)
1
1
H()
Frequency
response simply
determined by
running around
the unit circle
 Corresponds to
the Nyquist
Frequency (fs/2)
50
40
30
20
1
0
10
| H ( z ) |
1
2
H ( z )  1  2
-1
Sean Danaher
00
0.2
0.4
0.6
frequency (Nyquist=1)
0.8
1
Notch Filter
h( z ) 

5
( z  e )( z  e
i

5
i

5
1.4
)
i

5
( z  0.95e )( z  0.95e )
y[n]=x[n]-1.6180x[n-1]+x[n-2]
+1.5371y[n-1]-0.9025y[n-2]
Fs=1000Hz
1.2
1
|H()|
i
0.8
1
0.6
0
0.4
-1
0.2
0
0
/5
-1
Signal buried in 100Hz
100
1.5
200
300
frequency
0
1
400
500
Recovered signal
1
0.5
0
-0.5
-1
-1.5
0
Sean Danaher
0.05
0.1
0.15
Time
0.2
0.25
Spectral Analysis
Using the Fourier Transform is seldom
the best way to get a spectrum. Normally
methods based around
• Autocorrelation (AC)
•Linear Prediction are used
First 500 of 2048 data points
4
3
1
0
Autocorrelation
1
-1
-2
-3
-4
0
100
200
300
Sample number
400
500
0.18
Amplitude(Relative)
0
-100
-50
0
Lag
50
100
150
Linear Prediction
0.14
a0 y [n]  a1y [n  1]  a2 y [n  2]....  e[n]
0.12
0.1
All pole filter driven by white
noise. Need to choose the order
with care. But can now
reproduce the spectrum
0.08
0.06
0.04
0.02
0
0.5
-0.5
-150
fft
ac
lpc
0.16
AC method
Use FT of AC
To get PSD
Amplitude (relative)
Amplitude
2
0
0.2
Sean
0.4
0.6
Frequency
(relative)
Danaher
0.8
1
Spectral estimation example
pitch
gain
64k
•8 bit 8kHz
voiced/unvoiced
switch
vocal tract
filter
•Split speech into
frame 240 samples
long
synthesized
speech
•Use LPC10 to
estimate spectrum
noise
•Reconstruct
2.4k
8
spectrum of frame 48
7
1.5
1
5
h(t)
|h()|
6
933
4
0.5
3
0
2
-0.5
1
00
Sean Danaher
Impulse Response
2
1000
2000
3000
Frequency (Hz)
4000
-1
0
10
20
30
Sample [n]
40
50
Conclusions
• SP Techniques Useful
• Thanks for listening
• Questions
Sean Danaher