(57)【要約】  【課題】  極が集中している周波数帯域の分解能を必要なレベルにまで向上させる周波数ワープ・パラメータの最適化を行った白色化フィルタを自動構成する技法を提供すること。  【解決手段】  ワープ・パラメータλをもつオールパス・フィルタを通過させた信号を白色化フィルタが処理するようにし、その後、-λのワープ・パラメータのオールパス・フィルタにより周波数軸を戻すようにシステムを構成して、最適なλを決定することで白色化フィルタを最適化する。そのためにまず、λ = 0でのAR次数pを自動推定する。ワープ・パラメータλをもつオールパス・フィルタを通過させた離散フーリエ変換スペクトルの値と、ワープ・パラメータλをもつオールパス・フィルタを通過させた離散ARスペクトルの値を用いるスペクトル距離dλを計算する。そして、スペクトル距離dλを最小にするλがワープ・パラメータとして設定され、ワープ・パラメータの最適化を行った白色化フィルタが構成される。【選択図】  図５
2014030074.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータの処理により、観測信号から白色化フィルタを構成する方法であって、
ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号でのAR次数pを自動推定するステップと、
前記推定された次数pに基づき、ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号における、DFTスペクトルと、ARスペクトルを計算するステップと、
λの値を変化させながら下記のステップ(a)から(c)までを実行するステップと、
(a) 前記λの値を用いて、オールパスフィルタリングよる周波数軸変換を行なうステップ、
(b) 前記フィルタリングした信号のDFTスペクトルと、ARスペクトルを計算するステップ、
(c) 前記λ = 0を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、今回のλの値を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、ARスペクトルの二乗和誤差によりスペクトル距離dλを計算するステップ、
前記スペクトル距離dλを最小にするλを選ぶステップと、
前記選ばれたλを用いた全域通過フィルタ出力信号からARパラメータを推定し白色化フィルタを構成するステップを有する、
白色化フィルタの構成方法。
【請求項２】
前記DFTスペクトルの成分の和を計算するステップが、前記DFTスペクトルの成分を降順にソートし、上位の所定の個数の成分の和を計算する、請求項１に記載の白色化フィルタの構成方法。
【請求項３】
前記スペクトル距離dλが、max(d0,dλ)/min(d0,dλ)（ここで、d0は、λ=0の場合のdλの値)という項を掛けられている、請求項１に記載の白色化フィルタの構成方法。
【請求項４】
前記AR次数pを自動推定するステップが、AICに基づく、請求項１に記載の白色化フィルタの構成方法。
【請求項５】
コンピュータの処理により、観測信号から白色化フィルタを構成するプログラムであって、
前記コンピュータに、
ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号でのAR次数pを自動推定するステップと、
前記推定された次数pに基づき、ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号における、DFTスペクトルと、ARスペクトルを計算するステップと、
λの値を変化させながら下記のステップ(a)から(c)までを実行するステップと、
(a) 前記λの値を用いて、オールパスフィルタリングよる周波数軸変換を行なうステップ、
(b) 前記フィルタリングした信号のDFTスペクトルと、ARスペクトルを計算するステップ、
(c) 前記λ = 0を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、今回のλの値を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、ARスペクトルの二乗和誤差によりスペクトル距離dλを計算するステップ、
前記スペクトル距離dλを最小にするλを選ぶステップと、
前記選ばれたλを用いた全域通過フィルタ出力信号からARパラメータを推定し白色化フィルタを構成するステップを実行させる、
白色化フィルタの構成プログラム。
【請求項６】
前記DFTスペクトルの成分の和を計算するステップが、前記DFTスペクトルの成分を降順にソートし、上位の所定の個数の成分の和を計算する、請求項５に記載の白色化フィルタの構成プログラム。
【請求項７】
前記スペクトル距離dλが、max(d0,dλ)/min(d0,dλ)（ここで、d0は、λ=0の場合のdλの値)という項を掛けられている、請求項５に記載の白色化フィルタの構成プログラム。
【請求項８】
前記AR次数pを自動推定するステップが、AICに基づく、請求項５に記載の白色化フィルタの構成プログラム。
【請求項９】
コンピュータに、
請求項１で選ばれたλで推定したARパラメータを用いて逆フィルタリングするステップと、
-λを用いたオールパスフィルタリングを実行して周波数軸を元に戻すステップを実行させる、
白色化フィルタ。
【請求項１０】
コンピュータの処理により、観測信号から白色化フィルタを構成するシステムであって、
ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号でのAR次数pを自動推定する手段と、
前記推定された次数pに基づき、ワープ・パラメータλ = 0を用いた全域通過フィルタ出力信号における、DFTスペクトルと、ARスペクトルを計算する手段と、
λの値を変化させながら下記のステップ(a)から(c)までを実行する手段と、
(a) 前記λの値を用いて、オールパスフィルタリングよる周波数軸変換を行なうステップ、
(b) 前記フィルタリングした信号のDFTスペクトルと、ARスペクトルを計算するステップ、
(c) 前記λ = 0を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、今回のλの値を用いた全域通過フィルタ出力信号で計算されたDFTスペクトルの成分の和と、ARスペクトルの二乗和誤差によりスペクトル距離dλを計算するステップ、
前記スペクトル距離dλを最小にするλを選ぶ手段と、
前記選ばれたλを用いた全域通過フィルタ出力信号からARパラメータを推定し白色化フィルタを構成する手段を有する、
白色化フィルタの構成システム。
【請求項１１】
前記DFTスペクトルの成分の和を計算するステップが、前記DFTスペクトルの成分を降順にソートし、上位の所定の個数の成分の和を計算する、請求項１０に記載の白色化フィルタの構成システム。
【請求項１２】
前記スペクトル距離dλが、max(d0,dλ)/min(d0,dλ)（ここで、d0は、λ=0の場合のdλの値)という項を掛けられている、請求項１０に記載の白色化フィルタの構成システム。
【請求項１３】
前記AR次数pを自動推定するステップが、AICに基づく、請求項１０に記載の白色化フィルタの構成システム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、有色ノイズを除去するための白色化フィルタの構成技法に関するものである。
【背景技術】
【０００２】
従来より、高精度の極推定を行う手法は、NMR(Nuclear Magnetic Resonance spectroscopy: 核磁気共鳴分光法)におけるDNA構造解析等で用いられている様に（例えば、村田, 久保田, “AR過程に基づく高精度な極推定法とNMR構造解析への応用,” 信学論A, Vol.J91-A, No.8, pp.721-781, Aug. 2008.などを参照）、現在も多くの応用分野が存在する。
【０００３】
近年、ナノ衛星と呼ばれる超小型人工衛星を、従来の衛星よりもはるかに低いコストで打ち上げ・運用していく計画が進行中であるが、このナノ衛星における姿勢制御信号に対する擾乱成分解析も、その応用分野の一つとなってきている。
【０００４】
ナノ衛星は、これまでに無い小型人工衛星であるため、姿勢制御が観測精度に重要な役割を持つこととなり、かつ、その制御は適応的に行う必要がある。また従来よりも飛行高度が低いため、自然環境からの外乱要因が従来と異なり、残留磁気トルクの影響が最も大きくなる等といった特徴がある（例えば、T.Inamori,N.Sako, S.Nakasuka, “Attitude Control System for the Nano-Astrometry Satellite “Nano-JASMINE”,” Inter. Journal of Aircraft Eng. And Aerospace Techn, pp.221-228, 2011.などを参照）。これは、閉ループ制御系が極めて短時間の適応系となり、従来の人工衛星知見では想定できない擾乱がプラント入力として入り込み、観測結果の擾乱要因となることを意味する。また一般に、人工衛星にかかる外部擾乱は低周波数帯域の信号であることが多いことも問題を難しくしている。
【０００５】
従来は、例えば特開２０００−２２４０７７号公報や特開２００８−１９９６０２号公報に記載されているように、こうした擾乱要因を全体で白色ガウス性信号と仮定して制御・観測を行ってきた。しかし、実際には上記のとおり、複数の未知周波数信号（マルチトーン信号）がノイズに埋もれた低SN比信号とも考えられるものである。そこで、低SN比マルチトーン信号をモデル化し逆フィルタリングによって白色化、あるいは無視できるレベルの擾乱に変換することができれば制御・観測精度の改善が可能と考えられるため、そうしたモデル化のための効果的手法が必要とされている。
【０００６】
そのためには、AR(Auto-Regressive: 自己回帰)モデルとして観測信号系を同定し、その逆フィルタをかけることが有効と考えられる。非特許文献１のような高精度推定の場合、擾乱には弱いことが多く、一方、正弦波の個数推定であれば負のSN比であっても固有値分布などに基づき推定する手法が提案されているが（例えば、辻、大森、佐野, “ARスペクトル推定における正弦波信号の個数判定法,” 信学論A, Vol.J74-A, No.9, pp.1374-1384, Sept. 1991.を参照）、スペクトル包絡を構成できるパラメータ推定とは異なる。
【０００７】
そこで、上記のナノ衛星の制御で必要となるのは、極スペクトル包絡を崩すことなく、逆フィルタが構成できるレベルの低SN比、すなわち、数dB程度の観測信号を基に、極分布に偏りがある推定を自動的に高精度に行うことが課題となっていた。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２０００−２２４０７７号公報
【特許文献２】特開２００８−１９９６０２号公報
【非特許文献】
【０００９】
【非特許文献１】村田, 久保田, “AR過程に基づく高精度な極推定法とNMR構造解析への応用,” 信学論A, Vol.J91-A, No.8, pp.721-781, Aug. 2008.
【非特許文献２】T.Inamori,N.Sako, S.Nakasuka, “Attitude Control System for the Nano-Astrometry Satellite “Nano-JASMINE”", Inter. Journal of Aircraft Eng. And Aerospace Techn, pp.221-228, 2011.
【非特許文献３】辻、大森、佐野, “ARスペクトル推定における正弦波信号の個数判定法,” 信学論A, Vol.J74-A, No.9, pp.1374-1384, Sept. 1991.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従って、この発明の目的は、極が集中している周波数帯域の分解能を必要なレベルにまで向上させる周波数ワープ・パラメータの最適化を行った白色化フィルタを自動構成する技法を提供することにある。
【００１１】
この発明の他の目的は、低SN比マルチトーン信号において、擾乱成分の影響を少なくするために、スペクトルピーク存在帯域の分解能向上と、原信号（低SN比マルチトーン信号）のフーリエパワースペクトルを許容範囲内で保存させることを両立させつつ、ARシステム同定を行うことで、従来の次数最適化基準に基づく同定よりも誤差が少ない（白色化）逆フィルタを構成するための技法を提供することにある。
【課題を解決するための手段】
【００１２】
この発明は、ワープ・パラメータλをもつオールパス・フィルタを通過させた信号を白色化フィルタが処理するようにし、その後、ワープ・パラメータ  -λのオールパス・フィルタにより周波数軸を戻すようにシステムを構成して、最適なλを決定することで白色化フィルタを最適化するものである。
【００１３】
このためまず、この発明のシステムは、λ = 0でのAR次数pを自動推定する。ここで次数推定は低SN比マルチトーン信号(原信号)に対し、好適には、AIC(Akaike Information Criterion)によって自動推定するが、適宜他の次数推定基準を使用してもよい。
【００１４】
この発明のシステムは、ワープ・パラメータλをもつオールパス・フィルタを通過させた離散フーリエ変換スペクトルの値と、ワープ・パラメータλをもつオールパス・フィルタを通過させた離散ARスペクトルの値を用いるスペクトル距離dλを計算する。
【００１５】
ワープ・パラメータλは、各λ近辺で局所最適値を持つので、この発明のシステムは、ワープさせない原信号でのスペクトル距離d0以下となることで非反証されたものを最初に抽出する。最初のλ候補のセットは、|λ|<1の範囲内を等間隔にセットしても良いし、適宜集中的な候補をセットしても良い。このとき、近傍探索法は直接探索以外の既存のものであっても構わない。
【００１６】
次に、この発明のシステムは、αトリムド・スペクトル距離dλを計算することで、低SN比マルチトーン信号のAR分析分解能を自動的に決定し、最適なワープ・パラメータλを決定する。ここで、αトリムド（処理）とは、離散スペクトルが規格化周波数0から2πの2N-pointで計算されているとすると、0からπまでのN-pointのスペクトル・ビンを大きさでソートし、連続したα個のビンを選択する非線形フィルタリングを意味する。
【００１７】
この発明のシステムはまず、離散フーリエ変換(DFT)スペクトルと離散ARスペクトルについてαトリムドMSE(平均二乗誤差)を計算する。上記スペクトル距離のみでは、ワープさせたDFTスペクトルが原信号のものから著しく乖離している場合でもαトリムドMSE距離は小さくなる場合が存在するので（逆方向のワープで極をつぶしてしまう場合など）、ここで、好適にはペナルティ効果として原信号のDFTスペクトルのピークパワーとのαトリムド比を乗じる。
【００１８】
候補となったλ周辺の局所最適値を直接探索する場合の距離基準も、このペナルティ効果を伴うαトリムドMSE距離を用いる。
【００１９】
この最適ワープ・パラメータ決定後は、この発明のシステムは、そのλによるオールパスフィルタリングを行い、ARパラメータ推定を行った後、その推定パラメータによる逆フィルタリングを行い、−λによるオールパスフィルタリングによってワープした分解能を元に戻す。
【発明の効果】
【００２０】
この発明によれば、低SN比マルチトーン信号に対して、小さい誤差で動作し得る白色化フィルタが自動構成される。
【図面の簡単な説明】
【００２１】
【図１】人工衛星の制御システムの概要ブロック図である。
【図２】本発明を実施するためのコンピュータ・ハードウェア構成のブロック図である。
【図３】ワープ・パラメータλを決定する処理の機能ブロック図である。
【図４】ワープ・パラメータλを決定する処理のフローチャートを示す図である。
【図５】白色化フィルタの処理を示す図である。
【発明を実施するための形態】
【００２２】
以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。
【００２３】
図１は、人工衛星の制御システムの全体の概要ブロック図を示す図である。図１において、人工衛星１０２は、アクチュエータ１０４によって駆動制御される。アクチュエータ１０４は、姿勢制御に関するシステムと、推進に関するシステム等に関与する。
【００２４】
人工衛星１０２は、地球重力場の歪み、月や太陽からの引力、太陽風、大気中の空気分子やイオンなど、様々な擾乱要因を受けながら軌道上を飛ぶ。このような擾乱要因は、図１では、環境雑音１０６として示されている。
【００２５】
このような環境雑音１０６によって軌道や姿勢が乱されるので、人工衛星１０２は、センサにより姿勢制御情報などの情報を取得し、これらの観測信号を、アクチュエータ１０４にフィードバックすることにより、サーボ・システムを構成する。
【００２６】
その際、人工衛星１０２のセンサの観測信号には、環境雑音１０６に起因するノイズが含まれているので、観測信号は、環境雑音キャンセラ１０８を介して人工衛星１０２にフィードバックされる。環境雑音キャンセラ１０８は、観測信号中のノイズを除去するために、例えばカルマン・フィルタを使用する。同様に、観測信号を入力して、姿勢制御のためアクチュエータ１０４に制御信号を入力する観測姿勢制御器１１０も、観測信号中のノイズを除去するために、例えばカルマン・フィルタを使用する。
【００２７】
従来は、環境雑音キャンセラ１０８や観測姿勢制御器１１０で使用されるカルマン・フィルタは、観測信号に含まれるノイズが白色雑音であると想定して設計されており、比較的大型の人工衛星の場合はそれであまり問題にならなかった。
【００２８】
しかし、小型人工衛星の場合は、既述のように、飛行高度が低いので外乱要因が大型の人工衛星と異なったり、閉ループ制御系が極めて短時間の応答となる、などの相違点があり、観測信号に含まれるノイズが白色雑音であると想定すると、適切な制御が困難になることがある。
【００２９】
そこで、本願発明は、このような低SN比の観測信号の雑音成分を白色化するフィルタ１１２を提供するものである。
【００３０】
図２は、上記システムを実施するためのコンピュータ・ハードウェア構成の一例の概要ブロック図である。
【００３１】
図２において、システム・バス２０２には、ＣＰＵ２０４と、主記憶（ＲＡＭ）２０６と、フラッシュ・メモリ２０８が接続されている。システム・バス２０２にはさらに、ＰＣＩ、ＵＳＢなどの所定のインターフェース・カード２１０を介して、プラント２１２及び各種センサ（図示しない）が接続されている。
【００３２】
ＣＰＵ２０４は、Intel(R) 386EX、PowerPC、ARM(R)アーキテクチャのプロセッサなどの組み込みシステムに適合する任意のプロセッサを用いることができる。
【００３３】
フラッシュ・メモリ２０８には、組み込み用Linux(R)などのオペレーティング・システム、図１に示す環境雑音キャンセラ１０８、観測姿勢制御器１１０などの制御プログラム、本発明に係る白色化フィルタ及びその構成モジュール１１２、その他全体の制御モジュールなどが保存されており、これらのモジュールは、システムの起動時、あるいは必要時に主記憶２０６にロードされて実行される。
【００３４】
プラント２１２は、人工衛星本体１０２、アクチュエータ１０４、その他のセンサの総称であり、インターフェース２１０を介して、環境雑音キャンセラ１０８、観測姿勢制御器１１０などからの制御信号を入力して、駆動される。一方、ＣＰＵ２０４によって動作される制御プログラムは、インターフェース２１０を介して、プラント２１２から観測信号を受け取る。
【００３５】
次に、図３は、本発明に係る白色化フィルタ及びその構成モジュール１１２の機能構成を詳細に示すブロック図である。
【００３６】
図３において、メイン・プログラム３０２は、白色化フィルタ及びその構成モジュール１１２の全体を機能を制御する。パラメータ３０６は、フラッシュ・メモリ２０８、または図示しないが所定のＲＯＭなどに格納され、計算に使用されるパラメータを含む。
【００３７】
メイン・プログラム３０２は、信号入力ルーチン３０４、AR次数自動推定ルーチン３０８、DFTスペクトル計算ルーチン３１０、ARスペクトル計算ルーチン３１２、スペクトル距離dλ計算ルーチン３１４、ワープ・パラメータλ探索ルーチン３１６、オールパス・フィルタ３１８、白色化フィルタ３２０、及び信号出力ルーチン３２２を適宜呼び出して処理を行う。
【００３８】
これらのプログラムあるいはルーチンは、C、C++、アセンブラ、Java(R)などの任意の既知のプログラミング言語により作成することができる。
【００３９】
ＡＲ次数自動推定ルーチン３０８は、Akaike Information Criterion (AIC)によりAR次数を自動推定する。それは具体的には、次の計算による。すなわち、p次のARモデルを次のように記述する。
【数１】
2014030074.tif 000003
ここで、aiは、ARモデルの係数であり、enは、誤差項である。これに対する対数尤度は、次のように記述される。
【数２】
2014030074.tif 000004
ここで、θはARパラメータ・ベクトル、Nは観測データの数、
【数３】
2014030074.tif 000005
は、ARモデル予測誤差の二乗平均をあらわす。この式は、次の場合に最大値をとる:
【数４】
2014030074.tif 000006
その場合、式は次のようになる。
【数５】
2014030074.tif 000007
定数部は次数推定のAIC比較に無関係となるので、下記のように記述することができる。
【数６】
2014030074.tif 000008
この式を最小にするpを求めることによって、AR次数pが自動推定される。
【００４０】
なお、AR次数推定の方法として、AIC以外に、最終予測誤差規範(Final Prediction Error: FPE)を使用することができる。あるいは、観測信号の共分散行列を求め、その固有値を降べきの順に並べたときに、固有値の変化がほぼなくなる順番を推定次数とする、という手法も使用することができる。
【００４１】
さて、観測信号xnに対して、オールパス・フィルタ３１８を通した信号をxn,λとする。
すなわち、xn,λ = A(z-1,λ)xn
ここで、
【数７】
2014030074.tif 000009
【００４２】
すると、DFT(離散フーリエ変換)計算ルーチン３１０は、Nポイント離散フーリエ変換によるパワースペクトルを次式で計算する。
【数８】
2014030074.tif 000010
【００４３】
同様にして、次式のp次(p≦N) ARモデルが推定されたとき、
【数９】
2014030074.tif 000011
オールパス・フィルタ３１８を通した信号も同様にARモデルで以下のように表現することができる。
【数１０】
2014030074.tif 000012
【００４４】
すると、AR計算ルーチン３１２は、Nポイント・パワースペクトルを次式で計算する。
【数１１】
2014030074.tif 000013
ここで、
【数１２】
2014030074.tif 000014
は、予測誤差の二乗平均をあらわす。
【００４５】
dλ計算ルーチン３１４は、DFTスペクトル計算ルーチン３１０及びARスペクトル計算ルーチン３１２の計算結果を利用して、下記の式で、スペクトル距離dλを計算する。
【数１３】
2014030074.tif 000015
ここで、
【数１４】
2014030074.tif 000016
また、スペクトル距離dλの式の右辺の最初の、max()/min()の項と、次の1/αの項は、ペナルティ項として、逆方向のワープで極をつぶしてしまう場合を回避する働きをする。
【００４６】
【数１５】
2014030074.tif 000017
は、A(z-1,λ)のオールパス・フィルタ３１８を通した信号について、DFTスペクトル計算ルーチン３１０で計算した2N点DFTスペクトルを規格化周波数0からπの間に相当するN個について大きさの順にソートしたi番目の成分をあらわす。
【００４７】
【数１６】
2014030074.tif 000018
は、A(z-1,λ)のオールパス・フィルタ３１８を通した信号にp次システム同定を行い、ARスペクトル計算ルーチン３１２で計算した、AR係数の2N点ARスペクトルを規格化周波数0からπの間に相当するN個について大きさの順にソートしたi番目の成分をあらわす。
【００４８】
また、和の下位項をトリムするためのαの値は厳密ではないが、離散フーリエ変換のポイント数の２割程度にセットされる。例えば、離散フーリエ変換のポイント数Nは、好適には128、256、512のような２のべき乗の値が選ばれ、ポイント数が256であるなら、αの値は50程度にセットされる。αの値や、離散フーリエ変換のポイント数Nなどの値は、予めパラメータ３０６として保存されている。
【００４９】
図３に戻って、λ探索ルーチン３１６は、与えられたλの与えられた近傍の範囲で、一番スペクトル距離dλが小さくなるλの値を探索する機能をもつ。
【００５０】
白色化フィルタ３２０は、与えられたλで推定したARパラメータを用いて、逆フィルタリング、すなわち白色化を行う機能をもつ。白色化フィルタ３２０の処理は、図５を参照して、後で詳しく説明する。
【００５１】
信号出力ルーチン３２２は、白色化フィルタ３２０によって逆フィルタされ、-λのオールパス・フィルタリングされた信号を、図１に示す観測姿勢制御器１１０、あるいは環境雑音キャンセラ１０８などに提供する機能をもつ。
【００５２】
次に図４のフローチャートを参照して、最適ワープ・パラメータλを決定する処理について説明する。この処理は基本的に、メイン・プログラム３０２が、図３に示す様々なルーチンを適宜呼び出すことにより実行される。
【００５３】
図４において、メイン・プログラム３０２は、ステップ４０２で、信号入力ルーチン３０４により所定個数の観測信号を読み取ることにより、AR次数自動推定ルーチン３０６を呼び出して、λ=0でのAR次数pを自動推定する。
【００５４】
次にメイン・プログラム３０２は、ステップ４０４で、DFTスペクトル計算ルーチン３１０及びARスペクトル計算ルーチン３１２を呼び出して、観測信号のDFTスペクトルと、ARスペクトルを計算する。
【００５５】
次にメイン・プログラム３０２は、ステップ４０６で、スペクトル距離dλ計算ルーチン３１４を呼び出し、λ = 0のときのスペクトル距離dλの値である、d0を計算する。
【００５６】
次にメイン・プログラム３０２は、ステップ４０８で、λの初期候補をセットする。ここでのλの初期候補とは例えば、λを-1.0から1.0まで等分に分割して、-1.0, -0.9, -0.8, -0.7, ..., 0.7, 0.8,0.9, 1.0のような実数の組である。この分割数は細かいほど最適なλが計算されるが、細かすぎると計算時間が増加するので、適当な分割数が選ばれる。
【００５７】
次にメイン・プログラム３０２は、ステップ４１０で、λの初期候補から値を選ぶ。例えば最初は、λ = -1.0を選ぶ。
【００５８】
次にメイン・プログラム３０２は、ステップ４１２で、対象のλを用いて、オールパス・フィルタ３１８を呼び出し、観測信号のオールパス・フィルタリングによる周波数軸変換を行う。
【００５９】
次にメイン・プログラム３０２は、ステップ４１４で、DFTスペクトル計算ルーチン３１０及びARスペクトル計算ルーチン３１２を呼び出して、フィルタリングした観測信号のDFTスペクトルと、ARスペクトルを計算する。
【００６０】
次にメイン・プログラム３０２は、ステップ４１６で、dλ計算ルーチン３１４を呼び出し、スペクトル距離dλを計算する。
【００６１】
次にメイン・プログラム３０２は、ステップ４１８で、初期セットのλの候補が最後、すなわち上記の例では1.0に達したかどうかを判断し、そうでないなら、ステップ４１０に戻り、候補から次のλの値を選ぶ。ここで、次のλの値とは、前回のλの値が-0.9であったなら、-0.8を選ぶような具合である。
【００６２】
こうして初期セットの最後のλの値まで達すると、メイン・プログラム３０２は、ステップ４２０に進み、初期セットの最後のλのうち、初期候補の各々で計算したスペクトル距離dλが、d0以下になるようなλを選択候補として選びセットする。
【００６３】
次にメイン・プログラム３０２は、ステップ４２２で、候補からλを選び、λ探索ルーチン３１６を呼び出して、δをある小さい実数として、λの近傍[λ-δ,λ+δ]の範囲で、λを探索する。
【００６４】
次にメイン・プログラム３０２は、ステップ４２４で、探索されたλを用いて、オールパス・フィルタ３１８を呼び出し、観測信号のオールパス・フィルタリングによる周波数軸変換を行う。
【００６５】
次にメイン・プログラム３０２は、ステップ４２６で、DFTスペクトル計算ルーチン３１０及びARスペクトル計算ルーチン３１２を呼び出して、フィルタリングした観測信号のDFTスペクトルと、ARスペクトルを計算する。
【００６６】
次にメイン・プログラム３０２は、ステップ４２８で、dλ計算ルーチン３１４を呼び出し、スペクトル距離dλを計算する。
【００６７】
次にメイン・プログラム３０２は、ステップ４３０で、選択したλの直接探索が収束したかどうか判断する。ここで収束したとは例えば、探索したλについて、スペクトル距離dλが所定の閾値以下になった、あるいは選択したλのスペクトル距離dλに対して、探索したλのスペクトル距離dλが所定の値より減ったなどの基準を意味する。
【００６８】
選択したλの直接探索が収束していないと判断すると、処理はステップ４２４に戻り、周辺探索を続けて、別の探索したλについて、ステップ４２４の計算を続ける。ただ、ステップ４３０で所定時間あるいは所定回数の探索で収束しない場合は、そのλの直接探索は打ち切るようにする。
【００６９】
ステップ４３０で、選択したλの直接探索が収束したと判断すると、メイン・プログラム３０２は、ステップ４３２に戻り、選択セットのλの候補が最後かどうか判断する。そしてまだ最後でないならステップ４２２に戻って別の選択セットのλを選び、選択セットのλの候補が最後なら、ステップ４３４で、各λ近辺の最適値の中のうちで、最小スペクトル距離dλを生成するλを最適ワープ・パラメータと決定して処理を終わる。
【００７０】
次に図５を参照して、最適ワープ・パラメータλを用いた白色化フィルタ３２０の動作について説明する。すなわち、図５の参照番号５０２に示すように、図４のフローチャートで示す処理の結果、最適ワープ・パラメータλが決定されると、これは、最適なλで推定したARパラメータを用いて逆フィルタリング(白色化)を行うステップ５０４を実行する白色化フィルタ３２０と、周波数軸を元に戻すため-λを用いたオールパス・フィルタリングを行うステップ５０６を実行するオールパス・フィルタ３１８によって利用される。
【００７１】
より具体的には、観測信号をxnとする。すると、白色化フィルタ３２０は、周波数軸をワープさせた信号を次式により生成する。
【数１７】
2014030074.tif 000019
これをxnwと置くと、白色化フィルタ３２０は、このxnwにARモデルを仮定して、次式となるARパラメータを推定する。
【数１８】
2014030074.tif 000020
【００７２】
この伝達関数表現は次のようになる。
【数１９】
2014030074.tif 000021
よって逆フィルタリング(白色化)とは、以下の操作を意味する。
【数２０】
2014030074.tif 000022
すなわち、enwだけが出力されるようにフィルタリングがなされる。enwが白色雑音であるというのがARモデルの定義なので、白色化がなされたことになる。以上がステップ５０４で実行される処理である。
【００７３】
ステップ５０６では、下記のように、周波数軸を元に戻すオールパス・フィルタリングが行われる。
【数２１】
2014030074.tif 000023
【００７４】
以上、人工衛星の観測信号の白色化の実施例に基づき本発明を説明してきたが、本発明はこの実施例に限定されず、任意の分野の白色化フィルタに適用可能である。特に低SN比の状況で、本発明の白色化フィルタは、有利な効果を奏する。
【符号の説明】
【００７５】
１１２  雑音成分の白色化ブロック
２０４  ＣＰＵ
２０６  ＲＡＭ
３０８  AR次数自動推定ルーチン
３１０  DFTスペクトル計算ルーチン
３１２  ARスペクトル計算ルーチン
３１４  dλ計算ルーチン
３１６  オールパス・フィルタ
３１８  白色化フィルタ
【図１】
2014030074.tif 000024
【図２】
2014030074.tif 000025
【図３】
2014030074.tif 000026
【図４】
2014030074.tif 000027
【図５】
2014030074.tif 000028