1

Lucas Parra, CCNY

Cit y College of New Yor k
BME I5100: Biomedical Signal Processing
Stochastic Processes
https://blackboard-doorway.cuny.edu/
[email protected]
Some links from previous course: http://newton.bme.columbia.edu/~lparra/class/
Lucas C. Parra
Biomedical Engineering Department
City College of New York
CCNY
2
Lucas Parra, CCNY
Cit y College of New Yor k
Schedule
Week 1: Introduction
Linear, stationary, normal - the stuff biology is not made of.
Week 1-5: Linear systems
Impulse response
Moving Average and Auto Regressive filters
Convolution
Discrete Fourier transform and z-transform
Sampling
Week 6-7: Analog signal processing
Operational amplifier
Analog filtering
Week 8-11: Random variables and stochastic processes
Random variables
Moments and Cumulants
Multivariate distributions, Principal Components
Stochastic processes, linear prediction, AR modeling
Week 12-14: Examples of biomedical signal processing
Harmonic analysis - estimation circadian rhythm and speech
Linear discrimination - detection of evoked responses in EEG/MEG
Hidden Markov Models and Kalman Filter- identification and filtering
3
Lucas Parra, CCNY
Cit y College of New Yor k
A stochastic process is a random quantity that changes over
time. More specifically it is a set of random variables sorted
in time X[t]. Sampling a random process repeatedly gives and
ensemble:
Examples: Spike train, position of a particle moving randomly,
number of bacteria in a culture, thermal noise in an amplifier,
etc.
Stochastic Processes
e
1
e
3
e
2
S
x
1
[t]
x
2
[t]
x
3
[t]
4
Lucas Parra, CCNY
Cit y College of New Yor k
Stochastic Process - Poisson Process
Pr ( N (t +At )N (t )=1)=+At
Pr ( N (t )=k )=
(+t )
k
e
+t
k !
Pr ( N (t )=0)=e
+t
Poisson process arises when we count the number of occurrences
N(t) of independent random events where the likelihood of
occurrence of an event is proportional to the time interval t.
One can proof by induction that the count follows a Poisson
distribution with =t
Interesting is the first induction step. For k=0 using definition and
limes t 0 it is easy to see that
Which means that the ISI is exponentially distributed.
5
Lucas Parra, CCNY
Cit y College of New Yor k
Stochastic Process - Shot Noise
X (t )=

i=1
N
h(t )6(t t
i
)=

i=1
N
h(t t
i
)
Shot noise we call Poisson process with events at times t
i

convolved with a IR.

The expected value of X(t) is
Event or spike rate can be estimated by convolving a spike train
with a sufficiently long IR. This in effect is what neurons do,
where the IR is given by the post-synaptic potential (PSP).
Sometimes (t
i
-t
i-1
)
-1
is considered "instantaneous" spike rate (t).
E| X (t )=+(t )

h(u) du
6
Lucas Parra, CCNY
Cit y College of New Yor k
A stochastic process x(t) is said to be a Gaussian process if
x=[x(t
1
),x(t
2
),...,x(t
n
)]
T
are jointly Gaussian distributed
In general the covariance matrix and mean depend on the
time instances t
1
, t
2
, ... , t
n
.
Classic example is Brownian motion with
Stochastic Process - Gaussian Process
p( x)=
1
.
(2n)
n
2
exp
|

1
2
( xj)
T
2
1
( xj)

r
ij
=Bmin(t
i
, t
j
)
7
Lucas Parra, CCNY
Cit y College of New Yor k
If p(x(t
1
),x(t
2
),...,x(t
n
)) does not depend on times t
1
, t
2
, ... , t
n
the
process is called stationary.
Any stochastic process for which the mean and covariance are
independent of time is called Wide Sense Stationary (WSS).
A WSS process is called ergodic if the sample average
converges to the ensemble average. With samples x[k] = x(t
k
)
this reads:
Stochastic Process - Wide Sense Stationary
E| f ( x(t ))=lim
N -
1
N

k=1
N
f ( x| k )
8
Lucas Parra, CCNY
Cit y College of New Yor k
If we have equal distant samples, x[k] = x(kT), the second
moments correspond to the auto-correlation
For WSS stochastic processes:
The mean is independent of time
The Auto-correlation
depends only on time differences
it is symmetric,
the zero lag is the power
upper bounded by the power
Stochastic Process -Auto-correlation
r
x
| k , l =E| x| k x

|l
r
x
| k , l =r
x
| kl
m
x
=E| x| k
r
x
| k =r
x

|k
r
x
|0r
x
| k
r
x
|0=E|x| k
2

r
x
| k =E| x| n+k x

| n
9
Lucas Parra, CCNY
Cit y College of New Yor k
Assuming ergodic WSS a direct estimate of the auto-
correlation is the (unbiased) sample auto-correlation
Note that correlation is convolution with opposite sign. It can
be computed fast with the FFT. Use xcorr()
Stochastic Process - Auto-correlation
r
x
| k =
1
Nk

l =1
Nk
x|l +k x

|l
x*
x
k N
10
Lucas Parra, CCNY
Cit y College of New Yor k
Examples
Stochastic Process - Auto-correlation
11
Lucas Parra, CCNY
Cit y College of New Yor k
For WSS the auto-correlation matrix becomes Hermitian
Toeplitz
Important WSS process is (Gaussian) uncorrelated noise
>> randn()
Stochastic Process -Auto-correlation
R
xx
=
|
r
x
(0) r
x

(1) r
x

(2) r
x

( p)
r
x
(1) r
x
(0) r
x

(1) r
x

( p1)
r
x
(2) r
x
(1) r
x
(0) r
x

( p2)

r
x
( p) r
x
( p+1) r
x
( p+2) r
x
(0)

R
xx
=c
2
I m
x
=0
12
Lucas Parra, CCNY
Cit y College of New Yor k
Oscillations of the correlation are best analyzed in the
frequency domain, which leads to the Power Spectrum
One can show that P
x
(e
j
) is real, even and positive.
The auto-correlation can be recovered with the inverse
Fourier transform
Stochastic Process - Power spectrum
P
x
(e
j o
)=

k=

r
x
| k e
jk o
r
x
| k =
1
2n

n
n
d oP
x
(e
j o
)e
jk o
13
Lucas Parra, CCNY
Cit y College of New Yor k
In particular, the total power is given by
the power spectrum is sometimes called spectral density
because it is positive and the signal power can always be
normalized to r(0) =(2)
-1
.
Example: Uncorrelated noise has a constant power spectrum
Hence it is also called white noise.
Stochastic Process - Power spectrum
r | k =E| x| n+k x

| n=c
2
6( k )
P
x
(e
j o
)=

k=

c
2
6( k ) e
jk o
=c
2
r
x
|0=E|x| k
2
=
1
2n

n
n
d oP
x
(e
j o
)
14
Lucas Parra, CCNY
Cit y College of New Yor k
The power spectrum captures the spectral content of the
random sequence. It can be estimate directly from the Fourier
transform of the data:
This estimate improves in the mean with increasing N.
Remember that when estimating auto-correlation and power
spectrum we implicitly assume WSS processes!
Stochastic Process - Power spectrum

P
x
(e
j o
)=
1
N
X (e
j o
)
2
X (e
j o
)=

k=0
N1
x| k e
j ok
P
x
(e
j o
)=lim
N -
E|

P
x
(e
j o
)
15
Lucas Parra, CCNY
Cit y College of New Yor k
Unfortunately, the direct estimate is inconsistent, i.e. its
variance does not converge to 0 for increasing N.
A classic heuristic, called the periodogram, is to smooth
neighboring frequencies: Compute Fourier transform for
window of size N/K and average over K windows.
Stochastic Process - Power spectrum
16
Lucas Parra, CCNY
Cit y College of New Yor k
Examples: >> psd(x);
Note 1/f spectrum, noise for small N.
Stochastic Process - Periodogram
17
Lucas Parra, CCNY
Cit y College of New Yor k
The effect of filtering on the statistics of a stochastic process:
Filtering corresponds to convolution of the autocorrelation
Stochastic Process - Filtering
y| n=

k=

h| k x| nk
r
y
| k =

l =

m=

h|l r
x
| ml +k h

(m)
=r
x
| k h| k h

|k
y[n]
x[n]

h[k]
r
x
[n]
h[k]
h
*
[-k]
r
y
[n]
r
yx
[k]

18
Lucas Parra, CCNY
Cit y College of New Yor k
With the convolution theorem of the Fourier transform we find
for the power spectrum
Or in the z-domain
Stochastic Process - Filtering
P
y
(e
j o
)=P
x
(e
j o
) H (e
j o
) H

(e
j o
)=P
x
(e
j o
)H (e
j o
)
2
P
yx
(e
j
)
H(e
j
) H
*
(e
j
)
P
y
(e
j
)
P
x
(e
j
)
|H(e
j
)|
2
P
x
(e
j
) P
x
(e
j
)
|H(e
j
)|
2

P
y
( z)=P
x
( z) H ( z) H

(1/ z

)
19
Lucas Parra, CCNY
Cit y College of New Yor k
Filtering shapes the power spectrum
Stochastic Process - Filtering
abs(fft(h,N)).^2
20
Lucas Parra, CCNY
Cit y College of New Yor k
"The power spectrum gives the spectral content of the data."
To see that consider the power of a signal after filtering with a
narrow bandpass filter around
0
.
Stochastic Process - Spectral Content
E|y| n
2
=
1
2n

n
n
d oH (e
j o
)
2
P
x
(e
j o
)
=
1
2n

o
0
Ao/ 2
o
0
+Ao/ 2
d oP
x
(e
j o
)

Ao
2n
P
x
(e
j o
0
)

0

H(e
j
)
21
Lucas Parra, CCNY
Cit y College of New Yor k
Colored noise process x[n] is filtered white noise v[n]
If h[n] is causal and invertible, x[n] is called a regular process
and v[n] the corresponding innovation process with power
v
2
.
Its power spectrum is continuous and given by
The inverse filter 1/H(z) is called the whitening filter of x[n].
Stochastic Process - Regular process
x[n]
v[n]

h[n]
P
x
(e
j
)=
v
2
|H(e
j
)|
2

P
v
(e
j
) =
v
2

|H(e
j
)|
2

v[n] x[n]
1/H(z)
P
x
(e
j o
)=c
v
2
H (e
j o
)
2
22
Lucas Parra, CCNY
Cit y College of New Yor k
A process x[n] is called predictive or singular if it can be
predicted exactly from a linear combination of its past:
Examples are sinusoids or sums thereof:

The spectrum of a predictive process has only discrete lines:
Stochastic Process - Predictive process
x| n=

k=1

a| k x| nk
P
x
(e
j o
)=

k=1

o
k
6(oo
k
)
x| n=

k=1

A
k
e
i no
k
+
k
23
Lucas Parra, CCNY
Cit y College of New Yor k
The Wold decomposition theorem tells us that any WSS
process x[n] can be decomposed into a sum of a regular
process x
r
[n] and predictive process x
p
[n]:

where the processes are orthogonal:
Hence the spectrum of any WSS is composed of a continuous
component with invertible H(z) plus discrete lines:
Stochastic Process - Wold decomposition
x| n=x
p
| n+x
r
| n
E| x
p
| n x
r

| m=0
P
x
(e
j o
)=H (e
j o
)
2
+

k=1

o
k
6(oo
k
)
24
Lucas Parra, CCNY
Cit y College of New Yor k
Note that for orthogonal processes the power spectrum is
additive.
If
with
Then
Stochastic Process - Wold decomposition
z | n=x| n+y| n
E| x| n y

| m=0
P
z
(e
j o
)=P
x
(e
j o
)+P
y
(e
j o
)
25
Lucas Parra, CCNY
Cit y College of New Yor k
Spectrogram: Power spectrum with strong spectral
component can be estimate on short sequences, and hence,
followed as it develops over time.
Example: Speech >> specgram(x);
Assumption: WSS within each small window
Stochastic Process -Spectrogram
Note harmonic
components
26
Lucas Parra, CCNY
Cit y College of New Yor k
An important regular process is the ARMA process: white
noise v[n] filtered with a rational transfer function
Example:
Stochastic Process - ARMA process
H ( z)=
B( z)
A( z)
x[n]
a[k]
v[n]

h[k]
b[k]
H ( z)=
1+0.9025 z
2
10.5562 z
1
+0.81 z
2
27
Lucas Parra, CCNY
Cit y College of New Yor k
ARMA modeling means to find a[k], and b[k] that approximate
an observed P
x
(e
j
) for a given x[n]
Note that now the input is not given. Instead only statistical
assumptions on the input are know (white WSS process).
Stochastic Process - ARMA modeling
x[n]

a[k]
v[n]

h[k]

b[k]
P
x
(e
j o
)=
B(e
j o
)
2
A(e
j o
)
2
x| n=

l =1
P
a|l x| nl +

l =0
Q
b|l v| nl
28
Lucas Parra, CCNY
Cit y College of New Yor k
We derive now statistical conditions for a[k] and b[k]. Multiply
the difference equation by x
*
[n -k] and take the expectation:
with cross-correlation:
Using
we obtain (non-linear) conditions for a[k] and b[k]
Stochastic Process - ARMA modeling

l =0
P
a|l r
x
| kl =

l =0
Q
b|l r
vx
| kl
r
vx
| kl =E| v| nl x

| nk
x| n=h| nv| n
E| v| n v

| m=6(nm)

l =0
P
a|l r
x
| kl =

l =0
Q
b|l h

|l k
29
Lucas Parra, CCNY
Cit y College of New Yor k
These conditions are called the Yule-Walker equations

In matrix form and using Hermitian symmetry of the auto-
correlation we can write for lags 0...L:
Since c[k] depend on h[k] these equations are non-linear in a[k]
and b[k]. They simplify for Q=0, i.e. AR model only.
Stochastic Process - ARMA modeling

l =0
P
a|l r
x
| kl =

l =0
Qk
b|l +k h

|l c| k
|
r
x
|0 r
x

|1 r
x

| 2 r
x

| P
r
x
|1 r
x
|0 r
x

|1 r
x

| P1

r
x
|Q r
x
|Q1 r
x
|Q2 r
x
|QP
r
x
|Q+1 r
x
|Q r
x
|Q1 r
x
|QP+1

r
x
| L r
x
| L1 r
x
| L2 r
x
| LP

|
1
a|1
a| 2

a| P

=
|
c|0
c|1

c|Q
0

0

30
Lucas Parra, CCNY
Cit y College of New Yor k
For Q=0 the Yule-Walker equations simplify to

In matrix form and using Hermitian symmetry of the auto-
correlation we obtain the AR normal equations:
Which are a set of linear equations for the unknown a[k].
Stochastic Process - AR modeling

l =0
P
a|l r
x
| kl =b|0
2
6(k )
|
r
x
|0 r
x

|1 r
x

| 2 r
x

| P
r
x
|1 r
x
|0 r
x

|1 r
x

| P1
r
x
| 2 r
x
|1 r
x
|0 r
x

| P2

r
x
| P r
x
| P1 r
x
| P2 r
x
|0

|
1
a|1
a| 2

a| P

=
|
b|0
2
0
0

0

31
Lucas Parra, CCNY
Cit y College of New Yor k
The AR normal equations can be rewritten as:
Linear equations with Toeplitz Hermitian form can be solved
efficiently with the Levinson-Durbin recursion.
Computes the auto-correlations and solves this linear equation:
>> a = lpc(x,P);
Stochastic Process - AR modeling
|
r
x
|0 r
x

|1 r
x

| P1
r
x
|1 r
x
|0 r
x

| P2

r
x
| P1 r
x
| P2 r
x
|0

|
a|1
a| 2

a| P

=
|
r
x
|1
r
x
| 2

r
x
| P

R
xx
a=r
x
32
Lucas Parra, CCNY
Cit y College of New Yor k
The goal of linear prediction is to predict sample x[n] from its
past P samples:
And a[k] are called the linear prediction coefficients (LPC).
The prediction error is given by
The LPC that minimize the expected error
satisfy the same linear equations as the AR model parameters:
which we obtain by taking derivatives with respect to a
*
[k].
Stochastic Process - Linear prediction
e| n=x| n x| n
a=argmin E|e| n
2

l =0
P
a|l r
x
| kl =c6(k )
x| n=

k=1
P
a| k x| nk
33
Lucas Parra, CCNY
Cit y College of New Yor k
In matlab all this is done with:
>> a = lpc(x,P);
The error process e[n] represents the "new" or "innovative"
part over the linear prediction, hence the name innovation
process.
One can show that the innovation in linear prediction is white
and uncorrelated to the prediction,
Stochastic Process - LPC and Wold Decomp.
E| e| n e

| m=6(nm)
E| e| n x

| m=0
34
Lucas Parra, CCNY
Cit y College of New Yor k
Assignment 9:

Pick either speech or heart rate

>> load hrs.mat; x=hrs;
>> [x,fs]=wavread('speech.wav');x=x(5501:6012);

Generate an AR model by solving normal equations.

Plot an estimated power spectrum together with the AR

model spectrum.

Show the innovation process and its power spectrum.

Compute the LPC coefficient using lpc(x,P)

Compare your AR results with LPC.

Show the linear prediction estimation error.

Compute the signal to noise ratio (SNR) of your LPC model:

Show SNR in dB as a function of model order P.

Stochastic Process - AR and LPC
SNR=

n
x| n
2

n
e| n
2