Improving Evoked Potential Extraction
Michael Vinther, s973971 mv@logicnet.dk Ørsted, DTU November
8 Supervisor: Kaj-Åge Henneberg
## 1 Table of contents3.1
Constructing a Wiener filter 3.2 Epoch
optimization by Wiener filtering 4.3 Subspace method
and Wiener filtering combined Appendix 1 Correlation
estimates Appendix 2 Wiener filter
generator function Appendix 3 Functions for noise
reduction by PCA ## 2 IntroductionWhen obtaining EEG evoked potentials from scalp electrodes, background activity and other noise is added to the signal. The power of these unwanted components is typically greater than the power of the evoked potentials, so ways of reducing the background noise are necessary for analysis of the data. The background activity is caused by brain processes that may be unrelated to the experiment being performed, and which are uncorrelated with the evoked potentials. Other noise is typically EMG (movement) artifacts. The frequency content of the EMG artifacts lies in a higher band than the evoked potentials and can be removed with a lowpass filter. On the other hand, the spectrum of the background EEG overlaps that of the evoked potential, and conventional frequency-selective filters cannot be used. A simple and also widely used way of removing background activity and noise is by averaging many epochs of the same evoked potential. By averaging a number of recordings with the same stimulus, the other activities will cancel out and only evoked potentials remains because they are phase-locked. In general the signal to noise ratio is improved proportional to the square root of the number of epochs averaged. Conventional averaging typically requires more than 100 epochs, so methods requiring fewer epochs are sought after in order to save time and bother the patients/subjects at little as possible. Aunon et al. (1981) summarizes several approaches for improved averaging, including filtering and compensation for latency jitter. Cichocki et al. (2001) discuss the use Wiener filters and subspace methods. Here, I will first derive the computation of Wiener filters, and then give examples showing that applying optimized Wiener filters to each epoch before averaging can reduce the noise. Next subspace methods for noise reduction are described both computed via principal component analysis (PCA) and singular value decomposition (SVD). Finally, Wiener filtering and subspace methods are combined as proposed by Cichocki et al. (2001). Subspace methods require relative high signal-to-noise ratio (SNR). As evoked potential recordings often have SNR less than zero, it seems reasonable to initially reduce the noise with a Wiener filter. EEG recordings with responses to click sounds from a group of test subjects were used as example data. The stimulus was presented to each subject 120 times, and epochs where noise forced the A/D converter to saturation were removed. The sampling was done at 1 kHz. This data set was made available to me by Sidse Arnfred, dept. of psychiatry, Bispebjerg Hospital. ## 3 Wiener filteringA
Wiener filter is a digital filter, which is designed to minimize the mean
square difference between the filtered output and some desired signal. It is
sometimes called a minimum mean square error (MMSE) filter. This is illustrated
in Figure 1, where the difference is the signal In this report I will use matrix-vector
notation: Vectors are written in bold and matrices in bold capitals (
Figure 1: Model for noise reduction problem. This
kind of filter will reduce the noise, as the output ## 3.1 Constructing a Wiener filterApplying
the filter Where
|| From general optimization theory we have
that the solution to the problem ( ## 3.1.1 ExampleThe
following is a simple example of how The sequence of input values is longer than the output sequence. Often this is achieved by zero-padding when applying a FIR filter to some data sequence. ## 3.1.2 Filter delayAn
ordinary linear phase FIR filter introduces a delay equal to half the filter
order because the coefficients are symmetric, and the output sample corresponds
in time to the input sample multiplied by the middle filter coefficient. For a
Wiener filter, the delay introduced is decided by the time shift between In my implementation of this method, ## 3.1.3 Effective computationFinding
QR factorization of a matrix This is actually what happens in Matlab when solving an overdetermined system by writing `h = X\d`
as done in my implementation. Note that if the number of filter
coefficients is set too high, the rank of the matrix system may be
insufficient. This only happens if it is possible to construct a “perfect”
filter where the MSE is zero. In that case the number of filter coefficients
should be set equal to the rank of A Matlab implementation, which returns ## 3.1.4 ExampleA
simple test of the method is shown in Figure
2. Here, the source signal
Figure 2: Example use of Wiener filter to remove low frequency noise. The filter output in the bottom plot is almost a perfect reconstruction of the source. In
this example a good results is achieved with only four filter coefficients,
i.e. Excluding the first and last few samples,
MSE was computed to 2.3523·10
## 3.2 Epoch optimization by Wiener filteringAs
mentioned, some approximation to the true evoked potential is needed to
construct a Wiener filter for the EEG epochs. Cichocki et al. (2001) compute a
new filter for each epoch setting Figure 4 demonstrates the effect of the filtering. It shows both the conventional average and the average of the filtered epochs from a Wiener filter with 50 coefficients. It is not possible to compare the two averages with the true evoked potential to see which is best, as the true response is not available. But the conventional average contains much more high frequency noise, and EEG usually does not have that much high frequency contents, so it seems reasonable to assume that the average of the filtered epochs is the best.
Figure 4: Comparing conventional average and average of MMSE filtered epochs. The example is the evoked response to stimulus 1, Cz electrode on subject 102. 118 epochs were used. Epoch
averaging presupposes that the evoked potentials in all epochs have the same
latency. In this kind of experiments small variations can occur, and different
approaches for correcting latency jitter exists. MMSE optimized filters should
be able to compensate for time differences in an interval corresponding to the
filter length, because a filter introducing the necessary delay will result in
lower MSE when comparing the filter output to the average As the noise changes, the filters constructed for each epoch are quite different. The amplitude response of some of the filters used for the demonstration above is shown in Figure 5. They are obviously not just common band pass or band stop filters. Still the filters primarily remove frequency contents greater than 15 Hz, which is visible in the power specters in Figure 6. The histogram in Figure 7 indicates the filters ability to compensate for latency jitter: The median of the filter group delay computed in the frequency interval 0-50 Hz varies around zero.
Figure 5: Amplitude response of Wiener filters for four different epochs.
A better way of evaluating the method is to make two estimations of the evoked potential from smaller subsets of the epochs and compare these. In Figure 8 the conventional average and the average of MMSE filtered epochs is shown for two 20-epoch subsets from the 118 epochs used above. The EP might change a little over time, and to compensate for that I excluded the first 20 epochs and selected two sets of every second epoch starting from number 21 and number 22. This way any slow change in the responses over time will not influence the comparison of the subsets, nor will the subjects getting used to the sound in the first few minutes. The two conventional averages shown to the left are quite different, especially the peaks near sample 25 and sample 400. There are also differences in the filtered averages, but they are less obvious, and these two computed from only 20 epochs seems to be almost as good as the conventional average of all 118 epochs.
Figure 8: Comparing conventional average and average of MMSE filtered epochs. The example is the evoked response to stimulus 1, Cz electrode on subject 102. 20 epochs were used. Now,
to get an idea of how many epochs are necessary to achieve results comparable
to those obtained from conventional averaging of 118 epochs, I produced some
plots like the one in Figure
9. Here, the mean value of the squared difference
between normalized versions of the conventional average The
epochs used for computing i, _{}should be a better estimate than _{}, and MSEis almost constant._{i } This happens at about
## 4 Subspace methodsIn my first study of evoked potentials[3], spatial PCA was used to reduce the number of channels to analyze. A similar technique with the purpose of reducing noise in the epochs will be examined here and combined with Wiener filtering. ## 4.1 Temporal PCAIn data/channel reduction, further analysis is performed of some of the principal components instead of the recording. In this application, the data is reconstructed from a subset of the components. The signals obtained from the electrodes are treated independently, and instead of separation in space the different time epochs are used as input components for the PCA. It means that for 100 obtained epochs, the covariance matrix is 100 x 100 and the output is 100 components. If an eigenvalue in the analysis is large, many input components have contributed to the corresponding principal component. For 100 components, all principal components with a power contribution of less than 1% to the total power (eigenvalue less than 1 for normalized covariance matrix) can be considered as noise if the same signal is known to exists in all input components. The evoked potential should be present in all epochs, so actually there is only one source in a signal obtained from an electrode and only one principal component should be preserved. The Matlab function pcareduce in Appendix
3 can discard a number of principal components and then
reconstruct the epochs. The principal component ## 4.2 Alternative approach: SVDThe
same kind of noise reduction can be achieved using singular value
decomposition. This method is denoted the subspace method, because it divides
the data into signal and noise subspaces. The epochs should be arranged as
columns in a matrix
The columns of The signal part is extracted by computing
the projection of the observed signals where
the columns in Again, for evoked potentials, and
## 4.3 Subspace method and Wiener filtering combinedCichocki et al. (2001) propose combining the subspace method with Wiener filtering. The subspace method performs best for relative high SNR because this gives better separation between the signal and noise spaces, but EEG recordings typically have low SNR. As mentioned, a Wiener filter cannot remove noise in the frequency bands where desired signal is. But it can increase the SNR, and thus give a better starting point to compute the subspaces. The proposed method is to do singular value
decomposition of the filtered epochs – that is constructing I
did some tests with this method similar to the test described for the Wiener
filter in section 3.2, and they showed very little difference from just
averaging the filter output as in Figure
8. The output was not exactly the same, but I was not
able to find any improvement. The only thing the subspace method does when ## 5 ConclusionWiener filters both reduce background noise thereby improving the signal-to-noise ratio and compensate for latency jitter in the epochs. It was demonstrated that applying Wiener filters to the recordings can considerably improve the extracted evoked potentials, and in that way reduce the number of epochs required. In my tests the results achieved with only 20 epochs seemed as good as the results from conventional averaging of more than 100 unfiltered epochs. The subspace methods had little effect on the noise; both used alone and applied after Wiener filtering. An explanation could be that the noise was almost equally powerful in all epochs, which would mean that the weighting of the epochs was almost equal. In that case it is close to a conventional average. Another way of evaluating the methods might be to use artificial data simulating brain evoked potentials with background noise, where the SNR is known. That could give a quantitative estimate of the improvement in SNR from both the Wiener filtering and the subspace methods. ## 6 LiteratureAunon, Jorge I., McGillem, Clare D. and Childers, Donald G.
CRC Critical Reviews in Bioengineering July, 1981 Cichocki, Gharieb and Hoya
Lab. for advanced Brain Signal Processing Brain Science Institute, RIKEN, 2001 Madsen, Kaj and Nielsen, Hans Bruun
IMM, DTU, 2002-8-13 http://www.imm.dtu.dk/courses/02611/ Proakis, John G. and Manolakis, Dimitris G.
Third edition Prentice Hall 1996 Vinther, Michael
Ørsted, DTU, 2002-6-7 http://logicnet.dk/reports/ Widrow, Bernard and Stearns D. Samuel
Prentice Hall 1985 Appendix 1 Correlation estimates The autocorrelation matrix is[4] and the crosscorrelation is where
Appendix 2 Wiener filter generator function
` a = floor((n-1)/2);`
`x = x(:)';` `d = d(:);`
` error('x and d must have same length');`
`r = size(d,1)-n+1;`
` error('Too few known samples for specified filter order.');`
`X = zeros(r,n);`
` X(i,:) = x(i:i+n-1);`
`h = X\d(n-a:n-a-1+r);` ` `
Simple function to apply a FIR filter to a data sequence. The output vector is of same length as the input, and the filter is not mirrored, as it is when conv is used for filtering.
` a = floor((length(h)-1)/2);`
`x = conv(s,flipud(h(:)));` `x = x(1+a:a+length(s)); ` Appendix 3 Functions for noise reduction by PCA
[F,D,T] = pca(S);
error('n must be greater than
0') P = D./sum(D) pn = n; n = 1; p = P(1); n = n+1; p = p+P(n);
F(:,n+1: F = (T\F')';
R = cov(S); [U,D,V] = svd(R); D = diag(D)'; T = diag(sqrt(1./D))*V'; F = (T*S')'; |