nRF52833で測距(7) - IQデータの処理方法(4)

ちょっと本気出して測距を考えてみたい。IQデータの料理方法についてまとめておきたい。、、、の第4回目。MUSICやる。(タイトルにnRF52833ってもう要らんなー、、、が、今さら消すわけにもいかない)

CaponまでやったらMUSICもやってみる(<--Geminiのアドバイスによる(´ㅂ`; ))。で、例によってMUSICとはどういうものかGeminiに解説してもらったのでこちらをどうぞσ(￣∇￣; )。ちな、Geminiが結論のところに書いてくれるとよかったんだけど、MUSICスペクトルは信号そのものの強さではなく数学的な直交性(離れ具合)がプロットされているので、MUSICスペクトラムの「ピークの高さ」と「信号の電力（強さ）」には直接的な相関関係は無いです。(しかし、これ考えた人ほんとすごいね、発想が。どういう思考経路でここに至るのか見当もつかん。)
で、MUSICスペクトラムは
\[ P_{\text{MU}}(t) = \frac{\mathbf{a}^H(t) \mathbf{a}(t)}{\mathbf{a}^H(t) \mathbf{E}_n \mathbf{E}_n^H \mathbf{a}(t)} \tag{1} \] ってなっていて、$\mathbf{R}_{xx}$の固有値分解がこうなっていて、
\[ \mathbf{R}_{xx} = \mathbf{U} \mathbf{\Lambda} \mathbf{U}^H = \sum_{m=0}^{M-1} \lambda_m \mathbf{e}_m \mathbf{e}_m^H \tag{2} \] (添え字を0始まりに変更)
で、$M$個の固有値$\lambda_m$を昇順に並べたときに、どっかからどっかまでが信号で、どっかからどっかまでがノイズって(ある意味勝手に)決める。で、ノイズに該当するとした固有ベクトルを連結して$\mathbf{E}_n$とするわけだ
\[ \mathbf{E}_n = (\mathbf{e}_{0}, \mathbf{e}_{1}, \dots, \mathbf{e}_{M-1-K}) \tag{3} \] 、、、numpyで固有値計算すると、がっちゃんこなって固有ベクトルが(行列で)出てくるから連結するというよりは取り出すってかんじになる。
で、$a(t)$はCaponのときと同じステアリングベクトル
\[ \mathbf{a}(t) = \begin{pmatrix} 1 \\ e^{-j 2 \pi \Delta f t} \\ \vdots \\ e^{-j 2 \pi (N-1)\Delta f t} \end{pmatrix} \tag{4} \]
いちおうこれでMUSICスペクトラムが計算できるのでやってみる。
ble_music.py

結果はこうなる

たしかにMUSICってよさそうに見える。ここまで来たら一旦データ採りまくって検証したいけど、狭い我が家では2mくらいが精いっぱいT_T。公園でやったら怪しいし、会社でやったら会社の資産だから掲載できないし、、、過疎ってる実家付近ならやり放題なんだけどなー、、、ちな、今回使っているモジュールは技適マークあるので違法電波にはならないよー


物理的制約のある中、無理難題をうまく解決するのがエンジニアリングだと思っている。元エンジニアでも管理職になるとエンジニアじゃなくなる人いるんだね。残念なことです。

nRF52833で測距(6) - IQデータの処理方法(3)

ちょっと本気出して測距を考えてみたい。IQデータの料理方法についてまとめておきたい。、、、の第3回目。相関係数の前処理について。空間平均法やる。

信号源が一つ（または反射による相関が極めて強い）場合、そのままでは相関行列のランクが不足するという数学的課題を解決するために前回Tamed Caponってのをやったけど、空間平均法もやってみる。空間平均法は行列のサイズを擬似的に小さくして平均化することで「ランク（情報の深み）」として復元する手法らしい。

空間平均法のアルゴリズム
全体のデータ数をN、サブアレイのサイズをMとして、M個の要素を持つ部分ベクトル$\mathbf{x}_i$を作成する。
全体データは
\[ \mathbf{x}=\begin{pmatrix} x_1,x_2,...,x_{N-1} \end{pmatrix}^T \] 以下のように部分ベクトルを作っていく \[ \mathbf{x}_{0}=\begin{pmatrix} x_0,x_1,...,x_{M-1} \end{pmatrix}^T \] \[ \mathbf{x}_{1}=\begin{pmatrix} x_1,x_2,...,x_{M+1-1} \end{pmatrix}^T \] \[ \mathbf{x}_{2}=\begin{pmatrix} x_2,x_3,...,x_{M+2-1} \end{pmatrix}^T \] \[ \dots \] \[ \mathbf{x}_{i}=\begin{pmatrix} x_i,x_{i+1},...,x_{M+i-1} \end{pmatrix}^T \] \[ \dots \] \[ \mathbf{x}_{L-1}=\begin{pmatrix} x_{L-1},x_{L-1+1},...,x_{M+L-1-1} \end{pmatrix}^T \] で \[ M+L-1-1=N-1 \] なので \[ L=NｰM+1 \] で、各サブセットごとに$M \times M$の相関行列$\mathbf{R}_{{xx}_i}$を算出して、加算平均すると最終的な相関行列$\mathbf{\bar{R}}_{xx}$が得られるらしい。すなわち
\[ \mathbf{R}_{{xx}_i} = \mathbf{x}_i \mathbf{x}_i^H \tag{1} \] \[ \mathbf{\bar{R}}_{xx} = \frac{1}{L} \sum_{i=0}^{L-1} \mathbf{R}_{{xx}_i} \tag{2} \] なんでこれでランクが上がるのかをGeminiに聞いてみた
なぜこれで「ランク」が上がるのか
無線測距のデータがランク $1$ になってしまうのは、各周波数成分の位相関係が完全に固定（コヒーレント）されているからです。空間平均法で窓をずらしていくと、各サブセット $\mathbf{x}_i$ にとって「基準となる位相」が $i$ ごとに変化します。この「位相のずれ」を含んだ状態で平均をとることで、数学的には「信号源が統計的に独立であるかのように振る舞わせる」ことができ、行列のランクが回復（最大 $M$ まで）します。
、、、らしい(^∀^;)
で、これでもランクが足りない場合、「フォワード・バックワード空間平均（FBSS: Forward-Backward Spatial Smoothing）」というのをやるらしい。
Forwardとは(上に書いてあるのと同じ)、
\[ \mathbf{x}_{fi}=\begin{pmatrix} x_i,x_{i+1},...,x_{M+i-1} \end{pmatrix}^T \tag{3} \] Backwardとは、 \[ \mathbf{x}_{bi}=\begin{pmatrix} x_{M+i-1}^*,...,x_{i+1}^*,x_i^* \end{pmatrix}^T \tag{4} \] これにより新たに別の独立したデータを取得したようになるらしい。で、
\[ \mathbf{R}_{xxi} = \frac{\mathbf{x}_{fi}\mathbf{x}_{fi}^H + \mathbf{x}_{bi}\mathbf{x}_{bi}^H}{2} \tag{5} \] として、(2)の計算をすればいいってことらしい。
で、これをどう実装するかというと、上に記載したアルゴリズム通りに実装してもいいんだけど、オリジナルの$\mathbf{R}_{xx}$も欲しいってこともある。んで、オリジナルの$\mathbf{R}_{xx}$から部分行列を切り出し切り出しして足して足して最後に割って空間平均された相関行列$\mathbf{\bar{R}}_{xx}$を得ようと思う(参考サイト : 菊間信良著　科学技術出版社 アレーアンテナによる適応信号処理 FORTRAN, MATLAB ソースリスト)。
で、作成したコードがこちら
ble_spatial_rxx_capon2.py

解析に使ってるデータはこちら(今まで使いまわしてきたやつ...いちおう)
previous_log.txt

で、結果

いちおうこういう結果が出て、BeamformerよりCaponの方がまあ確かに鋭くはなったんだけど、空間平均のサブアレー要素数とか、Tamedの係数とかで、だらっとさせることもできるし、鋭さを求めると、別のピークが現れてきたり、めっちゃぶっ飛んだり、、、どういう設定にするかで結果が変わりやすいんだけど、じゃあその設定を何を根拠に調整していけばいいのかがわからない。実機でチューニングする類のやつなのかな( ˘•ω•˘ ).｡oஇ


実は彼は天才で今週末でスパッと戦争を止められる秘策を持っていてほしい願望

nRF52833で測距(5) - IQデータの処理方法(2)

ちょっと本気出して測距を考えてみたい。IQデータの料理方法についてまとめておきたい。、、、の第２回目。

BeamformerができるんだったらCaponもできるのでは？って考えちゃう。んで、Caponってどんなんだったっけ？ってのは、、、Geminiに調べてもらった。前回はIQのフーリエ逆変換をどう理解するかっていう少し人間っぽい思考があった(まあわかってやってたんだけど)けど、今回はその必要がない(?)ので一気にGeminiに解説してもらう。
ってんで、こちら。 んで要点だけ抜き出すと
こちら、ステアリングベクトル(というらしい)
\[ \mathbf{a}(t) = \begin{pmatrix} 1 \\ e^{-j 2 \pi \Delta f t} \\ \vdots \\ e^{-j 2 \pi (N-1)\Delta f t} \end{pmatrix} \tag{1} \] ウェイトベクトル$\mathbf w$によって電力$P$は
\[ P=\mathbf w^H \mathbf R_{xx} \mathbf w \tag{2} \] って表されるんでした。見たい距離での出力を一定にしたうえで、$P$を最小化するってのがCapon法。それを数学的に書くと
\[ \min_{\mathbf{w}} \mathbf{w}^H \mathbf{R}_{xx} \mathbf{w} \quad \text{subject to} \quad \mathbf{w}^H \mathbf{a}(t) = 1 \tag{3} \] となるらしい。
これを解いて、Capon法のスペクトルは
\[ P_{CP}(t) = \frac{1}{\mathbf{a}(t)^H \mathbf{R}_{xx}^{-1} \mathbf{a}(t)} \tag{4} \] となるらしい。
(ステアリングベクトルって「距離によって位相差がこうなる」ってのを示していて、ステアリングベクトルの距離(すなわち$t$)をスイープすることで出力を最大化するってのがBeamformerでした。)
ところで、このCapon法は実は困ったことがあって、相関行列$\mathbf{R}_{xx}$の逆行列が存在しないってこと。今回IQデータのスナップショット数は1なので$\mathbf{R}_{xx} = \mathbf{x} \mathbf{x}^H$において$rnk(\mathbf{R}_{xx})=1$らしくて、これだと逆行列は存在しない。さらに測距において、LocalとRemoteのIQデータから得られる信号は、単一の経路（あるいは固定された反射路）を通っていて、強い相関(コヒーレンシー)を持っていて、いくらスナップショットを重ねても、位相関係が固定されているので、情報は重なったまま（コヒーレントなまま）で、行列のランクが上がらないってのも問題。これを解決するには

1. 空間平均法=サブアレー的な処理をして信号間の相関を崩す
2. 複数シーケンス=独立した複数のデータ取得で様々な要因を取り込んで信号間の相関を崩す
3. 対角成分付加法=数学的には固有値の底上げ

ていうアイデアがあるらしいんだけど、3が簡単そうなのでやってみる。1が正攻法っぽいけど、3は確実に逆行列を安定して得られるという点で良い方法らしい。ちな、「Tamed Capon」とか「Diagonal Loadingを適用したRobust Capon Beamforming」とかいうらしい。

ble_tamedcapon.py

あれ？するどくなってる？なってないような、、、σ(￣∇￣; )
やっぱり空間平均法もやってみるかー


ほんとまじで戦争やめて！
○○○○が○○○○○の首をとって○○○に謝ればいいのに

nRF52833で測距(4) - IQデータの処理方法(1)

ちょっと本気出して測距を考えてみたい。IQデータの料理方法についてまとめておきたい。

nRF52833で測距(2)では周波数vs位相の傾きから推定する方法で距離を推定しました(本当はSDKの出力との差をどう埋めるかもやらんといかんのだけど、、、)。そのほかにも方法があるみたいなので、調べてみる。
まずはIQデータをから得られる各周波数での位相差は以下のように表せるんでした。 \[ x_{RB,RX}(t) \cdot x_{IB,RX}(t)=e^{-j ( \omega_I + \omega_R ) \frac{D}{c} } \tag{1} \] ここで、もはや(1)において$t$は変数ではなく、$f$を変数とすることから、 \[ x(f)=x_{RB,RX}(t) \cdot x_{IB,RX}(t) \tag{2} \] として、 \[ \omega_I=\omega_R=2 \pi f \tag{3} \] とすると、 \[ x(f)=e^{-j\frac{4 \pi fD}{c}} \tag{4} \] ここで、 \[ f=f_0+n\Delta f \tag{5} \] \[ n=0,1,2, \dots ,N-1 \tag{6} \] とし(4)を \[ x_n=e^{-j\frac{4 \pi D}{c}(f_0+n\Delta f)} \tag{7} \] としたうえで、$x_n$をフーリエ逆変換することを考える。
このとき時間ドメインでの時間の最大値$T_{max}$について \[ T_{max}=\frac{1}{\Delta f} \tag{8} \] 時間の刻み$T_{step}$について \[ T_{step}=\frac{1}{N \Delta f} \tag{9} \] これより \[ \frac{k}{N}=\Delta f t \tag{10} \] フーリエ逆変換の式 \[ X_k = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{j \frac{2\pi}{N} nk} \tag{11} \] \[ k = 0, 1,2, \dots, N-1 \tag{12} \] に(7)(10)を代入すると、 \[ X(t) = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{j (2 \pi n \Delta f t)} \tag{13} \] 整理して \[ X(t)=\frac{1}{N} e^{-j \frac{4\pi D f_0}{c} } \sum_{n=0}^{N-1} e^{j \left\{ 2\pi n \Delta f \left( t - \frac{2D}{c} \right) \right\} } \tag{14} \] この(14)において、 \[ t=\frac{2D}{c} \tag{15} \] (15)以外の$t$だと$\sum$において、周期関数の位相が$n$によってばらばらになる。すなわち(15)の場合にのみ各$f_n$の位相ずれが「補正」され、ばらばらだった波が重なり最大振幅になることを示している。すなわちフーリエ逆変換後のピークを与える$t$が距離による遅延時間である。ん？あたりまえやん。まあ数式にしたかったのよ。
ところで、(13)を見てみると、$\frac{1}{N}$を除いた部分 \[ y(t) = \sum_{n=0}^{N-1} x_n e^{j (2 \pi n \Delta f t)} \tag{16} \] においては、受信波形にウェイトを与えているように見ることもできる。ウェイトを$w$($w$は$t$の関数)として受信電力を最大化する$w$を求めるすなわち$t$を求めることが距離を推定することといえるが、これってアレーアンテナを使ったAoA(Angle of Arrival)のBeamformer法と同じ(AoAは複数のアンテナの幾何学的配置の差で受信波形に位相差が出ることを利用して発信源の方向(角度)を推定する。今回のものは複数の周波数によって受信波形に位相差が出ること利用して発信源の距離を推定する。)。数式的に全く同じことだってわかる。
(16)を少し書き換えてみる。 \[ \mathbf{x}=\begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{pmatrix} \tag{17} \] \[ \mathbf{w}=\begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_{N-1} \end{pmatrix} \tag{18} \] (この辺の式を転置で書いている文献が多いんだけど、当方おつむが弱くてイメージするのにひと呼吸必要になってしまうので、行数が増えてしまうけど素直に書いちゃう)
として、 \[ y(t)=\mathbf w^H \mathbf x \tag{19} \] であると考える。
(19)は \[ y(t)=\mathbf x^T \mathbf w^*=\mathbf x^H \mathbf w \tag{20} \] と書き換えることができるため、得られる電力$P$を \[ P=|y(t)|^2=\mathbf w^H \mathbf x \mathbf x^H \mathbf w = \mathbf w^H \mathbf R_{xx} \mathbf w \tag{21} \] \[ \mathbf R_{xx}=\mathbf x \mathbf x^H \tag{22} \] と表せる。$\mathbf R_{xx}$は$\mathbf x$に$\mathbf x$の複素共役転置を掛けたもので、相関行列と言われていて、各周波数ごとの位相の関係性を表している。 ちなみに、nRF52833からは各条件で1つのIQデータしか取得できないので、(21)(22)の表現になっているが、各条件で複数のIQデータが取得できる場合は、$R_{xx}$を平均すればいい。(式をじっくり見ればわかる)
(21)において$P$を最大化する$\mathbf w$($\mathbf w$は$t$の関数)を探索することが、信号源との距離によって発生する$t$を求めることになる。が、$\mathbf w$のエネルギーを一定に保つという条件を加える必要があるので$P$をウェイト$\mathbf w$のエネルギー$\mathbf w^H \mathbf w$で割って正規化する。これを$P_{BF}$とすると、 \[ P_{BF}(t)=\frac{\mathbf w^H \mathbf R_{xx} \mathbf w}{\mathbf w^H \mathbf w} \tag{23} \] $\mathbf w$について考えてみる。(16)より \[ y(t)= \begin{pmatrix} e^{j0\Delta \omega t} & e^{j1\Delta \omega t} & e^{j2\Delta \omega t} & \ldots & e^{j(N-1)\Delta \omega t} \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{N-1} \end{pmatrix} \tag{24} \] よって、 \[ \mathbf w= \begin{pmatrix} e^{-j0\Delta \omega t} \\ e^{-j1\Delta \omega t} \\ e^{-j2\Delta \omega t} \\ \vdots \\ e^{-j(N-1)\Delta \omega t} \end{pmatrix} \tag{25} \] これをもって$\mathbf w$を$t$で算出しながら実測から算出した$\mathbf R_{xx}$を使って$P$を計算していって、$P$のピークを探索すれば$t$が推定できるってこと。なので、それをやるコードはこちら(てかGeminiに出してもらってちょこっと修正)
ble_beamformer.py

これを実行した結果、

2.54mってなってる。SDKのifftでは0.88mなのでこりゃまた異なってる。

ちなみに、データを何も考えずにフーリエ逆変換するとこうなる。
ble_simple_ifft.py

2.11mってなってる。ちなみに単純にフーリエ逆変換すると(9)より \[ T_{step}=14.08ns \tag{26} \] \[ 2D_{step}=T_{step} \cdot c \tag{27} \] より \[ D_{step}=2.11m \tag{28} \] これってSDKのifftも単純なフーリエ逆変換じゃないってことよね(;￣Д￣)

ちなみに(11)を使ってフーリエ逆変換すると、こうなる
ble_ifft.py

ふむ。Beamformerと同じだよね。


合併を持ちかけて株が下がる。合併を取り下げて株が下がる。もはやコント。

C言語ソースコードをブログに埋め込む

ソースコードをブログにきれいに埋め込みたい

もうかなり長いことブログを書いているのに最適な方法が見つからない。Google BloggerもJimdoも同じコードで同じように表示されるようにっていう課題もある。
今まで試したもの

Jimdoの「文章」コンテンツ ==> Google Bloggerで使えない

srctohtml ==> 概ね満足してたけどCopyボタンがない

CarbonでiFrameを取得 ==> コードが長いと欠落

Geminiにコードそのものを食わせてhtmlで出力させる ==> Geminiがコードそのものを勝手に咀嚼したりする

で、ひねり出したアイデアが、carbonのiFrameをJavascriptで取得してhtml形式で出力するってんだけど、このJavascriptをGeminiに生成させるって超手抜き。

で、できたコードがこちら


これをブラウザで開いて、テクストボックスにコードを入力してボタンを押すとhtmlになる。このhtmlをコピーしてブログ内にペーストする。
しばらくはこれでいってみようと思う。