使用滤波器组的FBCCA算法
1、FBCCA算法原理
著称为光驱动响应(photic driving response)的SSVEP信号,可以通过刺激频率及其谐波频率下的正弦波形来表征。更确切的说,SSVEP信号由与刺激频率相同、谐波(基波的整数倍)和亚谐波的大脑反应组成。
除了基本的频率成分外,谐波频率成分也可以为频率识别提供有用的信息。然而,基于CCA的频率识别算法并没有很好的利用谐波频率成分。一部分研究人员尝试采用不同的谐波成分作为CCA算法中的模板信号,但对最终的分类精度并没有起到太大的作用。
为了提取嵌入在谐波频率成分中的有效信息,我们需要一种更有效的算法,即CCA算法的改进版——FBCCA算法。
FBCCA(Filter Bank Canonical Corrlation Analysis)算法将SSVEP信号分解成不同频带的多个子成分,即用一个带通滤波器组(An array of band-pass filters)将输入的脑电信号分解成N个子频带成分,再在这N个子成分上进行CCA算法。之后把N个子频带上的得到的相关系数做一个加权平均,对应于每一个刺激频率 f k , k = 1 , 2 , . . . , S f_k,k=1,2,...,S fk,k=1,2,...,S得到一个整体的相关系数值。最终从这S个相关系数值中选出最大的相关系数值,其对应频率即作为最终的识别结果。
2、FBCCA算法过程推导
FBCCA算法包含以下3个主要步骤:
(1). 滤波器组分析(Filter Bank Analysis)
(2). SSVEP子频带成分与正余弦参考信号的CCA算法(CCA between SSVEP sub-band components and sinusoidal reference signals)
(3). 目标识别(target identification)
2.1 滤波器组分析
首先,我们使用零相位的I型切比雪夫滤波器(其他滤波器也可以)将原始EEG信号分解成 X S B 1 , X S B 2 , . . . , X S B N X_{SB_1},X_{SB_2},..., X_{SB_N} XSB1,XSB2,...,XSBN N个子频带。
2.2 子频带成分的CCA
将标准CCA算法分别应用至每个子频带成分,得到每个子频带成分与预定义好的参考信号(对应于所有刺激频率 Y f k , k = 1 , 2 , . . . , S Y_{f_k},k=1,2,...,S Yfk,k=1,2,...,S,S为刺激数目)之间的相关系数值。
其中,第
k
k
k个参考信号对应的相关系数值由向量
ρ
k
⃗
\vec{ρ_k}
ρk表示,其包含N个子频带对应的相关系数值。
ρ
k
=
[
ρ
k
1
ρ
k
2
⋮
p
k
N
]
=
[
ρ
(
X
S
B
1
T
W
X
(
X
S
B
1
Y
f
k
)
,
Y
T
W
Y
(
X
S
B
1
Y
f
k
)
)
ρ
(
X
S
B
2
T
W
X
(
X
S
B
2
Y
f
k
)
,
Y
T
W
Y
(
X
S
B
2
Y
f
k
)
)
⋮
ρ
(
X
S
B
N
T
W
X
(
X
S
B
N
Y
f
k
)
,
Y
T
W
Y
(
X
S
B
N
Y
f
k
)
)
]
\begin{aligned} ρ_k&=\left[ \begin{matrix} ρ_{k_1} \\ ρ_{k_2} \\ \vdots \\ p_{k_N} \end{matrix} \right]\\ \\ &=\left[ \begin{matrix} ρ(X_{SB_1}^TW_X(X_{SB_1}Y_{f_k}),Y^TW_Y(X_{SB_1}Y_{f_k})) \\ \\ ρ(X_{SB_2}^TW_X(X_{SB_2}Y_{f_k}),Y^TW_Y(X_{SB_2}Y_{f_k})) \\ \vdots \\ ρ(X_{SB_N}^TW_X(X_{SB_N}Y_{f_k}),Y^TW_Y(X_{SB_N}Y_{f_k})) \end{matrix} \right]\\ \end{aligned}
ρk=⎣⎢⎢⎢⎡ρk1ρk2⋮pkN⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡ρ(XSB1TWX(XSB1Yfk),YTWY(XSB1Yfk))ρ(XSB2TWX(XSB2Yfk),YTWY(XSB2Yfk))⋮ρ(XSBNTWX(XSBNYfk),YTWY(XSBNYfk))⎦⎥⎥⎥⎥⎥⎤
2.3 目标识别
将
ρ
k
ρ_k
ρk中N个子频带成分
ρ
k
1
,
ρ
k
2
,
.
.
.
,
ρ
k
N
ρ_k^1,ρ_k^2,...,ρ_k^N
ρk1,ρk2,...,ρkN 做一个加权平方和融合,即
ρ
k
~
=
∑
n
=
1
N
w
(
n
)
∗
(
p
n
k
)
2
\tilde{ρ_k}=\sum_{n=1}^{N}w(n)*(p_n^k)^2
ρk~=n=1∑Nw(n)∗(pnk)2
其中,加权函数
w
(
n
)
w(n)
w(n)的定义为:
w
(
n
)
=
n
−
a
+
b
,
n
∈
[
1
,
N
]
w(n)=n^{-a} + b,n∈[1,N]
w(n)=n−a+b,n∈[1,N]
其中,a和b均为常量,它们的取值由分类器性能表现至最佳情况时所定。在实际的应用中,我们可以在离线实验中通过grid search(网格搜索)来确定a,b。
最终,我们得到了S个刺激频率对应的S个加权相关系数 ρ k ~ , k = 1 , 2 , . . . , S \tilde{ρ_k},k=1,2,...,S ρk~,k=1,2,...,S,取其中最大的 p k ~ \tilde{p_k} pk~对应的刺激频率,即为我们识别的目标频率。
f t a r g e t = max f i { ρ ~ 1 , ρ ~ 2 , . . . , ρ ~ S } , i ∈ [ 1 , S ] f_{target}={\underset{f_i} {\operatorname{max}}}\{ {\tildeρ_{1},\tildeρ_{2},...,\tildeρ_{S}} \},i∈[1,S] ftarget=fimax{ρ~1,ρ~2,...,ρ~S},i∈[1,S]
实现代码: https://github.com/YuDongPan/Canonical_Classifier
参考文献:
- Xiaogang Chen, Yijun Wang, Shangkai Gao, Tzyy-Ping Jung and Xiaorong Gao,"Filter bank canonical correlation analysis for implementing a high-speed SSVEP-based brain–computer interface,"Journal of Neural Engineering, Volume 12, Number 4