JP5985079B2 - スペクトル推定装置及びスペクトル推定方法 - Google Patents
スペクトル推定装置及びスペクトル推定方法 Download PDFInfo
- Publication number
- JP5985079B2 JP5985079B2 JP2015553383A JP2015553383A JP5985079B2 JP 5985079 B2 JP5985079 B2 JP 5985079B2 JP 2015553383 A JP2015553383 A JP 2015553383A JP 2015553383 A JP2015553383 A JP 2015553383A JP 5985079 B2 JP5985079 B2 JP 5985079B2
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- spectrum
- vector
- component
- setting
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
- G01S3/74—Multi-channel systems specially adapted for direction-finding, i.e. having a single antenna system capable of giving simultaneous indications of the directions of different signals
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/02—Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
- G01S13/06—Systems determining position data of a target
- G01S13/42—Simultaneous measurement of distance and other co-ordinates
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Variable-Direction Aerials And Aerial Arrays (AREA)
- Radar Systems Or Details Thereof (AREA)
Description
この発明は、レーダ装置に搭載されているアンテナ素子の受信信号から、例えば、波源の位置や方向を推定するスペクトル推定装置及びスペクトル推定方法に関するものである。
波源の位置や方向を推定する到来波推定装置（空間方向に対する電力スペクトルを推定するスペクトル推定装置）では、波源からの電波を受信する複数のアンテナ素子に対して、推定可能な方位角の範囲と刻み幅を設定する。
例えば、角度点数がＮ、アンテナの素子数がＬ、時間サンプル数がＭであるとき、アンテナ素子毎に時間サンプル数分の信号を受信して構成されるＬ×Ｍの観測行列ｙは、雑音成分ｎが含まれていることを踏まえ、アンテナ素子毎に角度点数分の走査ベクトルを並べて構成されるＬ×Ｎの辞書行列Ａによって、下記の式（１）のように表される。
式（１）において、Ｘは角度毎に時間サンプル数分の要素を有するＮ×Ｍの到来波行列であり、空間方向に対する電力スペクトルに相当する情報である。したがって、到来波推定処理は、既知の行列である観測行列ｙ及び辞書行列Ａから、未知行列である到来波行列Ｘの推定を行う処理となる。
例えば、角度点数がＮ、アンテナの素子数がＬ、時間サンプル数がＭであるとき、アンテナ素子毎に時間サンプル数分の信号を受信して構成されるＬ×Ｍの観測行列ｙは、雑音成分ｎが含まれていることを踏まえ、アンテナ素子毎に角度点数分の走査ベクトルを並べて構成されるＬ×Ｎの辞書行列Ａによって、下記の式（１）のように表される。
非特許文献１，２には、最適化アルゴリズムを用いて、未知行列である到来波行列Ｘを求める方法が開示されている。
まず、到来波行列Ｘを評価する目的関数ｆ（Ｘ）を下記の式（２）のように設定する。
また、λとｐは解の精度に寄与する計算用パラメータであり、０＜ｐ≦１の範囲となるように、それぞれ任意の値が設定される。
この目的関数ｆ（Ｘ）から、目的関数ｆ（Ｘ）の勾配方向へと探索を進める最適化アルゴリズムを利用するため、式（２）を近似することで、微分可能な下記の式（３）を目的関数ｆ（Ｘ）として使用する。
ここで、εはゼロ除算を回避するための計算用パラメータであり、任意の値が設定される。
まず、到来波行列Ｘを評価する目的関数ｆ（Ｘ）を下記の式（２）のように設定する。
この目的関数ｆ（Ｘ）から、目的関数ｆ（Ｘ）の勾配方向へと探索を進める最適化アルゴリズムを利用するため、式（２）を近似することで、微分可能な下記の式（３）を目的関数ｆ（Ｘ）として使用する。
上記の式（３）を微分することで、下記の式（４）に示す勾配▽ｆ（Ｘ）を得る。
ここで、ｄｉａｇ｛Ｓ｝は長さＮのベクトルＳを対角成分とするＮ×Ｎの対角行列であり、ＡＨはＮ×Ｎの行列Ａの複素共役転置行列である。
数学的な最適化アルゴリズムを導入した場合、未知の到来波行列Ｘを繰り返し修整していく処理を行う。
ｋ回目の反復時のパラメータを（ｋ）と表記すると、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新しい解）へと修整する式は、下記の式（１０）のように表される。
ここで、αは探索方向に対するステップ幅であり、下記の式（１１）から得られる値である。
ｋ回目の反復時のパラメータを（ｋ）と表記すると、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新しい解）へと修整する式は、下記の式（１０）のように表される。
また、ｄは探索方向を表すＮ×Ｍの行列であり、最適化アルゴリズムの種類に応じて計算法が異なる。最急降下法では、下記の式（１２）に示すように、目的関数ｆ（Ｘ）に対する勾配▽ｆ（Ｘ）を、そのまま探索方向ｄとして扱う。
非線形共役勾配法では、下記の式（１３）に示すように、目的関数ｆ（Ｘ）に対する勾配▽ｆ（Ｘ）に加えて、共役方向βを用いる。
共役方向βの計算法は多様に存在するが、例えば、下記の式（１４）で示すように、Fletcher-Reevesの公式を用いた計算を行う。
準ニュートン法では、目的関数ｆ（Ｘ）に対する勾配▽ｆ（Ｘ）と、目的関数ｆ（Ｘ）の二階微分であるヘッセ行列▽２ｆ（Ｘ）を近似した行列Ｈの逆行列との積を計算する。
式（１５）における行列Ｈは、式（５）のＨに相当する。αを１に設定し、式（１０）に対して、式（１５）と式（４）を代入すると、ヘッセ近似行列Ｈを用いた準ニュートン法は、以下のように展開され、連立方程式（１６）の反復計算となる。
この方法は、以上で示したような最適化アルゴリズムに基づき、以下の手順（１）〜（６）で未知行列である到来波行列Ｘを徐々に修整（反復計算）して、最終的な到来波行列Ｘの推定結果として、波源の方位に対応する要素以外をゼロ成分とした疎行列を得るものである。
（１）式（６）及び式（７）より、行列Ｈ０及び行列ｂを計算する。また、未知行列である到来波行列Ｘに任意の初期値を設定する。
（２）式（９）より、Ｎ×Ｍの到来波行列Ｘを時間サンプル数Ｍの方向に積み上げ、長さＮの到来波ベクトルｘを生成する。
（３）式（８）に対して、生成した到来波ベクトルｘを代入して、対角成分ベクトルＳを計算する。また、式（５）に対して、計算した対角成分ベクトルＳを代入して、行列Ｈ０の対角成分をベクトルＳに設定した行列Ｈを得る。
（１）式（６）及び式（７）より、行列Ｈ０及び行列ｂを計算する。また、未知行列である到来波行列Ｘに任意の初期値を設定する。
（２）式（９）より、Ｎ×Ｍの到来波行列Ｘを時間サンプル数Ｍの方向に積み上げ、長さＮの到来波ベクトルｘを生成する。
（３）式（８）に対して、生成した到来波ベクトルｘを代入して、対角成分ベクトルＳを計算する。また、式（５）に対して、計算した対角成分ベクトルＳを代入して、行列Ｈ０の対角成分をベクトルＳに設定した行列Ｈを得る。
（４）設定した行列Ｈを最適化アルゴリズムに入力する。
連立方程式の計算を伴う準ニュートン法を用いる場合は、式（１６）で与えられた連立方程式Ｈ（ｋ）Ｘ（ｋ＋１）＝ｂを解くことで、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新しい解）に更新する。
一方、連立方程式の計算を必要としない最急降下法や非線形共役勾配法を用いる場合は、まず、式（４）に基づき、行列積Ｈ（ｋ）Ｘ（ｋ）の計算を経て勾配▽ｆ（Ｘ（ｋ））を得る。
次に、勾配▽ｆ（Ｘ（ｋ））を式（１２）又は式（１３）に入力して探索方向行列ｄ（ｋ）を生成する。
最後に、生成した探索方向行列ｄ（ｋ）を式（１０）に入力することで、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新たな解）に更新する。
連立方程式の計算を伴う準ニュートン法を用いる場合は、式（１６）で与えられた連立方程式Ｈ（ｋ）Ｘ（ｋ＋１）＝ｂを解くことで、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新しい解）に更新する。
一方、連立方程式の計算を必要としない最急降下法や非線形共役勾配法を用いる場合は、まず、式（４）に基づき、行列積Ｈ（ｋ）Ｘ（ｋ）の計算を経て勾配▽ｆ（Ｘ（ｋ））を得る。
次に、勾配▽ｆ（Ｘ（ｋ））を式（１２）又は式（１３）に入力して探索方向行列ｄ（ｋ）を生成する。
最後に、生成した探索方向行列ｄ（ｋ）を式（１０）に入力することで、到来波行列Ｘ（ｋ）をＸ（ｋ＋１）（新たな解）に更新する。
（５）（２）と同様にして、到来波行列Ｘ（ｋ＋１）から到来波ベクトルｘ（ｋ＋１）を生成する。
（６）更新後の到来波行列Ｘ（ｋ＋１）が最適解に収束したか否かを、例えば、下記の判定式（１７）を用いて判定し、未だ収束していなければ、（３）の処理に戻る。
下記の式（１７）では、今回の反復計算で生成された到来波ベクトルｘ（ｋ＋１）と、前回の反復計算で生成された到来波ベクトルｘ（ｋ）との差分ベクトルのノルムを、前回の反復計算で生成された到来波ベクトルｘ（ｋ）のノルムで除算し、その除算結果が、予め設定された閾値δｔｈｒｅｓｈｏｌｄより小さければ、更新後の到来波行列Ｘ（ｋ＋１）が最適解に収束したと判定する。
（６）更新後の到来波行列Ｘ（ｋ＋１）が最適解に収束したか否かを、例えば、下記の判定式（１７）を用いて判定し、未だ収束していなければ、（３）の処理に戻る。
下記の式（１７）では、今回の反復計算で生成された到来波ベクトルｘ（ｋ＋１）と、前回の反復計算で生成された到来波ベクトルｘ（ｋ）との差分ベクトルのノルムを、前回の反復計算で生成された到来波ベクトルｘ（ｋ）のノルムで除算し、その除算結果が、予め設定された閾値δｔｈｒｅｓｈｏｌｄより小さければ、更新後の到来波行列Ｘ（ｋ＋１）が最適解に収束したと判定する。
Malioutov，D.，Cetin，M.，and Willsky，A.S.: A Sparse Signal Reconstruction Perspective for Source Localization With Sensor Arrays，IEEE Trans.Signal Processing，Vol.53，No.8，pp.3010-3022，2005.
Cetin，M.，Malioutov，D.，and Willsky，A. S.: A variational technique for source localization based on a sparse signal reconstruction perspective，in Proc. 2002 IEEE International Conference on Acoustics，Speech，and Signal Processing（ICASSP），Vol.3，pp.2965-2968，2002.
従来のスペクトル推定装置は以上のように構成されているので、反復計算（３）〜（６）の中に、演算負荷が高い到来波行列Ｘの修整操作（（４）における連立方程式ＨＸ＝ｂの計算、または、行列積ＨＸの計算）が含まれている。このため、角度点数が大きく、行列Ｈ及び到来波行列Ｘの行列サイズが大きくなると、演算量が爆発的に増加し、到来波行列Ｘの推定結果が得られるまでに多くの時間を要してしまう課題があった。
この発明は上記のような課題を解決するためになされたもので、角度点数が大きく、到来波行列Ｘの行列サイズが大きくても、短時間で到来波行列Ｘを推定することができるスペクトル推定装置及びスペクトル推定方法を得ることを目的とする。
この発明に係るスペクトル推定装置は、１以上のアンテナ素子の走査ベクトルからなる辞書行列の複素共役転置と辞書行列との積の行列Ｈと、１以上のアンテナ素子の受信信号からなる観測行列と辞書行列の複素共役転置との積の行列ｂと、受信信号に含まれているスペクトルに関する未知行列であるスペクトル行列とを設定する行列設定手段と、そのスペクトル行列を積み上げることでスペクトルベクトルを生成するベクトル生成手段と、ベクトル生成手段により生成されたスペクトルベクトルに含まれている各要素のうち、疎成分の要素が集まるように各要素を並び替える並び替え手段と、並び替え手段による並び替え後のスペクトルベクトルから定数成分が集まった設定用ベクトルを算出し、その設定用ベクトルによって行列Ｈの対角成分を設定する対角成分設定手段と、対角成分設定手段により設定用ベクトルが設定された行列Ｈと行列ｂを用いて、そのスペクトル行列を更新する行列更新手段とを設け、収束判定手段が、行列更新手段により更新されたスペクトル行列が収束するまで、スペクトルベクトルの生成からスペクトル行列の更新までの処理を繰り返し実施させるようにしたものである。
この発明によれば、行列設定手段により設定されたスペクトル行列を積み上げることでスペクトルベクトルを生成するベクトル生成手段と、ベクトル生成手段により生成されたスペクトルベクトルに含まれている各要素のうち、疎成分の要素が集まるように各要素を並び替える並び替え手段と、並び替え手段による並び替え後のスペクトルベクトルから定数成分が集まった設定用ベクトルを算出し、その設定用ベクトルによって行列Ｈの対角成分を設定する対角成分設定手段と、対角成分設定手段により設定用ベクトルが設定された行列Ｈと行列ｂを用いて、そのスペクトル行列を更新する行列更新手段とを設け、スペクトル行列の疎成分に対する演算結果を再利用するように構成したものであるため、角度点数が大きく、疎成分を多量に含むスペクトル行列の行列サイズが大きくても、短時間でスペクトル行列を推定することができる効果がある。
実施の形態１．
この実施の形態１では、スペクトル推定装置が、空間方向に対する電力スペクトルである到来波（アンテナ素子の受信信号に含まれている到来波）を推定する例を説明する。
図１はこの発明の実施の形態１によるスペクトル推定装置を示す構成図である。
図１において、行列設定部１はＬ本のアンテナ素子の走査ベクトルが角度点数Ｎだけ並べられたＬ×Ｎの辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０と、辞書行列Ａの複素共役転置ＡＨとＬ本のアンテナ素子の受信信号が時間サンプル数Ｍだけ並べられたＬ×Ｍの観測行列ｙとの積の行列である行列ｂと、Ｌ本のアンテナ素子の受信信号に含まれている到来波に関する未知行列であるＮ×Ｍの到来波行列（スペクトル行列）Ｘとを設定する処理を実施する。Ｎは電波を受信する方位角の点数である。なお、行列設定部１は行列設定手段を構成している。
この実施の形態１では、スペクトル推定装置が、空間方向に対する電力スペクトルである到来波（アンテナ素子の受信信号に含まれている到来波）を推定する例を説明する。
図１はこの発明の実施の形態１によるスペクトル推定装置を示す構成図である。
図１において、行列設定部１はＬ本のアンテナ素子の走査ベクトルが角度点数Ｎだけ並べられたＬ×Ｎの辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０と、辞書行列Ａの複素共役転置ＡＨとＬ本のアンテナ素子の受信信号が時間サンプル数Ｍだけ並べられたＬ×Ｍの観測行列ｙとの積の行列である行列ｂと、Ｌ本のアンテナ素子の受信信号に含まれている到来波に関する未知行列であるＮ×Ｍの到来波行列（スペクトル行列）Ｘとを設定する処理を実施する。Ｎは電波を受信する方位角の点数である。なお、行列設定部１は行列設定手段を構成している。
到来波ベクトル生成部２は行列設定部１により設定された到来波行列Ｘを時間サンプル数Ｍだけ積み上げることで、スペクトルベクトルである到来波ベクトルｘを生成する処理を実施する。なお、到来波ベクトル生成部２はベクトル生成手段を構成している。
要素並び替え部３は到来波ベクトル生成部２により生成された到来波ベクトルｘに含まれている各要素ｘｋのうち、疎成分の要素を抽出するように各要素ｘｋを並び替える処理を実施する。要素並び替え部３は並び替え手段を構成している。
要素並び替え部３は到来波ベクトル生成部２により生成された到来波ベクトルｘに含まれている各要素ｘｋのうち、疎成分の要素を抽出するように各要素ｘｋを並び替える処理を実施する。要素並び替え部３は並び替え手段を構成している。
対角成分設定部４は要素並び替え部３により各要素ｘｋが並び替えられた到来波ベクトルｘから定数成分が集まった設定用ベクトルＳを算出するとともに、その設定用ベクトルＳを行列設定部１により設定された初期行列Ｈ０の対角成分に設定し、設定後の初期行列Ｈ０を行列Ｈとして出力する処理を実施する。なお、対角成分設定部４は対角成分設定手段を構成している。
到来波行列更新部５は対角成分設定部４により設定用ベクトルＳが対角成分に設定された行列Ｈと行列設定部１により設定された行列ｂを用いて、到来波行列Ｘを更新する処理を実施する。なお、到来波行列更新部５は行列更新手段を構成している。
到来波行列更新部５は対角成分設定部４により設定用ベクトルＳが対角成分に設定された行列Ｈと行列設定部１により設定された行列ｂを用いて、到来波行列Ｘを更新する処理を実施する。なお、到来波行列更新部５は行列更新手段を構成している。
収束判定部６は到来波行列更新部５により更新された到来波行列Ｘが収束しているか否かを判定し、未だ収束していなければ、到来波ベクトル生成部２、要素並び替え部３、対角成分設定部４及び到来波行列更新部５に対して反復計算指令を出力することで、到来波ベクトルｘの生成から到来波行列Ｘの更新までの処理を繰り返し実施させる処理を実施する。
到来波ベクトル復元部７は収束判定部６により到来波行列Ｘが収束していると判定された場合、要素並び替え部３により並び替えられている到来波ベクトルｘの各要素ｘｋを元の位置に戻してから、当該到来波ベクトルｘを出力する処理を実施する。
なお、収束判定部６及び到来波ベクトル復元部７から収束判定手段が構成されている。
到来波ベクトル復元部７は収束判定部６により到来波行列Ｘが収束していると判定された場合、要素並び替え部３により並び替えられている到来波ベクトルｘの各要素ｘｋを元の位置に戻してから、当該到来波ベクトルｘを出力する処理を実施する。
なお、収束判定部６及び到来波ベクトル復元部７から収束判定手段が構成されている。
図１の例では、スペクトル推定装置の構成要素である行列設定部１、到来波ベクトル生成部２、要素並び替え部３、対角成分設定部４、到来波行列更新部５、収束判定部６及び到来波ベクトル復元部７のそれぞれが専用のハードウェア（例えば、ＣＰＵを実装している半導体集積回路、あるいは、ワンチップマイコンなど）で構成されているものを想定しているが、スペクトル推定装置がコンピュータで構成されていてもよい。
スペクトル推定装置をコンピュータで構成する場合、行列設定部１、到来波ベクトル生成部２、要素並び替え部３、対角成分設定部４、到来波行列更新部５、収束判定部６及び到来波ベクトル復元部７の処理内容を記述しているプログラムをコンピュータのメモリに格納し、当該コンピュータのＣＰＵが当該メモリに格納されているプログラムを実行するようにすればよい。
スペクトル推定装置をコンピュータで構成する場合、行列設定部１、到来波ベクトル生成部２、要素並び替え部３、対角成分設定部４、到来波行列更新部５、収束判定部６及び到来波ベクトル復元部７の処理内容を記述しているプログラムをコンピュータのメモリに格納し、当該コンピュータのＣＰＵが当該メモリに格納されているプログラムを実行するようにすればよい。
次に動作について説明する。
この実施の形態１では、到来波行列更新部５の手段として、連立方程式の計算を必要とする最適化アルゴリズムを利用する。この最適化アルゴリズムを説明する便宜上、「準ニュートン法」の操作手順について説明するが、準ニュートン法と同様に連立方程式の計算を必要とする他の最適化アルゴリズム（例えば「ニュートン法」）を使用するようにしてもよい。また、説明の便宜上、連立方程式の求解手法として、行列Ｈを二種類の三角行列に分解する「ＬＵ分解」を用いた方法を例として挙げるが、ＬＵ分解同様に、行列Ｈの対角成分を左上隅から右下隅へと順に選択し、未選択の対角成分を含む部分行列を更新していく他の方法（例えば、「コレスキー分解」や「修正コレスキー分解」を用いた方法）を使用するようにしてもよい。
この実施の形態１では、推定対象の波源から到来している信号データの観測環境として、電波を受信する方位角の点数がＮ、アンテナの素子数がＬ、時間サンプル数がＭであり、Ｌ本のアンテナ素子の受信信号が時間サンプル数Ｍだけ並べられた行列がＬ×Ｍの観測行列ｙ、Ｌ本のアンテナ素子による角度点数Ｎに対する走査ベクトルが並べられた行列がＬ×Ｎの辞書行列Ａであるとして説明する。
この実施の形態１では、到来波行列更新部５の手段として、連立方程式の計算を必要とする最適化アルゴリズムを利用する。この最適化アルゴリズムを説明する便宜上、「準ニュートン法」の操作手順について説明するが、準ニュートン法と同様に連立方程式の計算を必要とする他の最適化アルゴリズム（例えば「ニュートン法」）を使用するようにしてもよい。また、説明の便宜上、連立方程式の求解手法として、行列Ｈを二種類の三角行列に分解する「ＬＵ分解」を用いた方法を例として挙げるが、ＬＵ分解同様に、行列Ｈの対角成分を左上隅から右下隅へと順に選択し、未選択の対角成分を含む部分行列を更新していく他の方法（例えば、「コレスキー分解」や「修正コレスキー分解」を用いた方法）を使用するようにしてもよい。
この実施の形態１では、推定対象の波源から到来している信号データの観測環境として、電波を受信する方位角の点数がＮ、アンテナの素子数がＬ、時間サンプル数がＭであり、Ｌ本のアンテナ素子の受信信号が時間サンプル数Ｍだけ並べられた行列がＬ×Ｍの観測行列ｙ、Ｌ本のアンテナ素子による角度点数Ｎに対する走査ベクトルが並べられた行列がＬ×Ｎの辞書行列Ａであるとして説明する。
この実施の形態１は、ベクトルの疎成分に対する演算結果を引き継ぐ形で、連立方程式の計算に基づき未知行列の修整を行う方法に関するものであるが、この実施の形態１でのスペクトル推定方法の利点を明らかにするため、最初に、非特許文献１，２に開示されているスペクトル推定方法を説明し、このスペクトル推定方法との違いを示す形で、この実施の形態１でのスペクトル推定方法を説明する。
図２はこの発明の実施の形態１によるスペクトル推定装置の処理内容（スペクトル推定方法）を示すフローチャートである。
また、図３は非特許文献１，２に開示されている従来のスペクトル推定方法を示すフローチャートである。
図２はこの発明の実施の形態１によるスペクトル推定装置の処理内容（スペクトル推定方法）を示すフローチャートである。
また、図３は非特許文献１，２に開示されている従来のスペクトル推定方法を示すフローチャートである。
［従来のスペクトル推定方法］
図４は従来のスペクトル推定方法での処理フローを図解している説明図である。
まず、図１の行列設定部１に相当する行列設定部（図示せず）は、辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０（図４中、符号９８のＨ０）と、観測行列ｙと辞書行列Ａの複素共役転置ＡＨとの積の行列である行列ｂ（図４中、符号１００のｂ）と、到来波に関する未知行列であるＮ×Ｍの到来波行列Ｘ（図４中、符号１０１のＸの初期値）とを設定する（図３のステップＳＴ１）。
なお、初期行列Ｈ０は、連立方程式ＨＸ＝ｂにおける左辺の行列Ｈの初期値であり、未知行列である到来波行列Ｘの初期値は、任意の値で構わない。例えば、任意のスペクトル推定方法による推定結果や、全ての成分を任意の乱数生成方法で生成した乱数値でも構わない。
図４は従来のスペクトル推定方法での処理フローを図解している説明図である。
まず、図１の行列設定部１に相当する行列設定部（図示せず）は、辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０（図４中、符号９８のＨ０）と、観測行列ｙと辞書行列Ａの複素共役転置ＡＨとの積の行列である行列ｂ（図４中、符号１００のｂ）と、到来波に関する未知行列であるＮ×Ｍの到来波行列Ｘ（図４中、符号１０１のＸの初期値）とを設定する（図３のステップＳＴ１）。
なお、初期行列Ｈ０は、連立方程式ＨＸ＝ｂにおける左辺の行列Ｈの初期値であり、未知行列である到来波行列Ｘの初期値は、任意の値で構わない。例えば、任意のスペクトル推定方法による推定結果や、全ての成分を任意の乱数生成方法で生成した乱数値でも構わない。
図１の到来波ベクトル生成部２に相当する到来波ベクトル生成部（図示せず）は、行列設定部が到来波行列Ｘを設定すると、上記の式（９）に示すように、その到来波行列Ｘ（図４中、符号１０１のＸ）を時間サンプル数Ｍだけ積み上げることで、到来波ベクトルｘ（図４中、符号１０２のｘ）を生成する（図３のステップＳＴ４）。
図１の対角成分設定部４に相当する対角成分設定部（図示せず）は、到来波ベクトル生成部が到来波ベクトルｘを生成すると、上記の式（８）に示すように、その到来波ベクトルｘから設定用ベクトルＳ（図４中、符号１０３のＳ）を算出する（図３のステップＳＴ７）。
対角成分設定部は、設定用ベクトルＳを算出すると、式（５）に示すように、その設定用ベクトルＳを行列設定部により設定された初期行列Ｈ０（図４中、符号９８のＨ０）の対角成分に設定し、設定後の初期行列Ｈ０を行列Ｈ（図４中、符号９９のＨ）として出力する（図３のステップＳＴ８）。
即ち、対角成分設定部は、全ての成分が０である行列の対角成分に設定用ベクトルＳを設定することで、図４中、符号１０４の疎行列を定義し、その疎行列を行列設定部により設定された初期行列Ｈ０に足し合わせることで、図４中、符号９９の行列Ｈを生成する。
即ち、対角成分設定部は、全ての成分が０である行列の対角成分に設定用ベクトルＳを設定することで、図４中、符号１０４の疎行列を定義し、その疎行列を行列設定部により設定された初期行列Ｈ０に足し合わせることで、図４中、符号９９の行列Ｈを生成する。
到来波行列更新部５に相当する到来波行列更新部（図示せず）は、対角成分設定部が行列Ｈの対角成分を設定すると、式（１６）に示すように、連立方程式ＨＸ＝ｂを解くことで到来波行列Ｘ（図４中、符号１０１のＸ）を求め、その到来波行列Ｘを更新後の到来波行列Ｘ（新しい解）に設定する（図３のステップＳＴ１０）。
以下、到来波行列更新部の処理内容を具体的に説明する。
以下、到来波行列更新部の処理内容を具体的に説明する。
到来波行列更新部は、対角成分設定部が行列Ｈの対角成分を設定すると、下記の式（１８）に示すように、その行列Ｈを下三角行列Ｌ（上三角成分がゼロの行列）と上三角行列Ｕ（下三角成分がゼロの行列）との積（Ｈ＝ＬＵ）の形に分解する（ＬＵ分解操作）。
ＬＵ分解では、行列Ｈの左上隅の成分Ｈ１，１を除数として選択し、 同じ列の成分（Ｈ２，１，Ｈ３，１…ＨＮ，１）を除算して行列Ｈを更新する。ここで、下記に示すように除算後のＨ２，１，Ｈ３，１…ＨＮ，１を（Ｎ―１）行１列の行列、Ｈ１，１と同じ行の成分Ｈ１，２，Ｈ１，３…Ｈ１，Ｎを１行（Ｎ−１）列の行列とする。
この時点で、対角成分Ｈ１，１と上記の行列成分は、以降の操作で更新されない成分として確定する。このとき、上記の（Ｎ−１）行１列の行列と１行（Ｎ−１）列の行列との積をとった下記の（Ｎ−１）行（Ｎ−１）列の行列を、行列Ｈの残る成分値（部分行列）に対するマイナス値として、行列Ｈを更新する。
これと同様の操作を、行列Ｈの左上隅から右下隅へ順に、対角成分を選択して繰り返すことで、下三角行列Ｌと上三角行列Ｕをひとまとめに格納した行列Ｈが生成される。
到来波行列更新部は、行列Ｈを下三角行列Ｌと上三角行列Ｕとの積をＨ＝ＬＵの形に分解した後の連立方程式ＬＵＸ＝ｂを、ｔ＝ＵＸなる行列ｔを定義することで、下記の式（１９）に示す連立方程式へと変形する。
ここで、下三角行列Ｌは上三角成分がゼロであるため、未知行列ｔの成分値は、下記の式（２０）で表されるように、代入操作（前進代入）によって徐々に求まっていくものである。
この結果、得られた行列ｔを使用して、同様の代入操作をＵＸ＝ｔに対して行えば、連立方程式の解である到来波行列Ｘが求まる。
ただし、Ｕは上三角行列であるため、代入操作は逆順（後進代入）となる。
ただし、Ｕは上三角行列であるため、代入操作は逆順（後進代入）となる。
図１の収束判定部６に相当する収束判定部（図示せず）は、到来波行列更新部が到来波行列Ｘを更新する（到来波行列Ｘに対する修整操作を行う）と、式（１７）に示すように、更新後の到来波行列Ｘが収束しているか否かを判定する（図３のステップＳＴ１２）。
以上のように構成されている従来のスペクトル推定方法では、ステップＳＴ４で得られる到来波ベクトルｘが反復の度に疎ベクトルへと修整されていく特徴を有している。
図５は到来波ベクトルｘが疎ベクトルに修整されていく様子を示す説明図である。
反復１周目の到来波ベクトルｘ（図５中、符号１２０のｘ）は、非疎成分（斜線部）で埋められており、この到来波ベクトルｘから設定用ベクトルＳ（図５中、符号１２１のＳ）が算出され、この設定用ベクトルＳが対角成分（図５中、符号１２３の対角成分）となる行列Ｈ（図５中、符号１２２のＨ）が生成される。この行列Ｈから、ステップＳＴ１０での連立方程式の求解操作を経て、反復２周目の処理に移行する。
図５は到来波ベクトルｘが疎ベクトルに修整されていく様子を示す説明図である。
反復１周目の到来波ベクトルｘ（図５中、符号１２０のｘ）は、非疎成分（斜線部）で埋められており、この到来波ベクトルｘから設定用ベクトルＳ（図５中、符号１２１のＳ）が算出され、この設定用ベクトルＳが対角成分（図５中、符号１２３の対角成分）となる行列Ｈ（図５中、符号１２２のＨ）が生成される。この行列Ｈから、ステップＳＴ１０での連立方程式の求解操作を経て、反復２周目の処理に移行する。
反復２周目以降の処理では、ステップＳＴ４で到来波ベクトルｘ（図５中、符号１２４のｘ）を生成するが、この到来波ベクトルｘは、疎成分を含むベクトルに修整されている。
この到来波ベクトルｘから算出される設定用ベクトルＳ（図５中、符号１２５のＳ）には、疎成分に対する演算結果である定数が格納される。疎成分へと修整されたベクトル成分は、以降の反復計算中も常に疎であるものとして扱えば、行列Ｈ（図５中、符号１２６のＨ）の対角成分（図５中、符号１２７の対角成分）には、以降の反復計算でも常に定数となる成分が出現する特徴がある。
この反復処理が終了するとき、到来波ベクトルｘ（図５中、符号１２８のｘ）は、ほとんど全ての成分が疎で占められた疎ベクトルになる。
この到来波ベクトルｘから算出される設定用ベクトルＳ（図５中、符号１２５のＳ）には、疎成分に対する演算結果である定数が格納される。疎成分へと修整されたベクトル成分は、以降の反復計算中も常に疎であるものとして扱えば、行列Ｈ（図５中、符号１２６のＨ）の対角成分（図５中、符号１２７の対角成分）には、以降の反復計算でも常に定数となる成分が出現する特徴がある。
この反復処理が終了するとき、到来波ベクトルｘ（図５中、符号１２８のｘ）は、ほとんど全ての成分が疎で占められた疎ベクトルになる。
［実施の形態１のスペクトル推定方法］
図６は実施の形態１のスペクトル推定方法での処理フローを図解している説明図である。
まず、行列設定部１は、辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０（図６中、符号１４０のＨ０）と、観測行列ｙと辞書行列Ａの複素共役転置ＡＨとの積の行列である行列ｂ（図６中、符号１４３のｂ）と、到来波に関する未知行列であるＮ×Ｍの到来波行列Ｘ（図６中、符号１４５のＸの初期値）とを設定する（図２のステップＳＴ１）。
なお、初期行列Ｈ０は、連立方程式ＨＸ＝ｂにおける左辺の行列Ｈの初期値であり、未知行列である到来波行列Ｘの初期値は、任意の値で構わない。例えば、任意のスペクトル推定方法による推定結果や、全ての成分を任意の乱数生成方法で生成した乱数値でも構わない。
図６は実施の形態１のスペクトル推定方法での処理フローを図解している説明図である。
まず、行列設定部１は、辞書行列Ａの複素共役転置ＡＨと辞書行列Ａとの積の行列である行列Ｈの初期行列Ｈ０（図６中、符号１４０のＨ０）と、観測行列ｙと辞書行列Ａの複素共役転置ＡＨとの積の行列である行列ｂ（図６中、符号１４３のｂ）と、到来波に関する未知行列であるＮ×Ｍの到来波行列Ｘ（図６中、符号１４５のＸの初期値）とを設定する（図２のステップＳＴ１）。
なお、初期行列Ｈ０は、連立方程式ＨＸ＝ｂにおける左辺の行列Ｈの初期値であり、未知行列である到来波行列Ｘの初期値は、任意の値で構わない。例えば、任意のスペクトル推定方法による推定結果や、全ての成分を任意の乱数生成方法で生成した乱数値でも構わない。
行列設定部１は、各種の行列を設定すると、疎成分に対する部分演算結果の格納領域Ｐ（図６中、符号１４２のＰ）を初期化して、データ領域を確保する（図２のステップＳＴ２）。
また、行列設定部１は、到来波ベクトルｘの並べ替え情報（後述する要素並び替え部３により各要素が並び替えられた到来波ベクトルｘの各要素の並び順を示す情報）の格納領域Ｑ（図６中、符号１４４のＱ）を初期化して、データ領域を確保する（図２のステップＳＴ３）。
また、行列設定部１は、到来波ベクトルｘの並べ替え情報（後述する要素並び替え部３により各要素が並び替えられた到来波ベクトルｘの各要素の並び順を示す情報）の格納領域Ｑ（図６中、符号１４４のＱ）を初期化して、データ領域を確保する（図２のステップＳＴ３）。
到来波ベクトル生成部２は、行列設定部１が到来波行列Ｘを設定すると、上記の式（９）に示すように、その到来波行列Ｘ（図６中、符号１４５のＸ）を時間サンプル数Ｍだけ積み上げることで、到来波ベクトルｘ（図６中、符号１４６のｘ）を生成する（図２のステップＳＴ４）。
要素並び替え部３は、到来波ベクトル生成部２が到来波ベクトルｘを生成すると、その到来波ベクトルｘに含まれている要素ｘｋ毎に、当該要素ｘｋが疎であるか否かを判定し、疎成分の要素が左端に集まるように各要素ｘｋを並び替える処理を実施する。
ここで、図７は到来波ベクトルｘが疎ベクトルに修整されていく様子を示す説明図である。詳細は後述するが、反復２周目において、到来波ベクトルｘに含まれている疎成分の要素が左端に集められている。
要素並び替え部３は、各要素ｘｋを並び替える処理を実施すると、並び替え後の各要素ｘｋの並び順を示す並べ替え情報を格納領域Ｑ（図６中、符号１４４のＱ）に格納する（図２のステップＳＴ５）。
また、要素並び替え部３は、到来波ベクトルｘに含まれている各要素ｘｋの並び替えを行うと、並べ替えた要素に対する演算結果が変わらないようにするため、行列設定部１により生成された行列ｂ（図６中、符号１４３のｂ）の各要素の並び順が、到来波ベクトルｘの各要素ｘｋの並び順と一致するように並び替える（図２のステップＳＴ６）。
ここで、図７は到来波ベクトルｘが疎ベクトルに修整されていく様子を示す説明図である。詳細は後述するが、反復２周目において、到来波ベクトルｘに含まれている疎成分の要素が左端に集められている。
要素並び替え部３は、各要素ｘｋを並び替える処理を実施すると、並び替え後の各要素ｘｋの並び順を示す並べ替え情報を格納領域Ｑ（図６中、符号１４４のＱ）に格納する（図２のステップＳＴ５）。
また、要素並び替え部３は、到来波ベクトルｘに含まれている各要素ｘｋの並び替えを行うと、並べ替えた要素に対する演算結果が変わらないようにするため、行列設定部１により生成された行列ｂ（図６中、符号１４３のｂ）の各要素の並び順が、到来波ベクトルｘの各要素ｘｋの並び順と一致するように並び替える（図２のステップＳＴ６）。
対角成分設定部４は、要素並び替え部３が到来波ベクトルｘの各要素ｘｋを並び替えると、上記の式（８）に示すように、その到来波ベクトルｘから設定用ベクトルＳ（図６中、符号１４７のＳ）を算出する（図２のステップＳＴ７）。
対角成分設定部４は、設定用ベクトルＳを算出すると、式（５）に示すように、その設定用ベクトルＳを行列設定部１により設定された初期行列Ｈ０（図６中、符号１４０のＨ０）の対角成分に設定し、設定後の初期行列Ｈ０を行列Ｈ（図６中、符号１４１のＨ）として出力する（図２のステップＳＴ８）。
即ち、対角成分設定部４は、全ての成分が０である疎行列の対角成分に設定用ベクトルＳを設定することで、図６中、符号１４８の疎行列を定義し、その疎行列を行列設定部１により設定された初期行列Ｈ０に足し合わせることで、図６中、符号１４１の行列Ｈを生成する。
即ち、対角成分設定部４は、全ての成分が０である疎行列の対角成分に設定用ベクトルＳを設定することで、図６中、符号１４８の疎行列を定義し、その疎行列を行列設定部１により設定された初期行列Ｈ０に足し合わせることで、図６中、符号１４１の行列Ｈを生成する。
到来波行列更新部５は、対角成分設定部４が行列Ｈの対角成分を設定すると、格納領域Ｐから疎成分に対する部分演算結果をロードし（図２のステップＳＴ９）、その部分演算結果を引き継ぐ形で、上記の式（１６）に示すように、連立方程式ＨＸ＝ｂを解くことで到来波行列Ｘ（図６中、符号１４５のＸ）を求め、その到来波行列Ｘを更新後の到来波行列Ｘ（新しい解）に設定する（図３のステップＳＴ１０）。
また、到来波行列更新部５は、疎成分に対する連立方程式の求解操作を部分演算結果として格納領域Ｐに格納する（図３のステップＳＴ１１）。
なお、疎成分に対する部分演算結果を引き継ぐ形で連立方程式を解いているので、到来波ベクトルｘに疎成分が多く現れるほど、より多くの演算結果を引き継ぐことができるため、連立方程式の求解操作に対する演算量が削減される。
以下、到来波行列更新部の処理内容を具体的に説明する。
また、到来波行列更新部５は、疎成分に対する連立方程式の求解操作を部分演算結果として格納領域Ｐに格納する（図３のステップＳＴ１１）。
なお、疎成分に対する部分演算結果を引き継ぐ形で連立方程式を解いているので、到来波ベクトルｘに疎成分が多く現れるほど、より多くの演算結果を引き継ぐことができるため、連立方程式の求解操作に対する演算量が削減される。
以下、到来波行列更新部の処理内容を具体的に説明する。
到来波行列更新部５は、対角成分設定部４が行列Ｈの対角成分を設定すると、下記の式（２１）に示すように、その行列Ｈを下三角行列Ｌ（上三角成分がゼロの行列）と上三角行列Ｕ（下三角行列がゼロの行列）との積（Ｈ＝ＬＵ）の形に分解する（ＬＵ分解操作）。
通常のＬＵ分解では、行列Ｈの左上隅の成分Ｈ１，１を除数として選択し、同じ列の成分（Ｈ２，１，Ｈ３，１…ＨＮ，１）を除算して行列Ｈを更新する。ここで、下記に示すように除算後のＨ２，１，Ｈ３，１…ＨＮ，１を（Ｎ―１）行１列の行列、Ｈ１，１と同じ行の成分Ｈ１，２，Ｈ１，３…Ｈ１，Ｎを１行（Ｎ−１）列の行列とする。
この時点で、対角成分Ｈ１，１と上記の行列成分は、以降の操作で更新されない成分として確定する。
このとき、上記の（Ｎ−１）行１列の行列と１行（Ｎ−１）列の行列との積をとった下記の（Ｎ−１）行（Ｎ−１）列の行列を、行列Ｈの残る成分値に対するマイナス値として、行列Ｈを更新する。
このとき、上記の（Ｎ−１）行１列の行列と１行（Ｎ−１）列の行列との積をとった下記の（Ｎ−１）行（Ｎ−１）列の行列を、行列Ｈの残る成分値に対するマイナス値として、行列Ｈを更新する。
Ｈ１，1が到来波ベクトルｘの疎成分に対応した対角成分であり、Ｈ１，1及びＨ１，1と同じ行と列の成分が格納領域Ｐに格納されていない場合は、上記の手順で値を計算し、確定したＨ１，1及びＨ１，1と同じ行と列の成分を格納領域Ｐに格納する。同様に、（Ｎ−１）行（Ｎ−１）列のマイナス値行列を上記の手順で計算し、格納領域Ｐに格納する。
一方、Ｈ１，1が到来波ベクトルｘの疎成分に対応した対角成分であり、Ｈ１，1及びＨ１，1と同じ列と行の成分が格納領域Ｐに格納されている場合は、上記の手順でＨ１，1及びＨ１，1と同じ行と列の成分を計算する代わりに、格納領域Ｐから対応する成分値を行列Ｈに代入する。
残る行列Ｈの対角成分に対しても、同様の操作を右下隅の方向へ順に繰り返し、格納領域Ｐから代入操作ができない対角成分（到来波ベクトルｘの非疎成分に対応した対角成分）に遭遇した時点で、行列Ｈの残る成分値（部分行列）を、格納領域Ｐに格納したマイナス値行列を適用して更新する。その後、残る対角成分に対して通常のＬＵ分解操作を進めていくことで、丁度、下三角行列Ｌと上三角行列Ｕをひとまとめに格納した行列Ｈが生成される。
一方、Ｈ１，1が到来波ベクトルｘの疎成分に対応した対角成分であり、Ｈ１，1及びＨ１，1と同じ列と行の成分が格納領域Ｐに格納されている場合は、上記の手順でＨ１，1及びＨ１，1と同じ行と列の成分を計算する代わりに、格納領域Ｐから対応する成分値を行列Ｈに代入する。
残る行列Ｈの対角成分に対しても、同様の操作を右下隅の方向へ順に繰り返し、格納領域Ｐから代入操作ができない対角成分（到来波ベクトルｘの非疎成分に対応した対角成分）に遭遇した時点で、行列Ｈの残る成分値（部分行列）を、格納領域Ｐに格納したマイナス値行列を適用して更新する。その後、残る対角成分に対して通常のＬＵ分解操作を進めていくことで、丁度、下三角行列Ｌと上三角行列Ｕをひとまとめに格納した行列Ｈが生成される。
到来波行列更新部５は、行列Ｈを下三角行列Ｌと上三角行列Ｕとの積（Ｈ＝ＬＵ）の形に分解した後の連立方程式ＬＵＸ＝ｂを、ｔ＝ＵＸなる行列ｔを定義することで、下記の式（２２）に示す連立方程式へと変形する。
ここで、下三角行列Ｌは上三角成分がゼロであるため、未知行列ｔの成分値は、下記の式（２３）で表されるように、代入操作（前進代入）によって徐々に求まっていくものである。
この結果、得られた行列ｔを使用して、同様の代入操作をＵＸ＝ｔに対して行えば、連立方程式の解である到来波行列Ｘが求まる。
ただし、Ｕは上三角行列であるため、代入操作は逆順（後進代入）となる。
ただし、Ｕは上三角行列であるため、代入操作は逆順（後進代入）となる。
収束判定部６は、到来波行列更新部５が到来波行列Ｘを更新（到来波行列Ｘに対する修整操作）すると、上記の式（１７）に示すように、更新後の到来波行列Ｘが収束しているか否かを判定する（図２のステップＳＴ１２）。
収束判定部６は、到来波行列Ｘが収束していると判定すれば、到来波行列Ｘに対する修整操作を終了するが、未だ到来波行列Ｘが収束していないと判定すれば、到来波ベクトル生成部２、対角成分設定部４及び到来波行列更新部５に対して反復計算指令を出力することで、到来波ベクトルｘの生成から到来波行列Ｘの更新までの処理を繰り返し実施させる。
収束判定部６は、到来波行列Ｘが収束していると判定すれば、到来波行列Ｘに対する修整操作を終了するが、未だ到来波行列Ｘが収束していないと判定すれば、到来波ベクトル生成部２、対角成分設定部４及び到来波行列更新部５に対して反復計算指令を出力することで、到来波ベクトルｘの生成から到来波行列Ｘの更新までの処理を繰り返し実施させる。
到来波ベクトル復元部７は、収束判定部６により到来波行列Ｘが収束していると判定された場合（到来波行列Ｘに対する修整操作が終了した場合）、格納領域Ｑ（図６中、符号１４４のＱ）に格納されている並べ替え情報を参照して、要素並び替え部３により並び替えられている到来波ベクトルｘの各要素ｘｋを元の位置に戻してから、当該到来波ベクトルｘを出力する。
ここで、図７を参照しながら、到来波ベクトルｘが疎ベクトルに修整されていく様子を説明する。
反復１周目の処理は、従来のスペクトル推定方法と同様であり、到来波ベクトルｘ（図７中、符号１２９のｘ）は、非疎成分（斜線部）で埋められており、この到来波ベクトルｘから設定用ベクトルＳ（図７中、符号１３０のＳ）が算出され、この設定用ベクトルＳが対角成分（図７中、符号１３２の対角成分）となる行列Ｈ（図７中、符号１３１のＨ）が生成される。
この行列Ｈから、ステップＳＴ１０での連立方程式の求解操作を経て、反復２周目の処理に移行する。
反復１周目の処理は、従来のスペクトル推定方法と同様であり、到来波ベクトルｘ（図７中、符号１２９のｘ）は、非疎成分（斜線部）で埋められており、この到来波ベクトルｘから設定用ベクトルＳ（図７中、符号１３０のＳ）が算出され、この設定用ベクトルＳが対角成分（図７中、符号１３２の対角成分）となる行列Ｈ（図７中、符号１３１のＨ）が生成される。
この行列Ｈから、ステップＳＴ１０での連立方程式の求解操作を経て、反復２周目の処理に移行する。
反復２周目以降の処理では、到来波ベクトルｘ（図７中、符号１３３のｘ）に含まれている疎成分の要素が左端に集まるように、各要素ｘｋの並び替えが行われる。到来波ベクトルｘの各要素ｘｋの並び替えに伴って、設定用ベクトルＳ（図７中、符号１３４のＳ）の定数成分の位置も左端に集まるため、行列Ｈ（図７中、符号１３５のＨ）の対角成分（図７中、符号１３６の対角成分）の定数部が左上隅に集まるようになる。連立方程式ＨＸ＝ｂの求解操作では、対角成分を左上隅から右下隅へと順に選択し、選択した対角成分と、その対角成分と同じ行と列の成分の値を計算により確定させ、確定した成分値を用いて残る部分行列を更新する。このため、行列Ｈの対角成分の定数部を左上隅に集めておくことで、対角成分の定数部に対する反復２週目以降の演算は、格納領域Ｐから確定済みの演算結果を代入する操作に置き換えることができる。
このような並べ替えを含む形での最適化計算（到来波行列Ｘに対する修整操作）が収束すると、右端に非疎成分（斜線部）を有する到来波ベクトルｘ（図７中、符号１３７のｘ）が得られる。
この到来波ベクトルｘの各要素ｘｋを並べ替え以前の順序に再整列することで、従来のスペクトル推定方法で得られる到来波ベクトルｘ（図５中、符号１２８のｘ）と同様の到来波ベクトルｘ（図７中、符号１３８のｘ）を得ることができる。
このような並べ替えを含む形での最適化計算（到来波行列Ｘに対する修整操作）が収束すると、右端に非疎成分（斜線部）を有する到来波ベクトルｘ（図７中、符号１３７のｘ）が得られる。
この到来波ベクトルｘの各要素ｘｋを並べ替え以前の順序に再整列することで、従来のスペクトル推定方法で得られる到来波ベクトルｘ（図５中、符号１２８のｘ）と同様の到来波ベクトルｘ（図７中、符号１３８のｘ）を得ることができる。
以下、到来波行列更新部５における到来波行列Ｘの更新処理の具体例を説明する。
ここでは、一例として、Ｎ＝３、下記の行列ＨをＬＵ分解する操作例を説明する。図８はＬＵ分解の操作例を示す説明図である。
なお、ＬＵ分解は、演算誤差の最小化を目的に、行や列を入れ替えるピボット操作を行うことが一般的であるが、説明を簡単にするため、ピボット操作を行わないＬＵ分解を使用するものとする。また、行列Ｈの左上隅の成分Ｈ１，１は、定数であるものとする。
ここでは、一例として、Ｎ＝３、下記の行列ＨをＬＵ分解する操作例を説明する。図８はＬＵ分解の操作例を示す説明図である。
まず、反復１周目のＬＵ分解では、行列Ｈ（図８中、符号１５０のＨ）の左上隅の成分Ｈ１，１を除数として選択し、同じ列の成分（Ｈ２，１，Ｈ３，１）を除算して行列Ｈ（図８中、符号１５１のＨ）へと更新する。
ここで、下記に示すように、Ｈ２，１，Ｈ３，１を２行１列の行列、Ｈ１，１と同じ行の成分Ｈ１，２，Ｈ１，３を１行２列の行列とする。
１行２列の行列＝（５ ７）
ここで、下記に示すように、Ｈ２，１，Ｈ３，１を２行１列の行列、Ｈ１，１と同じ行の成分Ｈ１，２，Ｈ１，３を１行２列の行列とする。
このとき、上記の２行１列の行列と１行２列の行列との積をとった下記の２行２列の行列は、行列Ｈ（図８中、符号１５１のＨ）の残る成分値に対するマイナス値となり、行列Ｈ（図８中、符号１５２のＨ）へと更新する。
これと同様の操作を、行列Ｈの左上隅から右下隅へ順に、対角成分を選択して繰り返すことで、丁度、下三角行列Ｌ（図８中、符号１５４のＬ）と上三角行列Ｕ（図８中、符号１５５のＵ）をひとまとめに格納した行列Ｈ（図８中、符号１５３のＨ）が生成される。
このとき、対角成分の定数を基準に、行と列の成分（図８中、符号１５６の成分）と、マイナス値の行列（図８中、符号１５７の行列）を部分演算結果として格納領域Ｐに格納する。
このとき、対角成分の定数を基準に、行と列の成分（図８中、符号１５６の成分）と、マイナス値の行列（図８中、符号１５７の行列）を部分演算結果として格納領域Ｐに格納する。
次に、反復２周目で、行列Ｈ（図８中、符号１５０のＨ）の対角成分が更新された下記の行列Ｈ（図８中、符号１５８のＨ）をＬＵ分解する操作について説明する。
反復１周目に引き続き対角成分Ｈ１，１は定数である。
格納領域Ｐに格納されている部分演算結果のうち、計算結果（図８中、符号１５６の計算結果）を行列Ｈ（図８中、符号１５８のＨ）に代入して、行列Ｈ（図８中、符号１５９のＨ）へと更新する。
その後、マイナス値の行列（図８中、符号１５７の行列）を行列Ｈ（図８中、符号１５９のＨ）に適用して、行列Ｈ（図８中、符号１６０のＨ）へと更新する。
この時点で、対角成分Ｈ１，１に関するＬＵ分解操作が完了しているため、次の対角成分Ｈ２，２から通常のＬＵ分解を行うことで、丁度、下三角行列Ｌ（図８中、符号１６２のＬ）と上三角行列Ｕ（図８中、符号１６３のＵ）をひとまとめに格納した行列Ｈ（図８中、符号１６１のＨ）が生成される。
その後、マイナス値の行列（図８中、符号１５７の行列）を行列Ｈ（図８中、符号１５９のＨ）に適用して、行列Ｈ（図８中、符号１６０のＨ）へと更新する。
この時点で、対角成分Ｈ１，１に関するＬＵ分解操作が完了しているため、次の対角成分Ｈ２，２から通常のＬＵ分解を行うことで、丁度、下三角行列Ｌ（図８中、符号１６２のＬ）と上三角行列Ｕ（図８中、符号１６３のＵ）をひとまとめに格納した行列Ｈ（図８中、符号１６１のＨ）が生成される。
更に、以上の手順で、部分演算結果を引き継ぐＬＵ分解を行うと、連立方程式を解くための代入操作まで省略することができる。
下記の式（２６）において、対角成分Ｈ１，１を基準とした行と列（Ｈ１，１〜Ｈ１，Ｎ、Ｈ１，１〜ＨＮ，１）の演算結果が定数として格納領域Ｐに格納されている場合、下三角行列Ｌの列成分Ｌ１，１〜ＬＮ，１と上三角行列Ｕの行成分Ｕ１，１〜Ｕ１，Ｎも定数となる。
下記の式（２６）において、対角成分Ｈ１，１を基準とした行と列（Ｈ１，１〜Ｈ１，Ｎ、Ｈ１，１〜ＨＮ，１）の演算結果が定数として格納領域Ｐに格納されている場合、下三角行列Ｌの列成分Ｌ１，１〜ＬＮ，１と上三角行列Ｕの行成分Ｕ１，１〜Ｕ１，Ｎも定数となる。
即ち、代入操作で解く下記の式（２７）の連立方程式において、代入操作で要する乗算が、定数部分については、乗算結果の代入操作で置き換えられる。つまり、部分演算結果として、この定数部分に対する乗算結果も併せて格納領域Ｐに格納しておくようにする。
以上で明らかなように、この実施の形態１によれば、行列設定部１により設定された到来波行列Ｘを積み上げることで到来波ベクトルｘを生成する到来波ベクトル生成部２と、到来波ベクトル生成部２により生成された到来波ベクトルｘに含まれている各要素ｘｋのうち、疎成分の要素が左端に集まるように各要素ｘｋを並び替える要素並び替え部３と、要素並び替え部３による並び替え後の到来波ベクトルｘから定数成分が左端に集まっている設定用ベクトルＳを算出し、その設定用ベクトルＳによって行列Ｈの対角成分を設定する対角成分設定部４と、対角成分設定部４により対角成分が設定された行列Ｈと行列ｂを用いて、その到来波行列Ｘを更新する到来波行列更新部５とを設けるように構成したので、角度点数Ｎが大きく、到来波行列Ｘの行列サイズが大きくても、短時間で到来波行列Ｘを推定することができる効果を奏する。
即ち、この実施の形態１によれば、最適化の進行に伴って、到来波ベクトルｘに疎成分が増えるほど、連立方程式の求解操作手順を省略することができる。この結果、反復１回毎の処理時間を徐々に短くすることができるため、従来より短時間で到来波行列Ｘの推定結果を得ることができる。
また、連立方程式を解く行列サイズを徐々に小さくすることができるため、到来波推定処理へ入力する行列サイズを大きくした場合でも、演算時間の爆発的な増加を抑制することができる。
また、連立方程式を解く行列サイズを徐々に小さくすることができるため、到来波推定処理へ入力する行列サイズを大きくした場合でも、演算時間の爆発的な増加を抑制することができる。
なお、この実施の形態１では、要素並び替え部３が、到来波ベクトル生成部２により生成された到来波ベクトルｘに含まれている各要素ｘｋのうち、疎成分の要素が左端に集まるように各要素ｘｋを並び替えるものを示したが、予め設定された閾値（例えば、１０−４）より小さな値を有する要素については、疎成分の要素に置き換えてから、疎成分の要素が左端に集まるように各要素ｘｋを並び替えるようにしてもよい。
また、この実施の形態１では、到来波行列更新部５が、式（１６）の連立方程式を解くことで到来波行列Ｘを更新するものを示したが、行列Ｈの行列サイズに応じて、到来波行列Ｘの更新演算の演算精度を切り替えるようにしてもよい。
即ち、この実施の形態１では、最適化の進行に伴って、扱う行列サイズ（行列Ｈのサイズ、到来波行列Ｘのサイズ）が縮小していく特徴を有しており、行列サイズが小さい場合には、精度の低い演算処理を行っても、演算誤差が小さくなるので、行列Ｈの行列サイズが一定の行列サイズ（例えば、１００×１００）以下となった時点で、例えば、倍精度の演算処理から単精度の演算処理に切り替えることで、高速な処理を実現するようにしてもよい。
即ち、この実施の形態１では、最適化の進行に伴って、扱う行列サイズ（行列Ｈのサイズ、到来波行列Ｘのサイズ）が縮小していく特徴を有しており、行列サイズが小さい場合には、精度の低い演算処理を行っても、演算誤差が小さくなるので、行列Ｈの行列サイズが一定の行列サイズ（例えば、１００×１００）以下となった時点で、例えば、倍精度の演算処理から単精度の演算処理に切り替えることで、高速な処理を実現するようにしてもよい。
また、この実施の形態１では、図２のステップＳＴ４〜ＳＴ１２で反復計算を行うものを示したが、到来波行列更新部５が、到来波行列Ｘの更新を並列処理で行うようにしてもよい。
即ち、逐次的な反復処理である未知行列（到来波行列Ｘ）の修整操作（行列Ｈの対角成分を更新して連立方程式を解く操作）に対して、反復１回当りの未知行列の修整を、複数のＣＰＵの演算単位で並列処理する手段を備え、未知行列を並列に修整する（連立方程式を並列に解く）ようにしてもよい。
あるいは、グラフィックス用プロセッサであるＧＰＵで並列処理する手段を備えて、未知行列を並列に修整する（連立方程式を並列に解く）ようにしてもよい。
即ち、逐次的な反復処理である未知行列（到来波行列Ｘ）の修整操作（行列Ｈの対角成分を更新して連立方程式を解く操作）に対して、反復１回当りの未知行列の修整を、複数のＣＰＵの演算単位で並列処理する手段を備え、未知行列を並列に修整する（連立方程式を並列に解く）ようにしてもよい。
あるいは、グラフィックス用プロセッサであるＧＰＵで並列処理する手段を備えて、未知行列を並列に修整する（連立方程式を並列に解く）ようにしてもよい。
実施の形態２．
上記実施の形態１では、到来波行列更新部５が連立方程式を解く最適化アルゴリズムで構成され、格納領域Ｐから疎成分に対する部分演算結果をロードし、その部分演算結果を引き継ぐ形で演算量を削減する方式を示したが、到来波行列更新部５が、対角成分設定部４により対角成分が設定された行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を連立方程式の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで連立方程式を解くようにしてもよい。
上記実施の形態１では、到来波行列更新部５が連立方程式を解く最適化アルゴリズムで構成され、格納領域Ｐから疎成分に対する部分演算結果をロードし、その部分演算結果を引き継ぐ形で演算量を削減する方式を示したが、到来波行列更新部５が、対角成分設定部４により対角成分が設定された行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を連立方程式の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで連立方程式を解くようにしてもよい。
図９はこの発明の実施の形態２によるスペクトル推定装置の処理内容（スペクトル推定方法）を示すフローチャートである。
この実施の形態２では、上記実施の形態１と比べると、行列設定部１が部分演算結果の格納領域Ｐを初期化する処理（ステップＳＴ２）と、到来波行列更新部５が格納領域Ｐから疎成分に対する部分演算結果をロードする処理（ステップＳＴ９）と、到来波行列更新部５が疎成分に対する連立方程式の求解操作を部分演算結果として格納領域Ｐに格納する処理（ステップＳＴ１１）とがなくなっている。
また、この実施の形態２では、到来波行列更新部５が、行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を連立方程式の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで連立方程式を解くようにしている（図９のステップＳＴ２０）。
この実施の形態２では、上記実施の形態１と比べると、行列設定部１が部分演算結果の格納領域Ｐを初期化する処理（ステップＳＴ２）と、到来波行列更新部５が格納領域Ｐから疎成分に対する部分演算結果をロードする処理（ステップＳＴ９）と、到来波行列更新部５が疎成分に対する連立方程式の求解操作を部分演算結果として格納領域Ｐに格納する処理（ステップＳＴ１１）とがなくなっている。
また、この実施の形態２では、到来波行列更新部５が、行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を連立方程式の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで連立方程式を解くようにしている（図９のステップＳＴ２０）。
図１０はこの実施の形態２における到来波行列更新部５の更新処理を示す説明図である。
図１０はＮ＝３で、行列Ｈ（図１０中、符号２００のＨ）が下記の場合を示している。
この例では、対角成分Ｈ１，１が定数であるものとして扱う。
図１０はＮ＝３で、行列Ｈ（図１０中、符号２００のＨ）が下記の場合を示している。
この実施の形態２では、到来波行列更新部５が到来波ベクトルｘの疎成分の要素に対応する到来波行列Ｘの定数成分を演算対象から除外するため、対角成分Ｈ１，１を基点とする行成分Ｈ１，１〜Ｈ１，Ｎと列成分Ｈ１，１〜ＨＮ，１を削除する。
これにより、行列Ｈ（図１０中、符号２００のＨ）は、行列Ｈ（図１０中、符号２０１のＨ）へと更新される。
到来波行列更新部５は、連立方程式内の行列Ｈが符号２０１の行列Ｈであるとして、ＬＵ分解やコレスキー分解などを用いた方法で、連立方程式を解くようにする。
これにより、行列Ｈ（図１０中、符号２００のＨ）は、行列Ｈ（図１０中、符号２０１のＨ）へと更新される。
なお、この実施の形態２で示している未知行列である到来波行列Ｘに対する修整操作を、任意の反復タイミングで、上記実施の形態１で示している修整操作に切り替えるようにしてもよい。
逆に、上記実施の形態１で示している修整操作を、任意の反復タイミングで、上記実施の形態２で示している修整操作に切り替えるようにしてもよい。
また、連立方程式を解く方法（ＬＵ分解やコレスキー分解など）は、任意の反復タイミングで切り替えても構わない。連立方程式を解いた後に、演算誤差を判定し、別の連立方程式を解く方法で、解き直しても構わない。
逆に、上記実施の形態１で示している修整操作を、任意の反復タイミングで、上記実施の形態２で示している修整操作に切り替えるようにしてもよい。
また、連立方程式を解く方法（ＬＵ分解やコレスキー分解など）は、任意の反復タイミングで切り替えても構わない。連立方程式を解いた後に、演算誤差を判定し、別の連立方程式を解く方法で、解き直しても構わない。
以上で明らかなように、この実施の形態２によれば、到来波行列更新部５が、対角成分設定部４により対角成分が設定された行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を連立方程式の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで連立方程式を解くように構成しているので、上記実施の形態１と同様の効果を奏する他に、最適化の進行に伴って、到来波ベクトルｘに疎成分が増えるほど、連立方程式を解く解の数を削減して、処理の高速化を図ることができる効果を奏する。
実施の形態３．
上記実施の形態１，２では、到来波行列更新部５が、連立方程式を解く最適化アルゴリズムで構成されているものを示したが、到来波行列更新部５の手段として、連立方程式の計算を必要としない最適化アルゴリズムを利用し、その際の行列積計算に対して疎成分に対応する演算を省略するようにしても構わない。
この実施の形態３では、連立方程式の計算を必要としない最適化アルゴリズムに対する実施例を説明する便宜上、「最急降下法」の操作手順を例に説明するが、最急降下法と同様に連立方程式の計算を必要としない他の最適化アルゴリズム（例えば、「非線形共役勾配法」やヘッセ逆行列を直接的に近似した「準ニュートン法」など）を使用するようにしてもよい。
上記実施の形態１，２では、到来波行列更新部５が、連立方程式を解く最適化アルゴリズムで構成されているものを示したが、到来波行列更新部５の手段として、連立方程式の計算を必要としない最適化アルゴリズムを利用し、その際の行列積計算に対して疎成分に対応する演算を省略するようにしても構わない。
この実施の形態３では、連立方程式の計算を必要としない最適化アルゴリズムに対する実施例を説明する便宜上、「最急降下法」の操作手順を例に説明するが、最急降下法と同様に連立方程式の計算を必要としない他の最適化アルゴリズム（例えば、「非線形共役勾配法」やヘッセ逆行列を直接的に近似した「準ニュートン法」など）を使用するようにしてもよい。
図１１はこの発明の実施の形態３によるスペクトル推定装置の処理内容（スペクトル推定方法）を示すフローチャートである。
この実施の形態３では、上記実施の形態２と比べると、到来波行列更新部５が、式（４）に基づき、到来波行列Ｘを評価する目的関数の勾配▽ｆ（Ｘ）を行列積ＨＸと行列ｂとの差分として計算し、その計算した勾配▽ｆ（Ｘ）から式（１０）、式（１１）、式（１２）に基づき、到来波行列Ｘを更新する形式に置き換わっている。
即ち、上記実施の形態２における図９のフローチャートのうち、連立方程式の計算による到来波行列Ｘの更新処理（ステップＳＴ２０）が、この実施の形態３における図１１のフローチャートでは、行列積ＨＸの計算を中心とする到来波行列Ｘの更新処理（ステップＳＴ３０）に置き換わっている。このステップＳＴ３０以外の処理内容については、上記実施の形態２と同様である。
また、この実施の形態３では、到来波行列更新部５が、行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を行列積ＨＸの演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで行列積を計算するようにしている。
この実施の形態３では、上記実施の形態２と比べると、到来波行列更新部５が、式（４）に基づき、到来波行列Ｘを評価する目的関数の勾配▽ｆ（Ｘ）を行列積ＨＸと行列ｂとの差分として計算し、その計算した勾配▽ｆ（Ｘ）から式（１０）、式（１１）、式（１２）に基づき、到来波行列Ｘを更新する形式に置き換わっている。
即ち、上記実施の形態２における図９のフローチャートのうち、連立方程式の計算による到来波行列Ｘの更新処理（ステップＳＴ２０）が、この実施の形態３における図１１のフローチャートでは、行列積ＨＸの計算を中心とする到来波行列Ｘの更新処理（ステップＳＴ３０）に置き換わっている。このステップＳＴ３０以外の処理内容については、上記実施の形態２と同様である。
また、この実施の形態３では、到来波行列更新部５が、行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を行列積ＨＸの演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで行列積を計算するようにしている。
図１２はこの発明の実施の形態３によるスペクトル推定装置での到来波推定の処理フローを図解している説明図である。
上記実施の形態２と同様の手順で、行列Ｈ（図１２の符号３２３）、行列Ｘ（図１２の符号３２０）及びベクトルｘ（図１２の符号３２１）は、到来波ベクトル生成部２、行列設定部１、要素並び替え部３によって、疎成分（図１２の符号３２０から符号３２３までの斜線領域）と非疎成分（図１２の符号３２０から符号３２３の白色領域）とがそれぞれ集合するように並べ替えが行われている。
上記実施の形態２と同様の手順で、行列Ｈ（図１２の符号３２３）、行列Ｘ（図１２の符号３２０）及びベクトルｘ（図１２の符号３２１）は、到来波ベクトル生成部２、行列設定部１、要素並び替え部３によって、疎成分（図１２の符号３２０から符号３２３までの斜線領域）と非疎成分（図１２の符号３２０から符号３２３の白色領域）とがそれぞれ集合するように並べ替えが行われている。
到来波行列Ｘの更新処理において、勾配▽ｆ（Ｘ）を求めるための行列積ＨＸの計算では、ｎ行ｍ列目の要素Ｃｎ，ｍが、下記の式（３０）に示すように、行列Ｈのｎ行目の行ベクトルと、行列Ｘのｍ列目の列ベクトルとの積和計算によって求められる。
この際、到来波行列Ｘにおいて疎となった領域は、以降の反復計算中も常に疎として扱うことで、その後の反復計算では、疎と判明した領域に対する値の修整処理が必要なくなる。このため、行列積ＨＸの計算領域のうち、行列Ｘの疎成分に対応する領域（図１２の符号３２４の白色領域）の計算については全て省略しても到来波行列Ｘの更新に支障がなく、演算量を削減した勾配▽ｆ（Ｘ）の計算が実施される。
更に、到来波行列Ｘの非疎成分に対応する領域についても、Ｘｉ，ｍの疎となった領域（図１２の符号３２０）に対するＨｎ，ｉとの積の結果はゼロであるため、式（３０）における積和計算のうち、行列Ｘの疎となった領域に対応する積和計算を省略して扱う。このため、到来波行列Ｘに疎成分が増加していくにつれて、対応する行列Ｈの計算不要な領域（図１２の符号３２３の白色領域）が増加し、演算量を削減した勾配▽ｆ（Ｘ）の計算が実施される。
更に、到来波行列Ｘの非疎成分に対応する領域についても、Ｘｉ，ｍの疎となった領域（図１２の符号３２０）に対するＨｎ，ｉとの積の結果はゼロであるため、式（３０）における積和計算のうち、行列Ｘの疎となった領域に対応する積和計算を省略して扱う。このため、到来波行列Ｘに疎成分が増加していくにつれて、対応する行列Ｈの計算不要な領域（図１２の符号３２３の白色領域）が増加し、演算量を削減した勾配▽ｆ（Ｘ）の計算が実施される。
以上で明らかなように、この実施の形態３によれば、到来波行列更新部５が、対角成分設定部４により対角成分が設定された行列Ｈに含まれている到来波ベクトルｘの疎成分の要素に対応する行列成分を行列積の演算対象から除外し、その疎成分の要素に対応する行列成分以外の行列成分のみで行列積計算を解くように構成しているので、連立方程式の計算を必要としない最適化アルゴリズムを利用した場合においても、上記実施の形態１，２と同様に、最適化の進行に伴って、到来波ベクトルｘに疎成分が増えるほど、行列積の演算量を削減して、処理の高速化を図ることができる効果を奏する。
実施の形態４．
上記実施の形態１〜３では、スペクトル推定装置が、空間方向に対する電力スペクトルである到来波（アンテナ素子の受信信号に含まれている到来波）を推定するものを示したが、ｙ＝ＡＸ＋ｎとしてモデル化できるスペクトル推定問題であれば、汎用的に当てはめて利用して、未知行列Ｘの推定を実施することができる。
この実施の形態４では、その一例として、周波数スペクトルである遅延時間の推定について説明する。
上記実施の形態１〜３では、スペクトル推定装置が、空間方向に対する電力スペクトルである到来波（アンテナ素子の受信信号に含まれている到来波）を推定するものを示したが、ｙ＝ＡＸ＋ｎとしてモデル化できるスペクトル推定問題であれば、汎用的に当てはめて利用して、未知行列Ｘの推定を実施することができる。
この実施の形態４では、その一例として、周波数スペクトルである遅延時間の推定について説明する。
遅延時間の推定では、時間サンプル数Ｍで観測した時間軸上の受信信号（時間波形）をＭ×１の観測行列ｙとし、これに対応する周波数軸上での時間周波数をＭ×１の行列Ｘとすることで、例えば、レーダから目標に対して放射したパルス波が反射して返ってくるまでの遅延時間を求める。このとき、辞書行列Ａは、下記の式（３１）で表される。
上記実施の形態１〜３では、未知行列Ｘが到来波行列であり、その到来波行列Ｘを積み上げることで到来波ベクトルｘが生成されるが、この実施の形態４では、未知行列Ｘが遅延時間行列であり、その遅延時間行列を積み上げることで遅延時間ベクトルｘが生成される。
また、辞書行列Ａ及び未知行列Ｘと観測行列ｙは、上記の式（１）によって関連付けられる。
ここで、未知行列Ｘを推定するために、上記実施の形態１〜３の観測条件に合うように、時間サンプル数Ｍを１に、角度点数Ｎを実施の形態４における時間サンプル数Ｍに、アンテナ素子数Ｌを実施の形態４における時間サンプル数Ｍに設定する。
これにより、上記実施の形態１〜３のいずれかに示した方法で、この実施の形態４における周波数軸上での時間周波数のスペクトルＸが求まる。
ここで、未知行列Ｘを推定するために、上記実施の形態１〜３の観測条件に合うように、時間サンプル数Ｍを１に、角度点数Ｎを実施の形態４における時間サンプル数Ｍに、アンテナ素子数Ｌを実施の形態４における時間サンプル数Ｍに設定する。
これにより、上記実施の形態１〜３のいずれかに示した方法で、この実施の形態４における周波数軸上での時間周波数のスペクトルＸが求まる。
なお、本願発明はその発明の範囲内において、各実施の形態の自由な組み合わせ、あるいは各実施の形態の任意の構成要素の変形、もしくは各実施の形態において任意の構成要素の省略が可能である。
この発明は、角度点数が大きく、スペクトル行列Ｘの行列サイズが大きくても、短時間でスペクトル行列Ｘを推定する必要があるスペクトル推定装置及びスペクトル推定方法に適している。
１ 行列設定部（行列設定手段）、２ 到来波ベクトル生成部（ベクトル生成手段）、３ 要素並び替え部（並び替え手段）、４ 対角成分設定部（対角成分設定手段）、５ 到来波行列更新部（行列更新手段）、６ 収束判定部（収束判定手段）、７ 到来波ベクトル復元部（収束判定手段）。
Claims (10)
- １以上のアンテナ素子の走査ベクトルからなる辞書行列の複素共役転置と前記辞書行列との積の行列Ｈと、前記アンテナ素子の受信信号からなる観測行列と前記辞書行列の複素共役転置との積の行列ｂと、前記受信信号に含まれているスペクトルに関する未知行列であるスペクトル行列とを設定する行列設定手段と、
前記スペクトル行列を積み上げることでスペクトルベクトルを生成するベクトル生成手段と、
前記ベクトル生成手段により生成されたスペクトルベクトルに含まれている各要素のうち、疎成分の要素が集まるように前記各要素を並び替える並び替え手段と、
前記並び替え手段による並び替え後のスペクトルベクトルから定数成分が集まった設定用ベクトルを算出し、前記設定用ベクトルによって前記行列Ｈの対角成分を設定する対角成分設定手段と、
前記対角成分設定手段により対角成分が設定された行列Ｈと前記行列ｂを用いて、前記スペクトル行列を更新する行列更新手段と、
前記行列更新手段により更新されたスペクトル行列が収束するまで、前記スペクトルベクトルの生成から前記スペクトル行列の更新までの処理を繰り返し実施させる収束判定手段と
を備えたスペクトル推定装置。 - 前記収束判定手段は、前記行列更新手段により更新されたスペクトル行列が収束している場合、前記並び替え手段により並び替えられているスペクトルベクトルの各要素を元の位置に戻してから、当該スペクトルベクトルを出力することを特徴とする請求項１記載のスペクトル推定装置。
- 前記行列更新手段は、前記ベクトル生成手段により生成されたスペクトルベクトルに含まれている各要素のうち、予め設定された閾値より小さな値を有する要素については、疎成分の要素に置き換えてから、疎成分の要素が集まるように前記各要素を並び替えることを特徴とする請求項１記載のスペクトル推定装置。
- 前記行列更新手段は、前記対角成分設定手段により対角成分が設定された行列Ｈと前記行列ｂを用いて、前記スペクトル行列を更新する際、前記行列Ｈの行列サイズに応じて、前記スペクトル行列の更新演算の演算精度を切り替えることを特徴とする請求項１記載のスペクトル推定装置。
- 前記行列更新手段は、前記対角成分設定手段により対角成分が設定された行列Ｈと前記スペクトル行列の積が前記行列ｂと一致する連立方程式を解くことで更新後の前記スペクトル行列を求め、前記スペクトルベクトルの疎成分の要素に対応する行列成分への演算結果を格納して、２回目以降の連立方程式を解く際に再利用することを特徴とする請求項１記載のスペクトル推定装置。
- 前記行列更新手段は、前記対角成分設定手段により対角成分が設定された行列Ｈに含まれている前記スペクトルベクトルの疎成分の要素に対応する行列成分を前記連立方程式の演算対象から除外し、前記疎成分の要素に対応する行列成分以外の行列成分のみで前記連立方程式を解くことを特徴とする請求項５記載のスペクトル推定装置。
- 前記行列更新手段は、前記連立方程式の求解をＬＵ分解又はコレスキー分解の操作で実現することを特徴とする請求項５記載のスペクトル推定装置。
- 前記行列更新手段は、スペクトル行列の更新を並列処理で行うことを特徴とする請求項１記載のスペクトル推定装置。
- 前記行列設定手段は、前記受信信号に含まれているスペクトルに関するスペクトル行列として、前記受信信号に含まれている空間方向に対する電力スペクトルである到来波に関する到来波行列を設定し、
前記ベクトル生成手段は、前記到来波行列を積み上げることで到来波ベクトルを生成し、
前記並び替え手段は、前記ベクトル生成手段により生成された到来波ベクトルに含まれている各要素のうち、疎成分の要素が集まるように前記各要素の並び替えを実施し、
前記対角成分設定手段は、前記並び替え手段による並び替え後の到来波ベクトルから定数成分が集まった設定用ベクトルを算出し、前記設定用ベクトルによって前記行列Ｈの対角成分を設定し、
前記行列更新手段は、前記対角成分設定手段により対角成分が設定された行列Ｈと前記行列ｂを用いて、前記到来波行列を更新し、
前記収束判定手段は、前記行列更新手段により更新された到来波行列が収束するまで、前記到来波ベクトルの生成から前記到来波行列の更新までの処理を繰り返し実施させることを特徴とする請求項１記載のスペクトル推定装置。 - 行列設定手段が、１以上のアンテナ素子の走査ベクトルからなる辞書行列の複素共役転置と前記辞書行列との積の行列Ｈと、前記アンテナ素子の受信信号からなる観測行列と前記辞書行列の複素共役転置との積の行列ｂと、前記受信信号に含まれているスペクトルに関する未知行列であるスペクトル行列とを設定する行列設定処理ステップと、
ベクトル生成手段が、前記スペクトル行列を積み上げることでスペクトルベクトルを生成するベクトル生成処理ステップと、
並び替え手段が、前記ベクトル生成処理ステップで生成されたスペクトルベクトルに含まれている各要素のうち、疎成分の要素が集まるように前記各要素を並び替える並び替え処理ステップと、
対角成分設定手段が、前記並び替え処理ステップでの並び替え後のスペクトルベクトルから定数成分が集まった設定用ベクトルを算出し、前記設定用ベクトルによって前記行列Ｈの対角成分を設定する対角成分設定処理ステップと、
行列更新手段が、前記対角成分設定処理ステップで対角成分が設定された行列Ｈと前記行列ｂを用いて、前記スペクトル行列を更新する行列更新処理ステップと、
収束判定手段が、前記行列更新処理ステップで更新されたスペクトル行列が収束するまで、前記スペクトルベクトルの生成から前記スペクトル行列の更新までの処理を繰り返し実施させる収束判定処理ステップと
を備えたスペクトル推定方法。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2013262632 | 2013-12-19 | ||
| JP2013262632 | 2013-12-19 | ||
| PCT/JP2014/059446 WO2015093068A1 (ja) | 2013-12-19 | 2014-03-31 | スペクトル推定装置及びスペクトル推定方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP5985079B2 true JP5985079B2 (ja) | 2016-09-06 |
| JPWO2015093068A1 JPWO2015093068A1 (ja) | 2017-03-16 |
Family
ID=53402429
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2015553383A Expired - Fee Related JP5985079B2 (ja) | 2013-12-19 | 2014-03-31 | スペクトル推定装置及びスペクトル推定方法 |
Country Status (2)
| Country | Link |
|---|---|
| JP (1) | JP5985079B2 (ja) |
| WO (1) | WO2015093068A1 (ja) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105911349B (zh) * | 2016-05-31 | 2019-01-11 | 清华大学 | 基于重排时频谱的线性扫频信号基本参数估算方法及装置 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101980043A (zh) * | 2010-09-15 | 2011-02-23 | 电子科技大学 | 一种抗接收机相位跳变的多干扰源测向方法 |
| JP2011163922A (ja) * | 2010-02-09 | 2011-08-25 | Toshiba Corp | 信号到来方向推定方法 |
| US20120268309A1 (en) * | 2011-02-10 | 2012-10-25 | The Arizona Board Of Regents On Behalf Of The University Of Arizona | Virtual aperture radar (var) imaging |
| JP2013174584A (ja) * | 2012-02-27 | 2013-09-05 | Mitsubishi Electric Corp | シーン内の反射物を再構成する方法 |
| JP2013234871A (ja) * | 2012-05-07 | 2013-11-21 | Mitsubishi Electric Corp | 観測装置 |
-
2014
- 2014-03-31 JP JP2015553383A patent/JP5985079B2/ja not_active Expired - Fee Related
- 2014-03-31 WO PCT/JP2014/059446 patent/WO2015093068A1/ja active Application Filing
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2011163922A (ja) * | 2010-02-09 | 2011-08-25 | Toshiba Corp | 信号到来方向推定方法 |
| CN101980043A (zh) * | 2010-09-15 | 2011-02-23 | 电子科技大学 | 一种抗接收机相位跳变的多干扰源测向方法 |
| US20120268309A1 (en) * | 2011-02-10 | 2012-10-25 | The Arizona Board Of Regents On Behalf Of The University Of Arizona | Virtual aperture radar (var) imaging |
| JP2013174584A (ja) * | 2012-02-27 | 2013-09-05 | Mitsubishi Electric Corp | シーン内の反射物を再構成する方法 |
| JP2013234871A (ja) * | 2012-05-07 | 2013-11-21 | Mitsubishi Electric Corp | 観測装置 |
Non-Patent Citations (4)
| Title |
|---|
| JPN6014024286; 後町将人、外２名: '"高分解能ＤｏＡ推定法に対する高速実装方式の検討"' ２０１４年電子情報通信学会総合大会講演論文集 , 20140304, ｐ．２６３, 一般社団法人電子情報通信学会 * |
| JPN6014024287; Ｄｍｉｔｒｙ Ｍａｌｉｏｕｔｏｖ， ｅｔ ａｌ．: '"Ａ ｓｐａｒｓｅ ｓｉｇｎａｌ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｐｅｒｓｐｅｃｔｉｖｅ ｆｏｒ ｓｏｕｒｃｅ ｌｏｃａｌｉｚａｔｉｏｎ ｗｉｔｈ ｓｅｎｓｏｒ ａｒｒａｙｓ"' Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ， ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｖｏｌ．５３，Ｎｏ．８, 200508, ｐ．３０１０-３０２２ * |
| JPN6014024288; Ｍｕｊｄａｔ Ｃｅｔｉｎ， ｅｔ ａｌ．: '"Ａ ｖａｒｉａｔｉｏｎａｌ ｔｅｃｈｎｉｑｕｅ ｆｏｒ ｓｏｕｒｃｅ ｌｏｃａｌｉｚａｔｉｏｎ ｂａｓｅｄ ｏｎ ａ ｓｐａｒｓｅ ｓｉｇｎａｌ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｐｅｒｓｐｅｃｔｉｖｅ' Ａｃｏｕｓｔｉｃｓ， Ｓｐｅｅｃｈ， ａｎｄ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ （ＩＣＡＳＳＰ）， ２００２ ＩＥＥＥ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｖｏｌ．３, 20020513, ｐ．２９６５-２９６８ * |
| JPN6016025471; Ａｍｉｒ Ｂｅｃｋ， ｅｔ ａｌ．: '"Ａ Ｆａｓｔ Ｉｔｅｒａｔｉｖｅ Ｓｈｒｉｎｋａｇｅ-Ｔｈｒｅｓｈｏｌｄｉｎｇ Ａｌｇｏｒｉｔｈｍ ｆｏｒ Ｌｉｎｅａｒ Ｉｎｖｅｒｓｅ Ｐｒｏｂｌｅｍｓ"' ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｉｍａｇｉｎｇ Ｓｃｉｅｎｃｅ Ｖｏｌ．２， Ｎｏ．１, 2009, ｐｐ．１８３-２０２ * |
Also Published As
| Publication number | Publication date |
|---|---|
| JPWO2015093068A1 (ja) | 2017-03-16 |
| WO2015093068A1 (ja) | 2015-06-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bai et al. | High-speed compressed sensing reconstruction on FPGA using OMP and AMP | |
| JP6615062B2 (ja) | 画像を処理する方法及びシステム | |
| Yang et al. | Compressed sensing and Cholesky decomposition on FPGAs and GPUs | |
| CN107450045B (zh) | 基于focuss二次加权算法的doa估计方法 | |
| CN110174658B (zh) | 基于秩一降维模型和矩阵补全的波达方向估计方法 | |
| CN110244272B (zh) | 基于秩一去噪模型的波达方向估计方法 | |
| JP6953287B2 (ja) | 音源探査装置、音源探査方法およびそのプログラム | |
| US20210103849A1 (en) | Preparing and converting gottesman-kitaev-preskill states | |
| Jiang et al. | An efficient FPGA-based direct linear solver | |
| KR102026958B1 (ko) | 레이더용 압축 센싱을 위한 새로운 분할 역변환 기법을 사용한 감소된 계산 복잡성의 오엠피 방법 | |
| Niu et al. | SPEC2: Spectral sparse CNN accelerator on FPGAs | |
| JP5985079B2 (ja) | スペクトル推定装置及びスペクトル推定方法 | |
| CN110941980A (zh) | 一种密集环境中基于压缩感知的多径时延估计方法及装置 | |
| Mohanty et al. | Design and performance analysis of fixed-point jacobi svd algorithm on reconfigurable system | |
| Amin-Nejad et al. | Floating point versus fixed point tradeoffs in FPGA Implementations of QR Decomposition Algorithm | |
| US10713332B2 (en) | Hardware accelerated linear system solver | |
| Eidelman et al. | On the fast reduction of a quasiseparable matrix to Hessenberg and tridiagonal forms | |
| Liu et al. | An FPGA-based MVDR beamformer using dichotomous coordinate descent iterations | |
| US7945061B1 (en) | Scalable architecture for subspace signal tracking | |
| Prabhu | Fixed point pipelined architecture for QR decomposition | |
| Siddhartha et al. | Long short-term memory for radio frequency spectral prediction and its real-time fpga implementation | |
| JP2008269329A (ja) | 連立一次方程式の解を反復的に決定する方法 | |
| US10489481B1 (en) | Efficient matrix property determination with pipelining and parallelism | |
| Lucius et al. | An algorithm for extremal eigenvectors computation of hermitian matrices and its FPGA implementation | |
| CN117192471B (zh) | 基于hls实现二维doa估计的方法、装置及存储介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20160705 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20160802 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 5985079 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| LAPS | Cancellation because of no payment of annual fees |