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と同じだよね。


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