JP2017092577A - 無線チャネル解析装置 - Google Patents

Abstract

【課題】ノイズが多い環境でもチャネル統計量を高精度に推定する。
【解決手段】無線チャネル解析装置は、無線通信路を介して受信された複数のパイロット信号を含む観測データを時間軸上の各スロットについて取得する観測データ取得部２２と、観測データ取得部２２が取得した観測データを適用した尤度の算定を含む粒子フィルタにより、無線通信路のチャネル統計量およびチャネル歪の確率分布をスロット毎に推定する解析処理部２４とを具備し、解析処理部２４は、解析対象のスロットの後方のスロットの観測データを適用したカルマン平滑化により解析対象のスロットにおける尤度を算定する。
【選択図】図１

Description

本発明は、無線通信路（チャネル）を解析する技術に関する。

時系列に観測される観測データを解析する技術が従来から提案されている。例えば非特許文献１には、粒子フィルタとカルマンフィルタとを併用することで、システムの状態とパラメータとを同時に推定する技術が開示されている。

Maria Paula Rios and Hedibert Freitas Lopes, "The Extended Liu and West Filter: Parameter Learning in Markov Switching Stochastic Volatility Models," State-Space Models, p. 23-61, Springer, 2013

ところで、無線通信路を利用した無線通信では、障害物等による反射や散乱に起因した相異なる複数の経路の電波が相互に干渉することで信号歪（以下「チャネル歪」という）が発生し得る。チャネル歪（フェージング）を補償するには、携帯電話機等の端末装置の移動速度や遅延スプレッド等のチャネル統計量を推定することが必要である。前述の非特許文献１の技術を利用すれば、粒子フィルタとカルマンフィルタとの併用によりチャネル歪とチャネル統計量とを同時に推定することが可能である。

しかし、非特許文献１の技術のように、粒子フィルタに適用される尤度をカルマンフィルタのみで算定する技術では、無線通信路にノイズが多い環境のもとで過剰適合（overfitting）が発生し易く、結果的にチャネル統計量を高精度に推定できないという問題がある。以上の事情を考慮して、本発明は、ノイズが多い環境でもチャネル統計量を高精度に推定することを目的とする。

以上の課題を解決するために、本発明の好適な態様に係る無線チャネル解析装置は、無線通信路を介して受信された複数のパイロット信号を含む観測データを時間軸上の各スロットについて取得する観測データ取得部と、前記観測データ取得部が取得した観測データを適用した尤度の算定を含む粒子フィルタにより、前記無線通信路のチャネル統計量およびチャネル歪の確率分布をスロット毎に推定する解析処理部とを具備し、前記解析処理部は、解析対象のスロットの後方のスロットの観測データを適用したカルマン平滑化により前記解析対象のスロットにおける前記尤度を算定する。

本発明によれば、ノイズが多い環境でもチャネル統計量を高精度に推定することが可能である。

本発明の第１実施形態に係る通信システムの構成図である。 パイロット信号の説明図である。 粒子フィルタの説明図である。 解析処理のフローチャートである。 尤度算定処理のフローチャートである。 尤度算定処理のフローチャートである。 ノイズが少ない環境で対比例により移動速度を推定した結果のグラフである。 ノイズが多い環境で対比例により移動速度を推定した結果のグラフである。 ノイズが多い環境で第１実施形態により移動速度を推定した結果のグラフである。 第１実施形態の効果の説明図である。 第２実施形態における初期値設定処理のフローチャートである。 第２実施形態における第１初期化処理の説明図である。 第２実施形態における第２初期化処理の説明図である。

＜第１実施形態＞
図１は、第１実施形態の通信システム１００を例示する構成図である。第１実施形態の通信システム１００は、音声通話やデータ通信等の移動体通信サービスを提供するシステムであり、相互に無線通信が可能な無線基地局１２と端末装置１４とを具備する。実際には複数の無線基地局１２と複数の端末装置１４とが通信システム１００に包含されるが、図１では１個の無線基地局１２と１個の端末装置１４とを便宜的に図示した。

各端末装置１４は、例えば携帯電話機やスマートフォン等の可搬型の通信装置（ユーザ装置：User Equipment）である。タブレット端末，パーソナルコンピュータ，ＵＭＰＣ（Ultra-Mobile Personal Computer），携帯型ゲーム装置等の通信端末も端末装置１４として利用され得る。無線基地局１２は、移動体通信網を構成する通信装置であり、端末装置１４との間で無線通信が可能である。例えば、３ＧＰＰ（Third Generation Partnership Project）のＬＴＥ（Long Term Evolution）に準拠した通信方式で各端末装置１４と無線通信するｅＮＢ（evolved Node B）が無線基地局１２として好適に利用される。具体的には、下りリンクの無線送信方式としてＯＦＤＭＡ（Orthogonal Frequency Division Multiple Access）が採用され、上りリンクの無線送信方式としてＳＣ-ＦＤＭＡ（Single-Carrier Frequency Division Multiple Access）が採用される。ただし、無線基地局１２と端末装置１４との間の通信方式は以上の例示に限定されない。

図１に例示される通り、第１実施形態の端末装置１４は、無線通信部２０と観測データ取得部２２と解析処理部２４とチャネル補償部２６とを具備する。端末装置１４の各要素は、例えば記憶装置に記憶されたプログラムを演算処理装置（ＣＰＵ）が実行することで実現される。ただし、例えばＦＰＧＡ（Field Programmable Gate Array）やＤＳＰ（Digital Signal Processor）等の電子回路で端末装置１４の一部または全部の機能を実現することも可能である。

無線通信部２０は、無線通信路（無線チャネル）を介して無線基地局１２と通信する要素であり、例えば電波を送受信するアンテナを含んで構成される。具体的には、無線通信部２０は、無線基地局１２から送信された下りリンク信号を受信するとともに、上りリンク信号を生成して無線基地局１２に送信する。

図２は、無線基地局１２が端末装置１４に送信する下りリンク信号（ＯＦＤＭ信号）の説明図である。図２に例示される通り、第１実施形態の下りリンク信号は、周波数軸上の相異なるサブキャリアと時間軸上の相異なるサブフレームとに対応する複数のリソースエレメントを包含する。無線基地局１２と端末装置１４との間の無線通信路の状態の推定のために参照されるパイロット信号（具体的には、ＣＲＳ：Cell-specific Reference Signal）が、時間-周波数領域内の複数のリソースエレメントに分散して配置される。

観測データ取得部２２は、無線通信部２０が受信した下りリンク信号から時間軸上のスロット毎に観測データｙ(t)を取得する。記号ｔは、観測データ取得部２２が観測データｙ(t)を取得する任意の１個のスロットを意味する。観測データｙ(t)は、無線通信路を介して受信された複数のパイロット信号を配列したベクトルである。例えば、１個のサブフレーム（例えば図２の０,７,……）に配置された複数のパイロット信号の配列が観測データｙ(t)として生成される。すなわち、任意の１個のスロットは、時間軸上の複数のサブフレームのうち観測データ取得部２２が観測データｙ(t)を取得したサブフレームに相当する。

無線基地局１２と端末装置１４との間の無線通信路では、建造物等の障害物による反射や散乱に起因した相異なる複数の経路の電波が相互に干渉することでチャネル歪（フェージング）が発生し得る。図１の解析処理部２４は、観測データ取得部２２が取得した観測データｙ(t)を利用して、無線通信路のチャネルの状態であるチャネル歪と当該チャネルに関する統計量であるチャネル統計量とをスロット毎に推定する。具体的には、第１実施形態の解析処理部２４は、観測データｙ(t)を適用した粒子フィルタ（Particle filter)によりチャネル歪とチャネル統計量とを同時かつ実時間的に推定する。チャネル統計量は、例えば、端末装置１４の移動速度（すなわちドップラー周波数）または遅延スプレッドである。第１実施形態では、１種類のチャネル統計量（具体的には移動速度）の推定を便宜的に想定する。以上の説明から理解される通り、第１実施形態の端末装置１４は、下りリンクの無線通信路のチャネル歪およびチャネル統計量を解析する装置（無線チャネル解析装置）として機能する。

図１のチャネル補償部２６は、解析処理部２４による推定結果に応じたチャネル補償を実行する。具体的には、チャネル補償部２６は、制御チャネルおよびデータチャネルの振幅および位相を解析処理部２４による解析結果に応じてサブキャリア毎に調整する。解析処理部２４による解析結果を利用したチャネル補償には、公知の技術が任意に採用され得る。

解析処理部２４による無線通信路の解析について以下に詳述する。第１実施形態では無線通信路の解析に状態空間モデルを利用する。一般的な状態空間モデルは、以下の数式(1)の観測方程式と数式(2)の状態方程式とで表現される。

数式(1)の記号θ(t)は、無線通信路のチャネル歪に相当する状態ベクトルを意味する。また、数式(1)の記号Ｆ(t)は各サブキャリアのパイロット信号の観測行列を意味し、数式(2)の記号Ｇ(t)は各サブキャリアのパイロット信号の遷移行列を意味する。数式(1)の記号ｖ(t)および数式(2)の記号ｗ(t)は、複素正規分布に従うノイズベクトルである。

第１実施形態では、無線通信路のチャネル歪およびチャネル統計量の同時推定に拡張状態（augmented state）を利用する。拡張状態は、線形・ガウス型の部分θ(t)と非線形・非ガウス型の部分ψ(t)とを包含する。線形・ガウス型の部分θ(t)は無線通信路のチャネル歪に相当し、非線形・非ガウス型の部分ψ(t)は、超パラメータとしてのチャネル統計量に相当する。

第１実施形態の解析処理部２４は、無線通信路のチャネル統計量ψ(t)およびチャネル歪θ(t)を粒子フィルタにより推定する。具体的には、粒子フィルタは、チャネル統計量ψ(t)およびチャネル歪θ(t)の確率分布を複数（Ｎ個）の粒子の配列で近似的に表現する。粒子フィルタに利用される粒子の個数Ｎは任意（例えば１０個〜１００個）である。図３に模式的に例示される通り、Ｎ個の粒子の各々は、実現値（Ralization）と加重値（Weight）との組合せで指定される。

第１実施形態の解析処理部２４は、以下の数式(3)で表現される第(t-1)番目のスロットのＮ個の粒子から、数式(4)で表現される第ｔ番目のスロットのＮ個の粒子を生成する。

数式(4)で表現される通り、第ｔ番目のスロットにおける第ｎ番目（ｎ＝１〜Ｎ）の１個の粒子は、実現値（ψ(t)⁽ⁿ⁾，ｍ(t)⁽ⁿ⁾，Ｃ(t)⁽ⁿ⁾）と加重値ω(t)⁽ⁿ⁾とで表現される。具体的には、第ｎ番目の粒子の実現値は、無線通信路のチャネル統計量ψ(t)⁽ⁿ⁾と、チャネル歪θ(t)の確率分布（例えば複素正規分布）を規定する平均（平均ベクトル）ｍ(t)⁽ⁿ⁾および共分散（共分散行列）Ｃ(t)⁽ⁿ⁾とを包含する。第１実施形態の解析処理部２４は、解析対象となる第ｔ番目のスロットの各粒子の実現値（ψ(t)⁽ⁿ⁾，ｍ(t)⁽ⁿ⁾，Ｃ(t)⁽ⁿ⁾）および加重値ω(t)⁽ⁿ⁾を、直前（第(t-1)番目）のスロットの各粒子の実現値（ψ(t-1)⁽ⁿ⁾，ｍ(t-1)⁽ⁿ⁾，Ｃ(t-1)⁽ⁿ⁾）および加重値ω(t-1)⁽ⁿ⁾を利用した処理（以下「解析処理」という）で算定する。

図４は、任意の第ｔ番目のスロットを解析対象として解析処理部２４が実行する解析処理のフローチャートである。時間軸上の複数のスロットの各々が順番に解析対象として選択されてスロット毎に図４の解析処理が実行される。

解析処理を開始すると、解析処理部２４は、超パラメータとしてのチャネル統計量ψ(t-1)⁽ⁿ⁾の粒子毎の移動平均μ⁽ⁿ⁾と共分散Γとを算定する（ＱA1）。具体的には、解析処理部２４は、直前のスロットまでの過去のチャネル統計量ψ(t-1)⁽ⁿ⁾の移動平均μ⁽ⁿ⁾を以下の数式(5)の演算で粒子毎に算定するとともに、Ｎ個の粒子にわたるチャネル統計量ψ(t-1)⁽ⁿ⁾の共分散Γを以下の数式(6)の演算で算定する。

数式(5)の記号αは忘却係数であり、１以下の正数に設定される。数式(5)の記号Ｅ_ω(t-1) ⁽ⁿ⁾［ψ(t-1)⁽ⁿ⁾］は、加重値ω(t-1)⁽ⁿ⁾を適用したチャネル統計量ψ(t-1)⁽ⁿ⁾の加重平均を意味する。また、数式(4)の記号Ｖar_ω(t-1) ⁽ⁿ⁾［ψ(t-1)⁽ⁿ⁾］は、加重値ω(t-1)⁽ⁿ⁾を適用したチャネル統計量ψ(t-1)⁽ⁿ⁾の共分散を意味する。

次に、解析処理部２４は、第(t-1)番目のスロットにおけるＮ個の粒子から補助インデックス（auxiliary index）ｋ_nを選別する（ＱA2）。第１実施形態の例示のように補助インデックスｋ_nの選別により粒子の縮退を抑制する粒子フィルタは、特に補助粒子フィルタ（auxiliary particle filter）として、例えばMichael K Pitt and Neil Shephard, "Filtering Via Simulation: Auxiliary Particle Filters," Journal of the American Statistical Association (American Statistical Association) 94 (446): 590-591に開示されている。

第(t-1)番目のスロットにおける第ｎ番目の粒子が補助インデックスｋ_nとして選別される確率Ｐ(k_n=n)は、以下の数式(7)で表現される。

数式(5)の記号(1:t-1)は、第１番目から第(t-1)番目のスロットの時系列を意味する。また、記号~は推定値を意味する。例えば、数式(7)の推定値ψ~(t)⁽ⁿ⁾は、チャネル統計量ψ(t)⁽ⁿ⁾の推定値である。第１実施形態では、数式(5)で算定された移動平均μ⁽ⁿ⁾を推定値ψ~(t)⁽ⁿ⁾として利用する（ψ~(t)⁽ⁿ⁾＝μ⁽ⁿ⁾）。

数式(7)の記号ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)は、第１番目から第(t-1)番目までのスロットの観測データｙ(1:t-1)とチャネル統計量ψ(t)⁽ⁿ⁾の推定値ψ~(t)⁽ⁿ⁾とが観測されたときの第ｔ番目のスロットの観測データｙ(t)の尤度（第１尤度の具体例）を意味する。数式(7)から理解される通り、確率Ｐ(k_n=n)は、加重値ω(t-1)(n)と尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)との積に相当する確率分布に従う。第１実施形態のステップＱA2は、尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する処理（以下「尤度算定処理」という）ＱBを包含する。

図５は、尤度算定処理ＱBのフローチャートである。尤度算定処理ＱBを開始すると、解析処理部２４は、解析対象である第ｔ番目のスロットが特定のスロット（以下「特定スロット」という）に該当するか否かを判定する（ＱB1）。第１実施形態の解析処理部２４は、スロットのＭ個分を周期として複数のスロットの時系列から特定スロットを順次に選択する。具体的には、解析処理部２４は、解析対象のスロットの番号ｔを定数Ｍで除算した剰余（t mod M）が０である場合に、当該スロットは特定スロットに該当すると判定する一方、剰余が０でない場合には、解析対象のスロットは特定スロットに該当しないと判定する。定数Ｍは、所定の正数（例えば５）に設定される。

解析対象のスロットが特定スロットに該当しない場合（ＱB1：NO）、解析処理部２４は、第ｔ番目のスロットまでの観測データｙ(t)を適用したカルマンフィルタ（Kalman filter）により第ｔ番目のスロットの尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する（ＱB2）。すなわち、尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)は以下の数式(8)で表現される。

数式(8)の記号hypは超パラメータ（hyperparameter）を意味し、記号priは事前分布（prior distribution）を意味し、記号datは観測データ（data）を意味する。すなわち、解析処理部２４は、解析対象のスロットにおけるチャネル統計量ψ(t)⁽ⁿ⁾の推定値ψ~(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおけるチャネル歪θ(t-1)の確率分布（ｍ(t-1)⁽ⁿ⁾，Ｃ(t-1)⁽ⁿ⁾）と、解析対象のスロットの観測データｙ(t)とに応じた数式(8)の尤度Ｌ(hyp,pri,dat)を、数式(7)の尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)として算定する。

他方、解析対象のスロットが特定スロットに該当する場合（ＱB1：YES）、解析処理部２４は、解析対象である第ｔ番目のスロットに対して時間的に後方のスロットの観測データｙ(t:t+λ)を適用したカルマン平滑化（Kalman smoother）により第ｔ番目のスロットの尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する（ＱB3）。第１実施形態のカルマン平滑化は、ラグλを所定の正数（例えば５）に設定した固定ラグ平滑化である。

具体的には、解析対象のスロットが特定スロットに該当する場合、解析処理部２４は、以下の数式(9)で表現される尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する。

すなわち、解析処理部２４は、解析対象のスロットにおけるチャネル統計量ψ(t)⁽ⁿ⁾の推定値ψ~(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおけるチャネル歪θ(t-1)の確率分布（ｓ~(t-1)⁽ⁿ⁾，Ｓ~(t-1)⁽ⁿ⁾）と、解析対象のスロットの観測データｙ(t)とに応じた数式(9)の尤度Ｌ(hyp,pri,dat)を、数式(7)の尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)として算定する。解析処理部２４は、第ｔ番目から第(t+λ)番目までの(λ+１)個のスロットの観測データｙ(t:t+λ)を適用した数式(10)のカルマン平滑化により、チャネル歪θ(t-1)の確率分布（平滑化分布）の平均ｓ~(t-1)⁽ⁿ⁾および共分散Ｓ~(t-1)⁽ⁿ⁾を算定する。

記号ＫＳはカルマン平滑化を意味する。数式(10)から理解される通り、解析処理部２４は、解析対象のスロットにおけるチャネル統計量ψ(t)⁽ⁿ⁾の推定値ψ~(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおけるチャネル歪θ(t-1)の確率分布（ｍ(t-1)⁽ⁿ⁾，Ｃ(t-1)⁽ⁿ⁾）と、第ｔ番目から第(t+λ)番目までの(λ+１)個のスロットの観測データｙ(t:t+λ)とを適用したカルマン平滑化により、数式(9)の平均ｓ~(t-1)⁽ⁿ⁾および共分散Ｓ~(t-1)⁽ⁿ⁾を算定する。

以上の説明から理解される通り、第１実施形態の解析処理部２４は、複数のスロットのうち特定スロット（t mod M＝0）についてはカルマン平滑化により尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する一方、特定スロット以外の各スロット（t mod M≠0）については、当該スロットまでの観測データｙ(t)を適用したカルマンフィルタにより尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する。特定スロット以外の各スロットではカルマンフィルタにより算定される尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)が、特定スロットではカルマン平滑化により算定される尤度ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)に修正される、とも換言され得る。以上が尤度算定処理ＱBの具体例である。

図４に例示した移動平均μ⁽ⁿ⁾および共分散Γの算定（ＱA1）と補助インデックスｋ_nの設定（ＱA2）とを以上の手順で実行すると、解析処理部２４は、粒子フィルタのＮ個の粒子の各々について、第１処理ＱC1と第２処理ＱC2と第３処理ＱC3とを実行する。各処理の具体例を以下に詳述する。

第１処理ＱC1は、超パラメータであるチャネル統計量ψ(t)⁽ⁿ⁾（すなわち拡張状態のうち非線形・非ガウス型の部分）を粒子の実現値として特定する処理である。具体的には、解析処理部２４は、補助インデックスｋ_nについて数式(5)の演算で算定された移動平均μ^(kn)と数式(6)の演算で算定された共分散Γとで規定される確率分布からチャネル統計量ψ(t)⁽ⁿ⁾を抽出する。チャネル統計量ψ(t)⁽ⁿ⁾の抽出対象となる確率分布の種類は任意であるが、例えばベータ分布が好適に採用され得る。

第２処理ＱC2は、解析対象である第ｔ番目のスロットにおけるチャネル歪θ(t)（すなわち拡張状態のうち線形・ガウス型の部分）の確率分布の平均ｓ(t)⁽ⁿ⁾および共分散Ｓ(t)⁽ⁿ⁾を実現値として算定する処理である。第２処理ＱC2で算定された平均ｓ(t)⁽ⁿ⁾および共分散Ｓ(t)⁽ⁿ⁾は、チャネル歪θ(t)の確率分布の平均ｍ(t)⁽ⁿ⁾および共分散Ｃ(t)⁽ⁿ⁾の実現値として直後のスロットで使用される（ｍ(t)⁽ⁿ⁾←ｓ(t)⁽ⁿ⁾，Ｃ(t)⁽ⁿ⁾←Ｓ(t)⁽ⁿ⁾）。第２処理ＱC2で算定される平均ｓ(t)⁽ⁿ⁾および共分散Ｓ(t)⁽ⁿ⁾の確率分布（例えば複素正規分布）を以下では「平滑化分布」と表記する。

第１実施形態の解析処理部２４は、以下の数式(11)で表現されるカルマン平滑化で平均ｓ(t)⁽ⁿ⁾および共分散Ｓ(t)⁽ⁿ⁾（すなわち平滑化分布）を算定する。第２処理ＱC2のカルマン平滑化は、前述のステップＱB3のカルマン平滑化と同様に、所定の正数（例えば５）をラグλとする固定ラグ平滑化である。

数式(11)から理解される通り、解析処理部２４は、第１処理ＱC1で算定したチャネル統計量ψ(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおけるチャネル歪θ(t-1)の確率分布（ｍ(t-1)^(kn)，Ｃ(t-1)^(kn)）と、第ｔ番目から第(t+λ)番目までの(λ+１)個のスロットの観測データｙ(t:t+λ)とを適用したカルマン平滑化によりチャネル歪θ(t)の平滑化分布（ｓ(t)⁽ⁿ⁾，Ｓ(t)⁽ⁿ⁾）を粒子の実現値（ｍ(t)⁽ⁿ⁾，Ｃ(t)⁽ⁿ⁾）として算定する。第１実施形態の解析処理部２４は、観測データ取得部２２が観測データｙ(t)を取得するスロットの全部について、第２処理ＱC2ではカルマン平滑化により平滑化分布を算定する。

第３処理ＱC3は、第ｎ番目の粒子の加重値ω(t)⁽ⁿ⁾を算定する処理である。具体的には、第１実施形態の解析処理部２４は、以下の数式(12)の演算で加重値ω(t)⁽ⁿ⁾を算定する。

数式(12)の分子の記号ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)は、第１番目から第(t-1)番目までのスロットの観測データｙ(1:t-1)とチャネル統計量ψ(t)⁽ⁿ⁾とが観測されたときの第ｔ番目のスロットの観測データｙ(t)の尤度（第２尤度の例示）を意味する。第１実施形態の第３処理ＱC3は、数式(12)の尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)を算定する尤度算定処理ＱDを包含する。

図６は、尤度算定処理ＱDのフローチャートである。尤度算定処理ＱDを開始すると、解析処理部２４は、図５の尤度算定処理ＱBと同様に、解析対象である第ｔ番目のスロットが特定スロットに該当するか否かを判定する（ＱD1）。具体的には、解析処理部２４は、解析対象のスロットの番号ｔを定数Ｍで除算した剰余が０である場合に、当該スロットが特定スロットに該当すると判定する一方、剰余が０でない場合には、解析対象のスロットが特定スロットに該当しないと判定する。

解析対象のスロットが特定スロットに該当しない場合（ＱD1：NO）、解析処理部２４は、第ｔ番目のスロットまでの観測データｙ(t)を適用したカルマンフィルタにより第ｔ番目のスロットの尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)を算定する（ＱD2）。すなわち、尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)は以下の数式(13)で表現される。

数式(13)から理解される通り、解析処理部２４は、第１処理ＱC1で算定したチャネル統計量ψ(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおけるチャネル歪θ(t-1)の確率分布（ｍ(t-1)^(kn)，Ｃ(t-1)^(kn)）と、解析対象のスロットの観測データｙ(t)とに応じた尤度Ｌ(hyp,pri,dat)を、数式(12)の尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)として算定する。

他方、解析対象のスロットが特定スロットに該当する場合（ＱD1：YES）、解析処理部２４は、第(t-1)番目のスロットについて第２処理ＱC2のカルマン平滑化で算定されたチャネル歪θ(t-1)の確率分布（ｓ(t-1)⁽ⁿ⁾，Ｓ(t-1)⁽ⁿ⁾）から尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)を算定する。具体的には、解析対象のスロットが特定スロットに該当する場合、解析処理部２４は、以下の数式(14)で表現される尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)を算定する。

すなわち、解析処理部２４は、第１処理ＱC1で算定したチャネル統計量ψ(t)⁽ⁿ⁾と、第(t-1)番目のスロットにおける第２処理ＱC2でカルマン平滑化により算定されたチャネル歪θ(t-1)の平滑化分布（ｓ(t-1)⁽ⁿ⁾，Ｓ(t-1)⁽ⁿ⁾）と、当該解析対象のスロットまでの観測データｙ(t)とに応じた尤度Ｌ(hyp,pri,dat)を、数式(12)の尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)として算定する。

以上の説明から理解される通り、第１実施形態の解析処理部２４は、複数のスロットのうち特定スロット（t mod M＝0）については、当該特定スロットの前方（第(t-1)番目）のスロットにおける第２処理ＱC2でカルマン平滑化により算定されたチャネル歪θ(t-1)の確率分布（ｓ(t-1)⁽ⁿ⁾，Ｓ(t-1)⁽ⁿ⁾）に応じた尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)から加重値ω(t)⁽ⁿ⁾を算定する一方、特定スロット以外のスロット（t mod M≠0）については、当該スロットまでの観測データｙ(t)を適用したカルマンフィルタで算定される尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)から加重値ω(t)⁽ⁿ⁾を算定する。特定スロット以外の各スロットではカルマンフィルタにより算定される尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)が、特定スロットではカルマン平滑化により算定される尤度ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)に修正される、とも換言され得る。以上が尤度算定処理ＱDの具体例である。

粒子フィルタのＮ個の粒子の各々について図４の第１処理ＱC1から第３処理ＱC3を反復すると、解析処理部２４は、第３処理ＱC3で粒子毎に算定した加重値ω(t)⁽ⁿ⁾を正規化する（ＱE1）。具体的には、以下の数式(15)で表現される通り、Ｎ個の粒子にわたる加重値ω(t)⁽¹⁾〜ω(t)^(N)の総和により個々の粒子の加重値ω(t)⁽ⁿ⁾を除算することで、Ｎ個の粒子にわたる総和が１となる加重値ω(t)⁽¹⁾〜ω(t)^(N)を算定する。

以上に説明した通り、第１処理ＱC1で算定されたチャネル統計量ψ(t)⁽ⁿ⁾と、第２処理ＱC2で算定されたチャネル歪θ(t)の確率分布（ｍ(t)⁽ⁿ⁾，Ｃ(t)⁽ⁿ⁾）とを実現値として、第３処理ＱC3の加重値ω(t)⁽ⁿ⁾で規定されるＮ個の粒子が生成される。解析処理部２４は、Ｎ個の粒子で表現される確率分布から解析対象のスロットのチャネル統計量ψ(t)を推定する（ＱE2）。具体的には、第１実施形態の解析処理部２４は、Ｎ個の粒子で表現される確率分布の代表値をチャネル統計量ψ(t)として算定する。例えば、Ｎ個の粒子にわたるチャネル統計量ψ(t)⁽¹⁾〜ψ(t)^(N)の平均値、または、加重値ω(t)⁽ⁿ⁾が最大となるチャネル統計量ψ(t)⁽ⁿ⁾が、第ｔ番目のチャネル統計量ψ(t)の推定結果として算定される。図１のチャネル補償部２６は、前述の通り、解析処理部２４がスロット毎に算定するチャネル統計量ψ(t)を利用したチャネル補償を実行する。

以上に説明した通り、第１実施形態では、無線通信路のチャネル統計量ψ(t)およびチャネル歪θ(t)の確率分布を推定する粒子フィルタにおいて、解析対象のスロットに対して時間的に後方のスロットの観測データｙ(t:t+λ)を適用したカルマン平滑化により尤度（ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)，ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)が算定される。したがって、以下に例示される通り、尤度の算定にカルマンフィルタのみを利用する構成（以下「対比例」という）と比較して、無線通信路にノイズが多い環境でもチャネル統計量ψ(t)を高精度に推定できるという利点がある。対比例は、非特許文献１の技術を無線通信路の解析に適用した構成である。

図７および図８は、チャネル統計量ψ(t)の一例である移動速度Ｖを対比例の構成で推定した結果のグラフである。複数回にわたる推定結果の平均と全体の９０％の推定値の分布とが図７および図８では併記されている。図７は、伝送データのビット毎のエネルギーＥbとノイズのエネルギーＮ0との比Ｅb/Ｎ0が２０［ｄＢ］である環境、すなわちノイズが少ない環境での推定結果である。他方、図８は、比Ｅb/Ｎ0が０［ｄＢ］である環境、すなわち図７と比較してノイズが多い環境での推定結果である。ドップラー周波数ｆDは１４０［Ｈｚ］であり、遅延スプレッドσTは１.０［μs］である。また、実際の移動速度Ｖ（正解）は７５.６［ｋｍ/ｈ］である。

図７から理解される通り、ノイズが少ない環境では、対比例の構成でも移動速度Ｖを高精度に推定することが可能である。しかし、対比例の構成では、ノイズが多い環境において推定精度が低下する、という傾向が図８から確認できる。推定精度の低下の原因としては、例えば推定処理での過剰適合が想定される。

図９は、図８と同様に比Ｅb/Ｎ0が０［ｄＢ］である環境（すなわちノイズが多い環境）のもとで第１実施形態の構成により移動速度Ｖを推定した結果のグラフである。図８に図示した対比例の推定結果が図９には併記されている。図７および図８と同様に、ドップラー周波数ｆDは１４０［Ｈｚ］であり、遅延スプレッドσTは１.０［μs］である。また、実際の移動速度Ｖ（正解）は７５.６［ｋｍ/ｈ］である。第１実施形態によれば、ノイズが多い環境でも、対比例と比較して移動速度Ｖを高精度に推定できることが、図９から確認できる。また、第１実施形態によれば、対比例と比較して推定結果が迅速に収束するという利点もある。

図１０は、移動速度Ｖの誤差の絶対値に関する累積分布（ＣＤＦ：Cmulative Distribution Function）を第１実施形態と対比例とについて併記したグラフである。図１０では、比Ｅb/Ｎ0の複数の数値（Ｅb/Ｎ0＝｛０,４,８,１２,１６,２０｝）とドップラー周波数ｆDの複数の数値（ｆD＝１４×｛２,４,６,８,１０｝［Ｈｚ］）と遅延スプレッドσT（σT＝１.０［μｓ］）との全通りの組合せが想定されている。第１実施形態によれば、対比例と比較して移動速度Ｖを高精度に推定できることが図１０からも確認できる。

ところで、時間軸上の全部のスロットについてカルマン平滑化により尤度（ｐ(y(t)|y(1:t-1),ψ(t)⁽ⁿ⁾)，ｐ(y(t)|y(1:t-1),ψ~(t)⁽ⁿ⁾)を算定する構成のようにカルマン平滑化による尤度の算定の頻度が過度に高い構成では、チャネル統計量ψ(t)の推定精度が却って低下する可能性がある。第１実施形態では、複数のスロットから選択された特定スロットについてはカルマン平滑化により尤度を算定する一方、特定スロット以外のスロットについてはカルマンフィルタにより尤度が算定されるから、カルマン平滑化の頻度の過多に起因した推定精度の低下を抑制してチャネル統計量ψ(t)を高精度に推定できるという利点がある。

＜第２実施形態＞
本発明の第２実施形態を説明する。なお、以下に例示する各形態において作用や機能が第１実施形態と同様である要素については、第１実施形態の説明で使用した符号を流用して各々の詳細な説明を適宜に省略する。

第１実施形態では、無線通信路のチャネル統計量ψ(t)が１種類の統計量（移動速度Ｖ）を含む場合を便宜的に例示したが、第２実施形態のチャネル統計量ψ(t)は、複数種の統計量を包含する。具体的には、第２実施形態のチャネル統計量ψ(t)は、移動速度Ｖ（第１統計量の例示）および遅延スプレッドσT（第２統計量の例示）の２種類の統計量を包含する。複数種の統計量を含むチャネル統計量ψ(t)は、第１実施形態と同様の構成および動作により推定することが可能である。したがって、第２実施形態においても第１実施形態と同様の効果が実現される。

ところで、図４の解析処理は、解析対象のスロットの直前のスロットにおけるＮ個の粒子を更新する処理であり、第１回目（すなわち初回）の解析処理では、無線通信路のチャネル統計量ψ(t)を特定の初期値に設定する必要がある。解析処理による推定精度は、第１回目の解析処理に適用されるチャネル統計量ψ(t)の初期値の適否に依存する。以上の事情を考慮して、第２実施形態の解析処理部２４は、図４の解析処理の開始前に図１１の初期値設定処理ＱFを実行することでチャネル統計量ψ(t)（移動速度Ｖ，遅延スプレッドσT）の初期値を設定する。観測データ取得部２２が最初に所定個（例えばλ個）の観測データｙ(t)を取得した段階で図１１の初期値設定処理ＱFが実行され、初期値設定処理ＱFで設定されたチャネル統計量ψ(t)の初期値が第１回目の解析処理に適用される。観測データ取得部２２が解析処理の開始前に取得する観測データｙ(t)を以下では「初期観測データ」と表記する。

図１１に例示される通り、初期値設定処理ＱFは、第１初期化処理ｑAと第２初期化処理ｑBとを包含する。第１初期化処理ｑAは、チャネル統計量ψ(t)のうち移動速度Ｖの初期値を設定する処理であり、第２初期化処理ｑBは、チャネル統計量ψ(t)のうち遅延スプレッドσTの初期値を設定する処理である。第１初期化処理ｑAの実行後に第２初期化処理ｑBが開始される。

第１初期化処理ｑAを開始すると、解析処理部２４は、遅延スプレッドσTを暫定的な所定値（例えば０）βに設定する（ｑA1）。そして、解析処理部２４は、所定値βに設定された遅延スプレッドσTのもとで移動速度Ｖを相異なる候補値ηAに変化させた複数の場合の各々について周辺尤度（marginal likelihood）ρAを算定する（ｑA2）。周辺尤度ρAは、所定個の初期観測データｙ(t)が観測されたという条件のもとで移動速度Ｖが候補値ηAである尤度を意味する。具体的には、解析処理部２４は、所定個の初期観測データｙ(t)を適用したカルマン平滑化により候補値ηA毎に周辺尤度ρAを算定する。

図１２は、移動速度Ｖ（候補値ηA）と周辺尤度ρAとの関係を例示するグラフである。図１２に例示される通り、移動速度Ｖの数値範囲ＲAを所定の間隔（刻み幅）で区画する各境界の数値が候補値ηAとして利用される。数値範囲ＲAは、移動速度Ｖが現実的にとり得る数値の範囲に設定される。図１２では、移動速度Ｖの下限値０［ｋｍ/ｈ］から上限値３００［ｋｍ/ｈ］までの数値範囲ＲAを５［ｋｍ/ｈ］の間隔で区画する複数の数値が候補値ηAとして例示されている。

図１２に例示される通り、周辺尤度ρAは、特定の候補値ηAの近傍で最大となる。第２実施形態の解析処理部２４は、周辺尤度ρAが最大となる候補値ηA（図１２の例示では５５［ｋｍ/ｈ］）を移動速度Ｖの初期値として選択する（ｑA3）。すなわち、複数の初期観測データｙ(t)から算定される周辺尤度ρAを最大化する候補値ηAが移動速度Ｖの初期値として採択される（ＭＡＰ（Maximum a posteriori）推定）。以上が第１初期化処理ｑAの具体例である。

第１初期化処理ｑAの実行後に第２初期化処理ｑBを開始すると、解析処理部２４は、第１初期化処理ｑAで選択した初期値に移動速度Ｖを設定する（ｑB1）。そして、解析処理部２４は、当該移動速度Ｖのもとで遅延スプレッドσTを相異なる候補値ηBに変化させた複数の場合の各々について周辺尤度ρBを算定する（ｑB2）。周辺尤度ρBは、所定個の初期観測データｙ(t)が観測されたという条件のもとで遅延スプレッドσTが候補値ηBである尤度を意味する。具体的には、解析処理部２４は、所定個の初期観測データｙ(t)を適用したカルマン平滑化により候補値ηB毎に周辺尤度ρBを算定する。

図１３は、遅延スプレッドσT（候補値ηB）と周辺尤度ρBとの関係を例示するグラフである。図１３に例示される通り、遅延スプレッドσTの数値範囲ＲBを所定の間隔で区画する各境界の数値が候補値ηBとして利用される。数値範囲ＲBは、遅延スプレッドσTが現実的にとり得る数値の範囲に設定される。図１３では、遅延スプレッドσTの下限値０.２［μｓ］から上限値５.０［μｓ］までの数値範囲ＲBを０.１［μｓ］の間隔で区画する複数の数値が候補値ηBとして例示されている。

図１３から理解される通り、周辺尤度ρBは、特定の候補値ηBの近傍で最大となる。解析処理部２４は、周辺尤度ρBが最大となる候補値ηB（図１３の例示では０.９［μｓ］）を遅延スプレッドσTの初期値として選択する（ｑB3）。すなわち、複数の初期観測データｙ(t)から算定される周辺尤度ρBを最大化する候補値ηBが遅延スプレッドσTの初期値として採択される。以上が第２初期化処理ｑBの具体例である。前述の通り、初期値設定処理ＱFの実行後に、第１初期化処理ｑAで設定された移動速度Ｖの初期値と第２初期化処理ｑBで設定された遅延スプレッドσTの初期値とを適用した図４の解析処理が実行される。

以上に説明した通り、第２実施形態では、所定値βの遅延スプレッドσTのもとで周辺尤度ρAを最大化する候補値ηAが移動速度Ｖの初期値として設定され、当該初期値の移動速度Ｖのもとで周辺尤度ρBを最大化する候補値ηBが遅延スプレッドσTの初期値として設定される。したがって、例えば移動速度Ｖや遅延スプレッドσTの初期値を乱数に設定する構成と比較して、解析処理によりチャネル統計量ψ(t)を高精度に推定することが可能である。

＜変形例＞
前述の形態は多様に変形され得る。具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様は、相互に矛盾しない範囲で適宜に併合され得る。

（１）前述の各形態では、無線基地局１２から端末装置１４までの下りリンクの無線通信路を解析処理の対象としたが、端末装置１４から無線基地局１２までの上りリンクの無線通信路についても、前述の各形態と同様の構成および方法によりチャネル歪およびチャネル統計量を推定することが可能である。上りリンクの無線通信路を解析する構成では、観測データ取得部２２および解析処理部２４（さらにはチャネル補償部２６）が無線基地局１２に搭載される。

（２）第２実施形態では、第１初期化処理ｑAで移動速度Ｖの初期値を設定するとともに第２初期化処理ｑBで遅延スプレッドσTの初期値を設定したが、第１初期化処理ｑAで初期化される第１統計量と第２初期化処理ｑBで初期化される第２統計量との組合せは以上の例示に限定されない。例えば、移動速度Ｖ以外のチャネル統計量（例えば遅延スプレッドσT）を第１統計量として第１初期化処理ｑAで初期値を設定し、遅延スプレッドσT以外のチャネル統計量（例えば移動速度Ｖ）を第２統計量として第２初期化処理ｑBで初期値を設定することも可能である。

（３）前述の各形態では、移動速度Ｖおよび遅延スプレッドσTを例示したが、無線チャネル解析装置が解析するチャネル統計量は以上の例示に限定されない。例えば、サブキャリア毎のチャネルゲインに関する平均パワーをチャネル統計量として推定することも可能である。

（４）前述の各形態で例示した観測データ取得部２２および解析処理部２４としてコンピュータを機能させるプログラムは、コンピュータが読取可能な記録媒体に格納された形態で提供されてコンピュータにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ-ＲＯＭ等の光学式記録媒体が好例であるが、半導体記録媒体や磁気記録媒体等の公知の任意の形式の記録媒体を包含し得る。なお、通信網を介した配信の形態で前述のプログラムを提供することも可能である。

１００……通信システム、１２……無線基地局、１４……端末装置、２０……無線通信部、２２……観測データ取得部、２４……解析処理部、２６……チャネル補償部。

Claims (6)

1. 無線通信路を介して受信された複数のパイロット信号を含む観測データを時間軸上の各スロットについて取得する観測データ取得部と、
   前記観測データ取得部が取得した観測データを適用した尤度の算定を含む粒子フィルタにより、前記無線通信路のチャネル統計量およびチャネル歪の確率分布をスロット毎に推定する解析処理部とを具備し、
   前記解析処理部は、解析対象のスロットの後方のスロットの観測データを適用したカルマン平滑化により前記解析対象のスロットにおける前記尤度を算定する
   無線チャネル解析装置。
2. 前記解析処理部は、前記チャネル統計量の推定値と、前記解析対象のスロットの前方のスロットにおける前記チャネル歪の確率分布と、前記解析対象のスロットの後方のスロットの観測データとを適用したカルマン平滑化により、前記解析対象のスロットにおける第１尤度を前記尤度として算定する
   請求項１の無線チャネル解析装置。
3. 前記解析処理部は、複数のスロットのうち特定スロットについてはカルマン平滑化により前記第１尤度を算定する一方、前記特定スロット以外のスロットについては、当該スロットまでの観測データを適用したカルマンフィルタにより前記第１尤度を算定する
   請求項１または請求項２の無線チャネル解析装置。
4. 前記解析処理部は、前記粒子フィルタの複数の粒子の各々について、前記スロット毎に、
   前記チャネル統計量を当該粒子の実現値として確率分布から算定する第１処理と、
   前記第１処理で算定した前記チャネル統計量と、前記解析対象のスロットの前方のスロットにおける前記チャネル歪の確率分布と、前記解析対象のスロットの後方のスロットの観測データとを適用したカルマン平滑化により、前記解析対象のスロットにおける前記チャネル歪の確率分布を当該粒子の実現値として算定する第２処理と、
   前記第１処理で算定した前記チャネル統計量と、前記解析対象のスロットの前方のスロットにおける前記第２処理で算定された前記チャネル歪の確率分布と、当該解析対象のスロットまでの観測データとに応じた第２尤度を前記尤度として算定し、前記第２尤度から当該粒子の加重値を算定する第３処理とを実行する
   請求項１から請求項３の何れかの無線チャネル解析装置。
5. 前記解析処理部は、
   前記第３処理において、複数のスロットのうち特定スロットについては、当該特定スロットの前方のスロットにおける前記第２処理で算定された前記チャネル歪の確率分布に応じた前記第２尤度から前記加重値を算定する一方、前記特定スロット以外のスロットについては、当該スロットまでの観測データを適用したカルマンフィルタにより算定される前記チャネル歪の確率分布に応じた前記第２尤度から前記加重値を算定する
   請求項４の無線チャネル解析装置。
6. 前記チャネル統計量は、第１統計量および第２統計量を含み、
   前記解析処理部は、
   前記第２統計量を所定値に設定するとともに前記第１統計量を相異なる第１候補値に設定した複数の場合の各々について、前記観測データ取得部が最初に取得した所定個の観測データを適用したカルマン平滑化により周辺尤度を算定し、当該周辺尤度が最大となる前記第１候補値を前記第１統計量の初期値として選択する第１初期化処理と、
   前記第１初期化処理で選択された初期値に前記第１統計量を設定するとともに前記第２統計量を相異なる第２候補値に設定した複数の場合の各々について、前記所定個の観測データを適用したカルマン平滑化により周辺尤度を算定し、当該周辺尤度が最大となる前記第２候補値を前記第２統計量の初期値として選択する第２初期化処理とを実行する
   請求項１から請求項５の何れかの無線チャネル解析装置。
   

Images