本發(fā)明涉及信號(hào)處理技術(shù)領(lǐng)域,特別涉及一種基于二維空間的時(shí)間序列樣本熵的計(jì)算方法及系統(tǒng)。
背景技術(shù):
在對(duì)時(shí)間序列信號(hào)進(jìn)行處理時(shí),時(shí)域的非線性分析方法是腦電信號(hào)處理的重要方面,包括相關(guān)維數(shù)、李雅普諾夫指數(shù)、小波熵、近似熵、樣本熵等都是比較常用的方法。振幅和周期是表達(dá)波形振動(dòng)的兩個(gè)重要變量,現(xiàn)有的腦電可以依據(jù)波形的振動(dòng)特性從時(shí)間序列中提取振幅和周期特征,并結(jié)合兩種特征進(jìn)行非線性分析。
樣本熵是對(duì)振動(dòng)時(shí)間序列的復(fù)雜度進(jìn)行分析的一種指標(biāo),表征了時(shí)間序列中新模式的生成概率。以往的樣本熵研究中,通常只將波形作為一維序列進(jìn)行分析,而不能對(duì)波形的振幅-周期綜合特性進(jìn)行分析。
技術(shù)實(shí)現(xiàn)要素:
鑒于上述問題,提出了本發(fā)明以便提供一種克服上述問題或者至少部分地解決上述問題的一種基于二維空間的時(shí)間序列樣本熵的計(jì)算方法及系統(tǒng)。
依據(jù)本發(fā)明的一個(gè)方面,提供了一種基于二維空間的時(shí)間序列樣本熵的計(jì)算方法,所述方法包括:
S1:提取時(shí)間序列的極值點(diǎn)和極值點(diǎn)的時(shí)間位置,以生成極值序列P={p(i)}及時(shí)間位置序列T={t(i)},其中,i=1,2,3,…,n;
S2:生成序列Z1=Y(jié)1-X1,去除序列Z1的首項(xiàng)和尾項(xiàng),以生成單調(diào)振幅序列A={a(k)},并生成序列Z2=Y(jié)2-X2,去除序列Z2的首項(xiàng)和尾項(xiàng),以生成單調(diào)周期序列C={c(k)},其中,k=1,2,3,…,n-1,X1={0,p(i)},Y1={p(i),0},X1={0,t(i)},Y1={t(i),0};
S3:將序列A和序列C組成向量序列O={(a(k),c(k))},其中k=1,2,3,…,n-1,振幅序列為A={a(k)},周期序列為C={c(k)},k=1,2,3,…,n-1;
S4:從序列A中提取n-m+1個(gè)m維向量A(u)={a(u),a(u+1),a(u+2),…,a(u+m-1)},從C中提取n-m+1個(gè)m維向量C(u)={c(u),c(u+1),c(u+2),…,c(u+m-1)},并由A(u)和C(u)組成向量序列Q(u)={(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))};
S5:計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離;
S6:計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)];
S7:設(shè)相似容限為R,統(tǒng)計(jì)出D[Qm(u),Qm(v)]<R的個(gè)數(shù)num{D[Qm(u),Qm(v)]<R},計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R);
S8:計(jì)算所有向量Bum(R)的平均值Bm(R);
S9:將維數(shù)由m加到m+1,重復(fù)上述步驟S4~S8,得到Bm+1(R);
S10:根據(jù)所述Bm(R)和Bm+1(R)計(jì)算樣本熵。
可選地,步驟S5中,通過下式計(jì)算計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,
JD((a(u),c(u)),(a(v),c(v)))=1-J((a(u),c(u)),(a(v),c(v)))
其中,JD((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,J((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系數(shù),所述Jaccard相似系數(shù)的計(jì)算公式為ΔB為容納(a(u),c(u))的最小矩形B(u)與容納(a(v),c(v))的最小矩形B(v)的重合面積。
可選地,步驟S6中,通過下式計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)],
D[Qm(u),Qm(v)]=max(JD((a(u+x),c(u+x)),(a(v+x),c(v+x))))
其中,0≤x≤m-1,1≤u≤n-m+1,1≤v≤n-m+1。
可選地,步驟S7中,通過下式計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R),
其中,1≤v≤n-m+1且v≠u。
可選地,步驟S10中,根據(jù)所述Bm(R)和Bm+1(R)通過下式計(jì)算樣本熵,
D2SEn(m,R,n)=-ln[Bm+1(R)/Bm(R)]
其中,D2SEn(m,R,n)為樣本熵。
依據(jù)本發(fā)明的另一個(gè)方面,提供了一種基于二維空間的時(shí)間序列樣本熵的計(jì)算系統(tǒng),所述系統(tǒng)包括:
數(shù)據(jù)提取單元,用于提取時(shí)間序列的極值點(diǎn)和極值點(diǎn)的時(shí)間位置,以生成極值序列P={p(i)}及時(shí)間位置序列T={t(i)},其中,i=1,2,3,…,n;
序列生成單元,用于生成序列Z1=Y(jié)1-X1,去除序列Z1的首項(xiàng)和尾項(xiàng),以生成單調(diào)振幅序列A={a(k)},并生成序列Z2=Y(jié)2-X2,去除序列Z2的首項(xiàng)和尾項(xiàng),以生成單調(diào)周期序列C={c(k)},其中,k=1,2,3,…,n-1,X1={0,p(i)},Y1={p(i),0},X1={0,t(i)},Y1={t(i),0};
序列組成單元,用于將序列A和序列C組成向量序列O={(a(k),c(k))},其中k=1,2,3,…,n-1,振幅序列為A={a(k)},周期序列為C={c(k)},k=1,2,3,…,n-1;
向量提取單元,用于從序列A中提取n-m+1個(gè)m維向量A(u)={a(u),a(u+1),a(u+2),…,a(u+m-1)},從C中提取n-m+1個(gè)m維向量C(u)={c(u),c(u+1),c(u+2),…,c(u+m-1)},并由A(u)和C(u)組成向量序列Q(u)={(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))};
元素計(jì)算單元,用于計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離;
向量計(jì)算單元,用于計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)];
比值計(jì)算單元,用于設(shè)相似容限為R,統(tǒng)計(jì)出D[Qm(u),Qm(v)]<R的個(gè)數(shù)num{D[Qm(u),Qm(v)]<R},計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R);
均值計(jì)算單元,用于計(jì)算所有向量Bum(R)的平均值Bm(R);
維數(shù)調(diào)整單元,用于將維數(shù)由m加到m+1,得到Bm+1(R);
樣本熵計(jì)算單元,用于根據(jù)所述Bm(R)和Bm+1(R)計(jì)算樣本熵。
可選地,所述元素計(jì)算單元通過下式計(jì)算計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,
JD((a(u),c(u)),(a(v),c(v)))=1-J((a(u),c(u)),(a(v),c(v)))
其中,JD((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,J((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系數(shù),所述Jaccard相似系數(shù)的計(jì)算公式為ΔB為容納(a(u),c(u))的最小矩形B(u)與容納(a(v),c(v))的最小矩形B(v)的重合面積。
可選地,所述向量計(jì)算單元通過下式計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)],
D[Qm(u),Qm(v)]=max(JD((a(u+x),c(u+x)),(a(v+x),c(v+x))))
其中,0≤x≤m-1,1≤u≤n-m+1,1≤v≤n-m+1。
可選地,所述比值計(jì)算單元通過下式計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R),
其中,1≤v≤n-m+1且v≠u。
可選地,所述樣本熵計(jì)算單元根據(jù)所述Bm(R)和Bm+1(R)通過下式計(jì)算樣本熵,
D2SEn(m,R,n)=-ln[Bm+1(R)/Bm(R)]
其中,D2SEn(m,R,n)為樣本熵。
本發(fā)明通過各步驟之間的配合,實(shí)現(xiàn)了樣本熵的計(jì)算,能夠應(yīng)用于腦電復(fù)雜度計(jì)算,也可應(yīng)用于其他存在局部極值點(diǎn)的振動(dòng)序列或波形的復(fù)雜度計(jì)算。振動(dòng)熵可作為對(duì)波形進(jìn)行模式識(shí)別的特征指標(biāo),對(duì)不同復(fù)雜度波形進(jìn)行分類。在信號(hào)處理時(shí),信號(hào)中混有白噪聲的振幅和周期是在一定范圍內(nèi)呈混沌的正態(tài)分布的,而信號(hào)的波形較為規(guī)律,所以本發(fā)明也可用于對(duì)信號(hào)中噪聲的識(shí)別和剔除。
附圖說明
圖1是本發(fā)明一種實(shí)施方式的基于二維空間的時(shí)間序列樣本熵的計(jì)算方法的流程圖;
圖2是樣本熵與相似容限R之間的關(guān)系圖(其中,星號(hào)表示符合S1<S2<S3<S4<S5的關(guān)系);
圖3是在相似容限取一定值時(shí),仿真波形的樣本熵示意圖;
圖4是本發(fā)明一種實(shí)施方式的基于二維空間的時(shí)間序列樣本熵的計(jì)算系統(tǒng)的結(jié)構(gòu)框圖。
具體實(shí)施方式
下面結(jié)合附圖和實(shí)施例,對(duì)本發(fā)明的具體實(shí)施方式作進(jìn)一步詳細(xì)描述。以下實(shí)施例用于說明本發(fā)明,但不用來限制本發(fā)明的范圍。
圖1是本發(fā)明一種實(shí)施方式的基于二維空間的時(shí)間序列樣本熵的計(jì)算方法的流程圖;參照圖1,所述方法包括:
S1:提取時(shí)間序列的極值點(diǎn)和極值點(diǎn)的時(shí)間位置,以生成極值序列P={p(i)}及時(shí)間位置序列T={t(i)},其中,i=1,2,3,…,n;
S2:生成序列Z1=Y(jié)1-X1,去除序列Z1的首項(xiàng)和尾項(xiàng),以生成單調(diào)振幅序列A={a(k)},并生成序列Z2=Y(jié)2-X2,去除序列Z2的首項(xiàng)和尾項(xiàng),以生成單調(diào)周期序列C={c(k)},其中,k=1,2,3,…,n-1,X1={0,p(i)},Y1={p(i),0},X1={0,t(i)},Y1={t(i),0};
S3:將序列A和序列C組成向量序列O={(a(k),c(k))},其中k=1,2,3,…,n-1,振幅序列為A={a(k)},周期序列為C={c(k)},k=1,2,3,…,n-1;
S4:從序列A中提取n-m+1個(gè)m維向量A(u)={a(u),a(u+1),a(u+2),…,a(u+m-1)},從C中提取n-m+1個(gè)m維向量C(u)={c(u),c(u+1),c(u+2),…,c(u+m-1)},并由A(u)和C(u)組成向量序列Q(u)={(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))};
需要說明的是,A(u)和C(u)為某段振動(dòng)序列中的兩個(gè)表示振動(dòng)特征的向量。
S5:計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離;
可理解的是,在振幅-周期的歐式空間中,(a(u),c(u))與(a(v),c(v))的Jaccard相似系數(shù)為:容納(a(u),c(u))的最小矩形B(u)與容納(a(v),c(v))的最小矩形B(v)的重合面積ΔB,與此時(shí)B(u)和B(v)共占空間面積的比值為
其中ΔB=min(a(u),a(v))×min(c(u),c(v)),是矩形B(u)和B(v)重合面積的最大值。因此,向量(a(u),c(u))與(a(v),c(v))的Jaccard距離計(jì)算公式為:
JD((a(u),c(u)),(a(v),c(v)))=1-J((a(u),c(u)),(a(v),c(v)))
S6:計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)];
其中,向量序列Q(u)與Q(v)的距離表示為:Q(u)與Q(v)序列對(duì)應(yīng)向量元素Jaccard距離的最大值。
D[Qm(u),Qm(v)]=max(JD((a(u+x),c(u+x)),(a(v+x),c(v+x))))
其中,0≤x≤m-1,1≤u≤n-m+1,1≤v≤n-m+1。
S7:設(shè)相似容限為R,統(tǒng)計(jì)出D[Qm(u),Qm(v)]<R的個(gè)數(shù)num{D[Qm(u),Qm(v)]<R},計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R);
在具體實(shí)現(xiàn)中,通過下式計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R),
其中,1≤v≤n-m+1且v≠u。
S8:計(jì)算所有向量Bum(R)的平均值Bm(R);
在具體實(shí)現(xiàn)中,通過下式計(jì)算所有向量Bum(R)的平均值Bm(R),
S9:將維數(shù)由m加到m+1,重復(fù)上述步驟S4~S8,得到Bm+1(R);
其中,
S10:根據(jù)所述Bm(R)和Bm+1(R)計(jì)算樣本熵。
在具體實(shí)現(xiàn)中,根據(jù)所述Bm(R)和Bm+1(R)通過下式計(jì)算樣本熵,
D2SEn(m,R,n)=-ln[Bm+1(R)/Bm(R)]
其中,D2SEn(m,R,n)為樣本熵。
可理解的是,為了便于對(duì)上述實(shí)施方式的方法進(jìn)行驗(yàn)證,本實(shí)施方式中,在步驟S1之前,可依據(jù)波形特點(diǎn)由正弦波構(gòu)造仿真波形,子波函數(shù)為:
其中,當(dāng)i=1,3,5,…,2n-1時(shí),t∈[0,π),當(dāng)i=2,4,6,…,2n時(shí),t∈[π,2π),Ti和Ai分別表示第i個(gè)正弦子波的周期和振幅。
再設(shè)時(shí)間序列S的子波數(shù)為n,根據(jù)條件隨機(jī)生成隨機(jī)序列A={A1,A2,A3,…,An}和周期序列T={T1,T2,T3,…,Tn},由A和T組成振動(dòng)向量序列O={(A1,T1),(A2,T2),…,(An,Tn)},然后由O中的各元素生成余弦子波序列,將序列中的子波收尾相連構(gòu)成仿真時(shí)間序列S。
然后,構(gòu)造五種波形S1、S2、S3、S4、S5進(jìn)行實(shí)驗(yàn)。波形的復(fù)雜度關(guān)系為S1<S2<S3<S4<S5。
計(jì)算上述五種波形的樣本熵,振動(dòng)熵的計(jì)算方法則按照上述步驟S1~S10進(jìn)行計(jì)算,當(dāng)然,在進(jìn)行計(jì)算前,先對(duì)上述模擬信號(hào)進(jìn)行濾波,將其分解為若干子頻段信號(hào)。
由圖2可知,D2SEn的值與相似容限R的取值有關(guān),D2SEn隨相似容限R的增大而減小,僅當(dāng)相似容限的取值在一定范圍內(nèi)時(shí),仿真波形D2SEn值的大小關(guān)系才與預(yù)設(shè)條件相同。
圖3指示在相似容限取一定值時(shí)仿真波形的二維樣本熵示意圖。分別對(duì)R=0.05、0.1、0.15時(shí)的D2SEn值進(jìn)行方差分析。結(jié)果表明,在上述三種條件下,D2SEn值在五種波形間的組效應(yīng)顯著(R=0.05:F=3245.084,P<0.00001;R=0.1:F=3245.084,P<0.00001;R=0.15:F=3923.475,P<0.00001)。組間兩兩比較時(shí),每兩組波形間的D2SEn值均存在顯著差異(P<0.00001)。由圖2可知,當(dāng)相似容限R的取值在一定范圍內(nèi)時(shí),D2SEn值可以有效反映波形的復(fù)雜度特征。
本實(shí)施方式的方法可應(yīng)用于腦電復(fù)雜度計(jì)算,也可應(yīng)用于其他存在局部極值點(diǎn)的振動(dòng)序列或波形的復(fù)雜度計(jì)算。振動(dòng)熵可作為對(duì)波形進(jìn)行模式識(shí)別的特征指標(biāo),對(duì)不同復(fù)雜度波形進(jìn)行分類。在信號(hào)處理時(shí),信號(hào)中混有白噪聲的振幅和周期是在一定范圍內(nèi)呈混沌的正態(tài)分布的,而信號(hào)的波形較為規(guī)律,所以本實(shí)施方式的方法也可用于對(duì)信號(hào)中噪聲的識(shí)別和剔除。
圖4是本發(fā)明一種實(shí)施方式的基于二維空間的時(shí)間序列樣本熵的計(jì)算系統(tǒng)的結(jié)構(gòu)框圖;參照圖4,所述系統(tǒng)包括:
數(shù)據(jù)提取單元401,用于提取時(shí)間序列的極值點(diǎn)和極值點(diǎn)的時(shí)間位置,以生成極值序列P={p(i)}及時(shí)間位置序列T={t(i)},其中,i=1,2,3,…,n;
序列生成單元402,用于生成序列Z1=Y(jié)1-X1,去除序列Z1的首項(xiàng)和尾項(xiàng),以生成單調(diào)振幅序列A={a(k)},并生成序列Z2=Y(jié)2-X2,去除序列Z2的首項(xiàng)和尾項(xiàng),以生成單調(diào)周期序列C={c(k)},其中,k=1,2,3,…,n-1,X1={0,p(i)},Y1={p(i),0},X1={0,t(i)},Y1={t(i),0};
序列組成單元403,用于將序列A和序列C組成向量序列O={(a(k),c(k))},其中k=1,2,3,…,n-1,振幅序列為A={a(k)},周期序列為C={c(k)},k=1,2,3,…,n-1;
向量提取單元404,用于從序列A中提取n-m+1個(gè)m維向量A(u)={a(u),a(u+1),a(u+2),…,a(u+m-1)},從C中提取n-m+1個(gè)m維向量C(u)={c(u),c(u+1),c(u+2),…,c(u+m-1)},并由A(u)和C(u)組成向量序列Q(u)={(a(u),c(u)),(a(u+1),c(u+1)),(a(u+2),c(u+2)),…,(a(u+m-1),c(u+m-1))};
元素計(jì)算單元405,用于計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離;
向量計(jì)算單元406,用于計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)];
比值計(jì)算單元407,用于設(shè)相似容限為R,統(tǒng)計(jì)出D[Qm(u),Qm(v)]<R的個(gè)數(shù)num{D[Qm(u),Qm(v)]<R},計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R);
均值計(jì)算單元408,用于計(jì)算所有向量Bum(R)的平均值Bm(R);
維數(shù)調(diào)整單元409,用于將維數(shù)由m加到m+1,得到Bm+1(R);
樣本熵計(jì)算單元410,用于根據(jù)所述Bm(R)和Bm+1(R)計(jì)算樣本熵。
可選地,所述元素計(jì)算單元通過下式計(jì)算計(jì)算向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,
JD((a(u),c(u)),(a(v),c(v)))=1-J((a(u),c(u)),(a(v),c(v)))
其中,JD((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard距離,J((a(u),c(u)),(a(v),c(v)))為向量元素(a(u),c(u))和(a(v),c(v))的Jaccard相似系數(shù),所述Jaccard相似系數(shù)的計(jì)算公式為ΔB為容納(a(u),c(u))的最小矩形B(u)與容納(a(v),c(v))的最小矩形B(v)的重合面積。
可選地,所述向量計(jì)算單元通過下式計(jì)算向量序列Q(u)與Q(v)的距離D[Qm(u),Qm(v)],
D[Qm(u),Qm(v)]=max(JD((a(u+x),c(u+x)),(a(v+x),c(v+x))))
其中,0≤x≤m-1,1≤u≤n-m+1,1≤v≤n-m+1。
可選地,所述比值計(jì)算單元通過下式計(jì)算num{D[Qm(u),Qm(v)]<R}與所有距離總數(shù)的比值Bum(R),
其中,1≤v≤n-m+1且v≠u。
可選地,所述樣本熵計(jì)算單元根據(jù)所述Bm(R)和Bm+1(R)通過下式計(jì)算樣本熵,
D2SEn(m,R,n)=-ln[Bm+1(R)/Bm(R)]
其中,D2SEn(m,R,n)為樣本熵。
以上實(shí)施方式僅用于說明本發(fā)明,而并非對(duì)本發(fā)明的限制,有關(guān)技術(shù)領(lǐng)域的普通技術(shù)人員,在不脫離本發(fā)明的精神和范圍的情況下,還可以做出各種變化和變型,因此所有等同的技術(shù)方案也屬于本發(fā)明的范疇,本發(fā)明的專利保護(hù)范圍應(yīng)由權(quán)利要求限定。