Computers in Biology and Medicine 56 (2015) 13–19
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Filtering multifocal VEP signals using Prony’s method A. Fernández a, L. de Santiago a,n, R. Blanco b, C. Pérez-Rico b, J.M. Rodríguez-Ascariz a, R. Barea a, J.M. Miguel-Jiménez a, J.R. García-Luque a, M. Ortiz del Castillo a, E.M. Sánchez-Morla c, L. Boquete a a
Department of Electronics, University of Alcalá, Plaza de S. Diego, s/n, 28801 Alcalá de Henares, Spain Department of Surgery, University of Alcalá, Plaza de S. Diego, s/n, 28801 Alcalá de Henares, Spain c Department of Psychiatry, University Hospital of Guadalajara, Guadalajara, Spain b
art ic l e i nf o
a b s t r a c t
Article history: Received 8 September 2014 Accepted 25 October 2014
Background: This paper describes use of Prony’s method as a filter applied to multifocal visual-evokedpotential (mfVEP) signals. Prony’s method can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. Method: By selecting Prony method parameters, a frequency filter has been developed which improves signal-to-noise ratio (SNR). Three different criteria were applied to data recorded from control subjects to produce three separate datasets: unfiltered raw data, data filtered using the traditional method (fast Fourier transform: FFT), and data filtered using Prony’s method. Results: Filtering using Prony’s method improved the signals’ original SNR by 44.52%, while the FFT filter improved the SNR by 33.56%. The extent to which signal can be separated from noise was analysed using receiver–operating–characteristic (ROC) curves. The area under the curve (AUC) was greater in the signals filtered using Prony’s method than in the original signals or in those filtered using the FFT. Conclusion: filtering using Prony’s method improves the quality of mfVEP signal pre-processing when compared with the original signals, or with those filtered using the FFT. & 2014 Elsevier Ltd. All rights reserved.
Keywords: mfVEP Prony’s method Signal-to-noise ratio ROC curve Biosignal processing
1. Introduction Since development of the multifocal visual-evoked-potential (mfVEP) recording technique [1], it has been possible to obtain signals that reflect electrical activity in the primary visual cortex in response to stimuli in various sectors of the visual field. Results obtained by analysing mfVEP signals show that it is a valid tool for assessment of visual function in patients with glaucoma [2], amblyopia [3], optic neuritis, multiple sclerosis [4,5] and other pathologies that affect the visual pathway. As is the case with other electrophysiological recordings, mfVEP signals are contaminated by noise. In contrast-reversing stimulus, mfVEP signals’ amplitude is in the order of nanovolts (although in the case of pattern-pulse stimulus [6] it can be in the order of microvolts), and may even be essentially zero in some locations in control subjects due to local folding of the cortex [7]. Sources of contamination may include the EEG’s alpha waves, other physiological signals emitted by the patient (muscle-tension-associated
n
Corresponding author. Tel.: þ 34 91 8856572. E-mail addresses: [email protected], [email protected] (L. de Santiago). http://dx.doi.org/10.1016/j.compbiomed.2014.10.023 0010-4825/& 2014 Elsevier Ltd. All rights reserved.
noise[2], ECG, EOG, etc.), or interference from the power grid or from other electronic equipment. Very few papers investigate the influence of mfVEP signal preprocessing techniques. In some cases, the signals only undergo the band-pass filtering set by the amplifiers: cut-off frequencies of 3–100 Hz [5,8,9], 1–100 Hz [10,11] and 1–20 Hz [12,13]. In other papers after the hardware filter a digital filter based on fast Fourier transform is applied, with different cut-off frequencies depending of the work: 1–30 Hz [14], 3–30 Hz [15,16] and 3–35 Hz [17,18]. Filtering is essential when performing patient diagnosis based on the results of signal classification. Prony’s method [19] can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. Prony’s method has been widely applied to electrophysiological signals [20–24]. It is a natural transformation for impulse responses since it uses damped sinusoids as a basis and therefore representation is efficient in terms of the number of coefficients required. Another advantage is that Prony modelling produces higher frequency resolution than simple FFT methods due to its reliance on autoregressive modelling [25]. It is remarkable that the transient nature of evoked potentials is more suited to modelling as a sum of damped sinusoids in time [26].
14
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
To the best of our knowledge, Prony’s method has not been used to filter mfVEP signals. The purpose of this paper is to evaluate the contribution that filtering in the Prony domain makes to improving the signal-to-noise ratio (SNR) of mfVEP signals obtained from healthy subjects.
2. Material and methods The study protocol was approved by the Institutional Review Boards of hospitals affiliated with the University of Alcalá and adhered to the tenets of the Declaration of Helsinki. All participants provided informed consent. The control group consisted of mfVEP signals obtained from 48 eyes in 24 healthy subjects ranging in age from 23 to 38 (mean age 30.20 77.55 years) – 10 males and 14 females – with normal neurological and ophthalmologic examination results. Two 7-min recordings were obtained from monocular stimulation of each eye using the VERIS 5.9 software (Electro-Diagnostic Imaging, San Mateo, USA). The authors employed the following procedure [27]: m-sequence visual stimulation (215–1 steps). Three channels of continuous VEP recordings were obtained with gold cup electrodes (impedance o2 kΩ), signal amplification (gain: 105, bandwidth 3–100 Hz), and a sample frequency of 1200 Hz. For the midline channel, the electrodes were placed 4 cm above the inion (active), at the inion (reference), and on the forehead (ground). For the other two channels, the same ground and reference electrodes were used, but the active electrodes were placed 1 cm above and 4 cm lateral to the inion on either side. Taking the difference between pairs of channels, three additional “derived” channels were obtained. All analyses were performed with custom made scripts written in Matlab R12 (The MathWorks, Natick, MA). The SNR was calculated using the noise window signal-to-noise ratio (nwSNR) method explained in detailed in [7]. Briefly, this method divides the trace X(k) (0–500 ms) into two different intervals—the signal window (45–150 ms), which contains the evoked potential response, and the noise window (325–430 ms), which essentially contains noise. The root mean square (RMS) values in the signal and noise windows are calculated as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi ms ∑150 X ðkÞ mean X ðkÞ45 150 ms k ¼ 45 ms RMSS ðXðkÞÞ ¼ ð1Þ N 45 150 ms N45–150 ms is the number of samples between 45 and 150 ms. To compute the RMSN value, the noise window’s time values are applied to an equation similar to (1). RMSN is the mean of the same channel of 60 sectors of the RMS taken from the noise window (325–430 ms). The signal-to-noise ratios in the signal window (SNRS) and in the noise window (SNRN) are respectively: SNRS ¼
RMSS RMSN
SNRN ¼
RMSN RMSN
ð2Þ
2.1. Filtering mfVEP signals using Prony’s method Prony’s method allows a set of sampled data to be decomposed into a linear combination of exponential functions of unknown amplitudes (Ak), damping factors (αk: s 1), angular velocity (ωk: radians), and phase angles (φk: radians) [19]. When considering N complex data samples y[1],…y[N], Prony wished to estimate y[n] by applying to the data the P-term complex exponential model given by (3). P
y½n ¼ ∑ Ak eðαk þ jωk Þðn 1ÞT S þ jφk k¼1
ð3Þ
In (3), n ¼1, 2,…,N, and TS is the sampling period. The classic Prony method solves (3) by using an auxiliary system of linear equations to decouple this non-linear equation. The outline of the algorithm is present in [28] and described in detail in [29]. Given that the length of the mfVEP signals is composed of N ¼600 samples, approximation by applying (3) obtains P ¼300 components. Therefore, the mfVEP signal is characterized by the parameters Ak, αk, ωk and φk, (k ¼1,…, P ¼300). As the principal information components of mfVEP signals are in low frequencies, only the lowest M frequencies are chosen for the filter output (XProny[n]), as represented in (4). Establishing the optimal value of M obtains a simpler representation, while reducing the noise is one of our major findings. The results section demonstrates the approach followed to establish the optimal value of M. M
X Prony ½n ¼ ∑ Ap eðαp þ jωp Þðn 1ÞT S þ jφp ; p¼1
ω1 oω2 o ⋯ o ωM
ð4Þ
Note that frequency (ωp/2π) corresponds with Prony’s frequency, but is not an exact match for Fourier’s frequency. For example, the frequency calculated with Prony’s method can have arbitrary values and the frequency range is different for each signal.
3. Results This paper compares the quality of the output of three different processing methods: (a) The signal resulting from averaging the two recordings (raw data). The raw data is 600 samples (or coefficients) long. (b) The signal resulting from filtering raw data (0–500 ms, 600 samples) using a digital low-pass filter. Although it is possible to use more advanced signal-filtering options, such as a Chebyshev type 1 filter [4,30] filtering using the FFT is the most commonly employed option when processing mfVEP signals [27]. Based on both our own studies and those of other authors [31], we concluded that the best SNRS value is obtained by applying a sharp cut-off [1–20 Hz] between using a discrete Fourier Transform. Under these conditions, the filtered signal (XFFT) is represented by 20 coefficients. After the signal is filtered, IFFT is applied to obtain a register in time domain again with 600 samples-500 ms. (c) The signal resulting from filtering raw data using Prony’s method. Thus, SNRS_Prony is the signal-to-noise ratio in the signal window of the trace filtered, and RMSERAW–Prony is the root mean squared error (RMSE) between the original recordings and the signal filtered using Prony’s method: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ∑150 X RAW ðkÞ X Prony ðkÞ k ¼ 45 RMSERAWProny ¼ V units ð5Þ N 45 150 ms The best M value (number of components, ordered by frequency, in the Prony domain) is defined as the M value that maximizes the SNRS_Prony/RMSERAW-Prony ratio. This value represents a compromise between the SNR value of the filtered signal and the distortion between the original signal and the filtered signal. According to this criterion, the best value for our database is M¼7 (28 coefficients in Prony domain). With these components, the signals are reconstructed into time domain and 600 samples500 ms length signals are obtained. Fig. 1 shows the mfVEP signal for channel 6, sector 2, right eye, of one of the subjects selected at random. It shows the signal without pre-processing (XRAW), filtered in the frequency domain (XFFT), and filtered using Prony’s method (XProny).
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Table 1 shows the principal results obtained after analysing the mfVEP signals. The upper part of Fig. 2 shows the mean value of the SNRS obtained for each of the 24 subjects (each point is the mean of 2eyes 60sectors 6channels SNR values). For all of the signals obtained (NT ¼24 2 60 6¼ 17,280 traces), the mean value of the SNRS before filtering is SNRS_RAW ¼2.75 70.73 (95% CI: 2.44– 3.10); when the FFT is used, the value is SNRS_FFT ¼3.69 70.9 (95% CI: 3.31–4.1); and if the Prony filter is used, the value is SNRS_Prony ¼4.04 71.00 (95% CI: 3.62–4.47). An extremely significant difference is observed between the two types of filtering (Paired t-test: P o 0.0001). If the SNR gains of both filters are defined as follows: GFFT
%
¼
ðSNRS_Prony SNRS_RAW Þ ðSNRS_FFT SNRS_RAW Þ 100 and GProny% ¼ 100 SNRS _RAW SNRS_RAW
then complete set of recordings GFFT% ¼ þ34.18% and GProny% ¼ þ 46.90%. A criterion widely used to filter out low-quality recordings is to set a threshold of SNRS o ¼0.23 log units [32], which is equivalent to SNRS ¼1.7 (dashed blue line in Fig. 2). In our database, 63.88% of the unfiltered recordings would be analysable. If FFT filtering were used, 77.11% would be analysable, and if Prony’s method were used, 82.78% of the sectors would be analysable. 3.1. Results for the best channel The method most widely used to process mfVEP signals consists of selecting from each sector the channel with the best signal-to-noise ratio in the signal window. In this case, the number 0.5
15
of recordings available in the database amounts to 24 2 60 traces. The lower part of Fig. 2 shows the mean value of the SNRS obtained for each of the 24 subjects, selecting for each sector the best channel (each point is the mean of 2eyes 60sectors SNR values). Among the unfiltered signals, the mean value is SNRS_RAW ¼4.29 7 1.26 (95% CI: 3.76–4.83); for the signals filtered using the FFT, it is SNRS_FFT ¼5.73 71.42 (95% CI: 5.12–6.33); and for the signals filtered using Prony’s method, it is SNRS_Prony ¼ 6.20 71.58 (95% CI: 5.53–6.86). For the complete set of recordings, considering the best channel, GFFT% ¼ þ 33.56% and GProny% ¼ þ44.52%. An extremely significant difference exists between the two types of filtering (Paired t-test p o0.0001). In our database, and considering the best channel for each sector, without filtering, 90.93% of the visual field sectors are analysable (SNRS 41.7); with FFT filtering, 96.77% are analysable; and with Prony’s method, 97.91% are analysable. One question that arises is whether Prony’s method is better than the FFT filter at distinguishing mfVEP signals from noise in normal subjects. To address this, the SNRS and SNRN values defined in (2) are used. Fig. 3 shows the frequency distribution of the values of SNRS (red) and SNRN (blue) for 2880 records (24 individuals 2 eyes 60 mfVEP recordings from the best channel). The noise and signal distributions overlap. Therefore, some of the mfVEP responses from the signal window cannot be distinguished from noise. The capacity of the various methods to distinguish between signal and noise has been evaluated using Receiver Operating Characteristic (ROC) analysis. An ROC curve shows sensitivity against the False Positive Rate (FPR¼1 specificity). For all of the thresholds available between min(SNRN) and max(SNRS), the sensitivity (fraction of cases in which SNRS 4Threshold) and specificity
Raw FFT Prony
0.4 0.3
Level uV
0.2 0.1 0 -0.1 -0.2 -0.3 -0.4
0
50
100
150
200
250
300
350
400
450
500
Time ms
Fig. 1. Example of a waveform of mfVEP signal along the time in order to compare the effects of each method. The black trace is a raw signal, the blue trace is a signal filtered with a FFT filter and the red trace represents a signal processed with Prony method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. SNRS values for the 24 subjects. Each point represents the average SNR for each patient with each method. The upper part represents the values for All the channels and the lower part just for the Best Channel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 1 Comparison of results between 3 signal types.
All channels Best channels
RAW data
FFT filter
Prony filter
SNRS ¼2.75 70.73 (rank: 1.39–4.05); analysable sectors: 63.88% SNRS ¼4.29 71.26 (rank: 1.98–6.32); analysable sectors: 90.93% AUCRAW ¼0.9909 7 0.0191 (rank: 0.9069–0.9999) thresholds: 0–16.66
SNRS ¼ 3.69 7 0.9 (rank: 1.46–5.06) analysable sectors: 77.11% SNRS ¼ 5.737 1.42 (rank: 2.15–7.68) analysable sectors: 96.77% AUCFFT ¼ 0.9915 70.0183 (rank: 0.9083–0.9999) thresholds: 0–18.19
SNRS ¼4.04 71.00 (rank: 1.57–5.84). Analysable sectors: 82.78% SNRS ¼6.20 7 1.58 (rank: 2.30–8.78) analysable sectors: 97.91% AUCProny ¼ 0.9952 7 0.0125 (rank: 0.9376–0.9999) thresholds: 0–20.32
16
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Fig. 3. Distribution of SNRS and SNRN. X-axis represents all the possible SNR values for each case and the Y-axis represents the frequency of each value. Red traces correspond to SNR values of the signal window and blue traces correspond to SNR values from noise window. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
(fraction of cases in which SNRN oThreshold) values are calculated according to (6) and (7): Sensitivity ¼
casesðSNRS 4 ThresholdÞ casesðSNRS 4 ThresholdÞ þ casesðSNRS o ThresholdÞ ð6Þ
Specificity ¼
casesðSNRN o ThresholdÞ casesðSNRN 4ThresholdÞ þ casesðSNRN o ThresholdÞ ð7Þ
Fig. 4(a) shows the ROC curves for all of the subjects before filtering; Fig. 4(b) shows them after FFT filtering; Fig. 4(c) shows them after filtering using Prony’s method; and Fig. 4(d) shows the mean value of the ROC curve for the three signal types. Calculating an AUC value for each subject obtains the following results: AUCRAW ¼0.9909 70.0191 (95%CI: 0.9837–0.999), AUCFFT ¼0.9915 70.0183 (95%CI: 0.9838–0.9993) and AUCProny ¼ 0.9952 7 0.0125 (95%CI: 0.9901–1.000), with the mfVEP-Prony showing the largest area. The difference between AUCRAW and AUCFFT was not statistically significant (P¼0.081, Wilcoxon matched-pairs signed rank test). However, there was a difference between AUCFFT and AUCProny (P¼0.001) and between AUCRAW and AUCProny (P¼ 0.000).
4. Discussion The principal conclusion of this paper is that filtering healthy subjects’ mfVEP signals using Prony’s method improves the SNRS to a greater extent than filtering them in the frequency domain. As a consequence, this obtains a greater number of analysable visualfield sectors, making this method better at distinguishing mfVEP signals from noise than the standard filter method. This paper demonstrates that M ¼7 components of the Prony filter are sufficient to achieve a good approximation of the original signals in healthy subjects. The RAW signal is represented by 600 coefficients; the signal filtered using the FFT is represented by 20 coefficients; and the
signal filtered using Prony’s method is represented by 28¼ 7 4 coefficients, which is slightly higher than the FFT representation. If the quality of the mfVEP signal is quantified by its SNRS (analysing the six channels available per sector), for the overall database use of Prony’s method achieves an improvement of GProny% ¼ þ46.90%, while the Fourier filter achieves an improvement of GFFT% ¼ þ34.18%. Therefore, the advantage of using Prony’s method over the Fourier filter can be quantified as þ12.72%. These results mean that by using Prony’s method the percentage of traces that exceed the minimum value of SNRS 41.7 is 82.78% of visual field sectors, while using the FFT filter produces a figure of 77.11%. If for each sector the channel with the best SNRS (best channel) is selected, use of Prony’s method achieves a gain of GProny% ¼ þ44.52%, while use of the Fourier filter achieves a gain of GFFT% ¼ þ33.56%. In this case, with Prony’s method the proportion of sectors with signals that exceed the SNRS 41.7 threshold stands at 97.91%, while with the FFT filter the proportion stands at 96.77%. As is to be expected, the SNRS is greater when the best channel is selected. However, the gain obtained by selecting the best channel is smaller with both types of filtering than when 6 channels from each sector are considered. This effect is undoubtedly due to the fact that the signals with the best SNR in each channel are close to their upper quality limit and, as a consequence, the percentage gain is smaller. Analysis based on ROC curves shows that Prony’s method reduced the overlap between signal and noise distribution. Based on AUC criteria, Prony’s method obtains the best results (AUCProny ¼0.9952 7 0.0125) when compared against values of AUCFFT ¼0.9915 70.0183 and AUCRAW ¼0.9909 70.0191. No significant difference exists between AUCRAW and AUCFFT (the traditional processing method). However, a difference does exist between AUCFFT and AUCProny and between AUCRAW and AUCProny. The authors believe that the moderate margins for improvement in signal quality and in the number of analysable sectors according to their SNR is due to the fact that they have been recorded in healthy subjects. Therefore, they are generally of good quality and, consequently, the margin for improvement is very small. The next step would be to confirm that applying Prony’s
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
RAW
FFT 0.9 Sensitivity
Sensitivity
0.9 0.8 0.7
0
0.1
0.2 0.3 1-Specificity
Sensitivity
0.8 0.7 0.6 0
0.1
0.2 0.3 1-Specificity
0
0.4
0.1
0.2 0.3 1-Specificity
0.4
AVERAGE
1
0.9 Sensitivity
0.7
0.5
0.4
PRONY
0.5
0.8
0.6
0.6 0.5
17
0.9 0.8 0.7 0.6
0
0.1 0.2 1-Specificity
0.3
Fig. 4. (a) ROC curve for unfiltered signals. (b) ROC curve for signals filtered using FFT. (c) ROC curve for signals filtered using Prony’s method. In plots (a), (b) and (c) each trace represents a patient. (d) Representation of mean value.
method to patients with various pathologies achieves a greater gain in the parameters that define signal quality. Several papers confirm that the SNR of this study’s database is appropriate. For instance, [27] obtain the following formula for the SNR according to subject age: SNRS ¼ 5.1 (0.01 years), valid for the best channels of each sector and applied to data filtered using the Fourier transform (cut-off at 35 Hz). In their database, SNRS ¼5.59 (mean age¼49.0), a similar value to the one obtained in our signals (SNRS ¼5.73) after filtering them with the FFT (cut-off at 20 Hz). In [33] research the advantages of a prediction method based on 40 virtual electrodes. Using a database of 15 healthy participants, the standard method (bandwidth of 1 Hz–2 kHz) obtains an SNRS of 3.877 1.08, while with the prediction method the SNRS is 4.36 71.03, representing a gain of þ12.66%. Using Prony’s method, and when analysing the best channels, we obtained an improvement in the SNRS (GProny% ¼ þ44.52%) greater than that obtained using the prediction method. ROC analysis has been used previously to evaluate the overlap between SNR distribution in the signal and noise windows. For instance, Meigen and Kramer i [34] calculate the AUC of a sample of 30 subjects, positioning the four electrodes (A, B, C, D) within a 4-cm radius of the inion. Our original data (AUCRAW ¼0.9909 7 0.0191) exceed the best value obtained by Meigen and Kramer (0.9647 0.005), making it difficult to establish a direct comparison. In their paper, and using a single channel, they obtain a best value of AUC ¼0.93370.005 for A–B. In our case (AUCProny ¼ 0.9952 7 0.0125), we exceed that value. This is undoubtedly due to the different position of the electrodes and to the fact that we selected the best channel of the six available, while Meigen and Kramer obtained all their data from a single channel.
Ishikawa et al. [35] study a population of healthy Japanese subjects, obtaining an AUC of 0.983 7 0.027 for 54 left eyes and one of 0.9837 0.017 for 56 right eyes. Various channels are combined after obtaining them from electrodes positioned 4 cm above, left and right of the inion and 2.5 cm below the inion (ground electrode on the left earlobe). This value is comparable to that obtained in this paper using the FFT (AUCFFT ¼ 0.99157 0.0183). Nakamura et al. [36] study patients with glaucoma, obtaining for 62 control subjects a mean SNRS value of 3.09. They also analyse ROC curves, obtaining an AUC value of 0.98 70.02. This result can be compared with the results achieved in this paper using FFT filtering (AUCFFT ¼ 0.99157 0.0183), although the cut-off frequencies are different (35 Hz and 20 Hz) in our study. Our results exceed those obtained in [18] when using a PCA method to filter the signals. By selecting four principal components, they obtain a value of SNR ¼5, which is lower than the one obtained in this paper using Prony filtering. Several lines of research merit further investigation to ascertain if it is possible to improve on the results achieved here. For example, in the Prony filter (Eq. (4)), it would be interesting to see if better results are obtained by setting a cut-off frequency and selecting components below that value rather than by selecting a fixed number of components. Other mfVEP signal-filtering options also merit analysis, e.g. wavelet methods, classic principle components/ single value decomposition and hierarchical decomposition. [30]. In conclusion, we consider that this mfVEP-signal filtering method, although limited by the small sample size and by the fact that the research was carried out on healthy subjects, has great potential for use in diagnosing and monitoring visual pathway pathologies in patients with demyelinating and ischaemic diseases.
18
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Conflict of interest statement The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Summary Due to development of the multifocal visual-evoked-potential (mfVEP) recording technique, it has been possible to obtain signals that reflect electrical activity in the primary visual cortex in response to stimuli in various sectors of the visual field. As is the case with other electrophysiological recordings, mfVEP signals are contaminated by noise. The amplitude of this type of signals is in the order of nanovolts–microvolts, and may even be essentially zero in some locations in control subjects due to local folding of the cortex. Sources of contamination may include the EEG’s alpha waves, other physiological signals emitted by the patient (muscletension-associated noise, ECG, EOG, etc.), or interference from the power grid or from other electronic equipment. Filtering is essential when performing patient diagnosis based on the results of signal classification. Very few papers investigate the influence of mfVEP signal filter techniques. This paper describes use of Prony’s method as a filter applied to multifocal visual-evokedpotential (mfVEP) signals. Prony’s method can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. To the best of our knowledge, Prony’s method has not been used to filter mfVEP signals. By selecting Prony method parameters, a frequency filter has been developed which improves signal-to-noise ratio (SNR). Three different criteria were applied to data recorded from subjects (48 eyes in 24 healthy subjects ranging in age from 23 to 38) to produce three separate datasets: unfiltered raw data, data filtered using the traditional method (fast Fourier transform: FFT), and data filtered using Prony’s method. Filtering using Prony’s method improved the signals’ original SNR by 44.52%, while the FFT filter improved the SNR by 33.56%. The extent to which signal can be separated from noise was analysed using receiver–operating– characteristic (ROC) curves. The area under the curve (AUC) was greater in the signals filtered using Prony’s method than in the original signals or in those filtered using the FFT. The principal conclusion of this paper is that filtering healthy subjects’ mfVEP signals using Prony’s method improves the SNRS to a greater extent than filtering them in the frequency domain. As a consequence, this obtains a greater number of analysable visualfield sectors, making this method better at distinguishing mfVEP signals from noise than the standard filter method. The authors believe that the moderate margins for improvement in signal quality and in the number of analysable sectors according to their SNR is due to the fact that they have been recorded in healthy subjects, so the next step would be to confirm that applying Prony’s method to patients with various pathologies achieves a greater gain in the parameters that define signal quality.
Acknowledgments This research has been supported by Spain’s Ministerio de Ciencia e Innovación under the “Advanced analysis of multifocal
ERG and visual evoked potentials applied to the diagnosis of optic neuropathies” project (ref. TEC2011-26066) and by grants FIS PI11/ 00533 and RETICS RD12/0034/0006. We wish to thank Dr. D.C. Hood for his generosity with the software used to analyze the mfVEP data.
References [1] H. a Baseler, E.E. Sutter, S. a Klein, T. Carney, The topography of visual evoked response properties across the visual field, Electroencephalogr. Clin. Neurophysiol. 90 (1) (1994) 65–81. [2] D.C. Hood, V.C. Greenstein, Multifocal VEP and ganglion cell damage: applications and limitations for the study of glaucoma,, Prog. Retin. Eye Res. 22 (2) (2003) 201–251. [3] W. Zhang, K. Zhao, Multifocal VEP difference between early- and late-onset strabismus amblyopia, Doc. Ophthalmol. 110 (2-3) (2005) 173–180. [4] R. Ruseckaite, T. Maddess, G. Danta, C.J. Lueck, A.C. James, Sparse multifocal stimuli for the detection of multiple sclerosis, Ann. Neurol. 57 (6) (2005) 904–913. [5] L.K. Grover, D.C. Hood, Q. Ghadiali, T.M. Grippo, A.S. Wenick, V.C. Greenstein, M.M. Behrens, J.G. Odel, A comparison of multifocal and conventional visual evoked potential techniques in patients with optic neuritis/multiple sclerosis, Doc. Ophthalmol. 117 (2) (2008) 121–128. [6] A.C. James, The pattern-pulse multifocal visual evoked potential,, Invest. Ophthalmol. Vis. Sci. 44 (2) (2003) 879–890. [7] X. Zhang, D.C. Hood, C.S. Chen, J.E. Hong, A signal-to-noise analysis of multifocal VEP responses: an objective definition for poor records, Doc. Ophthalmol. 104 (3) (2002) 287–302. [8] M. Bengtsson, S. Andreasson, G. Andersson, Multifocal visual evoked potentials —a method study of responses from small sectors of the visual field, Clin. Neurophysiol. 116 (8) (2005) 1975–1983. [9] E. Gutiérrez-Díaz, C. Pérez-Rico, M.J. de Atauri, E. Mencía-Gutiérrez, R. Blanco, Evaluation of the visual function in obstructive sleep apnea syndrome patients and normal–tension glaucoma by means of the multifocal visual evoked potentials, Graefes Arch. Clin. Exp. Ophthalmol. 250 (11) (2012) 1681–1688. [10] V. Parisi, L. Ziccardi, G. Stifano, L. Montrone, G. Gallinaro, B. Falsini, Impact of regional retinal responses on cortical visually evoked responses: multifocal ERGs and VEPs in the retinitis pigmentosa model, Clin. Neurophysiol. 121 (3) (2010) 380–385. [11] Y. Shimada, M. Horiguchi, A. Nakamura, Spatial and temporal properties of interocular timing differences in multifocal visual evoked potentials, Vision Res. 45 (3) (2005) 365–371. [12] A. Klistorner, H. Arvind, R. Garrick, S.L. Graham, M. Paine, C. Yiannikas, Interrelationship of optical coherence tomography and multifocal visualevoked potentials after optic neuritis, Invest. Ophthalmol. Vis. Sci. 51 (5) (2010) 2770–2777. [13] P. Sriram, A. Klistorner, H. Arvind, S.L. Graham, Reproducibility of multifocal VEP latency using different stimulus presentations, Doc. Ophthalmol. 125 (1) (2012) 43–49. [14] S. Park, S.H. Park, J.H. Chang, Y.-H. Ohn, Study for analysis of the multifocal visual evoked potential, Korean J. Ophthalmol. 25 (5) (2011) 334–340. [15] (a) G.C. Herbik, J. Hölzl, Reupsch, M.B. Hoffmann, Differential effects of optic media opacities on mfERGs and mfVEPs, Clin. Neurophysiol. 124 (6) (2013) 1225–1231. [16] C. Pieh, M.B. Hoffmann, M. Bach, The influence of defocus on multifocal visual evoked potentials, Graefes Arch. Clin. Exp. Ophthalmol. 243 (1) (2005) 38–42. [17] D.C. Hood, X. Zhang, J.E. Hong, C.S. Chen, Quantifying the benefits of additional channels of multifocal VEP recording, Doc. Ophthalmol. 104 (3) (2002) 303–320. [18] X. Zhang, D.C. Hood, Increasing the sensitivity of the multifocal visual evoked potential (mfVEP) technique: incorporating information from higher order kernels using a principal component analysis method, Doc. Ophthalmol. 108 (3) (2004) 211–222. [19] Prony, Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques, et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alkool, á différentes températures, J. l’école Polytech 1 (22) (1795) 24–76. [20] M. Hansson, T. Gänsler, G. Salomonsson, Estimation of single event-related potentials utilizing the Prony method, IEEE Trans. Biomed. Eng. 43 (10) (1996) 973–981. [21] V. Garoosi, B.H. Jansen, Development and evaluation of the piecewise Prony method for evoked potential analysis, IEEE Trans. Biomed. Eng. 47 (12) (2000) 1549–1554. [22] S.W. Chen, A two-stage discrimination of cardiac arrhythmias using a total least squares-based Prony modeling algorithm, IEEE Trans. Biomed. Eng. 47 (10) (2000) 1317–1327. [23] B.H. Jansen, A. Hegde, N.N. Boutros, Contribution of different EEG frequencies to auditory evoked potential abnormalities in schizophrenia, Clin. Neurophysiol. 115 (3) (2004) 523–533. [24] V. Sajith, A. Sukeshkumar, J. Jacob, Identification of diabetic neuropathic conditions of simulated ulnar nerve response using Prony’s method,, in: H.K. Kim, S.-I. Ao, M.A. Amouzegar, B.B. Rieger (Eds.), IAENG Transactions on
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
[25] [26]
[27] [28]
[29] [30]
[31]
Engineering Technologies, vol. 247, Springer, Dordrecht, the Netherlands, 2014, pp. 683–695. P. Lander, D. Jones, E. Berbari, R. Lazzara, Time-frequency structure of the high-resolution ECG,, J. Electrocardiol. 27 (207-12) (1994). T. Demiralp, a Ademoglu, Y. Istefanopulos, H.O. Gülçür, Analysis of eventrelated potentials (ERP) by damped sinusoids, Biol. Cybern. 78 (6) (1998) 487–493. B. Fortune, X. Zhang, D.C. Hood, S. Demirel, C. a Johnson, Normative ranges and specificity of the multifocal VEP, Doc. Ophthalmol. 109 (1) (2004) 87–100. T. Lobos, J. Rezmer, P. Janik, H. Amarís, M. Alonso, C. Álvarez, Application of wavelets and Prony method for disturbance detection in fixed speed wind farms,, Int. J. Electr. Power Energy Syst 31 (9) (2009) 429–436. S. Lawrence, Marple, Digital Spectral Analysis with Applications, Prentice-Hall, Inc, Englewood Cliffs, New Jersey, 1987. T. Maddess, A.C. James, R. Ruseckaite, E.A. Bowman, Hierarchical decomposition of dichoptic multifocal visual evoked potentials, Vis. Neurosci 23 (5) (2006) 703–712. J.A. Macfarlane, The Multifocal Visual Evoked Cortical Potential in Visual Field Mapping: A Methodological Study, University of Glasgow, 2008.
19
[32] D.C. Hood, X. Zhang, C. Rodarte, E.B. Yang, N. Ohri, B. Fortune, C. a Johnson, Determining abnormal interocular latencies of multifocal visual evoked potentials, Doc. Ophthalmol. 109 (2) (2004) 177–187. [33] B. a E. Mazinani, T.D. Waberski, A.W. a Weinberger, P. Walter, G.F. Roessler, Improving the quality of multifocal visual evoked potential results by calculating multiple virtual channels, Jpn. J. Ophthalmol. 55 (4) (2011) 396–400. [34] T. Meigen, M. Krämer, Optimizing electrode positions and analysis strategies for multifocal VEP recordings by ROC analysis, Vis. Res. 47 (11) (2007) 1445–1454. [35] K. Ishikawa, T. Nagai, Y. Yamada, A. Negi, M. Nakamura, Optimal conditions for multifocal VEP recording for normal Japanese population established by receiver operating characteristic analysis, Doc. Ophthalmol. 122 (1) (2011) 29–37. [36] M. Nakamura, K. Ishikawa, T. Nagai, A. Negi, Receiver–operating characteristic analysis of multifocal VEPs to diagnose and quantify glaucomatous functional damage, Doc. Ophthalmol. 123 (2) (2011) 93–108.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Filtering multifocal VEP signals using Prony’s method A. Fernández a, L. de Santiago a,n, R. Blanco b, C. Pérez-Rico b, J.M. Rodríguez-Ascariz a, R. Barea a, J.M. Miguel-Jiménez a, J.R. García-Luque a, M. Ortiz del Castillo a, E.M. Sánchez-Morla c, L. Boquete a a
Department of Electronics, University of Alcalá, Plaza de S. Diego, s/n, 28801 Alcalá de Henares, Spain Department of Surgery, University of Alcalá, Plaza de S. Diego, s/n, 28801 Alcalá de Henares, Spain c Department of Psychiatry, University Hospital of Guadalajara, Guadalajara, Spain b
art ic l e i nf o
a b s t r a c t
Article history: Received 8 September 2014 Accepted 25 October 2014
Background: This paper describes use of Prony’s method as a filter applied to multifocal visual-evokedpotential (mfVEP) signals. Prony’s method can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. Method: By selecting Prony method parameters, a frequency filter has been developed which improves signal-to-noise ratio (SNR). Three different criteria were applied to data recorded from control subjects to produce three separate datasets: unfiltered raw data, data filtered using the traditional method (fast Fourier transform: FFT), and data filtered using Prony’s method. Results: Filtering using Prony’s method improved the signals’ original SNR by 44.52%, while the FFT filter improved the SNR by 33.56%. The extent to which signal can be separated from noise was analysed using receiver–operating–characteristic (ROC) curves. The area under the curve (AUC) was greater in the signals filtered using Prony’s method than in the original signals or in those filtered using the FFT. Conclusion: filtering using Prony’s method improves the quality of mfVEP signal pre-processing when compared with the original signals, or with those filtered using the FFT. & 2014 Elsevier Ltd. All rights reserved.
Keywords: mfVEP Prony’s method Signal-to-noise ratio ROC curve Biosignal processing
1. Introduction Since development of the multifocal visual-evoked-potential (mfVEP) recording technique [1], it has been possible to obtain signals that reflect electrical activity in the primary visual cortex in response to stimuli in various sectors of the visual field. Results obtained by analysing mfVEP signals show that it is a valid tool for assessment of visual function in patients with glaucoma [2], amblyopia [3], optic neuritis, multiple sclerosis [4,5] and other pathologies that affect the visual pathway. As is the case with other electrophysiological recordings, mfVEP signals are contaminated by noise. In contrast-reversing stimulus, mfVEP signals’ amplitude is in the order of nanovolts (although in the case of pattern-pulse stimulus [6] it can be in the order of microvolts), and may even be essentially zero in some locations in control subjects due to local folding of the cortex [7]. Sources of contamination may include the EEG’s alpha waves, other physiological signals emitted by the patient (muscle-tension-associated
n
Corresponding author. Tel.: þ 34 91 8856572. E-mail addresses: [email protected], [email protected] (L. de Santiago). http://dx.doi.org/10.1016/j.compbiomed.2014.10.023 0010-4825/& 2014 Elsevier Ltd. All rights reserved.
noise[2], ECG, EOG, etc.), or interference from the power grid or from other electronic equipment. Very few papers investigate the influence of mfVEP signal preprocessing techniques. In some cases, the signals only undergo the band-pass filtering set by the amplifiers: cut-off frequencies of 3–100 Hz [5,8,9], 1–100 Hz [10,11] and 1–20 Hz [12,13]. In other papers after the hardware filter a digital filter based on fast Fourier transform is applied, with different cut-off frequencies depending of the work: 1–30 Hz [14], 3–30 Hz [15,16] and 3–35 Hz [17,18]. Filtering is essential when performing patient diagnosis based on the results of signal classification. Prony’s method [19] can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. Prony’s method has been widely applied to electrophysiological signals [20–24]. It is a natural transformation for impulse responses since it uses damped sinusoids as a basis and therefore representation is efficient in terms of the number of coefficients required. Another advantage is that Prony modelling produces higher frequency resolution than simple FFT methods due to its reliance on autoregressive modelling [25]. It is remarkable that the transient nature of evoked potentials is more suited to modelling as a sum of damped sinusoids in time [26].
14
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
To the best of our knowledge, Prony’s method has not been used to filter mfVEP signals. The purpose of this paper is to evaluate the contribution that filtering in the Prony domain makes to improving the signal-to-noise ratio (SNR) of mfVEP signals obtained from healthy subjects.
2. Material and methods The study protocol was approved by the Institutional Review Boards of hospitals affiliated with the University of Alcalá and adhered to the tenets of the Declaration of Helsinki. All participants provided informed consent. The control group consisted of mfVEP signals obtained from 48 eyes in 24 healthy subjects ranging in age from 23 to 38 (mean age 30.20 77.55 years) – 10 males and 14 females – with normal neurological and ophthalmologic examination results. Two 7-min recordings were obtained from monocular stimulation of each eye using the VERIS 5.9 software (Electro-Diagnostic Imaging, San Mateo, USA). The authors employed the following procedure [27]: m-sequence visual stimulation (215–1 steps). Three channels of continuous VEP recordings were obtained with gold cup electrodes (impedance o2 kΩ), signal amplification (gain: 105, bandwidth 3–100 Hz), and a sample frequency of 1200 Hz. For the midline channel, the electrodes were placed 4 cm above the inion (active), at the inion (reference), and on the forehead (ground). For the other two channels, the same ground and reference electrodes were used, but the active electrodes were placed 1 cm above and 4 cm lateral to the inion on either side. Taking the difference between pairs of channels, three additional “derived” channels were obtained. All analyses were performed with custom made scripts written in Matlab R12 (The MathWorks, Natick, MA). The SNR was calculated using the noise window signal-to-noise ratio (nwSNR) method explained in detailed in [7]. Briefly, this method divides the trace X(k) (0–500 ms) into two different intervals—the signal window (45–150 ms), which contains the evoked potential response, and the noise window (325–430 ms), which essentially contains noise. The root mean square (RMS) values in the signal and noise windows are calculated as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi ms ∑150 X ðkÞ mean X ðkÞ45 150 ms k ¼ 45 ms RMSS ðXðkÞÞ ¼ ð1Þ N 45 150 ms N45–150 ms is the number of samples between 45 and 150 ms. To compute the RMSN value, the noise window’s time values are applied to an equation similar to (1). RMSN is the mean of the same channel of 60 sectors of the RMS taken from the noise window (325–430 ms). The signal-to-noise ratios in the signal window (SNRS) and in the noise window (SNRN) are respectively: SNRS ¼
RMSS RMSN
SNRN ¼
RMSN RMSN
ð2Þ
2.1. Filtering mfVEP signals using Prony’s method Prony’s method allows a set of sampled data to be decomposed into a linear combination of exponential functions of unknown amplitudes (Ak), damping factors (αk: s 1), angular velocity (ωk: radians), and phase angles (φk: radians) [19]. When considering N complex data samples y[1],…y[N], Prony wished to estimate y[n] by applying to the data the P-term complex exponential model given by (3). P
y½n ¼ ∑ Ak eðαk þ jωk Þðn 1ÞT S þ jφk k¼1
ð3Þ
In (3), n ¼1, 2,…,N, and TS is the sampling period. The classic Prony method solves (3) by using an auxiliary system of linear equations to decouple this non-linear equation. The outline of the algorithm is present in [28] and described in detail in [29]. Given that the length of the mfVEP signals is composed of N ¼600 samples, approximation by applying (3) obtains P ¼300 components. Therefore, the mfVEP signal is characterized by the parameters Ak, αk, ωk and φk, (k ¼1,…, P ¼300). As the principal information components of mfVEP signals are in low frequencies, only the lowest M frequencies are chosen for the filter output (XProny[n]), as represented in (4). Establishing the optimal value of M obtains a simpler representation, while reducing the noise is one of our major findings. The results section demonstrates the approach followed to establish the optimal value of M. M
X Prony ½n ¼ ∑ Ap eðαp þ jωp Þðn 1ÞT S þ jφp ; p¼1
ω1 oω2 o ⋯ o ωM
ð4Þ
Note that frequency (ωp/2π) corresponds with Prony’s frequency, but is not an exact match for Fourier’s frequency. For example, the frequency calculated with Prony’s method can have arbitrary values and the frequency range is different for each signal.
3. Results This paper compares the quality of the output of three different processing methods: (a) The signal resulting from averaging the two recordings (raw data). The raw data is 600 samples (or coefficients) long. (b) The signal resulting from filtering raw data (0–500 ms, 600 samples) using a digital low-pass filter. Although it is possible to use more advanced signal-filtering options, such as a Chebyshev type 1 filter [4,30] filtering using the FFT is the most commonly employed option when processing mfVEP signals [27]. Based on both our own studies and those of other authors [31], we concluded that the best SNRS value is obtained by applying a sharp cut-off [1–20 Hz] between using a discrete Fourier Transform. Under these conditions, the filtered signal (XFFT) is represented by 20 coefficients. After the signal is filtered, IFFT is applied to obtain a register in time domain again with 600 samples-500 ms. (c) The signal resulting from filtering raw data using Prony’s method. Thus, SNRS_Prony is the signal-to-noise ratio in the signal window of the trace filtered, and RMSERAW–Prony is the root mean squared error (RMSE) between the original recordings and the signal filtered using Prony’s method: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ∑150 X RAW ðkÞ X Prony ðkÞ k ¼ 45 RMSERAWProny ¼ V units ð5Þ N 45 150 ms The best M value (number of components, ordered by frequency, in the Prony domain) is defined as the M value that maximizes the SNRS_Prony/RMSERAW-Prony ratio. This value represents a compromise between the SNR value of the filtered signal and the distortion between the original signal and the filtered signal. According to this criterion, the best value for our database is M¼7 (28 coefficients in Prony domain). With these components, the signals are reconstructed into time domain and 600 samples500 ms length signals are obtained. Fig. 1 shows the mfVEP signal for channel 6, sector 2, right eye, of one of the subjects selected at random. It shows the signal without pre-processing (XRAW), filtered in the frequency domain (XFFT), and filtered using Prony’s method (XProny).
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Table 1 shows the principal results obtained after analysing the mfVEP signals. The upper part of Fig. 2 shows the mean value of the SNRS obtained for each of the 24 subjects (each point is the mean of 2eyes 60sectors 6channels SNR values). For all of the signals obtained (NT ¼24 2 60 6¼ 17,280 traces), the mean value of the SNRS before filtering is SNRS_RAW ¼2.75 70.73 (95% CI: 2.44– 3.10); when the FFT is used, the value is SNRS_FFT ¼3.69 70.9 (95% CI: 3.31–4.1); and if the Prony filter is used, the value is SNRS_Prony ¼4.04 71.00 (95% CI: 3.62–4.47). An extremely significant difference is observed between the two types of filtering (Paired t-test: P o 0.0001). If the SNR gains of both filters are defined as follows: GFFT
%
¼
ðSNRS_Prony SNRS_RAW Þ ðSNRS_FFT SNRS_RAW Þ 100 and GProny% ¼ 100 SNRS _RAW SNRS_RAW
then complete set of recordings GFFT% ¼ þ34.18% and GProny% ¼ þ 46.90%. A criterion widely used to filter out low-quality recordings is to set a threshold of SNRS o ¼0.23 log units [32], which is equivalent to SNRS ¼1.7 (dashed blue line in Fig. 2). In our database, 63.88% of the unfiltered recordings would be analysable. If FFT filtering were used, 77.11% would be analysable, and if Prony’s method were used, 82.78% of the sectors would be analysable. 3.1. Results for the best channel The method most widely used to process mfVEP signals consists of selecting from each sector the channel with the best signal-to-noise ratio in the signal window. In this case, the number 0.5
15
of recordings available in the database amounts to 24 2 60 traces. The lower part of Fig. 2 shows the mean value of the SNRS obtained for each of the 24 subjects, selecting for each sector the best channel (each point is the mean of 2eyes 60sectors SNR values). Among the unfiltered signals, the mean value is SNRS_RAW ¼4.29 7 1.26 (95% CI: 3.76–4.83); for the signals filtered using the FFT, it is SNRS_FFT ¼5.73 71.42 (95% CI: 5.12–6.33); and for the signals filtered using Prony’s method, it is SNRS_Prony ¼ 6.20 71.58 (95% CI: 5.53–6.86). For the complete set of recordings, considering the best channel, GFFT% ¼ þ 33.56% and GProny% ¼ þ44.52%. An extremely significant difference exists between the two types of filtering (Paired t-test p o0.0001). In our database, and considering the best channel for each sector, without filtering, 90.93% of the visual field sectors are analysable (SNRS 41.7); with FFT filtering, 96.77% are analysable; and with Prony’s method, 97.91% are analysable. One question that arises is whether Prony’s method is better than the FFT filter at distinguishing mfVEP signals from noise in normal subjects. To address this, the SNRS and SNRN values defined in (2) are used. Fig. 3 shows the frequency distribution of the values of SNRS (red) and SNRN (blue) for 2880 records (24 individuals 2 eyes 60 mfVEP recordings from the best channel). The noise and signal distributions overlap. Therefore, some of the mfVEP responses from the signal window cannot be distinguished from noise. The capacity of the various methods to distinguish between signal and noise has been evaluated using Receiver Operating Characteristic (ROC) analysis. An ROC curve shows sensitivity against the False Positive Rate (FPR¼1 specificity). For all of the thresholds available between min(SNRN) and max(SNRS), the sensitivity (fraction of cases in which SNRS 4Threshold) and specificity
Raw FFT Prony
0.4 0.3
Level uV
0.2 0.1 0 -0.1 -0.2 -0.3 -0.4
0
50
100
150
200
250
300
350
400
450
500
Time ms
Fig. 1. Example of a waveform of mfVEP signal along the time in order to compare the effects of each method. The black trace is a raw signal, the blue trace is a signal filtered with a FFT filter and the red trace represents a signal processed with Prony method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. SNRS values for the 24 subjects. Each point represents the average SNR for each patient with each method. The upper part represents the values for All the channels and the lower part just for the Best Channel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 1 Comparison of results between 3 signal types.
All channels Best channels
RAW data
FFT filter
Prony filter
SNRS ¼2.75 70.73 (rank: 1.39–4.05); analysable sectors: 63.88% SNRS ¼4.29 71.26 (rank: 1.98–6.32); analysable sectors: 90.93% AUCRAW ¼0.9909 7 0.0191 (rank: 0.9069–0.9999) thresholds: 0–16.66
SNRS ¼ 3.69 7 0.9 (rank: 1.46–5.06) analysable sectors: 77.11% SNRS ¼ 5.737 1.42 (rank: 2.15–7.68) analysable sectors: 96.77% AUCFFT ¼ 0.9915 70.0183 (rank: 0.9083–0.9999) thresholds: 0–18.19
SNRS ¼4.04 71.00 (rank: 1.57–5.84). Analysable sectors: 82.78% SNRS ¼6.20 7 1.58 (rank: 2.30–8.78) analysable sectors: 97.91% AUCProny ¼ 0.9952 7 0.0125 (rank: 0.9376–0.9999) thresholds: 0–20.32
16
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Fig. 3. Distribution of SNRS and SNRN. X-axis represents all the possible SNR values for each case and the Y-axis represents the frequency of each value. Red traces correspond to SNR values of the signal window and blue traces correspond to SNR values from noise window. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
(fraction of cases in which SNRN oThreshold) values are calculated according to (6) and (7): Sensitivity ¼
casesðSNRS 4 ThresholdÞ casesðSNRS 4 ThresholdÞ þ casesðSNRS o ThresholdÞ ð6Þ
Specificity ¼
casesðSNRN o ThresholdÞ casesðSNRN 4ThresholdÞ þ casesðSNRN o ThresholdÞ ð7Þ
Fig. 4(a) shows the ROC curves for all of the subjects before filtering; Fig. 4(b) shows them after FFT filtering; Fig. 4(c) shows them after filtering using Prony’s method; and Fig. 4(d) shows the mean value of the ROC curve for the three signal types. Calculating an AUC value for each subject obtains the following results: AUCRAW ¼0.9909 70.0191 (95%CI: 0.9837–0.999), AUCFFT ¼0.9915 70.0183 (95%CI: 0.9838–0.9993) and AUCProny ¼ 0.9952 7 0.0125 (95%CI: 0.9901–1.000), with the mfVEP-Prony showing the largest area. The difference between AUCRAW and AUCFFT was not statistically significant (P¼0.081, Wilcoxon matched-pairs signed rank test). However, there was a difference between AUCFFT and AUCProny (P¼0.001) and between AUCRAW and AUCProny (P¼ 0.000).
4. Discussion The principal conclusion of this paper is that filtering healthy subjects’ mfVEP signals using Prony’s method improves the SNRS to a greater extent than filtering them in the frequency domain. As a consequence, this obtains a greater number of analysable visualfield sectors, making this method better at distinguishing mfVEP signals from noise than the standard filter method. This paper demonstrates that M ¼7 components of the Prony filter are sufficient to achieve a good approximation of the original signals in healthy subjects. The RAW signal is represented by 600 coefficients; the signal filtered using the FFT is represented by 20 coefficients; and the
signal filtered using Prony’s method is represented by 28¼ 7 4 coefficients, which is slightly higher than the FFT representation. If the quality of the mfVEP signal is quantified by its SNRS (analysing the six channels available per sector), for the overall database use of Prony’s method achieves an improvement of GProny% ¼ þ46.90%, while the Fourier filter achieves an improvement of GFFT% ¼ þ34.18%. Therefore, the advantage of using Prony’s method over the Fourier filter can be quantified as þ12.72%. These results mean that by using Prony’s method the percentage of traces that exceed the minimum value of SNRS 41.7 is 82.78% of visual field sectors, while using the FFT filter produces a figure of 77.11%. If for each sector the channel with the best SNRS (best channel) is selected, use of Prony’s method achieves a gain of GProny% ¼ þ44.52%, while use of the Fourier filter achieves a gain of GFFT% ¼ þ33.56%. In this case, with Prony’s method the proportion of sectors with signals that exceed the SNRS 41.7 threshold stands at 97.91%, while with the FFT filter the proportion stands at 96.77%. As is to be expected, the SNRS is greater when the best channel is selected. However, the gain obtained by selecting the best channel is smaller with both types of filtering than when 6 channels from each sector are considered. This effect is undoubtedly due to the fact that the signals with the best SNR in each channel are close to their upper quality limit and, as a consequence, the percentage gain is smaller. Analysis based on ROC curves shows that Prony’s method reduced the overlap between signal and noise distribution. Based on AUC criteria, Prony’s method obtains the best results (AUCProny ¼0.9952 7 0.0125) when compared against values of AUCFFT ¼0.9915 70.0183 and AUCRAW ¼0.9909 70.0191. No significant difference exists between AUCRAW and AUCFFT (the traditional processing method). However, a difference does exist between AUCFFT and AUCProny and between AUCRAW and AUCProny. The authors believe that the moderate margins for improvement in signal quality and in the number of analysable sectors according to their SNR is due to the fact that they have been recorded in healthy subjects. Therefore, they are generally of good quality and, consequently, the margin for improvement is very small. The next step would be to confirm that applying Prony’s
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
RAW
FFT 0.9 Sensitivity
Sensitivity
0.9 0.8 0.7
0
0.1
0.2 0.3 1-Specificity
Sensitivity
0.8 0.7 0.6 0
0.1
0.2 0.3 1-Specificity
0
0.4
0.1
0.2 0.3 1-Specificity
0.4
AVERAGE
1
0.9 Sensitivity
0.7
0.5
0.4
PRONY
0.5
0.8
0.6
0.6 0.5
17
0.9 0.8 0.7 0.6
0
0.1 0.2 1-Specificity
0.3
Fig. 4. (a) ROC curve for unfiltered signals. (b) ROC curve for signals filtered using FFT. (c) ROC curve for signals filtered using Prony’s method. In plots (a), (b) and (c) each trace represents a patient. (d) Representation of mean value.
method to patients with various pathologies achieves a greater gain in the parameters that define signal quality. Several papers confirm that the SNR of this study’s database is appropriate. For instance, [27] obtain the following formula for the SNR according to subject age: SNRS ¼ 5.1 (0.01 years), valid for the best channels of each sector and applied to data filtered using the Fourier transform (cut-off at 35 Hz). In their database, SNRS ¼5.59 (mean age¼49.0), a similar value to the one obtained in our signals (SNRS ¼5.73) after filtering them with the FFT (cut-off at 20 Hz). In [33] research the advantages of a prediction method based on 40 virtual electrodes. Using a database of 15 healthy participants, the standard method (bandwidth of 1 Hz–2 kHz) obtains an SNRS of 3.877 1.08, while with the prediction method the SNRS is 4.36 71.03, representing a gain of þ12.66%. Using Prony’s method, and when analysing the best channels, we obtained an improvement in the SNRS (GProny% ¼ þ44.52%) greater than that obtained using the prediction method. ROC analysis has been used previously to evaluate the overlap between SNR distribution in the signal and noise windows. For instance, Meigen and Kramer i [34] calculate the AUC of a sample of 30 subjects, positioning the four electrodes (A, B, C, D) within a 4-cm radius of the inion. Our original data (AUCRAW ¼0.9909 7 0.0191) exceed the best value obtained by Meigen and Kramer (0.9647 0.005), making it difficult to establish a direct comparison. In their paper, and using a single channel, they obtain a best value of AUC ¼0.93370.005 for A–B. In our case (AUCProny ¼ 0.9952 7 0.0125), we exceed that value. This is undoubtedly due to the different position of the electrodes and to the fact that we selected the best channel of the six available, while Meigen and Kramer obtained all their data from a single channel.
Ishikawa et al. [35] study a population of healthy Japanese subjects, obtaining an AUC of 0.983 7 0.027 for 54 left eyes and one of 0.9837 0.017 for 56 right eyes. Various channels are combined after obtaining them from electrodes positioned 4 cm above, left and right of the inion and 2.5 cm below the inion (ground electrode on the left earlobe). This value is comparable to that obtained in this paper using the FFT (AUCFFT ¼ 0.99157 0.0183). Nakamura et al. [36] study patients with glaucoma, obtaining for 62 control subjects a mean SNRS value of 3.09. They also analyse ROC curves, obtaining an AUC value of 0.98 70.02. This result can be compared with the results achieved in this paper using FFT filtering (AUCFFT ¼ 0.99157 0.0183), although the cut-off frequencies are different (35 Hz and 20 Hz) in our study. Our results exceed those obtained in [18] when using a PCA method to filter the signals. By selecting four principal components, they obtain a value of SNR ¼5, which is lower than the one obtained in this paper using Prony filtering. Several lines of research merit further investigation to ascertain if it is possible to improve on the results achieved here. For example, in the Prony filter (Eq. (4)), it would be interesting to see if better results are obtained by setting a cut-off frequency and selecting components below that value rather than by selecting a fixed number of components. Other mfVEP signal-filtering options also merit analysis, e.g. wavelet methods, classic principle components/ single value decomposition and hierarchical decomposition. [30]. In conclusion, we consider that this mfVEP-signal filtering method, although limited by the small sample size and by the fact that the research was carried out on healthy subjects, has great potential for use in diagnosing and monitoring visual pathway pathologies in patients with demyelinating and ischaemic diseases.
18
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
Conflict of interest statement The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Summary Due to development of the multifocal visual-evoked-potential (mfVEP) recording technique, it has been possible to obtain signals that reflect electrical activity in the primary visual cortex in response to stimuli in various sectors of the visual field. As is the case with other electrophysiological recordings, mfVEP signals are contaminated by noise. The amplitude of this type of signals is in the order of nanovolts–microvolts, and may even be essentially zero in some locations in control subjects due to local folding of the cortex. Sources of contamination may include the EEG’s alpha waves, other physiological signals emitted by the patient (muscletension-associated noise, ECG, EOG, etc.), or interference from the power grid or from other electronic equipment. Filtering is essential when performing patient diagnosis based on the results of signal classification. Very few papers investigate the influence of mfVEP signal filter techniques. This paper describes use of Prony’s method as a filter applied to multifocal visual-evokedpotential (mfVEP) signals. Prony’s method can be viewed as an extension of Fourier analysis that allows a signal to be decomposed into a linear combination of functions with different amplitudes, damping factors, frequencies and phase angles. To the best of our knowledge, Prony’s method has not been used to filter mfVEP signals. By selecting Prony method parameters, a frequency filter has been developed which improves signal-to-noise ratio (SNR). Three different criteria were applied to data recorded from subjects (48 eyes in 24 healthy subjects ranging in age from 23 to 38) to produce three separate datasets: unfiltered raw data, data filtered using the traditional method (fast Fourier transform: FFT), and data filtered using Prony’s method. Filtering using Prony’s method improved the signals’ original SNR by 44.52%, while the FFT filter improved the SNR by 33.56%. The extent to which signal can be separated from noise was analysed using receiver–operating– characteristic (ROC) curves. The area under the curve (AUC) was greater in the signals filtered using Prony’s method than in the original signals or in those filtered using the FFT. The principal conclusion of this paper is that filtering healthy subjects’ mfVEP signals using Prony’s method improves the SNRS to a greater extent than filtering them in the frequency domain. As a consequence, this obtains a greater number of analysable visualfield sectors, making this method better at distinguishing mfVEP signals from noise than the standard filter method. The authors believe that the moderate margins for improvement in signal quality and in the number of analysable sectors according to their SNR is due to the fact that they have been recorded in healthy subjects, so the next step would be to confirm that applying Prony’s method to patients with various pathologies achieves a greater gain in the parameters that define signal quality.
Acknowledgments This research has been supported by Spain’s Ministerio de Ciencia e Innovación under the “Advanced analysis of multifocal
ERG and visual evoked potentials applied to the diagnosis of optic neuropathies” project (ref. TEC2011-26066) and by grants FIS PI11/ 00533 and RETICS RD12/0034/0006. We wish to thank Dr. D.C. Hood for his generosity with the software used to analyze the mfVEP data.
References [1] H. a Baseler, E.E. Sutter, S. a Klein, T. Carney, The topography of visual evoked response properties across the visual field, Electroencephalogr. Clin. Neurophysiol. 90 (1) (1994) 65–81. [2] D.C. Hood, V.C. Greenstein, Multifocal VEP and ganglion cell damage: applications and limitations for the study of glaucoma,, Prog. Retin. Eye Res. 22 (2) (2003) 201–251. [3] W. Zhang, K. Zhao, Multifocal VEP difference between early- and late-onset strabismus amblyopia, Doc. Ophthalmol. 110 (2-3) (2005) 173–180. [4] R. Ruseckaite, T. Maddess, G. Danta, C.J. Lueck, A.C. James, Sparse multifocal stimuli for the detection of multiple sclerosis, Ann. Neurol. 57 (6) (2005) 904–913. [5] L.K. Grover, D.C. Hood, Q. Ghadiali, T.M. Grippo, A.S. Wenick, V.C. Greenstein, M.M. Behrens, J.G. Odel, A comparison of multifocal and conventional visual evoked potential techniques in patients with optic neuritis/multiple sclerosis, Doc. Ophthalmol. 117 (2) (2008) 121–128. [6] A.C. James, The pattern-pulse multifocal visual evoked potential,, Invest. Ophthalmol. Vis. Sci. 44 (2) (2003) 879–890. [7] X. Zhang, D.C. Hood, C.S. Chen, J.E. Hong, A signal-to-noise analysis of multifocal VEP responses: an objective definition for poor records, Doc. Ophthalmol. 104 (3) (2002) 287–302. [8] M. Bengtsson, S. Andreasson, G. Andersson, Multifocal visual evoked potentials —a method study of responses from small sectors of the visual field, Clin. Neurophysiol. 116 (8) (2005) 1975–1983. [9] E. Gutiérrez-Díaz, C. Pérez-Rico, M.J. de Atauri, E. Mencía-Gutiérrez, R. Blanco, Evaluation of the visual function in obstructive sleep apnea syndrome patients and normal–tension glaucoma by means of the multifocal visual evoked potentials, Graefes Arch. Clin. Exp. Ophthalmol. 250 (11) (2012) 1681–1688. [10] V. Parisi, L. Ziccardi, G. Stifano, L. Montrone, G. Gallinaro, B. Falsini, Impact of regional retinal responses on cortical visually evoked responses: multifocal ERGs and VEPs in the retinitis pigmentosa model, Clin. Neurophysiol. 121 (3) (2010) 380–385. [11] Y. Shimada, M. Horiguchi, A. Nakamura, Spatial and temporal properties of interocular timing differences in multifocal visual evoked potentials, Vision Res. 45 (3) (2005) 365–371. [12] A. Klistorner, H. Arvind, R. Garrick, S.L. Graham, M. Paine, C. Yiannikas, Interrelationship of optical coherence tomography and multifocal visualevoked potentials after optic neuritis, Invest. Ophthalmol. Vis. Sci. 51 (5) (2010) 2770–2777. [13] P. Sriram, A. Klistorner, H. Arvind, S.L. Graham, Reproducibility of multifocal VEP latency using different stimulus presentations, Doc. Ophthalmol. 125 (1) (2012) 43–49. [14] S. Park, S.H. Park, J.H. Chang, Y.-H. Ohn, Study for analysis of the multifocal visual evoked potential, Korean J. Ophthalmol. 25 (5) (2011) 334–340. [15] (a) G.C. Herbik, J. Hölzl, Reupsch, M.B. Hoffmann, Differential effects of optic media opacities on mfERGs and mfVEPs, Clin. Neurophysiol. 124 (6) (2013) 1225–1231. [16] C. Pieh, M.B. Hoffmann, M. Bach, The influence of defocus on multifocal visual evoked potentials, Graefes Arch. Clin. Exp. Ophthalmol. 243 (1) (2005) 38–42. [17] D.C. Hood, X. Zhang, J.E. Hong, C.S. Chen, Quantifying the benefits of additional channels of multifocal VEP recording, Doc. Ophthalmol. 104 (3) (2002) 303–320. [18] X. Zhang, D.C. Hood, Increasing the sensitivity of the multifocal visual evoked potential (mfVEP) technique: incorporating information from higher order kernels using a principal component analysis method, Doc. Ophthalmol. 108 (3) (2004) 211–222. [19] Prony, Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques, et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alkool, á différentes températures, J. l’école Polytech 1 (22) (1795) 24–76. [20] M. Hansson, T. Gänsler, G. Salomonsson, Estimation of single event-related potentials utilizing the Prony method, IEEE Trans. Biomed. Eng. 43 (10) (1996) 973–981. [21] V. Garoosi, B.H. Jansen, Development and evaluation of the piecewise Prony method for evoked potential analysis, IEEE Trans. Biomed. Eng. 47 (12) (2000) 1549–1554. [22] S.W. Chen, A two-stage discrimination of cardiac arrhythmias using a total least squares-based Prony modeling algorithm, IEEE Trans. Biomed. Eng. 47 (10) (2000) 1317–1327. [23] B.H. Jansen, A. Hegde, N.N. Boutros, Contribution of different EEG frequencies to auditory evoked potential abnormalities in schizophrenia, Clin. Neurophysiol. 115 (3) (2004) 523–533. [24] V. Sajith, A. Sukeshkumar, J. Jacob, Identification of diabetic neuropathic conditions of simulated ulnar nerve response using Prony’s method,, in: H.K. Kim, S.-I. Ao, M.A. Amouzegar, B.B. Rieger (Eds.), IAENG Transactions on
A. Fernández et al. / Computers in Biology and Medicine 56 (2015) 13–19
[25] [26]
[27] [28]
[29] [30]
[31]
Engineering Technologies, vol. 247, Springer, Dordrecht, the Netherlands, 2014, pp. 683–695. P. Lander, D. Jones, E. Berbari, R. Lazzara, Time-frequency structure of the high-resolution ECG,, J. Electrocardiol. 27 (207-12) (1994). T. Demiralp, a Ademoglu, Y. Istefanopulos, H.O. Gülçür, Analysis of eventrelated potentials (ERP) by damped sinusoids, Biol. Cybern. 78 (6) (1998) 487–493. B. Fortune, X. Zhang, D.C. Hood, S. Demirel, C. a Johnson, Normative ranges and specificity of the multifocal VEP, Doc. Ophthalmol. 109 (1) (2004) 87–100. T. Lobos, J. Rezmer, P. Janik, H. Amarís, M. Alonso, C. Álvarez, Application of wavelets and Prony method for disturbance detection in fixed speed wind farms,, Int. J. Electr. Power Energy Syst 31 (9) (2009) 429–436. S. Lawrence, Marple, Digital Spectral Analysis with Applications, Prentice-Hall, Inc, Englewood Cliffs, New Jersey, 1987. T. Maddess, A.C. James, R. Ruseckaite, E.A. Bowman, Hierarchical decomposition of dichoptic multifocal visual evoked potentials, Vis. Neurosci 23 (5) (2006) 703–712. J.A. Macfarlane, The Multifocal Visual Evoked Cortical Potential in Visual Field Mapping: A Methodological Study, University of Glasgow, 2008.
19
[32] D.C. Hood, X. Zhang, C. Rodarte, E.B. Yang, N. Ohri, B. Fortune, C. a Johnson, Determining abnormal interocular latencies of multifocal visual evoked potentials, Doc. Ophthalmol. 109 (2) (2004) 177–187. [33] B. a E. Mazinani, T.D. Waberski, A.W. a Weinberger, P. Walter, G.F. Roessler, Improving the quality of multifocal visual evoked potential results by calculating multiple virtual channels, Jpn. J. Ophthalmol. 55 (4) (2011) 396–400. [34] T. Meigen, M. Krämer, Optimizing electrode positions and analysis strategies for multifocal VEP recordings by ROC analysis, Vis. Res. 47 (11) (2007) 1445–1454. [35] K. Ishikawa, T. Nagai, Y. Yamada, A. Negi, M. Nakamura, Optimal conditions for multifocal VEP recording for normal Japanese population established by receiver operating characteristic analysis, Doc. Ophthalmol. 122 (1) (2011) 29–37. [36] M. Nakamura, K. Ishikawa, T. Nagai, A. Negi, Receiver–operating characteristic analysis of multifocal VEPs to diagnose and quantify glaucomatous functional damage, Doc. Ophthalmol. 123 (2) (2011) 93–108.
