Next Article in Journal
Hybrid Adaptive Cubature Kalman Filter with Unknown Variance of Measurement Noise
Next Article in Special Issue
Riemannian Spatio-Temporal Features of Locomotion for Individual Recognition
Previous Article in Journal
Investigation on Perceptron Learning for Water Region Estimation Using Large-Scale Multispectral Images
Previous Article in Special Issue
Two-Way Affective Modeling for Hidden Movie Highlights’ Extraction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Novel Instantaneous Phase Detection Approach and Its Application in SSVEP-Based Brain-Computer Interfaces

1
School of Electrical and Information Engineering, Tianjin University, Tianjin 300072, China
2
College of Intelligence and Computing, Tianjin University, Tianjin 300072, China
*
Author to whom correspondence should be addressed.
Sensors 2018, 18(12), 4334; https://doi.org/10.3390/s18124334
Submission received: 17 October 2018 / Revised: 29 November 2018 / Accepted: 3 December 2018 / Published: 7 December 2018

Abstract

:
This paper proposes a novel phase estimator based on fully-traversed Discrete Fourier Transform (DFT) which takes all possible truncated DFT spectra into account such that it possesses two merits of ‘direct phase extraction’ (namely accurate instantaneous phase information can be extracted without any correction) and suppressing spectral leakage. This paper also proves that the proposed phase estimator complies with the 2-parameter joint estimation model rather than the conventional 3-parameter joint model. Numerical results verify the above two merits and demonstrate that the proposed estimator can extract phase information from noisy multi-tone signals. Finally, real data analysis shows that fully-traversed DFT can achieve a better classification on the phase of steady-state visual evoked potential (SSVEP) brain-computer interface (BCI) than the conventional DFT estimator does. Besides, the proposed phase estimator imposes no restrictions on the relationship between the sampling rates and the stimulus frequencies, thus it is capable of wider applications in phase-coded SSVEP BCIs, when compared with the existing estimators.

1. Introduction

Estimating the frequency, the phase and the amplitude of a signal is a standard classical problem of signal processing. In particular, phase estimation has been widely applied in synchronization in communication [1], power analysis [2], speech enhancement [3], GPS navigation [4] and much more. Until now, phase estimators based on direct Discrete Fourier Transform (DFT) are the mainstream [5].
In [5], Liguori pointed out that, the existing DFT-based phase estimators heavily rely on the result of frequency estimate. To illustrate this dependency, let us study the sampled version (the sampling rate is F s ) of a complex exponential signal x ( t ) = a 0 exp [ j ( 2 π f 0 t + θ 0 ) ] as
x ( n ) = x ( t ) t = n Δ t = a 0 e j ( 2 π f 0 n Δ t + θ 0 ) , n = 0 , , N 1 ,
where a 0 , f 0 , θ 0 are amplitude, frequency and phase respectively and Δ t is the sampling interval 1/ F s . Accordingly, the frequency resolution of the N-point DFT X ( k ) is Δ f = F s /N = 1/N Δ t. Assume the peak DFT bin is at k = k 0 , k 0 Z + ( Z + refers to the set of positive integer numbers). Then, it can be deduced that the ideal phase value of X ( k 0 ) is
θ 0 = arg X ( k 0 ) π δ ( N 1 ) / N .
The variable δ in (2) is the fractional frequency offset ( 0.5 δ < 0.5 ) with the value δ = f 0 / Δ f k 0 , which is closely related to two sampling cases (coherent sampling or noncoherent sampling). For the coherent sampling case, i.e., f 0 = k 0 Δ f , δ = 0 and thus the sequence { x ( n ) , n = 0 , 1 , , N − 1} exactly covers k 0 signal periods, it can be easily inferred from (2) that the phase θ 0 can be estimated by directly taking the phase at the peak bin. Nevertheless, in case of noncoherent sampling, i.e., f 0 = ( k 0 + δ ) Δ f , δ 0 and thus the sequence does not exactly contain integer times of signal periods, it can be inferred from (2) that the phase estimate θ 0 is related to both the peak index k 0 and the frequency offset δ .
Up to now, a lot of frequency estimators have been proposed to estimate the fractional offset δ . For example, Offelli proposed an energy-based approach [6] in which the frequency offset of some component can be derived from the energy summation of the windowed DFT spectral bins around the peak bin. Other approaches based on interpolated FFT were reported in [2,7,8,9], in which δ is calculated via interpolation between several successive high-amplitude DFT bins centered with peak bin. In recent years, a lot of high-accuracy interpolation-based estimators such as Provencher estimator [10], Jacobsen estimator [11], Candan estimator [12], and phase difference-based estimator [13] were proposed. However, in [5,14,15], Liguori emphasized that “the bias of frequency estimate will introduce the uncertainty propagation to the phase acquisition”. In other words, the error of frequency offset δ in (2) will inevitably give rise to the error of phase estimate θ 0 accordingly.
Besides, spectral leakage is also a non-ignorable factor of degrading DFT-based phase estimators. In essence, for the non-coherent sampling case, performance degradation of phase estimation actually results from DFT’s inherent spectral leakage effect [16]. Therefore, to enhance the accuracy of phase estimation, some approach capable of suppressing spectral leakage is expected to be developed.
This paper proposes a novel approach of correcting DFT spectra which can alleviate phase estimate’s dependence on frequency estimate (or rather, the dependence on frequency offset δ ). This improvement actually attributes to the property that the proposed fully-traversed DFT spectrum has a much slighter spectral leakage compared to original uncorrected DFT spectrum in non-coherent sampling case. To verify these advantages, we also apply this corrected DFT-based phase estimator into a phase-coded steady-state visual evoked potential-based brain-computer interface (SSVEP-based BCI).
Recently, for the purpose of increasing realizable targets, phase information is popularly integrated into frequency-coding in SSVEP-BCI, such as joint/mixed frequency and phase coding [17,18,19]. In [20], Lee introduced one scheme of phase coding in SSVEP-BCI system, in which 8 LEDs flickering at 31.25 Hz with the phase interval 45 were used to represent 8 cursor functions. However, the phase decoding in [20] was realized through detecting the maximum amplitude peak of the averaged SSVEP waveform, i.e., it was realized in time domain instead of frequency domain. In [21], Lee segmented one stimulus into reference epoch and phase-shift epoch, in which the phase differences between these two epochs were used to distinguish different targets. Shortly after that, in [22], Lee designed six SDFS (stepping delay flickering sequences) with stimulus frequency of 32 Hz, among which 6 different phase-tagged objects are represented by 6 distinct delay segments such that they can be detected by the normalized power of averaged responses. However, SSVEP-BCIs in [21,22] actually employ special stimulus sequences to facilitate phase decoding, which may be not suitable for the common situations of phase measurement. In [23], to increase the number of visual stimuli, Gao and Jia proposed a frequency and phase mixed coding scheme using multiple frequencies. Then, Gao proposed the phase constrained canonical correlation analysis (p-CCA) to better distinguish frequency-tagged stimuli from multichannel SSVEPs [24], in which the phase information was estimated from the FFT result of the SSVEP records over the occipital cortex. Following this, Gao et.al introduced phase coding in CCA to discriminate six 60 -interval targets [25]. However, the stimuli frequencies in [23,24,25] were actually deliberately chosen such that each stimuli frequency exactly equals the integer times of FFT resolution and thus the frequency offset is also exactly zero. Therefore, in SSVEP-BCIs, novel phase estimator capable of detecting the phases of stimuli with any frequency offset is expected to be developed.
This paper will demonstrate that the proposed phase corrected DFT-based estimator is in accordance with the above demand of SSVEP-BCIs. Experimental results show that, the proposed phase estimator is superior to the conventional estimator in SSVEP phase extraction.
The rest of this paper is organized as follows. Section 2 introduces the derivation of the proposed fully-traversed DFT. Section 3 elaborates the properties of the proposed fully-traversed DFT in the noiseless case. Section 4 gives accuracy analysis of this phase estimator in noisy case. Section 5 demonstrates the proposed DFT-based phase estimator’s superiority to conventional phase estimators in the phase-coded SSVEP-BCI with auto-calibration. Finally, we conclude with a summary of results in Section 6.

2. Derivation of the Fully-Traversed DFT Spectrum

2.1. Phase Property of DFT Spectrum

Consider a sampled signal { x ( 0 ) , x ( 1 ) , , x (N 1 ) } (or denoted as the vector x 0 ). As is known, the normalized conventional DFT of x ( n ) is
X ( k ) = 1 N n = 0 N 1 x ( n ) e j 2 π N k n , k = 0 , 1 , , N 1 .
Without loss of generality, assume that Δ ω = 2 π / N and ω 0 = β Δ ω = ( k 0 + δ ) Δ ω , k 0 Z + and 0.5 δ < 0.5 . Substituting the complex exponential sequence x ( n ) = exp(j ( ω 0 n + θ 0 ) ) into (3) yields
X ( k ) = 1 N n = 0 N 1 e j 2 π N β n + θ 0 e j 2 π N k n = e j θ 0 N n = 0 N 1 e j 2 π N ( β k ) n = e j θ 0 N n = 0 N 1 e j 2 π N ( k 0 k + δ ) n .
Using geometric series summation and Euler equation, one can further deduce (4) as
X ( k ) = e j θ 0 + N 1 N ( k 0 + δ k ) π · sin ( k 0 + δ k ) π N sin ( k 0 + δ k ) π / N , δ 0 e j θ 0 , δ = 0
Substituting the peak index k = k 0 into (4) and synthesize the two cases of δ 0 and δ = 0 , we can derive a uniform expression of the phase value at peak DFT bin as
ϕ X ( k 0 ) = θ 0 + N 1 / N · δ π
Obviously, (6) is in accordance with (1).
From (6), it can be found that, only for the coherent sampling case (the frequency offset δ = 0 ), the observed phase ϕ X ( k 0 ) at the peak DFT bin is an accurate estimate. In contrast to this, for the non-coherent sampling case ( δ 0 ), in order to estimate the phase, δ has to be estimated in advance. As a result, the phase estimator based on conventional DFT heavily depends on the frequency offset δ .

2.2. The Proposed Fully-Traversed DFT Spectrum

In this work, we aim to develop a phase estimator capable of removing the dependency on the frequency offset δ . To be specific, our goal is to construct a corrected DFT spectrum whose peak bin provides accurate phase information whether δ = 0 or not.
To derive the proposed fully-traversed DFT, it is necessary to study the relationship between the ideal Fourier transform and the conventional DFT. As known to us, the ideal Fourier transform X ( j ω ) of an infinite-length sequence { x ( n ) } = { , x ( N + 1 ) , , x ( 0 ) , , x ( N − 1 ) , } is
X ( j ω ) = n = + x ( n ) e j n ω
However, X ( j ω ) cannot be realized, since the calculation in (7) consumes innumerous samples and memory. Thus, in practical applications, X ( j ω ) is replaced by the conventional normalized DFT X ( k ) defined in (3). Furthermore, if we define a sequence x 0 ( n ) , which is truncated from the infinite-length sequence { x ( n ) , n }, i.e.,
x 0 ( n ) = x ( n ) , n [ 0 , N 1 ] 0 , others ,
then, apparently, combining (3) with (7), (8) yields a representation of X ( k ) as
X ( k ) = 1 N X 0 ( j ω ) ω = k 2 π / N = 1 N n = + x 0 ( n ) e j n k 2 π / N = 1 N n = 0 N 1 x ( n ) e j n k 2 π / N
Now, let us focus on the sample x ( 0 ) at which the ideal phase θ 0 is located. It can be inferred from (9) that conventional DFT only considers one truncation case. However, there are N truncated sequences x m ( m = 0 , 1 , , N 1 ) containing x ( 0 ) listed as
x 0 = { x ( 0 ) , x ( 1 ) , , x ( N 2 ) , x ( N 1 ) } , x 1 = { x ( 1 ) , x ( 0 ) , , x ( N 3 ) , x ( N 2 ) } , x N 1 = { x ( N + 1 ) , x ( N + 2 ) , , x ( 1 ) , x ( 0 ) } ,
where the elements of sequence x m = { x m ( n ) , n = 0 , , N − 1} are
x m ( n ) = x ( n m ) , n [ 0 , N 1 ] 0 , others
Similar to the discrete spectrum of the truncated sequence x 0 expressed in (8), a reasonable discrete spectrum for x m = { x ( m ) , x ( m + 1 ) , , x ( m + N − 1 ) } should be
X m ( k ) = 1 N n = m N 1 m x ( n ) e j n k 2 π / N , m , k = 0 , 1 , , N 1 .
Furthermore, aiming to rewrite X m ( k ) in the form of conventional DFT in which variable n ranges from 0 to N − 1, we should detach the series summation in (12) into two terms as
X m ( k ) = 1 N n = m 1 x ( n ) e j n k 2 π / N + 1 N n = 0 N m + 1 x ( n ) e j n k 2 π / N ,
The first summation in (13) can be denoted as
1 N n = m 1 x ( n ) e j n k 2 π / N n = n + N ̲ ̲ 1 N n = N m N 1 x ( n N ) e j n N k 2 π / N = 1 N n = N m N 1 x ( n N ) e j n k 2 π / N .
Hence, it is necessary to introduce N new truncated sequences x ˜ m ( m = 0 , 1 , , N 1 ) whose elements are
x ˜ m ( n ) = x ( n ) , n = 0 , , N m 1 x ( n N ) , n = N m , , N 1 .
In terms of (15), N sequences x ˜ 0 x ˜ N 1 are listed as
x ˜ 0 = { x ( 0 ) , x ( 1 ) , , x ( N 2 ) , x ( N 1 ) } , x ˜ 1 = { x ( 0 ) , x ( 1 ) , , x ( N 2 ) , x ( 1 ) } , x ˜ N 1 = { x ( 0 ) , x ( N + 1 ) , , x ( 2 ) , x ( 1 ) } .
Then, combining (13) with (15), we have
X m ( k ) = 1 N n = 0 N 1 x ˜ m ( n ) e j n k 2 π / N = D F T [ x ˜ m ( n ) ] , m = 0 , 1 , , N 1 .
Finally, averaging all the DFT results of x ˜ 0 x ˜ N 1 yields the proposed corrected spectrum, i.e.,
Y ( k ) = 1 N m = 0 N 1 X m ( k ) , k = 0 , 1 , , N 1 .
From (8)–(18), it can be noticed that, the center sample x ( 0 ) fully traverses all the possible starting positions of N truncated sequences x m ( m = 0 , 1 , , N 1 ) in (10). Moreover, all the DFT results (17) of these traversed sub-sequences are fully averaged in (18) to yield the final spectrum Y ( k ) in (18). Hence, this novel spectral analysis is named as Fully-traversed DFT, whose dataflow is summarized as:
step 1. 
Construct N N-length sub sequences x 0 x N 1 from the given ( 2 N 1 ) input samples x ( N + 1 ) , , x ( 0 ) , , x ( N 1 ) , as (10) lists;
step 2. 
For each index m , m = 0 , 1 , , N 1 , circularly left move m samples of x m = { x ( m ) , , x ( 0 ) , , x ( m + N − 1 ) } to generate N sequences x ˜ m = { x ( 0 ) , x ( 1 ) , , x ( N m + 1 ) , x ( m ) , , x ( 1 ) } .
step 3. 
Implement normalized DFT on each x ˜ m to obtain the discrete sub-spectrum X m ( k ) , m = 0 , , N − 1;
step 4. 
Average all N sub-spectra X m ( k ) to acquire the proposed phase corrected spectrum Y ( k ) .

3. Property of Phase Estimation in the Noiseless Circumstance

3.1. Single-Tone Case

Consider a single-tone signal x ( n ) = exp(j ( ω 0 n + θ 0 ) ) , ω 0 = β Δ ω = ( k 0 + δ ) 2 π / N . Since Y ( k ) is obtained by averaging X 0 ( k ) X N 1 ( k ) , substituting x ( n ) into X m ( k ) in (12) and (18) yields
Y ( k ) = 1 N m = 0 N 1 1 N n = m m + N 1 e j β 2 n π / N + θ 0 e j n k 2 π / N n = n + m ̲ ̲ e j θ 0 N 2 m = 0 N 1 n = 0 N 1 e j ( β k ) 2 ( n m ) π / N = e j θ 0 N 2 m = 0 N 1 e j ( β k ) 2 m π / N n = 0 N 1 e j ( β k ) 2 n π / N
Obviously, (19) can be rewritten as
Y ( k ) = e j θ 0 e j θ 0 N m = 0 N 1 e j 2 π N ( k 0 k + δ ) m · e j θ 0 N n = 0 N 1 e j 2 π N ( k 0 k + δ ) n
where the superscript ‘∗’ represents complex conjugate operation. Then, combining (20) with (4) yields
Y ( k ) = e j θ 0 X ( k ) X ( k ) = e j θ 0 X ( k ) 2 , k = 0 , 1 , , N 1 .
Taking the phase part of (21), we can obtain the phase spectrum φ Y ( k ) as
φ Y ( k ) = θ 0 , k = 0 , 1 , , N 1 .
Equation (22) shows that, for the signal x ( n ) = exp ( j ( ω 0 n + θ 0 ) ) , N + 1 n N 1 , the phase values at all N spectral bins (including the peak bin k = k 0 ) uniformly equal the ideal phase θ 0 (also refers to the instantaneous phase of the center sample x ( 0 ) ). In other words, the synthesized spectrum Y ( k ) directly provides the accurate phase estimate at any spectral bin, thus removing the conventional DFT-based phase estimator’s dependency on the frequency offset δ .
Furthermore, since Y ( k ) is obtained by linearly averaging the N sub DFT spectra as (18) shows, the proposed fully-traversed DFT spectrum is of linearity, also. Thus, for a single-tone signal with amplitude a 0 ( a 0 1 ) , the following holds
Y ( k ) = a 0 e j θ 0 sin 2 ( k 0 + δ k ) π N 2 sin 2 ( k 0 + δ k ) π / N , k = 0 , 1 , , N 1 ,
among which the peak spectral bin is
Y ( k 0 ) = a 0 e j θ 0 sin 2 ( δ π ) N 2 sin 2 ( δ π / N ) .

3.2. Multi-Tone Case

Consider a signal containing Q ( Q 2 ) components expressed as
x ( n ) = q = 1 Q a q e j ( ω q n + θ q ) = q = 1 Q a q e j ( k q + δ q ) Δ ω n + θ q , k q Z + , 0 < | δ q | < 0.5 .
Since the conventional DFT is a linear transform, combining (5) and (25), its multi-tone spectrum X ( k ) is composed of Q single-tone spectra, i.e.,
X ( k ) = q = 1 Q X q ( k ) = q = 1 Q a q e j θ q + N 1 N ( k q + δ q k ) sin ( k q + δ i k ) π N sin ( k q + δ q k ) π / N , k = 0 , , N 1
Similarly, since the proposed phase corrected DFT spectrum is also of linearity, combining (5) with (23), its multi-tone spectrum Y ( k ) equals the summation of Q single-tone spectra as
Y ( k ) = q = 1 Q Y q ( k ) = q = 1 Q a q e j θ q sin 2 ( k q + δ q k ) π N 2 sin 2 ( k q + δ q k ) π / N , k = 0 , , N 1
Furthermore, to measure the phase of the i-th tone, the peak bin at k = k i should be focused. Therefore, the conventional DFT spectrum X ( k i ) can be written as
X ( k i ) = X i ( k i ) + q i X q ( k i ) = a i e j ( θ i + N 1 N δ i ) sin ( δ i π ) N sin ( δ i π / N ) + q i a q e j θ q + N 1 N ( k q + δ q k i ) sin ( k q + δ q k i ) π N sin ( k q + δ q k i ) π / N .
Similarly, the proposed fully-traversed DFT spectrum Y ( k i ) can be expressed as
Y ( k i ) = Y i ( k i ) + q i Y q ( k i ) = a i e j θ i sin 2 ( δ i π ) N 2 sin 2 ( δ i π / N ) + q i a q e j θ q sin 2 ( k q + δ q k i ) π N 2 sin 2 ( k q + δ q k i ) π / N .
In (29), the first term is the expected i-th single-tone spectrum which directly provides the phase estimate, and the second item represents the interference of other tones. Obviously, for either the conventional DFT spectrum or the fully-traversed DFT spectrum, the accuracy of the i-th tone’s phase estimation depends on the intensity of interference. Particularly, this accuracy depends on the relative magnitude ratio between X q ( k i ) , q i , and X i ( k i ) . From (28) and (29), we have
X q ( k i ) X i ( k i ) = a q a i · sin ( ( k q + δ q k i ) π ) / sin ( ( k q + δ q k i ) π / N ) sin ( δ i π ) / sin ( δ i π / N ) Y q ( k i ) Y i ( k i ) = a q a i · sin 2 ( ( k q + δ q k i ) π ) / sin 2 ( ( k q + δ q k i ) π / N ) sin 2 ( δ i π ) / sin 2 ( δ i π / N )
From (30), we have
Y q ( k i ) / Y i ( k i ) X q ( k i ) / X i ( k i ) = sin ( ( k q + δ q k i ) π ) / sin ( ( k q + δ q k i ) π / N ) sin ( δ i π ) / sin ( δ i π / N )
Please note that 0 <   | δ i | , | δ q |   < 0.5 , | k q k i |   1 and thus | k q + δ q k i |   > 0.5 . Since sin ( δ π ) / sin ( δ π / N ) is a monotonously descending even function, we have
sin ( ( k q + δ q k i ) π ) sin ( ( k q + δ q k i ) π / N ) < sin ( δ i π ) sin ( δ i π / N )
Combining (31) with (32), we have
Y q ( k i ) Y i ( k i ) < X q ( k i ) X i ( k i ) , q i .
Equation (33) shows that, for the fully-traversed DFT spectrum, the relative magnitude ratio between any other tone and the tone of interest is smaller than that of the conventional DFT. In other words, compared to the conventional DFT spectrum, the fully-traversed DFT magnitude spectrum does better in suppressing spectral leakage and interferences, thus yielding a higher accuracy of phase estimation.

3.3. Simplified Dataflow of Phase Estimation

From the aforementioned 4-step procedure listed in Section II, one can find that the fully-traversed DFT experiences N times of DFT operation. Therefore, to reduce computation complexity, this procedure needs to be simplified.
According to the linearity of DFT, the average of N sub-spectra is equivalent to the DFT result of the averaged data. Therefore, if we average the N sub-sequences x ˜ 0 x ˜ N 1 in (15) and (16), then one new sequence { y ( n ) , n = 0 , 1 , , N − 1} can be constructed as
y ( n ) = 1 N m = 0 N 1 x ˜ m ( n ) = N n N x ( n ) + n N x ( n N ) , n = 0 , 1 , , N 1
Accordingly, the DFT result Y(k) of y(n) is
Y ( k ) = 1 N n = 0 N 1 N n N x ( n ) + n N x ( n N ) e j 2 π N n k , k = 0 , 1 , , N 1
Clearly, Equation (35) only involves one time of DFT operation, thus greatly reducing computation complexity. Accordingly, its simplified dataflow is illustrated in Figure 1 (take N = 4 as an example), from which a low-complexity procedure of multi-tone phase estimation can be summarized as follows.
Firstly, weight the input ( 2 N 1 )-length data sequence [ x ( N + 1 ) , , x ( 0 ) , , x ( N − 1 ) ] with one ( 2 N − 1)-length triangular window [ 1 / N , , ( N 1 ) / N , 1 , ( N 1 ) / N , , 1 / N ] ;
Secondly, (N − 1) weighted data pairs (in each pair, the two data are spaced with N samples) are individually summed up to generate (N − 1) data y ( 1 ) , , y ( N − 1), except the center sample x ( 0 ) due to y ( 0 ) = x ( 0 ) ;
Thirdly, implement normalized DFT on the data sequence [ y ( 0 ) , y ( 1 ) , , y ( N − 1)] to provide the final spectral result Y = [ Y ( 0 ) , Y ( 1 ) , , Y ( N − 1)].
Lastly, collect all the peak indices k 1 , , k Q of Y ( k ) . For each index k q , directly taking the phase value of Y ( k q ) provides its phase estimates θ ^ q .

4. Variance Analysis of Phase Estimation in Noisy Circumstances

4.1. CRLB for Conventional DFT Phase Estimator

Now we consider the phase estimation in noisy case, in which one random complex Gaussian process η ( n ) should be considered in (1), i.e.,
s ( n ) = x ( n ) + η ( n ) = a 0 e j ( ω 0 n + θ 0 ) + η ( n ) ,
where n = N + 1 , , N 1 , ω 0 = β Δ ω = ( k 0 + δ ) 2 π / N , and η ( n ) is a complex Gaussian variable with mean zero and variance σ 2 . As mentioned in (2), the conventional phase estimate is
θ ^ 0 = φ X ( k 0 ) 1 1 / N β k 0 π = φ X ( k 0 ) 1 1 / N π δ .
(37) indicates that, since the conventional DFT-based phase estimator relies on the frequency offset δ , the estimate error of frequency offset δ will propagate to the phase estimate. In fact, this dependency makes phase estimation obey a 3-parameter mathematical model parameterized with α = [ ω 0 , θ 0 , a 0 ] T . With regard to this model, previous studies [26,27] have derived a CRLB (Crammer-rao lower bound) for the variance of the phase estimate (consuming 2 N 1 samples) as
C R L B 3 ( θ 0 ) = 4 N 3 2 N ( 2 N 1 ) ρ ,
where ρ = a 0 2 / 2 σ 2 is the SNR (Signal to Noise Ratio). Constrained by CRLB 3 ( θ 0 ) in (38), the error variance of phase estimate obtained by any conventional DFT-based estimator cannot exceed this bound. This has been especially claimed for algorithms in [1,2,5,6,7,8,9,27,28,29].

4.2. CRLB for the Proposed Phase Estimator

Different from conventional 3-parameter model-based phase estimators, mathematical model for the proposed corrected DFT-based phase estimator can be simplified. The reasons are as three aspects.
Firstly, (29) implies that the proposed phase detector only requires roughly searching the peak spectral bin and then taking the phases directly, i.e., independent of frequency offset δ .
Secondly, as previously mentioned, the proposed fully-traversed DFT spectrum equals the average of N sub DFT spectra. This actually reflects the following mechanism: averaging N sub vectors plays the role of compensating the angles of N sub DFT spectra with each other, which leads the synthesized phase to automatically fall at the ideal phase value, whether the frequency offset δ = 0 or not.
Lastly, as previously mentioned, the proposed phase estimator can directly determine the ideal phase θ 0 , referring to the ‘instantaneous phase’ of the center sample among the 2 N 1 input samples x ( N + 1 ) x ( N 1 ) . Obviously, the position for the center sample is at n = 0 , which can be easily determined in advance. Thus, if we rewrite x ( n ) = a 0 e j ( ω 0 n + θ 0 ) as a 0 e j φ n , then, for the position n = 0 , the entire term φ n equals the ideal phase θ 0 , which also indicates that estimation of frequency ω 0 can be omitted.
For the above 3 reasons, the error variance of fully-traversed DFT phase estimator obeys a 2-parameter mathematical model with α = [ θ 0 , a 0 ] T , in which the estimate of frequency offset δ is bypassed. Now we deduce the CRLB of this 2-parameter model using the classical parameter estimation theory [27].
Since η ( n ) is a Gaussian noise, the joint probability density function (pdf) of the ( 2 N 1 ) -length observation sequence S = [ s N + 1 , , s 0 , , s N 1 ] T conditioned on the unknown vector α = [ θ 0 , a 0 ] T is given by [26]
f ( S | α ) = 1 2 π σ 2 2 N 1 2 exp 1 2 σ 2 n = N + 1 N 1 s n x n 2
From (39), we can derive a 2 × 2 Fisher information matrix J whose entries are
J i j = E In f ( S | α ) α i α j = 1 σ 2 n = N + 1 N 1 x n α i · x n α j , i , j = 1 , 2 .
Since x ( n ) = u ( n ) + j v ( n ) = a 0 cos ( ω 0 n + θ 0 ) + j a 0 sin ( ω 0 n + θ 0 ) , (40) can be further expressed as
J i j = 1 σ 2 n = N + 1 N 1 u n α i · u n α j + v n α i · v n α j , i , j = 1 , 2
with
u n α 1 = u n θ 0 = a 0 sin ( ω 0 n + θ 0 ) ,
u n α 2 = u n a 0 = cos ( ω 0 n + θ 0 ) ,
v n α 1 = v n θ 0 = a 0 cos ( ω 0 n + θ 0 ) ,
v n α 2 = v n a 0 = sin ( ω 0 n + θ 0 ) .
Substituting (42)∼(45) into (41) yields
J 11 = a 0 2 σ 2 n = N + 1 N 1 sin 2 ( ω 0 n + θ 0 ) + cos 2 ( ω 0 n + θ 0 ) = a 0 2 σ 2 ( 2 N 1 ) ,
J 12 = J 21 = 0 ,
J 22 = 1 σ 2 n = N + 1 N 1 sin 2 ( ω 0 n + θ 0 ) + cos 2 ( ω 0 n + θ 0 ) = 1 σ 2 ( 2 N 1 ) .
Thus the 2 × 2 Fisher information matrix J takes the following diagonal form
J = 1 σ 2 ( 2 N 1 ) a 0 2 0 0 2 N 1 ,
and has the inverse as
J 1 = σ 2 1 ( 2 N 1 ) a 0 2 0 0 1 ( 2 N 1 ) a 0 2 .
Thus, the CRLB for the error variance of the proposed phase estimator is
C R L B 2 ( θ 0 ) = σ 2 ( 2 N 1 ) a 0 2 = 1 2 ( 2 N 1 ) ρ
Since the two CRLBs in (38) and (51) share a same sample length ( 2 N 1 ) , their ratio is
C R L B 2 ( θ 0 ) C R L B 3 ( θ 0 ) = N ( 4 N 3 ) 1 4 .
Hence, the CRLB for 2-parameter joint estimation is only 25% of that for 3-parameter joint estimation.

4.3. Numerical Results

To further verify the superiority of the proposed phase estimator, simulations performed under various noisy conditions and different spectral orders are presented. The phase error variance of the proposed estimator was also compared with that of conventional DFT-based estimator (we choose the ratio method based on interpolated DFT [2]). Assume k 0 = 3 , N = 128 and θ 0 = 60 . Then, the specific signal based on (36) is
s ( n ) = a 0 e j ( 3 + δ ) × 2 π N n + π 3 + η ( n ) ,
where the frequency offset δ is specified with 3 values: 0.1 , 0.2 , 0.3 . Figure 2a–c gives the error variance curves versus SNRs for these 3 frequency offsets, respectively. For each SNR and δ case, 500 Monte Carlo trials were conducted.
As can be seen in each figure, the majority of phase corrected DFT’s error variance curve (marked in ‘∗’) lies below CRLB 3 , proving that the proposed phase estimator is independent of the conventional 3-parameter joint estimation model. Furthermore, the proposed estimator’s error variance curve is bounded by CRLB 2 curve, verifying the correctness of (51). Figure 2 also demonstrates that the error variances of the conventional interpolated DFT estimator are nearly one order of magnitude higher than that of the proposed estimator.

5. Applying Fully-Traversed DFT in Phase-Coded SSVEP-BCI

Recently, BCIs have become a very hot topic in neural engineering. A BCI detects an user’s ongoing brain activities and translates them into meaningful messages, which helps patients with severe motor disabilities to express their messages to external world [30]. In particular, BCI based steady-state visual evoked potential (SSVEP) has received much attention in bioengineering research due to its satisfactory performance [31]. To increase the number of recognizable targets, the phase-coded SSVEP-BCIs use phase information to encode subject’s visual intention. In this system, the phase-tagged visual stimuli are characterized with flashing at one frequency but different phases, resulting in that subjects’ SSVEPs also differ in phase features. Therefore, through extracting the phase information of SSVEP potentials, a computer is able to distinguish which flicker the subject desires to select.

5.1. Experiment Paradigm

In general, SSVEP is always elicited after some latency time L (or labeled as lag phase ‘ θ L ’), which actually corresponds to a phase difference between flicker’s phase ‘ θ F ’ and SSVEP’s phase ‘ θ S ’. Hence, the relationship between the phase difference ‘ θ L ’ and latency L is described by
θ S = θ F + θ L
θ L = L × 360 × f s q × 360
where f s , q denote the stimulus frequency and the integer cycles, respectively. Under normal condition, L is stable in a short period of time but differs in inter-subject such that θ L cannot be calculated in advance [21,23,24]. From (54), it can be found that the flicker’s phase is usually not equal to SSVEP’s phase (or θ F θ S ) due to θ L , which means that we cannot directly use θ s to identify which flicker the subject desires to select if θ L is unknown. As a result, an additional measure of phase calibration is necessary for phase-coded SSVEP-BCI to calibrate this error θ L in the detection algorithm.
Hence, we build up a BCI system which uses a half-field phase-tagged stimulus to evoke the SSVEP with two different frequency components f 1 , f 2 at the same time. This system does not adopt any phase calibration since it is able to identify the flicker by introducing the phase difference instead of the phase under the assumption f 1 f 2 (this frequency distinction makes a subject more sensitive to flickering stimuli than to those with the same frequency), i.e.,
θ S ( f 1 ) θ S ( f 2 ) = θ F ( f 1 ) θ F ( f 2 ) + θ L ( f 1 ) θ L ( f 2 )
If f 1 f 2 , then θ L ( f 1 ) θ L ( f 2 ) . Therefore, the difference between θ S ( f 1 ) and θ S ( f 2 ) approximately equals the difference between θ F ( f 1 ) and θ F ( f 2 ) . Hence, the lag phase difference can be removed.
In our phase-coded SSVEP-BCI system, a visual stimulator (ViewSonic, 22 inch, 120 Hz refresh rate, 1680 × 1050 screen resolution) presenting two phase-tagged flickers (with the size 12 cm × 8 cm each) was used to evoke subjects’ SSVEPs (Figure 3 and Table 1).
It should be noted that, in our SSVEP-BCI illustrated in Figure 3, the selected flickering frequencies f 1 and f 2 should be as close as possible (this helps to remove the possible jump change of multiple of 360 for the lag phase difference between θ L ( f 1 ) and θ L ( f 2 ) , which was solved in [17] by means of a exhaustive search procedure based on least-square fitting). Otherwise, their lag phase difference ( θ L ( f 1 ) θ L ( f 2 ) ) would not be removed. However, due to the fact that all the stimulus frequencies in our SSVEP-BCI are acquired by integer dividing a fixed LCD display refresh frequency F r = 120 Hz (see [32]), we can only obtain a limited number of flikering frequencies (they are 120/7 Hz, 120/8 Hz, 120/9 Hz, 120/10 Hz, 120/11 Hz) falling at the visual sensitive region (10 Hz, 20 Hz). Therefore, among these candidate flickering frequencies, the frequency pair with the minimum interval is ( f 1 , f 2 ) = (120/11 Hz, 120/10 Hz) = (10.9 Hz, 12 Hz). Obviously, in this case, the lag phase difference between θ L ( f 1 ) and θ L ( f 2 ) will not be removed as f 1 and f 2 are not close enough. Hence, this phase difference should be taken into account in order to achieve an accurate result. In practice, it can be roughly estimated according to their empirical results of the apparent latency of SSVEPs [24]. In our experiment, ( θ L ( f 1 ) θ L ( f 2 ) ) cam be roughly estimated as 36 .
Different from the well-known CCA (Canonical Correlation Analysis) method, which also uses phase information to enhance the classification accuracy of SSVEPs, our proposed scheme consumes lower hardware cost. As [17] pointed out, the CCA method needs multiple stimulus frequencies (6 frequencies were adopted) to remove the possible jump change of multiple of 360 for the lag phase difference. In contrast, our proposed scheme only employs 2 frequencies, thereby lowering the hard cost.
Three subjects (S1∼S3, two males and one female) were seated on a comfortable chair before the visual stimulator in an illuminated room. The subjects’ EEG signals were recorded by a g.USBamp EEG amplifier from 13 electrodes (PO 3 , PO 5 , PO 7 , POZ, PO 4 , PO 6 , PO 8 , P1, PZ, P2, O1, Oz, and O 2 ). Specify the sampling rate F s = 600 samples/s. This experiment consisted of 5 runs containing 10 trials each. Each trial lasted for 8 s. Subjects were instructed to focus on one of flickers according to the following paradigm: From 0 to 2 s a cue appeared indicating which flashing flicker was required to focus on; From 2 s to 8 s the subjects gazed at the specified flicker; Then the next trial started. The order of gazed-flickers was ‘1212121212’ in each run. Thus this dataset had 25 trials for each flicker. The whole experiment lasted about 30 minutes. During this experiment, the subjects’ EEG signals were recorded for offline analysis later.

5.2. Procedure of SSVEP Phase Extraction

We collected a total amount of 50 trials for each subject. Basically, the procedure of SSVEP phase extraction would contain the following steps:
step 1. 
Apply both conventional DFT and phase corrected DFT to extract the phase values θ S ( f 1 ) and θ S ( f 2 ) , respectively;
step 2. 
Substitute θ S ( f 1 ) θ S ( f 2 ) and the estimate ( θ L ( f 1 ) θ L ( f 2 ) ) = 36 into (56) to estimate the difference Δ θ ^ = θ F ( f 1 ) θ F ( f 2 ) ;
step 3. 
Use the estimate Δ θ ^ to identify the gazed-flicker. If Δ θ ^ is close to 0 (or 180) deg, then the gazed-flicker is judged as flicker 1 (or 2).
The judgement involved in step 3 can actually be extended to distinguish C targets ( C 2 ) by finding the maximum among C decision variables R k as
R k = cos ( Δ θ ^ φ k ) , k = 1 , , C
where φ k refers to the ideal phase difference ( θ F ( f 1 ) θ F ( f 2 ) ) , i.e., the ideal clustering center of pattern recognition. For the above 2-category recognition problem, we have φ 1 = 0 , φ 2 = 180 . From (57), it can be found if the detected phase difference Δ θ ^ φ k , then R k 1 . In other words, if R 1 is close to 1, the gazed-flicker is judged as flicker 1 and vice versa. Furthermore, (57) also allows to assume more clustering centers. As will be elaborated, introducing more assumed clustering centers helps to evaluate the performance of phase estimator.

5.3. Result of Offline Analysis

The classification rates of this 2-class experiment are listed in Table 2.
It can be found that the proposed phase estimator can achieve a higher classification rate than conventional DFT-based estimator does, which is entirely over 10%.
Table 3 lists not only these two estimators’ detected phase values ( θ s (12), θ s (10.9)) but also 4 columns of phase decision values R k ( j ) , k = 1 , 2 , 3 , 4 , calculated by (57). Here decision variable R k ( j ) corresponds to the kth assumed clustering center (4 assumed phase clustering centers are 0 , 180 , 90 , 270 , respectively) while the subject actually gazes at the flicker j. Therefore, for any flicker index j ( j = 1 , 2 ) , R j ( j ) (marked with shadow in Table 3) is close to 1 and tends to be the maximum among R k ( j ) , k = 1 , 2 , 3 , 4 , if the accuracy of selected phase estimator is high enough. The data for all trials involved in Table 3 are recorded within 4 s.
In Table 3, generally speaking, the detected phases θ S ( 10.9 ) by the proposed phase estimator are closer to the ideal phase ‘ 0 ’ than that of conventional DFT estimator. This result can be explained as follows: Since the sampling rate F s = 600 samples/s and the window length equals 4s, it follows that N = 2400 samples are recorded and thus DFT frequency resolution Δ f = F s / N = 0.25 Hz. Hence, the stimulus frequency f 1 = 10.9 Hz can also be written as f 1 = 43.6 Δ f = ( 44 0.4 ) Δ f , i.e., the frequency offset δ equals nonzero value 0.4 , indicating that severe spectral leakage causing large phase measurement error will arise in the DFT spectrum. In contrast to this, for our proposed phase estimator, due to the property of ‘direct phase extraction’, the phase measurement error caused by frequency offset is much smaller than that of the conventional DFT case, as Table 3 lists.
Different from the case of θ S ( 10.9 ) , both estimators’ detected phases θ S ( 12 ) are uniformly close to the ideal phase 0 (or 180 ) (Even both mean and standard deviation are similar). This is because f 2 = 12 Hz = 48 Δ f , i.e., the frequency offset δ = 0 , resulting in that spectral leakage will not appear in both conventional DFT spectra and corrected DFT spectra.
Hence, since SSVEP phase extraction in this experiment is based on the phase difference between θ S ( 10.9 ) and θ S ( 12 ) rather than either of them, the conventional DFT phase estimator is more likely to cause large errors than the proposed phase estimator. From Table 3, one can notice that, for the conventional DFT phase estimator, the phase differences between θ S ( 10.9 ) and θ S ( 12 ) are far away from 36 deg (see the results in 3-rd, 7-th and 11-th row, respectively). Accordingly, the expected maximum decision values R 1 ( 1 ) or R 2 ( 2 ) detected by conventional DFT estimator are entirely smaller than that detected by the proposed estimator, as Table 3 lists. This reflects the fact that the conventional DFT is not so good as the proposed fully-traversed DFT in extracting the phase information of SSVEP.
To further evaluate the performance of the fully-traversed DFT in this experiment, as previously mentioned, we assume that there were 4 flickers with stimulus phases 0 , 180 , 90 , 270 on screen. Hence, the gazed-flicker was identified by finding the maximum among R 1 , R 2 , R 3 , R 4 . Moreover, the average classification rate for these 4 assumed flickers are also listed in Table 4. In this simulated 4-category experiment, compared to the 2-category case in Table 2, the classification rate of the fully-traversed DFT decreases but still over 80%. In contrast, the classification rate of conventional DFT is only about 37%, as Table 4 lists. The underlying reason lies in that conventional DFT cannot accurately extract the phase information at 10.9 Hz due to spectral leakage, while the proposed phase corrected DFT still works well for its property of ‘direct phase extraction’.
In summary, fully-traversed DFT is well suitable for the phase-coded SSVEP-BCI, especially in our proposed duel-frequency stimulus-based system which uses phase difference to recognize the gazed flick. Due to the fact that fully-traversed DFT does well in extracting the phase information at any frequency offset, the proposed phase estimator can achieve a higher classification rate. More importantly, as listed in Table 3, we find that the corresponding phase’s standard deviation of the proposed estimator is smaller than that of conventional DFT estimator, although the frequency resolution of fully-traversed DFT is less than DFT. The reason is actually mentioned in Section 2, i.e., the mechanism that fully-traversed DFT is the average of N DFT sub-spectra (see formula (18)) also helps to reduce the averaged noise’s power.
Based on our preliminary results, we believe that fully-traversed DFT can be applied in the online system. Although the simulated 4-category experiment results show that the corresponding average classification rate is only 80%, the classification rate is surely to be further enhanced once the estimate of the lag phase difference is more accurate. In this experiment, we only roughly estimate it as 36 deg. In fact, we can use another two very close frequencies (for example, we can replace LCD stimuli with LED stimuli) such that the lag phase difference will be close to zero and thus we do not need to estimate it. Moreover, it was also clearly found that the phase variance of fully-traversed DFT estimator does not get worse than conventional DFT under our experiment condition.

6. Conclusions

A novel phase estimator based on corrected-phase DFT was proposed in this paper. Due to considering all possible truncated sequences containing the center sample, spectral leakage in corrected-phase DFT is greatly reduced and thus the instantaneous phase information of the center sample can be directly extracted. In addition, we have also proved that the process of corrected-phase DFT is equivalent to a streamline dataflow, and thus the proposed phase estimator has a relatively low computational complexity. Furthermore, the phase error variance of the proposed phase estimator follows a 2-parameterized mathematical model and thus it has a higher accuracy than the conventional DFT-based phase estimators in noisy circumstances.
We also applied the proposed phase estimators to phase-coded SSVEP-BCIs. In particular, compared to the conventional DFT-based estimators, our offline experiment results demonstrate that the fully-traversed DFT does better in extracting the phase information of phase-coded SSVEP-BCI. Moreover, the proposed phase estimator imposes on restrictions on the relationship between the sampling rates and the stimulus frequencies (i.e., non-synchronous sampling is allowed), it is of wider applications in phase-coded SSVEP BCIs than the existing estimators. Our future work is to improve our system design and apply fully-traversed DFT in the practical SSVEP-BCI online system.

Author Contributions

Conceptualization and Methodology, X.H. and J.X.; Writing-Original Draft Preparation, Z.W.; Writing-Review and Editing, J.X.; Software and Validation, J.X. and X.H.; Supervision, Z.W.

Funding

This research was funded by National Natural Science Foundation of China, grant number 61671012.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Rice, F.; Rice, M.; Cowley, B. A new bound and algorithm for Star 16-QAM carrier phase estimation. IEEE Trans. Commun. 2003, 51, 161–165. [Google Scholar] [CrossRef]
  2. Andria, G.; Savino, M.; Trotta, A. Windows and interpolation algorithms to improve electrical measurement accuracy. IEEE Trans. Instrum. Meas. 1989, 38, 856–863. [Google Scholar] [CrossRef]
  3. Abe, T.; Honda, M. Sinusoidal model based on instantaneous frequency attractors. IEEE Trans. Audio Speech Lang. Process. 2006, 14, 1292–1300. [Google Scholar] [CrossRef]
  4. Dach, R.; Schildknecht, T.; Springer, T.; Dudle, G.; Prost, L. Continuous time transfer using GPS carrier phase. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2002, 49, 1480–1490. [Google Scholar] [CrossRef] [PubMed]
  5. Liguori, C.; Paolillo, A.; Pignotti, A. Estimation of signal parameters in frequency domain in presence of harmonic interference: A comparative analysis. IEEE Trans. Instrum. Meas. 2006, 55, 562–569. [Google Scholar] [CrossRef]
  6. Offelli, C.; Petri, D. A frequency-domain procedure for accurate real–time signal parameter measurement. IEEE Trans. Instrum. Meas. 1990, 39, 363–368. [Google Scholar] [CrossRef]
  7. Offelli, C.; Petri, D. Interpolation techniques for real-time multifrequency waveform analysis. In Proceedings of the Conference Record, 6th IEEE, IMTC-89 Instrumentation and Measurement Technology Conference, Washington, DC, USA, 25–27 April 1989; pp. 325–331. [Google Scholar] [CrossRef]
  8. Schoukens, J.; Pintelon, R.; Van Hamme, H. The interpolated fast Fourier transform: A comparative study. IEEE Trans. Instrum. Meas. 1992, 41, 226–232. [Google Scholar] [CrossRef]
  9. Agrez, D. Weighted multipoint interpolated DFT to improve amplitude estimation of multifrequency signal. IEEE Trans. Instrum. Meas. 2002, 51, 287–292. [Google Scholar] [CrossRef]
  10. Provencher, S. Estimation of Complex Single-Tone Parameters in the DFT Domain. IEEE Trans. Signal Process. 2010, 58, 3879–3883. [Google Scholar] [CrossRef]
  11. Jacobsen, E.; Kootsookos, P. Fast, Accurate Frequency Estimators [DSP Tips Tricks]. IEEE Signal Process. Mag. 2007, 24, 123–125. [Google Scholar] [CrossRef]
  12. Candan, C. Analysis and Further Improvement of Fine Resolution Frequency Estimation Method From Three DFT Samples. IEEE Signal Process. Lett. 2013, 20, 913–916. [Google Scholar] [CrossRef]
  13. Huang, X.; Xia, X.G. A Fine Resolution Frequency Estimator Based on Double Sub-segment Phase Difference. IEEE Signal Process. Lett. 2015, 22, 1055–1059. [Google Scholar] [CrossRef]
  14. Betta, G.; Liguori, C.; Pietrosanto, A. Propagation of uncertainty in a discrete Fourier transform algorithm. Measurement 2000, 27, 231–239. [Google Scholar] [CrossRef]
  15. Novotny, M.; Slepicka, D.; Sedlacek, M. Uncertainty Analysis of the RMS Value and Phase in the Frequency Domain by Noncoherent Sampling. IEEE Trans. Instrum. Meas. 2007, 56, 983–989. [Google Scholar] [CrossRef]
  16. Agrez, D. Improving phase estimation with leakage minimization. IEEE Trans. Instrum. Meas. 2005, 54, 1347–1353. [Google Scholar] [CrossRef]
  17. Ke, L.; Wang, Y.; Gao, X. Time-frequency joint coding method for boosting information transfer rate in an SSVEP based BCI system. In Proceedings of the International Conference of the IEEE Engineering in Medicine and Biology Society, Orlando, FL, USA, 16–20 August 2016; pp. 5873–5876. [Google Scholar]
  18. Youssef, A.A.A.; Wittevrongel, B.; Van Hulle, M.M. Accurate Decoding of Short, Phase-Encoded SSVEPs. Sensors 2018, 18, 794. [Google Scholar] [CrossRef]
  19. Zhao, X.; Zhao, D.; Wang, X.; Hou, X. A SSVEP Stimuli Encoding Method Using Trinary Frequency-Shift Keying Encoded SSVEP (TFSK-SSVEP). Front. Hum. Neurosci. 2017, 11, 278. [Google Scholar] [CrossRef] [PubMed]
  20. Lee, P.L.; Sie, J.J.; Liu, Y.J.; Wu, C.H.; Lee, M.H.; Shu, C.H.; Li, P.H.; Sun, C.W.; Shyu, K.K. An SSVEP-Actuated Brain Computer Interface Using Phase-Tagged Flickering Sequences: A Cursor System. Ann. Biomed. Eng. 2010, 38, 2383–2397. [Google Scholar] [CrossRef]
  21. Wu, H.Y.; Lee, P.L.; Chang, H.C.; Hsieh, J.C. Accounting for Phase Drifts in SSVEP-Based BCIs by Means of Biphasic Stimulation. IEEE Trans. Biomed. Eng. 2011, 58, 1394–1402. [Google Scholar] [CrossRef]
  22. Chang, H.C.; Lee, P.L.; Lo, M.T.; Lee, I.H.; Yeh, T.K.; Chang, C.Y. Independence of Amplitude-Frequency and Phase Calibrations in an SSVEP-Based BCI Using Stepping Delay Flickering Sequences. IEEE Trans. Neural Syst. Rehabil. Eng. 2012, 20, 305–312. [Google Scholar] [CrossRef]
  23. Jia, C.; Gao, X.; Hong, B.; Gao, S. Frequency and Phase Mixed Coding in SSVEP-Based Brain–Computer Interface. IEEE Trans. Biomed. Eng. 2011, 58, 200–206. [Google Scholar] [CrossRef] [PubMed]
  24. Pan, J.; Gao, X.; Duan, F.; Yan, Z.; Gao, S. Enhancing the classification accuracy of steady-state visual evoked potential-based brain-computer interfaces using phase constrained canonical correlation analysis. J. Neural Eng. 2011, 8, 036027. [Google Scholar] [CrossRef] [PubMed]
  25. Li, Y.; Bin, G.; Gao, X.; Hong, B.; Gao, S. Analysis of phase coding SSVEP based on canonical correlation analysis (CCA). In Proceedings of the 2011 5th International IEEE/EMBS Conference on Neural Engineering (NER), Cancun, Mexico, 27 April–1 May 2011; pp. 368–371. [Google Scholar]
  26. Rife, D.; Boorstyn, R. Single tone parameter estimation from discrete-time observations. IEEE Trans. Inf. Theory 1974, 20, 591–598. [Google Scholar] [CrossRef]
  27. Kay, M. Fundamentals of Statistical signal processing, Volume 2: Detection theory. In Blind Equalization and System Identification; Springer: London, UK, 1998; pp. 83–182. [Google Scholar]
  28. Reisenfeld, S.; Aboutanios, E. A new algorithm for the estimation of the frequency of a complex exponential in additive Gaussian noise. IEEE Commun. Lett. 2003, 7, 549–551. [Google Scholar] [CrossRef]
  29. Zhu, L.; Song, X.; Li, H.; Ding, H. High accuracy estimation of multi-frequency signal parameters by improved phase linear regression. Signal Process. 2007, 85, 1066–1077. [Google Scholar] [CrossRef]
  30. Wolpaw, J.R.; Birbaumer, N.; McFarland, D.J. Brain computer interfaces for communication and control. Clin. Neurophysiol. 2002, 113, 767–791. [Google Scholar] [CrossRef]
  31. Wang, Y.; Gao, X.; Hong, B.; Jia, C.; Gao, S. Brain-Computer Interfaces Based on Visual Evoked Potentials. IEEE Eng. Med. Biol. Mag. 2008, 27, 64–71. [Google Scholar] [CrossRef]
  32. Wong, C.M.; Wang, B.; Wan, F.; Mak, P.U.; Mak, P.I.; Vai, M.I. An improved phase-tagged stimuli generation method in steady-state visual evoked potential based brain-computer interface. In Proceedings of the 2010 3rd International Conference on Biomedical Engineering and Informatics (BMEI), Yantai, China, 16–18 October 2010; Volume 2, pp. 745–749. [Google Scholar]
Figure 1. Simplified flow diagram of phase corrected DFT (N = 4).
Figure 1. Simplified flow diagram of phase corrected DFT (N = 4).
Sensors 18 04334 g001
Figure 2. Error variances of the phase estimator with (a) δ = 0.1 , (b) δ = 0.2 and (c) δ = 0.3 .
Figure 2. Error variances of the phase estimator with (a) δ = 0.1 , (b) δ = 0.2 and (c) δ = 0.3 .
Sensors 18 04334 g002
Figure 3. Visual stimulator presenting two half-field phase-tagged stimuli. (Left field of Flicker 1: 10.9 Hz and 0 deg, Right field of Flicker 1: 12 Hz and 0 deg, Left field of Flicker 2: 10.9 Hz and 0 deg, Right field of Flicker 2: 12 Hz and 180 deg).
Figure 3. Visual stimulator presenting two half-field phase-tagged stimuli. (Left field of Flicker 1: 10.9 Hz and 0 deg, Right field of Flicker 1: 12 Hz and 0 deg, Left field of Flicker 2: 10.9 Hz and 0 deg, Right field of Flicker 2: 12 Hz and 180 deg).
Sensors 18 04334 g003
Table 1. The parameters of two phase-tagged half-field flickers.
Table 1. The parameters of two phase-tagged half-field flickers.
Left-FieldRight-Field
Freq. (Hz)Phase (deg)Freq. (Hz)Phase (deg)
Flicker 110.90120
Flicker 210.9012180
Table 2. The classification rate of each subject under different window length.
Table 2. The classification rate of each subject under different window length.
SubjectMethodElectrodeWindow Length (sec)Average
3456
S1ProposedPOZ0.941.001.000.980.98
DFTPO70.920.880.780.780.84
S2ProposedO20.981.001.001.001.00
DFTP20.920.880.840.640.82
S3ProposedP10.800.960.920.940.91
DFTPZ0.920.940.880.740.87
AverageProposed0.910.990.970.970.96
DFT0.920.900.830.730.84
Table 3. The phase and phase feature R j of each subject (mean ± S.D.)
Table 3. The phase and phase feature R j of each subject (mean ± S.D.)
SubjectMethodFlicker j θ S (12) (deg) θ S (10.9) (deg) R 1 (j) R 2 (j) R 3 (j) R 4 (j)
S1Proposed1322.74 ± 29.36.84 ± 31.40.985 ± 0.030.133 ± 0.110.741 ± 0.10.651 ± 0.13
2126.52 ± 23.12.25 ± 27.20.253 ± 0.160.954 ± 0.050.56 ± 0.180.789 ± 0.18
DFT1329.22 ± 39.357.09 ± 38.00.857 ± 0.170.416 ± 0.260.896 ± 0.10.378 ± 0.22
2151.46 ± 31.949.97 ± 34.10.422 ± 0.220.865 ± 0.160.445 ± 0.230.848 ± 0.19
S2Proposed1334.22 ± 36.123.10 ± 35.50.983 ± 0.030.132 ± 0.120.771 ± 0.090.62 ± 0.12
2146.96 ± 31.017.47 ± 35.70.193 ± 0.130.972 ± 0.040.601 ± 0.160.774 ± 0.12
DFT1332.91 ± 32.272.68 ± 37.70.833 ± 0.110.511 ± 0.190.95 ± 0.060.25 ± 0.18
2148.24 ± 29.450.71 ± 28.80.421 ± 0.190.881 ± 0.120.362 ± 0.160.916 ± 0.07
S3Proposed115.04 ± 42.434.66 ± 43.20.899 ± 0.110.347 ± 0.250.538 ± 0.310.75 ± 0.24
2198.30 ± 26.024.81 ± 35.50.339 ± 0.210.913 ± 0.090.822 ± 0.190.471 ± 0.27
DFT19.71 ± 39.073.81 ± 36.80.899 ± 0.130.349 ± 0.240.794 ± 0.220.508 ± 0.27
2198.94 ± 24.581.04 ± 38.30.323 ± 0.260.892 ± 0.190.522 ± 0.220.814 ± 0.15
AverageProposed10.956 ± 0.060.204 ± 0.160.683 ± 0.170.674 ± 0.17
20.262 ± 0.170.946 ± 0.060.661 ± 0.180.678 ± 0.19
DFT10.863 ± 0.140.425 ± 0.230.88 ± 0.130.379 ± 0.22
20.389 ± 0.220.879 ± 0.160.443 ± 0.20.859 ± 0.14
Table 4. The classification rate of each subject under different window length (Simulated 4-classes experiment).
Table 4. The classification rate of each subject under different window length (Simulated 4-classes experiment).
SubjectMethodElectrodeWindow Length (sec)Average
3456
S1ProposedPOZ0.800.900.900.840.86
DFTPO70.670.430.180.180.37
S2ProposedO20.920.900.980.900.93
DFTP20.540.280.160.060.26
S3ProposedP10.540.580.680.740.64
DFTPZ0.560.620.460.300.49
AverageProposed0.750.790.850.830.81
DFT0.590.440.270.180.37

Share and Cite

MDPI and ACS Style

Huang, X.; Xu, J.; Wang, Z. A Novel Instantaneous Phase Detection Approach and Its Application in SSVEP-Based Brain-Computer Interfaces. Sensors 2018, 18, 4334. https://doi.org/10.3390/s18124334

AMA Style

Huang X, Xu J, Wang Z. A Novel Instantaneous Phase Detection Approach and Its Application in SSVEP-Based Brain-Computer Interfaces. Sensors. 2018; 18(12):4334. https://doi.org/10.3390/s18124334

Chicago/Turabian Style

Huang, Xiangdong, Jingwen Xu, and Zheng Wang. 2018. "A Novel Instantaneous Phase Detection Approach and Its Application in SSVEP-Based Brain-Computer Interfaces" Sensors 18, no. 12: 4334. https://doi.org/10.3390/s18124334

APA Style

Huang, X., Xu, J., & Wang, Z. (2018). A Novel Instantaneous Phase Detection Approach and Its Application in SSVEP-Based Brain-Computer Interfaces. Sensors, 18(12), 4334. https://doi.org/10.3390/s18124334

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop