(57)【要約】      （修正有）【課題】１つの原信号を複数のチャネルで観測した観測信号から、各チャネルの信号雑音比や原信号を推定する信号処理装置等を提供する。【解決手段】信号処理装置１０１において、入力受付部１０２はチャネルiの時刻tの観測信号xi(t)を受け付け、相関計算部１０４は、チャネルi,jの標本相関ri,jを計算し、初期値設定部１０５は、計算結果からチャネルsの原信号と雑音の混合比asの初期値として〔rs,tru,s/rt,u〕1/2を設定し、更新値設定部１０６は混合比を交互最小自乗法により混合比aiを更新し、収束したら、混合比出力部１０７は、各混合比を出力する。ベクトル計算部１０９は、得られた混合比と観測信号に基づいて重みベクトルを計算し、原信号推定部１１０は、計算された重みベクトルに基づいて観測信号の重み付き平均を求めて、これを原信号の推定信号とする。【選択図】図１
2008293163.tif 000002
【特許請求の範囲】
【請求項１】
時刻t (t = 1,…,T)における原信号g(t)を、n個のチャネルのうち、チャネルi (i = 1,…,n)で観測した観測信号xi(t)の入力を受け付ける入力受付部、
当該原信号g(t)について、
Σt=1T g(t)2/T = 1
と仮定し、雑音pi(t)が存在すると仮定して、当該観測信号xi(t)を、
aig(t) + pi(t)
とモデル化するための混合比aiの推定値が記憶される記憶部、
当該チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関
ri,j = rj,i = Σt=1T xi(t)xj(t)/〔T-1〕
を計算する相関計算部、
整数s (s = 1,…,n)のそれぞれについて、整数t,u (t = 1,…,n；u = 1,…,n；s≠t,t≠u,u≠s)を選択し、前記記憶部に記憶される混合比asの初期値として
〔rs,tru,s/rt,u〕1/2
を設定する初期値設定部、
前記記憶部に記憶される混合比a1,…,anが収束するまで、整数s (s = 1,…,n)のそれぞれについて、前記記憶部に記憶される混合比asの更新値として、
Σj≠s,j=1n ajrs,j/Σk≠s,k=1n ak2
を設定する更新値設定部、
前記記憶部に記憶される混合比a1,…,anが収束すると、前記記憶部に記憶される混合比a1,…,anを出力する混合比出力部
を備えることを特徴とする信号処理装置。
【請求項２】
請求項１に記載の信号処理装置であって、
前記初期値設定部は、整数s = 1,…,nのそれぞれに対して、チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関において、s≠i,s≠jを満たし、ri,jの値が最大となるi,jを選択し、当該i,jをそれぞれ整数t,uとして選択する
ことを特徴とする信号処理装置。
【請求項３】
請求項１に記載の信号処理装置であって、
前記初期値設定部は、整数s = 1,…,nのそれぞれに対して、整数t,uの組をランダムに所定の個数選択し、当該組のそれぞれについて計算される
〔rs,tru,s/rt,u〕1/2
の中央値を、当該混合比asの初期値とする
ことを特徴とする信号処理装置。
【請求項４】
請求項１から３のいずれか１項に記載の信号処理装置であって、
当該チャネルi (i=1,…,n)の信号雑音比ξiを、
ξi = ai2/〔Σt=1T xi(t)2/T - ai2〕
により計算する信号雑音比推定部
をさらに備えることを特徴とする信号処理装置。
【請求項５】
請求項１から４のいずれか１項に記載の信号処理装置であって、
前記出力された混合比から、チャネルiに含まれる共通信号のパワーPiを
Pi = ai2
により計算する信号強度推定部
をさらに備えることを特徴とする信号処理装置。
【請求項６】
請求項１から５のいずれか１項に記載の信号処理装置であって、
行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対する一般固有値問題
S w = λ(Q - B)w
の最大固有値解λに対応付けられる固有ベクトルw = [w1,…,wn]†を計算するベクトル計算部、
原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する原信号推定部
をさらに備えることを特徴とする信号処理装置。
【請求項７】
請求項１から５のいずれか１項に記載の信号処理装置であって、
行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対して、ベクトルw = [w1,…,wn]†を、
w = (S + Q - B)-1a
により計算するベクトル計算部、
原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する原信号推定部
をさらに備えることを特徴とする信号処理装置。
【請求項８】
入力受付部、記憶部、相関計算部、初期値設定部、更新値設定部、混合比出力部を有する信号処理装置にて実行される信号処理方法であって、
前記入力受付部が、時刻t (t = 1,…,T)における原信号g(t)を、n個のチャネルのうち、チャネルi (i = 1,…,n)で観測した観測信号xi(t)の入力を受け付ける入力受付工程を備え、
前記記憶部には、当該原信号g(t)について、
Σt=1T g(t)2/T = 1
と仮定し、雑音pi(t)が存在すると仮定して、当該観測信号xi(t)を、
aig(t) + pi(t)
とモデル化するための混合比aiの推定値が記憶され、
前記相関計算部が、当該チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関
ri,j = rj,i = Σt=1T xi(t)xj(t)/〔T-1〕
を計算する相関計算工程、
前記初期値設定部が、整数s (s = 1,…,n)のそれぞれについて、整数t,u (t = 1,…,n；u = 1,…,n；s≠t,t≠u,u≠s)を選択し、前記記憶部に記憶される混合比asの初期値として
〔rs,tru,s/rt,u〕1/2
を設定する初期値設定工程、
前記更新値設定部が、前記記憶部に記憶される混合比a1,…,anが収束するまで、整数s (s = 1,…,n)のそれぞれについて、前記記憶部に記憶される混合比asの更新値として、
Σj≠s,j=1n ajrs,j/Σk≠s,k=1n ak2
を設定する更新値設定工程、
前記混合比出力部が、前記記憶部に記憶される混合比a1,…,anが収束すると、前記記憶部に記憶される混合比a1,…,anを出力する混合比出力工程
を備えることを特徴とする信号処理方法。
【請求項９】
請求項８に記載の信号処理方法であって、
前記初期値設定工程では、整数s = 1,…,nのそれぞれに対して、チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関において、s≠i,s≠jを満たし、ri,jの値が最大となるi,jを選択し、当該i,jをそれぞれ整数t,uとして選択する
ことを特徴とする信号処理方法。
【請求項１０】
請求項８に記載の信号処理方法であって、
前記初期値設定工程では、整数s = 1,…,nのそれぞれに対して、整数t,uの組をランダムに所定の個数選択し、当該組のそれぞれについて計算される
〔rs,tru,s/rt,u〕1/2
の中央値を、当該混合比asの初期値とする
ことを特徴とする信号処理方法。
【請求項１１】
請求項８から１０のいずれか１項に記載の信号処理方法であって、前記信号処理装置は、信号雑音比推定部をさらに有し、
前記信号雑音比推定部が、当該チャネルi (i=1,…,n)の信号雑音比ξiを、
ξi = ai2/〔Σt=1T xi(t)2/T - ai2〕
により計算する信号雑音比推定工程
をさらに備えることを特徴とする信号処理装置。
【請求項１２】
請求項８から１１のいずれか１項に記載の信号処理方法であって、前記信号処理装置は、信号強度推定部をさらに有し、
前記信号強度推定部が、前記出力された混合比から、チャネルiに含まれる共通信号のパワーPiを
Pi = ai2
により計算する信号強度推定工程
をさらに備えることを特徴とする信号処理方法。
【請求項１３】
請求項８から１２のいずれか１項に記載の信号処理方法であって、前記信号処理装置は、ベクトル計算部、原信号推定部をさらに有し、
前記ベクトル計算部が、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対する一般固有値問題
S w = λ(Q - B)w
の最大固有値解λに対応付けられる固有ベクトルw = [w1,…,wn]†を計算するベクトル計算工程、
前記原信号推定部が、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する原信号推定工程
をさらに備えることを特徴とする信号処理方法。
【請求項１４】
請求項８から１２のいずれか１項に記載の信号処理方法であって、
前記信号処理装置は、ベクトル計算部、原信号推定部をさらに有し、
前記ベクトル計算部が、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対して、ベクトルw = [w1,…,wn]†を、
w = (S + Q - B)-1a
により計算するベクトル計算工程、
前記原信号推定部が、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する原信号推定工程
をさらに備えることを特徴とする信号処理方法。
【請求項１５】
コンピュータを、請求項１から７のいずれか１項に記載の信号処理装置の各部として機能させることと特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、１つの原信号を複数のチャネルで観測した観測信号から、各チャネルの信号雑音比や原信号を推定するのに好適な信号処理装置、信号処理方法、ならびに、これをコンピュータにて実現するプログラムに関する。
【背景技術】
【０００２】
従来から、通信情報技術、医療技術、天文気象技術環境技術等、さまざまな技術分野では、１つの信号を複数のセンサーで観測するマルチセンサーシステムが利用されている。
【０００３】
たとえば、脳波図（electroencephalogram；ＥＥＧ）を作成するためのＥＥＧデバイスでは、数百の電極を利用することもある。そして、ＥＥＧ信号には、極めて大きい雑音を含むことが多く、信号雑音比が−２０ｄＢ未満となることも多い。したがって、脳の活動を表すクリーンな信号を、雑音が含まれた複数の信号から抽出する必要がある。
【０００４】
一方で、下記の文献に開示されるように、信号分離の技術としては、独立成分分析（Independent Component Analysis；ＩＣＡ）が広く使われている。
【特許文献１】特願２００１−３３９３９６
【非特許文献１】C.E.Davila and M.S.Mobin，Weighted averaging of evoked potentials，IEEE transactions on Biomedical Engineering，Vol.39，No.4，pp.338-345，１９９２年
【０００５】
ここで、[特許文献１]では、複数の化学物質の量の変化の時系列を原信号と考え、独立成分分析を施すことにより、化学物質の分類を行う技術が提案されている。
【０００６】
しかしながら、ＩＣＡは、ガウス分布（Gaussian distribution）に近い雑音であるガウス信号（Gaussian signal)源が複数あるようなモデルにおいては適用ができない。
【０００７】
さらに、ＩＣＡは、観測対象となる信号が混合されるモデルを採用するため、上記のように、原信号の数が観測信号の数より多い状況下では、ＩＣＡをそのまま適用することはできない。
【０００８】
一方で、[非特許文献１]では、重み付き平均により原信号を求める手法が提案されている。しかしながら、センサーの数が多くならないと収束せず、センサー数が有限の場合には精度が低い、という悪い性質を持ち、性能も満足のできるものではない。
【発明の開示】
【発明が解決しようとする課題】
【０００９】
したがって、各センサーにて得られた測定信号を総合して、各センサーごとに原信号と雑音との混合比や信号雑音比を適切に推定したり、原信号そのものを推定する技術が強く望まれている。
【００１０】
本発明は、以上のような課題を解決するためのもので、１つの原信号を複数のチャネルで観測した観測信号から、各チャネルの信号雑音比や原信号を推定するのに好適な信号処理装置、信号処理方法、ならびに、これをコンピュータにて実現するプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明の第１の観点に係る信号処理装置は、入力受付部、記憶部、相関計算部、初期値設定部、更新値設定部、混合比出力部を備え、以下のように構成する。
【００１２】
すなわち、入力受付部は、時刻t (t = 1,…,T)における原信号g(t)を、n個のチャネルのうち、チャネルi (i = 1,…,n)で観測した観測信号xi(t)の入力を受け付ける。
【００１３】
一方、記憶部には、当該原信号g(t)について、
Σt=1T g(t)2/T = 1
と仮定し、雑音pi(t)が存在すると仮定して、当該観測信号xi(t)を、
aig(t) + pi(t)
とモデル化するための混合比aiの推定値が記憶される。
【００１４】
さらに、相関計算部は、当該チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関
ri,j = rj,i = Σt=1T xi(t)xj(t)/〔T-1〕
を計算する。
【００１５】
そして、初期値設定部は、整数s (s = 1,…,n)のそれぞれについて、整数t,u (t = 1,…,n；u = 1,…,n；s≠t,t≠u,u≠s)を選択し、記憶部に記憶される混合比asの初期値として
〔rs,tru,s/rt,u〕1/2
を設定する。
【００１６】
一方、更新値設定部は、記憶部に記憶される混合比a1,…,anが収束するまで、整数s (s = 1,…,n)のそれぞれについて、記憶部に記憶される混合比asの更新値として、
Σj≠s,j=1n ajrs,j/Σk≠s,k=1n ak2
を設定する。
【００１７】
さらに、混合比出力部は、記憶部に記憶される混合比a1,…,anが収束すると、記憶部に記憶される混合比a1,…,anを出力する。
【００１８】
また、本発明の信号処理装置において、初期値設定部は、整数s = 1,…,nのそれぞれに対して、チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関において、s≠i,s≠jを満たし、ri,jの値が最大となるi,jを選択し、当該i,jをそれぞれ整数t,uとして選択するように構成することができる。
【００１９】
また、本発明の信号処理装置において、初期値設定部は、整数s = 1,…,nのそれぞれに対して、整数t,uの組をランダムに所定の個数選択し、当該組のそれぞれについて計算される
〔rs,tru,s/rt,u〕1/2
の中央値を、当該混合比asの初期値とするように構成することができる。
【００２０】
また、本発明の信号処理装置は、信号雑音比推定部をさらに備え、信号雑音比推定部は、当該チャネルi (i=1,…,n)の信号雑音比ξiを、
ξi = ai2/〔Σt=1T xi(t)2/T - ai2〕
により計算するように構成することができる。
【００２１】
また、本発明の信号処理装置は、信号強度推定部をさらに備え、信号強度推定部は、出力された混合比から、チャネルiに含まれる共通信号のパワーPiを
Pi = ai2
により計算するように構成することができる。
【００２２】
また、本発明の信号処理装置は、ベクトル計算部、原信号推定部をさらに備え、以下のように構成することができる。
【００２３】
すなわち、ベクトル計算部は、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対する一般固有値問題
S w = λ(Q - B)w
の最大固有値解λに対応付けられる固有ベクトルw = [w1,…,wn]†を計算する。
【００２４】
一方、原信号推定部は、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する。
【００２５】
また、本発明の信号処理装置は、以下のように構成することができる。
【００２６】
すなわち、ベクトル計算部は、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対して、ベクトルw = [w1,…,wn]†を、
w = (S + Q - B)-1a
により計算する。
【００２７】
一方、原信号推定部は、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する。
【００２８】
本発明のその他の観点に係る信号処理方法は、入力受付部、記憶部、相関計算部、初期値設定部、更新値設定部、混合比出力部を有する信号処理装置にて実行され、入力受付工程、相関計算工程、初期値設定工程、更新値設定工程、混合比出力工程を備え、以下のように構成する。
【００２９】
すなわち、入力受付工程では、入力受付部が、時刻t (t = 1,…,T)における原信号g(t)を、n個のチャネルのうち、チャネルi (i = 1,…,n)で観測した観測信号xi(t)の入力を受け付ける。
【００３０】
一方、記憶部には、当該原信号g(t)について、
Σt=1T g(t)2/T = 1
と仮定し、雑音pi(t)が存在すると仮定して、当該観測信号xi(t)を、
aig(t) + pi(t)
とモデル化するための混合比aiの推定値が記憶される。
【００３１】
さらに、相関計算工程では、相関計算部が、当該チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関
ri,j = rj,i = Σt=1T xi(t)xj(t)/〔T-1〕
を計算する。
【００３２】
そして、初期値設定工程では、初期値設定部が、整数s (s = 1,…,n)のそれぞれについて、整数t,u (t = 1,…,n；u = 1,…,n；s≠t,t≠u,u≠s)を選択し、記憶部に記憶される混合比asの初期値として
〔rs,tru,s/rt,u〕1/2
を設定する。
【００３３】
一方、更新値設定工程では、更新値設定部が、記憶部に記憶される混合比a1,…,anが収束するまで、整数s (s = 1,…,n)のそれぞれについて、記憶部に記憶される混合比asの更新値として、
Σj≠s,j=1n ajrs,j/Σk≠s,k=1n ak2
を設定する。
【００３４】
さらに、混合比出力工程では、混合比出力部が、記憶部に記憶される混合比a1,…,anが収束すると、記憶部に記憶される混合比a1,…,anを出力する。
【００３５】
また、本発明の信号処理方法において、初期値設定工程では、整数s = 1,…,nのそれぞれに対して、チャネルi,j (i = 1,…,n；j = 1,…,n；i≠j)の標本相関において、s≠i,s≠jを満たし、ri,jの値が最大となるi,jを選択し、当該i,jをそれぞれ整数t,uとして選択するように構成することができる。
【００３６】
また、本発明の信号処理方法において、初期値設定工程では、整数s = 1,…,nのそれぞれに対して、整数t,uの組をランダムに所定の個数選択し、当該組のそれぞれについて計算される
〔rs,tru,s/rt,u〕1/2
の中央値を、当該混合比asの初期値とするように構成することができる。
【００３７】
また、本発明の信号処理方法において、信号処理装置は、信号雑音比推定部をさらに有し、信号雑音比推定部が、当該チャネルi (i=1,…,n)の信号雑音比ξiを、
ξi = ai2/〔Σt=1T xi(t)2/T - ai2〕
により計算する信号雑音比推定工程をさらに備えるように構成することができる。
【００３８】
また、本発明の信号処理装方法において、信号処理装置は、信号強度推定部をさらに有し、信号強度推定部が、出力された混合比から、チャネルiに含まれる共通信号のパワーPiを
Pi = ai2
により計算する信号強度推定工程をさらに備えるように構成することができる。
【００３９】
また、本発明の信号処理方法において、信号処理装置は、ベクトル計算部、原信号推定部をさらに有し、当該信号処理方法は、ベクトル計算工程原信号推定工程をさらに備え、以下のように構成することができる。
【００４０】
すなわち、ベクトル計算工程では、ベクトル計算部が、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対する一般固有値問題
S w = λ(Q - B)w
の最大固有値解λに対応付けられる固有ベクトルw = [w1,…,wn]†を計算する。
【００４１】
一方、原信号推定工程では、原信号推定部が、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する。
【００４２】
また、本発明の信号処理方法は、以下のように構成することができる。
【００４３】
すなわち、ベクトル計算工程では、ベクトル計算部が、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対して、ベクトルw = [w1,…,wn]†を、
w = (S + Q - B)-1a
により計算する。
【００４４】
一方、原信号推定工程では、原信号推定部が、原信号g(t)の推定信号h(t)を、
h(t) = Σi=1n wixi(t)
により計算する。
【００４５】
本発明のその他の観点に係るプログラムは、コンピュータを上記の信号処理装置の各部として機能させるように構成する。
【００４６】
また、本発明のプログラムは、コンパクトディスク、フレキシブルディスク、ハードディスク、光磁気ディスク、ディジタルビデオディスク、磁気テープ、半導体メモリ等のコンピュータ読取可能な情報記憶媒体に記録することができる。
【００４７】
上記プログラムは、プログラムが実行されるコンピュータやディジタル信号プロセッサとは独立して、コンピュータ通信網を介して配布・販売することができる。また、上記情報記憶媒体は、コンピュータやディジタル信号プロセッサとは独立して配布・販売することができる。
【発明の効果】
【００４８】
本発明によれば、１つの原信号を複数のチャネルで観測した観測信号から、各チャネルの信号雑音比や原信号を推定するのに好適な信号処理装置、信号処理方法、ならびに、これをコンピュータにて実現するプログラムを提供することができる。
【発明を実施するための最良の形態】
【００４９】
以下に本発明の実施形態を説明する。なお、以下に説明する実施形態は説明のためのものであり、本発明の範囲を制限するものではない。したがって、当業者であればこれらの各要素もしくは全要素をこれと均等なものに置換した実施形態を採用することが可能であるが、これらの実施形態も本発明の範囲に含まれる。
【実施例１】
【００５０】
図１は、本実施形態に係る信号処理装置の概要構成を示す模式図である。以下、本図を参照して説明する。
【００５１】
本図に示すように、本実施形態に係る信号処理装置１０１は、入力受付部１０２、記憶部１０３、相関計算部１０４、初期値設定部１０５、更新値設定部１０６、混合比出力部１０７を備える。
【００５２】
また、付加的な要素として、図示するように、信号雑音比推定部１０８、ベクトル計算部１０９、原信号推定部１１０を備えるように構成することができる。このほか、図示しない信号強度推定部を備えるように構成しても良い。
【００５３】
これらの各要素は、所定のプログラムをコンピュータ上で動作させることによって実現されるのが典型的であるが、それぞれの処理や演算は、専用ハードウェアの電気回路や電子回路、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等によって実現することも可能であり、これらの実施形態も本発明の範囲に含まれる。
【００５４】
以下、ＣＰＵ（Central Processing Unit）の制御化で、一時記憶領域としてＲＡＭ（Random Access Memory）を、不揮発性記憶領域としてＨＤ（Hard Disk）を使用して、コンピュータ上で信号処理装置１０１を実現する形態を主な実施例として説明する。
【００５５】
本発明では、時刻t (t = 1,…,T)における原信号g(t)を、n個のチャネルi (i = 1,…,n)で観測し、時刻tにおける観測信号xi(t)を得る。すなわち、xi(t)は、各センサーによる出力の時系列である。
【００５６】
また、チャネルiには、時刻tにおいて、雑音pi(t)があったものと考える。
【００５７】
また、
Σt=1T g(t)2/T = 1
と仮定する。この「１」は正規化のための単位と考えることができる。本手法においては、g(t)についての「定数倍」の冗長性を持っているため、このように仮定しても、一般性を失わない。
【００５８】
そして、各観測信号を、混合比aiを用いて、
xi(t) = aig(t) + pi(t)
とモデル化する。
【００５９】
図２は、本実施形態に係る信号処理装置にて実行される信号処理方法の制御の流れを示すフローチャートである。以下、本図を参照して説明する。
【００６０】
まず、ＣＰＵは、ＲＡＭ内に用意された各種の変数領域を確保する（ステップＳ２０１）。本実施形態において、記憶部１０３を構成するＲＡＭには、以下のような情報が記憶される。
（１）観測信号xi(t)。i,tを添字とする二次元配列にて実現するのが典型的である。
（２）混合比ai。aiを同定する途中の計算値、および同定されたaiの値を、iを添字とする一次元配列に記憶するのが典型的である。
（３）標本相関ri,jとその添字i,jの組。チャネルiとチャネルjの標本相関の計算値を保持するために使用するが、標本相関にはri,j = rj,iという性質があること、本実施形態では自己相関（i=jの場合）は利用しないこと、ri,jの大きさの順にi,jの組を選ぶ処理が含まれることから、これらの条件を利用して、i<jについてのiとjとri,jとの三つ組を要素として持つ一次元配列とするのが典型的である。
（４）推定されたチャネルiの信号雑音比ξi。iを添字とする一次元配列にて実現される。
（５）推定された原信号h(t)。tを添字とする一次元配列にて実現される。
（４）この他、繰り返し計算におけるカウンタ変数など、種々の一時的な値。以下では、理解を容易にするため、計算式全体を提示することとし、具体的な変数領域をどのように計算に用いるかについては、適宜説明を省略することとする。
【００６１】
なお、これらの値は、ＲＡＭではなくＨＤに記憶して適宜読み書きをすることとしても良い。また、これら以外の記憶装置を利用することも可能である。
【００６２】
ついで、入力受付部１０２は、コンピュータに接続されたセンサーから直接的に、あるいは、センサーの測定結果を保存したＨＤから間接的に、観測信号xi(t)の入力を受け付ける（ステップＳ２０２）。前者の場合には、ＣＰＵの制御の下コンピュータに接続されたインターフェースカードを経由してセンサーと通信することで、入力受付部１０２が実現される。後者の場合には、ＣＰＵの制御下で、ＨＤやＣＤ−ＲＯＭ（Compact Disk Read Only Memory）等に記録した測定結果を読み出すことで、入力受付部１０２が実現される。
【００６３】
本実施形態では、受け付けられた観測信号xi(t)の値は、上記のように、ＲＡＭ内に記憶される。
【００６４】
次に、相関計算部１０４が、整数i,j (i = 1,…,n；j = 1,…,n；i≠j)のそれぞれについて（ステップＳ２０３）、チャネルiとチャネルjの標本相関
ri,j = rj,i = Σt=1T xi(t)xj(t)/〔T-1〕
を計算して（ステップＳ２０４）、三組(i,j,ri,j)をＲＡＭに蓄積する（ステップＳ２０５）処理を繰り返す（ステップＳ２０６）。ここで、Σを用いた計算は単純なループ構造によって実現が可能である（以下同様）。
【００６５】
ri,jは、後述するように、チャネルiとチャネルjのアンサンブル相関の推定値に相当するものである。
【００６６】
上述のように、それぞれのi,jおよびri,jの組は、ＲＡＭ内に記憶される。
【００６７】
ついで、初期値設定部１０５は、上記の(i,j,ri,j)の三つ組の配列を、ri,jが大きい順にソートする（ステップＳ２０７）。
【００６８】
そして、整数s (s = 1,…,n)のそれぞれについて（ステップＳ２０８）、ステップＳ２０６においてソートされた三つ組の配列をri,jが大きい方から、i≠sかつj≠sであるものを走査検索して（ステップＳ２０９）、見つかったi,jをそれぞれ、整数t,uに代入し（ステップＳ２１０）、
〔rs,tru,s/rt,u〕1/2
を計算して（ステップＳ２１１）、この計算結果を混合比asの初期値として、ＲＡＭ内に代入する処理（ステップＳ２１２）を繰り返す（ステップＳ２１３）。
【００６９】
したがって、ＣＰＵは、ＲＡＭと共働して、初期値設定部１０５として機能する。
【００７０】
アンサンブル相関をE[・]と書くこととすると、チャネルi，jの観測値のアンサンブル相関は、
E[xi,xj] = E[aig+pi,ajg+pj]
【００７１】
と書くことができるが、雑音pi,piと原信号g(t)とは無相関であるから、i≠jについて、
E[xi,xj] = aiaj
が成立する。
【００７２】
一方、サンプル相関ri,jは、サンプル数Tが無限大に近付くと、サンプル相関E[xi,xj]に収束する。したがって、サンプル数が十分に大きければ、
aiaj = ri,j
とすることができる。
【００７３】
これから、互いに異なる整数s,t,uについて、
asatau = 〔rs,trt,uru,s〕1/2
であることがわかる。
【００７４】
これを、
atau = rt,u
で割れば、
as = 〔rs,tru,s/rt,u〕1/2
となる。
【００７５】
したがって、ステップＳ２１０において設定される値は、asの推定近似計算の初期値としてふさわしいものであることがわかる。
【００７６】
ここで、ステップＳ２０８において、rt,uとして、できるだけ値が大きいものを走査探索しているが、これは、rt,uが小さい値であると、推定の悪条件となるからである。
【００７７】
このほか、rt,uの選択の手法として、整数s = 1,…,nのそれぞれに対して、整数t,u (s≠t,t≠u,u≠s)の組をランダムに所定の個数選択し、当該組のそれぞれについて
〔rs,tru,s/rt,u〕1/2
を計算し、これらの計算結果の中央値を、当該混合比asの初期値とするように構成することも可能である。
【００７８】
上記の手法では、同じ組に対するi,jが繰り返し選択されることとなるが、当該i,jに何らかの測定誤差などが含まれることもありうる。本手法では、何個かサンプルを選んだときの中央値を用いるので、このような異常値を、より適切に避けることが可能である。もっとも、ほとんどの場合は上記の手法で十分である。
【００７９】
さて、更新値設定部１０６は、記憶部１０３に記憶される混合比a1,…,anが収束したか判断する（ステップＳ２１４）。
【００８０】
a1,…,anが更新されるごとに、当該更新における変化量の絶対値や自乗値を積算して、これが一定値未満となった場合に、収束した、と判定するのが典型的であり、したがって、ループの開始段階では、自動的に収束していないと判定されることになる。
【００８１】
収束していなければ（ステップＳ２１４；Ｎｏ）、以下の処理を繰り返す。
【００８２】
すなわち、整数s (s = 1,…,n)のそれぞれについて（ステップＳ２１５）、
Σj≠s,j=1n ajrs,j/Σk≠s,k=1n ak2
を計算して（ステップＳ２１６）、これを記憶部１０３に記憶される混合比asの更新値として、代入設定する（ステップＳ２１７）処理を繰り返して（ステップＳ２１８）、ステップＳ２１２に戻る。
【００８３】
a1,…,anは、n=3であれば一意な解が得られるが、n≧4の場合には、測定誤差やモデル化誤差、雑音等の影響で、一意な結果が得られないのが一般的である。
【００８４】
そこで、
Σi=1nΣj=1ni≠|aiaj - ri,j|2
を最小化する最適化問題に帰着させることとする。ステップＳ２１４〜ステップＳ２１８の繰り返しは、交互最小自乗法により、局所最適解を求めるものであるが、asの初期値として、上記のような適切なものを選んでいるため、大域解を求めることができる。
【００８５】
なお、各asを正の値とするような事前条件が設定されている場合には、ステップＳ２１４において得られた結果の絶対値をステップＳ２１５において設定する交互凸射影法を適用することも可能である。
【００８６】
したがって、ＣＰＵは、ＲＡＭ等と共働して、更新値設定部１０６として機能する。
【００８７】
記憶部１０３に記憶される混合比a1,…,anが収束すると（ステップＳ２１４；Ｙｅｓ）、混合比出力部１０７は、記憶部１０３に記憶される混合比a1,…,anを出力する（ステップＳ２１９）。
【００８８】
典型的には、ＣＰＵの制御の下、ＲＡＭに記憶された値をＨＤに書き出す。したがって、ＣＰＵは、ＲＡＭ等と共働して、混合比出力部１０７として機能する。
【００８９】
上記のようにして、観測信号のモデルに用いるasが求められたら、必要に応じて（以降の処理は省略可能である。）、各チャネルにおける信号雑音比（signal noise ratio；ＳＮ比）を求めることができる。
【００９０】
すなわち、信号雑音比推定部１０８は、整数i = 1,…,nのそれぞれについて（ステップＳ２２０）、当該チャネルiの信号雑音比ξiを、
ai2/〔Σt=1T xi(t)2/T - ai2〕
により計算して（ステップＳ２２１）、これをＲＡＭ内のξiの記憶域に代入する（ステップＳ２２２）処理を繰り返す（ステップＳ２２３）。また、計算が終わった後にξiをＨＤに書き出すこととしても良い。したがって、ＣＰＵは、ＲＡＭ等と共働して、信号雑音比推定部１０８として機能する。
【００９１】
このほか、各チャネルに含まれる共通信号のパワーを、出力された混合比a1,…,anから求めることも可能である。この場合、図示しない信号強度推定部は、チャネルiに含まれる共通信号のパワーPiを
Pi = ai2
により計算する。したがって、ＣＰＵは、ＲＡＭ等と共働して、信号強度推定部として機能する。
【００９２】
さらに、観測信号のモデルに用いるasが求められたら、必要に応じて（以降の処理は省略可能である。）、原信号g(t)の推定値h(t)を求めることができる。すなわち、観測信号xi(t)に対する重みwiとして適切なものを推定し、h(t)は、当該重みを用いた重み付き平均とするものである。
【００９３】
最適な重みを得るためには、評価関数が必要である。評価関数としては、以下のような２種類を考えることができる。
（１）信号雑音比J1を評価関数とする。具体的には、
J1 = E[(Σi=1n aiwig)2]/E[(Σi=1n wipi)2]
（２）原信号g(t)とのユークリッド距離J2を評価関数とする。具体的には、
J2 = E[(g - h)2]
【００９４】
ここで、手法（１）は一般固有値問題に帰着され、手法（２）は逆行列問題に帰着される。
【００９５】
そこで、ベクトル計算部１０９は、xi(t)およびasから、重みベクトルw = [w1,…,wn]†を計算する（ステップＳ２２４）。
【００９６】
具体的には、手法（１）では、行列
a = [a1,…,an]†；
S = a a†；
Q = diag(Σt=1T x1(t)2/T,…,Σt=1T xn(t)2/T)；
B = diag(a12,…,an2)
に対する一般固有値問題
S w = λ(Q - B)w
の最大固有値解λに対応付けられる固有ベクトル（|w| = 1）が、重みベクトルwとなる。
【００９７】
ここで、・†はベクトルや行列の転置を意味し、diag(・)は、指定された値列を対角要素としそれ以外の要素を０とする対角行列を意味する。
【００９８】
一般固有値問題の解法そのものは、線形代数における種々の公知技術をそのまま利用することができる。
【００９９】
また、手法（２）では、重みベクトルwは、
w = (S + Q - B)-1a
のように計算できる。手法（２）が適用できるのは、全チャネルにおいてE[xi2] > 0が成立する場合であり、これは、センサーから０以外の観測値が得られる場合に成立するから、ほとんどの場合に利用が可能である。
【０１００】
雑音を含まないセンサーがあると仮定した場合、Q - Bは特異となるが、実際には、すべてのセンサーが雑音を含むので、実用的なほとんどの場合Q - Bは正則となる。したがって、手法（１）においては０でない固有値が得られ、手法（２）においては逆行列を求めることができる。
【０１０１】
これらの手法においては、各行列やベクトルの値や一時的な値をＲＡＭに格納し、ＣＰＵにて計算を行うので、ＣＰＵがＲＡＭと共働して、ベクトル計算部１０９として機能することとなる。
【０１０２】
このようにして、重みベクトルwが得られたら、原信号推定部１１０は、原信号g(t)の推定信号h(t)を、重み付き平均により、
h(t) = Σi=1n wixi(t)
により計算して（ステップＳ２２５）、本処理を終了する。計算結果は、ＲＡＭに格納することとしても良いし、ＨＤに書き出すこととしても良い。
【０１０３】
したがって、ＣＰＵは、ＲＡＭ等と共働して、原信号推定部１１０として機能する。
【０１０４】
（実験結果）
以下では、上記実施形態の手法により行ったコンピュータシミュレーションによる実験とその結果について簡単に説明する。
【０１０５】
まず原信号g(t)は、平均０、分散１のガウス分布にしたがうものとし、各センサーにおける混合比aiは、区間[1,α]の一様分布の逆数にしたがうものとし、雑音pi(t)は、平均０、分散１のガウス分布にしたがうものとする。
【０１０６】
すなわち、雑音のエネルギーは各チャネルで固定であるが、各チャネルにおける原信号のエネルギーはランダムになる。また、センサーと信号源との距離を区間[1,α]の一様分布と考え、信号の強さは距離に反比例するものとしている。
【０１０７】
すると、aiの期待値はlog α/(α-1)となり、最大値は１である。
【０１０８】
以下では、
（１）手法（１）により得られた信号をh1、
（２）手法（２）により得られた信号をh2、
（３）ξiが最大となったチャネルの信号をh3、
（４）各チャネルの単純平均信号をh4
として説明する。
【０１０９】
そして、
（１）hiの信号雑音比、
（２）hiの自乗誤差の平均E[(g/|g| - hi/|hi|)2]、
（３）それぞれについて推定された混合比ベクトルの方向ベクトルと、真の混合比ベクトルの方向ベクトルと、の自乗誤差の平均
をシミュレーションによって求めた。実験においては、異なる乱数の種に対して行った１０００回の試行の平均を求めている。
【０１１０】
図３、図４は、各信号の信号雑音比を示すグラフである。以下、本図を参照して説明する。
【０１１１】
縦軸は、両図とも信号雑音比であり、横軸は、図３ではE[ai]であり、図４ではサンプル数Tである。
【０１１２】
両図を見ると、h1とh2とはほとんど同じ姿をしていることがわかる。図１を見ると、平均信号h4が-20dBであっても、h1とh2では-11dB程度であることがわかる。
【０１１３】
したがって、h1やh2のエネルギーは、h4の約８倍ということになる。
【０１１４】
図５、図６は、各信号の原信号との自乗誤差を示すグラフである。以下、本図を参照して説明する。
【０１１５】
縦軸は、両図とも自乗誤差の平均であり、横軸は、図５ではE[ai]であり、図６ではサンプル数Tである。
【０１１６】
前述の図と同様、h1とh2とはほとんど同じ姿をしており、全範囲において最良の推定が行われていることがわかる。
【０１１７】
図７、図８は、推定された混合比ベクトルと真の混合比ベクトルとの自乗誤差を示すグラフである。以下、本図を参照して説明する。
【０１１８】
縦軸は、両図とも自乗誤差の平均であり、横軸は、図７ではE[ai]であり、図８ではサンプル数Tである。
【０１１９】
本図を見ると、本発明による手法では、E[ai] < 0.02であるような混合比ベクトルについては、推定ができないことがある。しかしながら上述のように、推定された信号は、単純平均よりも大きな信号雑音比となっている。
【０１２０】
したがって、本願の手法では、混合比ベクトルの自乗誤差が大きい場合であっても、混合比が大きいチャネルが適切に検出できていることになる。
【０１２１】
また、図７、図８において、Conventionalとなっているものは、[非特許文献１]による手法である。また、LLSEは、混合比を求める際に対数を用いる他の手法である。また、before ALSは、aiの初期値をそのまま利用した場合（収束するまでの反復計算を行わない場合）である。
【０１２２】
本図を見てもわかる通り、[非特許文献１]による手法およびLLSEの性能は不十分であるが、本願の一部をなす技術before ALSは、これらよりも性能が良好と考えられる。
【０１２３】
本願による手法after ALSは、サンプル数が１万、E[ai]が０.０３より大きければ、すなわち、センサーの信号雑音比が−３０ｄＢより大きければ、十分な性能が得られていることがわかる。
【０１２４】
このほか、実際にＥＥＧデータからＳＳＶＥＰ（Steady State Visually Evoked Potential）を抽出する実験をおこなった。目的は、ある対象についての視覚刺激が与えられたときの応答をフラッシングパターンとして抽出し、これに対応付けられる周波数と、信号源とを見つけることである。
【０１２５】
ＳＳＶＥＰの応答は、脳の後頭部に現れることが多いが、個人差がかなり大きく、正確な点を同定するのは難しい。
【０１２６】
前処理として、バタワースハイパスフィルタ（５Ｈｚ）とノッチフィルタ（５０Ｈｚ）を用いてＤＣ成分と電源からの雑音を取り除いた。
【０１２７】
実験では、６４チャネルのＥＥＧシステムを利用し、その中の、後頭部に相当する２９チャネルを調べた。サンプリング周波数は２０４８Ｈｚであり、フラッシング周波数は１８.７Ｈｚである。
【０１２８】
そして、２つの対象について、各対象ごとに試行を５１回行った。
【０１２９】
以下に、得られた信号h1，…，h4、[非特許文献１]、による手法において、対象１と対象２のそれぞれを与えたときの、フラッシング周波数とその半分の周波数に対するパワーの、全信号パワーに対する割合を１００分率（対象１の場合／対象２の場合）で表すと、それぞれ以下のようになる。
（１）h1  …  ７.８７／８.３６
（２）h2  …  ７.８７／８.３６
（３）h3  …  ８.６９／９.７３
（４）h4  …  ６.３６／７.２７
（５）[非特許文献１]  …  ７.１３／８.０１
【０１３０】
これらから、h2とh4を比較したときの片側ｔ−テストのｐ値は、０.６／２９.７であり、h2と[非特許文献１]を比較したときの片側ｔ−テストのｐ値は、５.３／２８.３であることがわかる。これから、本願の手法が有効であることがわかる。
【０１３１】
図９は、推測された信号のパワー分布図である。以下、本図を参照して説明する。
【０１３２】
本図左側は対象１（Subject 1）、右側は対象２（Subject 2）に対するものであり、暗くなっている部分が、推定された信号のパワーの強さを表す。この実験によれば、対象１、２の両方について、右側（Ｒｉｇｈｔ）の後頭部が信号のソースとなっていることがわかる。
【０１３３】
視覚のセンターは脳の後半分に対応することがよく知られているが、本発明の手法による本実験の結果は、この定説によく合致している。
【産業上の利用可能性】
【０１３４】
上記のように、本発明によれば、１つの原信号を複数のチャネルで観測した観測信号から、各チャネルの信号雑音比や原信号を推定するのに好適な信号処理装置、信号処理方法、ならびに、これをコンピュータにて実現するプログラムを提供することができる。
【図面の簡単な説明】
【０１３５】
【図１】本実施形態に係る信号処理装置の概要構成を示す模式図である。
【図２】本実施形態に係る信号処理装置にて実行される信号処理方法の制御の流れを示すフローチャートである。
【図３】各信号の信号雑音比を示すグラフである。
【図４】各信号の信号雑音比を示すグラフである。
【図５】各信号の原信号との自乗誤差を示すグラフである。
【図６】各信号の原信号との自乗誤差を示すグラフである。
【図７】推定された混合比ベクトルと真の混合比ベクトルとの自乗誤差を示すグラフである。
【図８】推定された混合比ベクトルと真の混合比ベクトルとの自乗誤差を示すグラフである。
【図９】推測された信号のパワー分布図である。
【符号の説明】
【０１３６】
１０１  信号処理装置
１０２  入力受付部
１０３  記憶部
１０４  相関計算部
１０５  初期値設定部
１０６  更新値設定部
１０７  混合比出力部
１０８  信号雑音比推定部
１０９  ベクトル計算部
１１０  原信号推定部
【図１】
2008293163.tif 000003
【図２】
2008293163.tif 000004
【図３】
2008293163.tif 000005
【図４】
2008293163.tif 000006
【図５】
2008293163.tif 000007
【図６】
2008293163.tif 000008
【図７】
2008293163.tif 000009
【図８】
2008293163.tif 000010
【図９】
2008293163.tif 000011