Signal Processing for Medical Applications Frequency Domain - - PowerPoint PPT Presentation

Muthuraman Muthuraman

Christian-Albrechts-Universität zu Kiel Department of Neurology / Faculty of Engineering Digital Signal Processing and System Theory

Signal Processing for Medical Applications – Frequency Domain Analyses

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-2

Coherence

• The coherence analysis is an extensively used method to study the

correlations in frequency domain, between two simultaneously measured signals.

• Let and be two simultaneously recorded data sets of length .
• We estimate the short-time power spectra of , , and cross-spectrum,

which is the Fourier transform of the cross-correlation function of the signals and in each segment.

• Finally , we average the power spectra and the cross-spectrum across all the

segments and calculate the coherence as follows:

(4.1)

) (t x

) (t y

N

xx

s

yy

s

xy

s

) (t x

) (t y

) ( ) ( ) ( ) (

2

   

yy xx xy

s s s C 

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-3

Significance of Coherence

• The coherence spectrum represents the strength of correlation between two signals,

and .

• The confidence limit for coherence at the is given by

(4.2) = 0.99; and is the number of disjoint segements; hence the confidence limit is .

• If the signal length =10000; =1000; = ?;
• Frequency resolution If (i.e. the number of data points sampled per second)

is the sampling frequency, then the frequency resolution is .

• Thus, one should optimally choose the value of depending on the purpose of

analysis, to compromise between sensitivity and reliability. ) (t x

) (t y

 % 100

) 1 ( 1

) 1 ( 1



  

M L

C 



) 1 ( 1

01 . 1





M

M

N

D

L

C

s

f

D fs

D

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-4

Significance of Coherence

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-5

Dynamical Coherence

Lecture 4 – Different windows used for estimation

• The dynamical coherence analysis is done by estimating the coherence spectra

for a moving 30-second windows with an overlap of 28-seconds, resulting in an apparent time resolution of 2s.

• A model is created by coupling two AR2 processes . One

AR2 (V1) had narrow band characteristics and the other (V2) had broadband spectral characteristics .

• These two processes were simulated for a duration of 150 seconds at a sampling rate
• f 1,000 Hz. The narrow band AR2 was then band-pass filtered around its spectral

peak between 8 and 15 Hz and then combined by point-by-point summation with the broadband AR2 (V2) as follows: V=V2+0.2 V1.

• Independent white noises were added to V and V1 whose amounts were tuned so that
• verall coherence between V and V1 was around 0.1.

] [ ] 2 [ ] 1 [ ] [

2 1

n n y a n y a n y       ) 9753 . , 9691 . 1 (

2 1

   a a ) 36788 . , 37486 . (

2 1

   a a

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-6

Dynamical Coherence

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-7

Welch periodogram and Multitaper Method

• Estimation of coherence with these two methods which work on the

same principle Multitaper Method with single hanning taper Hanning window in Welch periodogram Method

Tapers Hanning Window

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-8

Increasing Time Resolution

• If we consider a finite length sample of a discrete time process , . Let us

assume a spectral representation for the process, (4.3)

• The Fourier transform of the data sequence is therefore given by

(4.4)

• The Welch periodogram method is capable of also analysing the signals using larger

time windows in the time domain which inturn gives a good estimation of the signal components.

• However, when shorter time windows need to be used, some disadvantages of this

method become evident when applied to non-linear signals.

• In this case for a stationary process, the spectrum is given by

(4.5)

) (t x

N t  , 2 , 1 

df ift f X t x ) 2 ( exp ) ( ) (

2 / 1 2 1









 



      

2 / 1 2 / 1 1

) ( ) , ( ) 2 ( exp ) ( ) ( ~ f d f X N f f K ift t x f x

N

 ) ( ~ f x

 

2

) ( ) ( f X E df f S 

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-9

Increasing Time Resolution

• A simple estimation of the spectrum (apart from the normalization constant) is
• btained by squaring the Fourier transform of the data sequence, i.e; .
• This suffers from two difficulties:
• Firstly, is not equal to , except when the data length is infinite, in which

case the kernel in equation (4.5) becomes a delta function. Rather it is related to by a convolution as given in equation (4.4).

• This problem is usually referred to as „bias“ corresponding to a mixing of information

from different frequencies of the underlying process due to a finite window length.

• Secondly, if the data are stochastic, then the squared Fourier transform of a time

series is an inconsitent estimator of the spectrum, because it does not converge to the „true“ spectrum when the data series tends to infinite length.

• Inorder to overcome all these disadvantages, the signals can be analysed with the

multitaper and the extended continous wavelet-transform method.

2

) ( ~ f x

) ( ~ f x ) ( f X ) ( f X

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-10

Multitaper Method

• If is the signal then the spectrum in this method is calculated by multiplying

the data with several orthogonal tapers (windows) (4.6) where is the Fourier transform

• f the data

(4.7) where are the orthogonal tapers.

• A particular choice of these taper functions, with optimal spectral concentration

properties, is given by the discrete prolate spheroidal sequences (DPSS). ) (t x

2 1

) ( ~ 1 ) (







K k k MT

X K S  

) ( ~ 

k

X         

  

dt t i t x X ) 2 exp( ) ( ) (   

) (t x





 

N t t t k

t i x k w X

1

) 2 exp( ) ( ) ( ~   

) , 2 , 1 )( ( K k k wt  

K

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-11

Multitaper Method

• Let be the DPSS of length and frequency bandwidth parameter
• Consider a sequence of length whose Fourier transform is given by

, we find the sequences so that the spectral amplitude is maximally concentrated in the interval , i.e. (4.8) is maximised.

• The maximisation problem leads to the matrix eigenvalue equation

(4.9)

• Eigen vector – Let be a square matrix, a non-zero vector is called a eigen

vector of if and only if there exists a number (real/complex) such that

) , , ( N W k wt

th

k

N

W

t

w

N



 

N t

t i w U

1

) 2 exp( ) (   

t

w

) ( U

 

W W, 

df f U

W W 2

) (





 

t t N t

w w t t t t W        

 



) ( ) ( 2 sin

C

A

A



. C AC  

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-12

Extended Continous Wavelet Transform Method

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-13

• If is the signal then the continuous wavelet transform can be written as

(4.10) where is the wavelet function and is the scaling factor.

• The wavelet used in this method is the Morlet wavelet which is defined as

(4.11)

• The relative bandwidth of this wavelet can be defined as

(4.12)

• By adjusting the ratio of which gives the flexibility in having a particular

frequency resolution at a particular frequency.

) (t x



        dt a t h t x a a CWTx  

*

) ( 1 ) , (

*

h

a

) 2 exp( ) exp( ) (

2

t jct t h   

c BWrel  2 2 

c 

Extended Continous Wavelet Transform Method

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-14

Model data

Lecture 4 – Different windows used for estimation

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-15

Partial Coherence

Lecture 4 – Different windows used for estimation

• Let and be three simultaneously measured signals of length .
• The partial coherence is estimating by first calculating the power spectra

and cross-spectra and in each of the disjoint windows. Finally we average these quantities across all the segments to get the estimate of the same.

• We estimate the partial coherence between the signals and as follows:

(4.13) where is a complex-valued function whose magnitude is called coherency between the two signals and .

• The confidence limit for the partial coherence at is .

) ( ), ( t y t x

) (t z

N

zz yy xx

S S S , ,

xz xy S

S ,

yz

S ) ( ), ( t y t x

) (t z

)) ( 1 ( )) ( 1 ( ) ( ) ( ) ( ) (

2

     

zy xz zy xz xy z xy

C C CY CY CY C    

) (

ij

CY

i

j

 % 100

) 2 ( 1

) 1 ( 1



 

M



Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-16

Partial Coherence

Lecture 4 – Different windows used for estimation

C2-F2/EMG C2-EMG/F2 F2-EMG/C2 F2-EMG C2-EMG C2-F2

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-17

Methods for finding direction of information flow

Lecture 4 – Different windows used for estimation Non-Parametric Methods Parametric Methods

(Based on modeling of system by linear VAR processes)

Partial Cross spectrum Granger Causality Index Partial Coherence Directed Transfer Function Partial Directed Coherence

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-18

Partial Directed Coherence

Lecture 4 – Different windows used for estimation

• Partial directed coherence (PDC) can be formulated as follows:

(4.14) (4.15) PDC: (4.16)

• Inorder to calculate PDC, we need to find the appropriate order ‚p‘ for the

underlying process.

• Optimum order can be found by Akaike Information Criterion (AIC).
• After finding the approprite ‚p‘, we need to find the coefficients for autoregressive

model which closely depicts the process.

) ( ) ( ) ( ) (

1

t r t x r a t x

p r

    





 

  

p r r i

e r a I A

1

) ( ) (









 k kj ij j i

A A

2

) ( ) ( ) (    

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-19

Partial Directed Coherence

Lecture 4 – Different windows used for estimation

• Application of PDC on model data with AR process of order p=5
• The five time series as the information flow between them as given in the figure

below:

X1 X2 X3 X4 X5

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-20

Partial Directed Coherence

Lecture 4 – Different windows used for estimation

X1 X2 X3 X4 X5

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-21

Lecture 4 – Different windows used for estimation

Granger Causality Index

• The principle of Granger causality states that if some series contains information

in past terms that helps in the prediction of series , then is said to cause .

• For predicting a value of using previous values of the series only, we get a

prediction error : (4.17)

• If we try to predict a value of using previous values of the series and

previous values of we get another predicton error : (4.18)

• If the variance of (after including series to the prediction) is lower than the

variance of we say that causes in the sense of Granger causality. ) (t y ) (t x

) (t y

) (t x ) (t x

p

X

e 



  

p i

t e j t X j A t X

1 ' 11

) ( ) ( ) ( ) (

) (t X

p X p

Y

1

e

) ( ) ( ) ( ) ( ) ( ) (

1 1 12 1 11

t e j t Y j A j t X j A t X

p j p j

    

 

 

1

e

Y

e

Y X

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-22

Lecture 4 – Different windows used for estimation

Granger Causality Index

• In the same way we can say that causes in the sense of Granger causality

when the variance of is reduced after including series in the prediction of series : (4.19)

• Granger causlity index is based directly on the defintion of causality, namely it

shows, if the information contributed by second channel improves the prediction of the first channel.

• The logarithm ratio of the residual variances for one and two-channel models is

computed: (4.20)

• This defintion can be extended to multichannel case by considering how the

inclusion of the given channels changes the residual variance ratios.

• GCI is an estimator in the time domain.

X Y

2

e

X Y

) ( ) ( ) ( ) ( ) ( ) (

2 1 21 1 22

t e j t Y j A j t X j A t Y

p j p j

    

 

 

) ( ln

1 2 1

e e GCI 



Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-23

Lecture 4 – Different windows used for estimation

Multivariate autoregressive Model

• Granger causality was defined for two channels, however, he later stated that the

causality principle holds only, if there are no other channels influencing the process.

• To account for the whole multivariate structure of a process of channels the

multichannel autoregressive model (MVAR) has to be considered. For the MVAR

• channel process :

(4.21) The model takes the form (4.22) where are vectors of size and the coefficients are - sized matrices

• Equation (4.22) can be easily transformed to describe relations in the frequency
• domain. After changing the sign of the and application of transform we get:

(4.23)

k k

) (t X

)) ( , ), ( ), ( ( ) (

2 1

t X t X t X t X

k

 





  

p i

t E j t X j A t X

1

) ( ) ( ) ( ) (

) (t E

k A k k 

A Z

) ( ) ( ) ( f X f A f E 

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-24

Lecture 4 – Different windows used for estimation

Multivariate autoregressive Model

• It can be derived as follows:

(4.24) (4.25)

• From the form of the above equations we can consider the model as a linear filter

with white noises on its input and the signals on its output. The matrix of filter coefficients is called the transfer matrix of the system.

• It contains information about all realtions between channels in the given set

including the phase relations between the signals.

) ( ) ( ) ( ) ( ) (

1

f E f H f E f A f X  

 1

) 2 ( exp ) ( ) (

 

           

p m

t f im m A f H 

) ( f E ) ( f X ) ( f H

Digital Signal Processing and System Theory| Signal Processing for Medical Applications | Introduction Slide I-25

Directed Transfer Function

Lecture 4 – Different windows used for estimation

• The formulation of the DTF is based on the properties of the transfer function
• f multivariate autoregressive process as:

(4.24)

• The DTF describes causal influence of channel on channel at frequency .
• The above equation defines a normalized version of DTF, which takes values from

0 to 1 producing a ration between the inflow of from channel to channel to all the inflows to channel .



 



k m im ij i j

f H f H f DTF

1 2 2 2

) ( ) ( ) (

j i f j i i