Transcript Presentation (PowerPoint File)

Fast Wavelet Estimation of Weak
Biosignals
By
Elvir Causevic
Department of Applied Mathematics
Yale University
Founder and President
Everest Biomedical Instruments
1
Overview





Introduction and Motivation
Human auditory system
Measurement of auditory function and difficulties in signal processing
Introduction to wavelets and conventional wavelet denoising
Novel wavelet denoising algorithm




Frame recombination
Denoising
Variable threshold selection
Estimation of rate of convergence
 Experimental results
 Future work
 Conclusion and summary
2
Introduction
 Overall goal
 Creation of a fast estimator of weak biosignals based on
wavelet signal processing. Application to auditory
brainstem responses (ABRs) and other evoked potentials
 Specific objectives
 Reduce the length of time to acquire a valid ABR signal.
 Allow ABR signal acquisition in a noisy environment.
 Key obstacles
 Very large amount of acoustical and electrical noise
present .
 Signals collected from ear and brain have very low SNR
and require long averaging times
3
Infant Hearing Screening
• Infant hearing screening is critically important in early
intervention of treating deafness.
• Hearing loss affects 3 in 1,000 infants: most commonly occurring
birth defect.
• 25,000 hearing impaired babies born annually in the U.S. alone.
• Lack of early detection often leads to permanent loss of ability to
acquire normal language skills.
• Early detection allows intervention that commonly results in
development of normal speech by school age.
• Intervention involves hearing aids, cochlear implants and
extensive parent and child education and training.
• 38 U.S. states mandate hearing screening, Europe, Australia,
Asia following closely.
4
Measurement of Hearing Function
 Auditory
Brainstem Response (ABR) neural test
– Response of the VIIIth nerve - auditory neuropathway to brain
VIIIth Nerve
5
Auditory Brainstem Response (ABR)
Signal Processing & Clinical Issues
for Infant Hearing Screening








Stimulus: 37 clicks per second, 65 dB SPL (30 dB nHL).
Response: scalp electrodes measure μV level signals.
Noise: completely buries the response (-35dB).
Pass: signal to noise ratio measure (called Fsp) greater than an
experimentally determined value (NIH Multicenter study).
With linear averaging, reliable results are obtained within ~15
minutes of averaging of ~ 4000-8000 frames at a single level.
We would like to test multiple levels (up to 10) , and with multiple
tone pips (vs. clicks). This test normally takes over an hour, in a
sound attenuated booth, manually administered by an expert.
Currently only a single level response is tested and only a
pass/fail result is provided, with over 5% false positive rate.
Substantial improvement in rate of signal averaging is required to
obtain a full diagnostic and reliable test.
6
Auditory Brainstem Response
example
7
8
Infant Hearing Screening
Space Limitations
Time Constraints
Patient Tracking
Acoustic Noise
Electrical Noise
9
Auditory Brainstem Response (ABR)
Signal Processing & Clinical Issues
G
ar
Qra
e
uickT
ph
need
ics
im
ed
e™
deco
toan
se
md
pre
ea
th
ssor
is picture.
0
-10
-20
-30
dBV
-40
-50
-60
-70
-80
-90
-100
100
1000
Frequencyin Hz
10000
Frequency domain characteristics of a typical
ABR click stimulus as measured in the ear using the ER-10C transducer
10
Auditory Brainstem Response (ABR)
Signal Processing & Clinical Issues
Typical single 512-sample frame with the final average ovelaid
(Subject 3; right ear; 65 dB click)
Typical ABR waveform with manually labeled peak latencies
(Subject 3; right ear; 65 dB click; 8,192 frame average, filtered)
20
0.5
0.4
15
0.3
10
0.2
peak V
Amplitude (mV)
Amplitude (mV)
5
0
0.1
0
-0.1
-5
-0.2
-10
-0.3
-15
-0.4
-20
-0.5
0
1
2
3
4
5
6
7
8
Latency after click presentation (ms)
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
Latency after click presentation (ms)
11
Linear Averaging

Linear averaging - sample mean estimate
1
Aˆ 
N
N 1
 x[n]
n 0
1
E{ Aˆ}  E 
N
n 0
1
var{ Aˆ}  var 
N



N 1
1
N 1
1
 x[n]  N  E{x[n]}  N NA  A
n 0
 1
x[n]  2

n 0
 N
N 1
N 1
 var{ x[n]} 
n 0
1
2
2
N


N2
N
Linear averaging increases the amplitude SNR by a factor of N1/2
Cramer Rao lower bound on variance
var( Aˆ ) 
1
c
1
 , where c 
 const .
2
  ln p ( x[n]; A)  N
  ln p ( x[n]; A) 
 N E
E


A2
A2




2
12
Linear Averaging
Typical Fsp comparison for ABR recordings
with 65 dB stimulus vs no stimulus
24
No stimulus
65 dB stimulus
21
Fsp value
18
15
12
9
6
3
0
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
Frame number
Comparison of Fsp values with and without stimulus presentation
13
Wavelet Basics


Traditional Fourier Transform
Representation of signals in orthonormal
basis using complex exponentials (real
and imaginary sinusoidal components).
N 1
X [k ]   x[n]W Nkn ,
n 0



where W Nkn  e  j ( 2 / N ) kn
Signal represented in frequency domain
by a one-dimensional sequence.
“Loses” time information.
Features like transients, drifts, trends, etc.
may be lost upon reconstruction.


Wavelet Transform
Representation of signals in unconditional
orthonormal basis using waveforms of
limited durations with average value of
zero.
N 1
j
C ( j , k )   x[n] j ,k n, where  j ,k n  2 2  (2  j n  k ).

n 0



Makes no assumption about length or
periodicity of signals.
Contains time information in coefficients
Signal can be fully reconstructed using
inverse transform, and local time features
are preserved.
14
Wavelet Transform
• Discrete wavelet transform (DWT)
C  ,    C ( j , k )   f n g j ,k n,
n
for   2 j ,   k 2 j , j  N , k  Z
(α = scale coefficient, β=translation coefficient)
Signal x[n]
LP filter with H
 HP filter with G
Hx, Gx
HHx HGx
HHHx HHGx
….
Final level
15
Example Wavelet Filters
LP Decomposition filter H
HP Decomposition filter G
1
0.5
0.5
0
0
-0.5
-0.5
0 1 2 3 4 5 6 7 8 9 10
-1
LP Reconstruction filter H'
0.5
0.5
0
0
-0.5
0
1
2
3
4
5
6
7
1
2
3
4
5
6
7
8
HP Reconstruction filter G'
1
-0.5
0
8
-1
0 1 2 3 4 5 6 7 8 9 10
16
Wavelet Decomposition Example
17
Conventional Wavelet Denoising

Conventional denoising
1. Perform wavelet transform.
2. Set coefficients |C(α,β)|<δ to zero, δ – threshold value.
These coefficients are more likely to represent noise than
signal.
3. Perform inverse wavelet transform.

Characteristics of conventional denoising
•
•
•
Assumes that signal is smooth and coherent, noise rough
and incoherent.
Operation is performed on a single frame of data.
Non-linear operation – reduces the coefficients differently
depending on their amplitude.
18
Conventional Wavelet Denoising

Why does wavelet denoising work?
• The underlying signal is smooth and coherent, while the
noise is rough and incoherent
• A function f(t) is smooth if
n  N 
d n f (t )
is a continuous function.
dt n
• A function f(t) is smooth to a degree d, if
d n f (t )
 0  n  d 
is a continuous function.
dt n
• Bandlimited functions are smooth
• Measured biologic functions are smooth (such as ABR)
19
Conventional Wavelet Denoising

Coherent vs. incoherent
• A signal is coherent if its energy is concentrated in
both time and frequency domains.
• A reasonable measure of coherence is the
percentage of wavelet coefficients required to
represent 99% of signal energy.
• An example well-concentrated signal may require
5% of coefficients to represent 99% of its energy.
• Completely incoherent noise requires 99% of
coefficients to represent 99% of its energy.
20
Conventional Wavelet Denoising
21
Conventional Wavelet Denoising
-20 dB
Noisy sinewave
Simple low pass filter
20
20
Conventional denoising
20
0
0
0
-20
-20
-20
-10 dB
10
0
0.5
1
0
0 dB
0.5
1
+10 dB
5
0.5
1
0
2
0.5
1
1
0
2
0.5
1
0.5
1
5
0
0.5
1
2
0
0.5
1
0
0.5
1
0
0.5
1
0
-2
0
0.5
1
2
0
-2
0
1
-5
0
0
-2
0.5
0
-2
0.5
0
-10
0
0
0
10
0
-5
0
-2
+20 dB
1
0
-5
2
0.5
-10
0
0
2
0
0
-10
5
10
-2
0
0.5
1
22
Novel Wavelet Denoising

Conventional denoising applied to weak biosignals
• Setting coefficients |C(α,β)|< δ to zero, effectively removes all
the coefficients, including the ones that represent the signal.
• SNR must be large (>20dB).

Novel Wavelet Denoising
• Take advantage of multiple frames of data available.
• Create new frames through recombination and denoising.
• Apply a different δk for each new set of recombined frames.
Proprietary confidential information
23
Tree Denoising

Create a tree
7.
Collect a set of N frames of original data [f1, f2, …, fN]
Take the first two frames of the signal, f1 and f2, and average together, f12=
(f1+f2)/2
Denoise this average f12 using a threshold δk , fd12=den(f12 ,δ1).
Linearly average together two more frames of the signal, f34 ,and denoise
that average, fd34=den(f34 ,δ1). Continue this process for all N frames
Create a new level of frames consisting of [fd12, fd34, …, fdN-1,N].
Linearly average each two adjacent new frames to create f1234=(fd12 +fd34),
and denoise that average to create fd1234=den(f1234 ,δ2).
Continue to apply in a tree like fashion.
8.
Apply a different δk for denoising frames at each new level .
1.
2.
3.
4.
5.
6.
Proprietary confidential information
24
Tree Denoising Graph
1. Original signal x[n] consisting of N=8 frames of data of n signal samples in each frame
f1
f2
f3
f4
f5
f6
f7
f8
N
2. Create a signal x1[n] at level k=1 by averaging frames of x[n] (f 12=(f 1+f 2)/2 and denoising by δ1
fd12
fd 34
fd56
fd 78
N/2
3. Create a signal x2[n] at level k=2 by averaging frames of x1[n] and denoising by δ2
fd1234
fd 5678
N/4
4. Create a signal x3[n] at level k=3 by averaging frames of x2[n] and denoising by δ3
fd12345678
N/8
Proprietary confidential information
25
Cyclic Shift Tree Denoising (CSTD)
1. Original signal x[n] consisting of N=8 frames of data of n signal samples in each frame
f1
f2
f3
f4
f5
f6
f7
f8
N
2. Create a signal x1[n] at level k=1 by averaging frames of x[n] (f 12=(f 1+f 2)/2), then cyclic shifting
frames x[n] to create new cyclic shift averages (f23=(f 2+f3)/2), and denoising
fd12
fd 34
fd56
fd 78
fd 23
fd45
N/2
fd 67
fd 81
N/2
3. Create a signal x2[n] at level k=2 by averaging frames of x1[n], then cyclic shifting frames x1[n]
to create new cyclic shift averages, and denoising
fd1234
fd 5678
fd3456
N/4
fd 7812
fd 2345
N/4
fd6781
fd 4567
N/4
fd 8123
N/4
4. Create a signal x3[n] at level k=3 by averaging frames of x2[n], then cyclic shifting frames x2[n]
to create new cyclic shift averages, and denoising
fd12345678
N/8
fd 56781234
N/8
fd34567812
N/8
Proprietary confidential information
fd 78123456
N/8
fd 23456781
N/8
fd67812345
N/8
fd 45678123
N/8
fd 81234567
N/8
26
Cyclic Shift Tree Denoising (CSTD)
Original signal
Denoise with δ1
k=1
 Denoise with δ2
k=2
…
 Denoise with δ3
…
…
Final level
 Denoise with δk
Proprietary confidential information
27
Frame Permutations
- Create new arrangements of original frames prior to CSTD
- xnew=(p*xold) mod N
- Increase total number of new frames by a factor of
0.5*N*log2(N)
p
1
3
5
7
Frame 0 Frame 1
0
1
0
3
0
5
0
7
Proprietary confidential information
Frame 2
2
6
2
6
Frame 3 Frame 4 Frame 5 Frame 6 Frame 7
3
4
5
6
7
1
4
7
2
5
7
4
1
6
3
5
4
3
2
1
28
Threshold Selection
1

k and k

k 2 and 1 ,

k2

1
log (k) and
,

log (k)

k   k
1
e and e k ,

1
 2 k and
,
k

2

 k
1
k 

, 0
 .
cos( 4 K ) and
 k
4K 4
cos(
)

4K

 
 
Proprietary confidential information
29
Estimated Rate of Convergence

Linear averaging - sample mean estimate
1
Aˆ 
N




N 1
 x[n]
n 0
var( Aˆ ) 
1
c
1
 , where c 
 const .
2
  ln p ( x[n]; A)  N
  ln p ( x[n]; A) 
 N E
E


A2
A2




2
CSTD Creates M=log2(N)*N new frames.
Permutations prior to CSTD create at most M=0.5*(N2 * log2(N) new
frames.
CSTD can improve the Cramer-Rao lower bound by at most a factor of
0.5*N*log2(N).
The new frames are not linearly dependent, but also not all statistically
independent.
30
Experimental Results
Noisy Sinewaves
Linear Average and CSTD for Sinewave data at -20 dB
512 frames with  =1
1
6
Variance of Linear Avgerage and CSTD for Sinewave at -20 dB
512 frames with 1 =1
Linear Avg.
CSTD
Original
5
Linear Average
CSTD
4
Estimator Variance
Magnitude (with plotting offset)
0
10
3
2
1
-1
10
0
-1
-2
0
2
4
6
8
Time (ms)
Proprietary confidential information
10
12
10
2
10
3
10
Number of frames averaged
31
Experimental Results
ABR Data
Linear Avgerage and CSTD for ABR Data (Subject 3)
512 frames with 1 =1
0
3
Variance of Linear Avgerage and CSTD for ABR Data (Subject 3)
512 frames with 1 =1
10
Linear Average
CSTD
Final Average
Linear Average
CSTD
2.5
Estimator Variance
Magnitude (with plotting offset)
2
1.5
1
-1
10
0.5
0
-2
-0.5
10
0
2
4
6
Latency after click presentation (ms)
8
10
2
12
3
10
10
Number of frames averaged
Frames
2
4
8
16
32
64
128
256
512
σdenoised δav eraged
5.2457 10.1841
2.5953
3.9456
2.6032
3.0031
1.5141
1.8796
0.5419
0.9337
0.077
0.2533
0.0448
0.1143
0.0379
0.075
0.0139
0.0357
Ratio
1.94
1.52
1.15
1.24
1.72
3.29
2.55
1.98
2.57
32
Experimental Results
ABR Data
Linear Average and CSTD for ABR Data (Subject 3)
128 frames with  1 =1
Linear Average and CSTD for ABR Data (Subject 3)
512 frames with  1 =1
Linear Average and CSTD for ABR Data (Subject 3)
256 frames with  1 =1
7
3
5
Linear Avg.
CSTD
Final Avg.
Linear Avg.
CSTD
Final Avg.
Linear Avg.
CSTD
Final Avg.
6
2.5
4
1.5
1
0.5
Magnitude (with plotting offset)
Magnitude (with plotting offset)
Magnitude (with plotting offset)
5
2
4
3
2
3
2
1
1
0
0
0
-1
-0.5
0
2
4
6
8
Latency after click presentation (ms)
10
12
-1
0
2
4
6
8
Latency after click presentation (ms)
10
12
0
2
4
6
8
10
12
Latency after click presentation (ms)
33
Experimental Results
AMLR Data
mV
Pa
Na
Pb
Nb
Time (ms)
(a)
(b)
(c)
(d)
(e)
Performance of CSTD algorithm compared to linear averaging 256 data frames. (a): Template of AMLR
34
evoked potential waveform from Spehlmann; (b): linear average of 8192 AMLR frames; (c): Single frame
consisting of AMLR model plus WGN; (d): Linear average of 256 frames; (e): Result of CSTD algorithm
The Final Product
35
Future Work & Other applications


Wavelet denoising using wavelet packets
EEG/EP Recording and Monitoring
• Use in ambulances and emergency rooms
• At-home patient monitoring

Depth of Anesthesia Monitoring
• Monitor brain stem and cortex activity during surgery
• Use in all operating rooms

Oto-toxic drug administration
• Certain strong antibiotics cause hearing loss - ototoxic
• Dosage can be monitored on-line
• Use in intensive care units
36
ED Bedside
in minutes
Non-patient care
Environment-hours
37
38
HLB PRELIMINARY CONCEPT
39
HLB PRELIMINARY CONCEPT
40
Thank you!
 Questions?
41
Experimental Results
Noisy Sinewaves
Frames
2
4
8
16
32
64
128
256
512
SNR
Linear Denoised Variance
Linear
Denoise
dB
variance variance
ratio
SNR (dB) SNR (dB) Improvement
23.6363
2.0081
11.77
-16.74
-6.03
10.71
12.6286
0.9706
13.01
-14.02
-2.87
11.14
5.967
0.6821
8.75
-10.76
-1.34
9.42
2.967
0.3341
8.88
-7.72
1.76
9.48
1.6339
0.161
10.15
-5.13
4.93
10.06
0.8
0.0827
9.67
-2.03
7.82
9.86
0.4052
0.0479
8.46
0.92
10.2
9.28
0.1902
0.0271
7.02
4.21
12.68
8.47
0.1023
0.0184
5.56
6.9
14.35
7.45
Proprietary confidential information
42
Example Wavelet Filters
An additional property of a basis is being unconditional. A basis {φn} is an unconditional basis for a
normed space if there is some constant C<∞ such that

c  
n 0
n
n
n
C

c 
n 0
n
n
for coefficients cn, and any sequence {εn} of zeros and ones. This means that if some coefficients cn are set to zero by
the sequence {εn}, the norm of the remaining series is always bounded. Sines and cosines are NOT unconditional
bases.
43