前回、パイロット信号が正しく取り出せないのはキャリア周波数オフセット(Carrier Frequency Offset, CFO)が原因ではないかと書きました。それ以後、 OFDM Baseband Receiver Design for Wireless Communicationsを参考に、 CFO, SCO(Sampling Clock Offset)を除去するためのPLLを実装してみたりもしましたが、どうにも思うとおりにはなりませんでした。

考え直してさらに調査した結果、DFTのデータを取り出すウィンドウが間違っていたということが判明しました。

これまでは、相関値が最大になる場所をDFTウィンドウの先頭としていましたが、 それだとCyclie Prefix(ガードインターバル)の先頭からOFDMシンボルを切り出してしまうことになります。 この領域は、信号に遅延波があるときに、他のシンボルからの干渉を避けるためのものなので、 本来DFTに使用するのは望ましくありません。

遅延波が全く無かったとしても、上記の文献の式(5.11)の通り、周波数領域の信号Zkは、送信信号Xk、チャンネル応答Hkに対して次のようになります。

ここに、TdはDFTウィンドウがCyclic Prefixに入り込んでいる時間です。

これまでのTs=1, Td=128, N=1024で計算してみると、e^{-j2pi128k/1024}= e^{-jpi*k/4}となります。つまり、隣り合ったサブキャリアの間で-pi/4だけ位相が本来の信号からずれていくことになります。

これまでの結果を振り返ってみると、サブキャリア12本ごとにpi位相がずれている様子でした。 計算すると、-pi/4*12=-3pi=pi (modulo 2pi)ですので、ちょうどこの結果と一致します。

それで、CPの開始からDFTウィンドウを設定するのではなく、CPの終わりの少し前からDFTウィンドウを設定することで、ほぼ所望の結果が得られるようになりました。

上はチャンネル補正をする前の、生のDFT後のデータです。 これだけバラけていると、QPSKであることも疑わしい感じですね。

SP自体は上のように、まだキャリアごとに、わずかに回転しています。このSPによって、残りのチャンネルは線形補間で推定します。推定結果でデータを補正すると、次になりました。

このように、CFO, SCOのためのPLLは一切使用せずに、相関値を算出したときのガンマの値による回転と、整数倍のチャンネル周波数補正、SPによるチャンネル補正のみで、QPSKの信号がほぼ復元できました。

流石に、この程度のコンステレーションでは64QAMは無理ですが、4象限しか区別しないQPSK程度であれば、これでも十分に思われます。また、TMCC, ACに相当するI軸上の点も見つけられます。

ここまで来れば、TMCCのデコードをはじめ、ワンセグのデコード処理も規格書に則って進めれば完成まで持っていけそうな気がしてきました。

前回までで、ワンセグの生データからのOFDMシンボル切り出し、パイロット信号を見つけるところまで実装しました。ところが、パイロット信号からチャンネル補正ができそうにない、というところで止まってしまいました。もう少し調査したところ、原因が分かってきました。 自分の理解の範囲でのまとめです。

OFDM Baseband Receiver Design for Wireless CommunicationsのChanter 5. Synchronizationの式(5.6)は次のようになっています。

これは、キャリア周波数のオフセットが、FFTした結果に及ぼす影響を記述した式です。 左辺Z_{i, k}は、(時系列に並んだ)i番目のシンボルの、k番目のサブキャリアの(周波数領域の)値を示しています。また、ε_I, ε_fがそれぞれキャリア周波数オフセットの整数部、小数部です。 Hはそれぞれのサブキャリアに対応するチャンネルの周波数応答です。

ごちゃごちゃしていますが、とにかく、ε_I, ε_fが0でなければ、Z = X Hとは記述できず、他のキャリア周波数の影響を受けてしまうことを示しています。

右辺の第一項を見てみると、大きく4つの掛け算になっていますが、最後の2つのe^{…}の部分は純粋な回転です。 これまでの実装で、小数部のキャリア周波数補正はできていると思われます。 そのため、ε_fはほとんど0と考えます。

したがって、振幅だけを見るならば、受信したk番目のサブキャリアのFFT結果は、送信側のk-ε_I番目のサブキャリアの値に、チャンネル応答を掛けたものになっています(sinの項はε_f->0で1になるため)。 ですから、パイロット信号も含めてXの振幅は既知なので、振幅について積算することで、 ε_Iを求めることができたのでした。

問題は3番目のeの項による回転です。 この式を見るとε_Iが0でなければ、時間と共に盛大に回転します。 (2019/2/20追記：これ自体は正しいですが、kには依存していませんので、 Scattered Pilotがキャリアによって回転量が異なることの説明にはなりません。 それより(5.6)右辺第2項が寄与しているのかも知れません)

これの逆数を求めて回転を補正することもできるかも知れませんが、 上記の式では右辺には第2項以降も存在します。 そのため、基本的にはε_I, ε_fを0に収束させるような制御をしたほうが良さそうです。

上記の本の5.2.3節にも、周波数領域の推定アルゴリズムには限度があるので、同期エラーはその限界内に来ていることを確かめないといけない、といった記述があります。そのためには、時間領域で補正することができる、と書かれています。

時間領域での補正方法は下記の図のようになります(Figure 5.13より)。

求められたε_I, ε_fをloop filterを通してNCOに与え、その周波数でサンプリングしたデータを回転させる、というものです。 これまでの実装ではε_fのみしか考慮していませんでしたが、ε_Iも考慮する必要がある、ということです。

以上を参考にして、PLLでキャリア周波数オフセットを0に制御できるか試してみたいと思います。

前回までで、ワンセグの生データからのOFDMシンボル切り出しあたりまで実装しました。今回はそれ以降の処理についてです。

前提知識

ワンセグも含めて、地上波デジタル放送は、ARIBにて規格書が配布されています。

標準規格（放送分野）一覧表のSTD-B31「地上デジタルテレビジョン放送の伝送方式」から入手可能です。私は会員でもないので、英語版を入手しています。

FFT以降の処理

前回の部分も少し修正して[11]については、次のようにしました。

重要なのはfftshiftを実施していることです。データ取得時のサンプリング周波数は、放送帯域の中心に設定しているので、(キャリア周波数オフセットが0であれば)fftshift後の配列中央の値が、放送帯域の中心キャリアになるはずです。

また、後でキャリア周波数の整数倍の補正をするために、周波数シフトをする必要があります。そのためのシフト量を保持するために、F_leftという変数を定義しています。シフト量が0であった場合に、FFT結果Fのうちキャリア番号0に対応する値を、初期値として設定しています。

パイロット信号を検出するためにPRBSが必要になるので、その定義と、TMCC、ACという特別なパイロット信号のキャリア番号を定義しています。PRBSの最初の値がキャリア0に、次の値がキャリア1に、…というように、432本のキャリアにそれぞれ0, 1のいずれかが対応します。SP(Scattered Pilot)信号は、このキャリア周波数毎の0, 1の値がBPSKにより、IQ空間上のそれぞれ(+4/3, 0), (-4/3, 0)にマップされます。これを上記では配列Wiで定義しています。

上記は、TMCCとACの場所を推定するための関数です。G2では、引数mで指定された分だけ、キャリア周波数を移動した場所をTMCC, ACとみなして、絶対値の積算をしています。

Gs = [G2(i, Wi) for i in range(-20, 20)]

にて、キャリア周波数を-20から20まで整数倍シフトして、G2の積算値が最大になる場所を求めます。最大値が得られるインデックスでF_leftを調整しています。これでキャリア番号0から431までの場所が求まったことになるはずです。上記のプロット結果は省略します。

次はSPの推定です。SP信号は12キャリアに1本あります。OFDMシンボル毎に場所が3キャリアずつ移動します。つまり、OFDMシンボル毎に、SPの場所はi12, i12+3, i12+6, i12+9, i*12,…のように巡回します。上記find_spでは、現OFDMシンボルでSPが何処に含まれているかを求めています。SPは、それ以外のデータキャリアに比べて4/3倍の絶対値を持っているので、上記のように積算した最大値を求めるという方法が使えます。

ここからがまだうまく行っていないと思われるのですが、上記では次のような出力が得られます(行が長くなりすぎるので、小数点以下は一部省略しています)。

[35.74547, 36.26827, 49.57394, 37.08745]
2
SPs:
6 1.0953 (-0.9832743533464458+0.482584084887267j) 1.3333
18 1.1773 (0.9411424714348882-0.7073685094875273j) 1.3333
30 1.2927 (0.9226357042956151-0.905485176200131j) -1.3333
42 1.3544 (-0.9451630023246448+0.9701050133998432j) -1.3333
54 1.3868 (-1.041607866189769+0.9155528236391836j) 1.3333
66 1.2944 (0.6906246945088386-1.0947704458085812j) 1.3333
78 1.0671 (0.6685409527954194-0.8316726977365281j) -1.3333
90 1.3859 (0.7967161344290873-1.1339926757265295j) 1.3333
102 1.1538 (-0.34105207413393135+1.102293053251151j) 1.3333
114 1.2401 (0.4762959068094629-1.144985106068488j) 1.3333
126 1.5607 (0.725557306133454-1.381834432391039j) -1.3333
138 1.2347 (-0.6547734652174335+1.0467712733713745j) -1.3333
150 1.3194 (-0.3329626554307+1.2767023965626096j) 1.3333
162 1.5988 (0.5091991640004263-1.5155690510396185j) 1.3333
174 1.3907 (-0.6562514301534865+1.2261640531390314j) 1.3333
186 1.6392 (-0.5092859972519473+1.5580715920605588j) -1.3333
198 1.5107 (0.26879055172568034-1.4866408415508643j) -1.3333
210 1.6382 (0.7524038499463859-1.4552230916819244j) 1.3333
222 1.3896 (0.44737449638796445-1.3156195604230028j) -1.3333
234 1.3833 (0.5356609612670982-1.2753992069341795j) 1.3333
246 1.559 (0.3497743406085074-1.5192501362727482j) -1.3333
258 1.1333 (-0.42997707286624176+1.048549775299747j) -1.3333
270 1.3702 (-0.2787479092644133+1.3415065163031206j) 1.3333
282 1.6204 (-0.5708941597565+1.5164865623719384j) -1.3333
294 1.4883 (0.6930350159378798-1.317100025661765j) -1.3333
306 1.5654 (-0.7261325631313592+1.3867477986830758j) -1.3333
318 1.5449 (0.4601928773417365-1.4747525850693373j) -1.3333
330 1.4912 (0.5746828035844906-1.3760673756943445j) 1.3333
342 1.4701 (0.825690594740095-1.2162618274644958j) -1.3333
354 1.4209 (0.24893414863314506-1.398884931842842j) 1.3333
366 1.3512 (-0.6292706725628592+1.1957269229318115j) 1.3333
378 1.4125 (-0.5413499265712505+1.3046784927444803j) -1.3333
390 1.4662 (0.7501030028129436-1.2598034492601253j) -1.3333
402 1.1753 (-0.7940964094918974+0.8664208985515854j) -1.3333
414 1.3589 (-0.7997872953740967+1.098570575241878j) 1.3333
426 1.0329 (-0.3727620060886928+0.96334271774669j) -1.3333

最初の行と次の行は、SPの場所と絶対値の積算値を表示しています。49は残りの3つの値36付近に比べると、ほぼ4/3倍となっているので、結果は正しそうに思われます。

ただ、それ以降、SPのキャリア番号と絶対値、その値、また、期待される値を各行表示していますが、 これがおかしいと思われます。例えば、キャリア番号6は、値が第2象限である一方、期待される値は(4/3, 0)です。次に、キャリア番号18では、値が第4象限である一方、こちらも期待される値は(4/3, 0)です。つまり、この結果が正しいのであれば、キャリア番号6と18とで、チャンネルの影響で回転した方向が逆になってしまっています。

しかも、よく観察すると、この結果は、キャリア番号6+24nと18+24nとでいつも成り立っているようです。つまり、回転方向がこれら2つの間でいつも逆です。

さらに、上記SPを実際にプロットした結果が次です。

この表示自体は、実数軸上にあるSPがどちらかに回転したと考えれば、それなりの結果に見えます。 ただ、上記の通り、回転方向がいつも逆になっているので、恐らくどこかに間違いがあるはずです。 現在のところここで躓いています。これだと、単純に残りのキャリアの推定を線形補間で行うと、絶対値が0近辺の値が出てきてしまうことになり、正しく補正できません。チャンネル推定をして補正を行う前に、上記の問題を解決する必要があります。

前回、GNU RadioとADALM PLUTOを使って、 ワンセグの生データを取得してファイルに保存しました。

いきなり余談(Jupyterを使うまで)

現在、Pythonにて上記のデータをデコードできるかを実験しています。 最初はCでプログラムを作成し始めました。でも、実装していると、 途中途中までで、思い通りに処理ができているかを確認したくなります。 そのために、グラフを作成する必要が出てきました。CImgなど、 C++で比較的に簡単に使えるライブラリも試しました。これは画像生成は簡単ですが、 グラフ生成は面倒に感じました。

それで、結局Jupyter Notebookを使ったほうが簡単だということで、Pythonに移行しました。アルゴリズムが固まったら、データを一気に処理するためにC,C++に戻る予定ですが。

ここから本題

キャプチャした生データをソフトウェアで処理して、MPEG TSを得るところまで持っていきたいと考えています。そのためには、まずなによりも、OFDMの復調処理が必要です。 ここがクリアできれば、残りはそれほど難しくないと予想しています。

以下、Notebookの画像を貼り付けていきます。

上記In [1]は、各種ライブラリを使用することを宣言しています。[2]でcapture.rawという、GNU Radioで保存したファイルを開いて、dataという変数に読み込みます。 ファイルの一番先頭は少し飛ばしています。

[3]はガードインターバルとシンボル長の定義です。 ADALM PLUTOでは、64/63MHzでサンプリングしています。1シンボルは1008[us]なので、 計算してみると、ちょうど1024というキリの良い個数になります(そうなるようにサンプリング周波数を設定している)。

(少なくとも日本の)ワンセグのガードインターバルはシンボル長の1/8なので、128になります。

[4]は、シンボルの境界を検出するための相関値を求める関数correlateの定義です。 この方法については、 “An Open and Free ISDB-T Full_Seg Receiver Implemented in GNU Radio” Federico Larroca et al.で最初に知りましたが、 “ML Estimation of Time and Frequency Offset in OFDM Systems” Van De Beek et al.がオリジナルのようです。

最初の論文の(1)式の計算になります。ざっくり説明すると、シンボル長だけ離れた、ガードインターバル長サンプルの相関(積和, 上記gamma)を計算しています。また、RHOは本来であればSNR/(1+SNR)で定義されるのですが、今回は定数としています。十分SNRが大きければ1に近づく値なので、 ひとまずこのように置いても問題なさそうです。

[5]では、実際に先頭をずらしながら相関値を計算していき、配列corに入れます。 correlateが返すのはタプルなので、[6]で分解してそれぞれの配列を得ます。 配列corrsの値が極大になっているところが、OFDMシンボルの区切りです。

[7]はグラフを見やすくするためのサイズ調整です(横幅を伸ばしている)。

[8]で相関値とガンマの値を時系列のグラフにしています。結果が次です。

上記Van De Beekの論文のFigure 4と同じような結果です。横軸はサンプルの番号です。赤い線が左の縦軸に対応しており、相関値を表します。青い線が右の縦軸に対応しており、ガンマの偏角を表します。

相関値がおよそ1024+128サンプル(OFDMシンボル長+ガードインターバル長)間隔で極大値を取り、 その時にガンマの偏角は大体同じ値(-2.0近辺)を示しています。

[9]は、シンボルごとの相関値の極大値(配列maxes)、 極大値を与える配列インデックス(配列indices)、その際のガンマの偏角を求めています。

ここまで、OFDMシンボル境界検出は上手く行っていますが、[10]以降はまだ試行錯誤しています。 [10]は、ガンマの偏角を用いて、サンプルデータを回転させてからFFTで周波数領域に移しています。 これは送信機と受信機のOFDMキャリア周波数の偏差(小数部のみ)を補正することに相当するはずです。

[11]では試しに、1つの相関のピーク(インデックス5530)を先頭として、 1OFDMシンボルをFFTしています。これを虚数平面にプロットしたのが次の画像です。

なんだかバラバラです。キャリア周波数の補正ができていれば、 あとはキャリアの位相差だけしか残らないと考えていました。 そうであれば、パイロット信号(Scatterd Pilot: BPSK, TMCC: DBPSK)を除くワンセグのデータはQPSKのため、 GNU Radioで上手く行っているときのコンステレーション画像：

を回転したような結果が得られるかと期待したのですが、全くうまく行っていないようです。

[12]は、FFTした結果の最大値を見つけています。計算が正しければ、これがパイロット信号になっているはずなのですが…

ひとまずここまでです。

• FFTで戻すデータの取り出し部分が間違っている?
• (OFDMシンボル期間をまたぐ)時間方向でパイロット信号と相関を取ってチャンネル補正なりをしないと(コンステレーションとしても)データが取り出せない?
• キャリア周波数が整数倍でずれていて、その補正が必要?

原因まで理解できていません。もう少し調査してみたいと思います。

2019/2/12追記

上記[10]が間違っていたようです。

rots = np.exp([-1j*arg*k/SYMBOL_LEN for k in range(SYMBOL_LEN)])

ではなく、

rots = np.exp([1j*arg*k/SYMBOL_LEN for k in range(SYMBOL_LEN)])

に直したところ(補正の回転方向が逆)、

上記程度にはなりました。少しはQPSKらしく見えてきました。これでTMCC検出やSP検出を行って等価処理をすれば、それなりの結果が得られそうです。

ダイアグラムを2つに分離して、ファイルを経由するように変更

これまでで、ADALM-PLUTOとGNU Radioを使って、 ワンセグ受信ができることは確認できました。 今回は、受信したデータを一旦ファイルに保存して、GNU Radioを使ってオフライン処理できるようにしました。 目的は、ADCのデータからMPEG TSのデコードまでを、GNU Radioを使用せず自作のプログラムで行うためです。

上記のようなダイアグラムを作成して実行すると、capture.rawというファイルが作成されました。 受信したデータのサンプリング周波数を調整して、LPFを通しただけのデータとなります。

次に、Pluto SDR Sourceの代わりにFile Sourceを置いて、OFDM Synchronization 1segの入力に繋ぎます。

上記を作成しました。こちらも実行すると、capture.rawに保存されたデータがオフライン処理され、 test_out.tsという名前でMPEG TSファイルが作成されました。

今後の目標としては、このcapture.rawを処理して、test_out.tsと同等の出力を得るプログラムを作成したいと考えています。

ちょっと苦労した点

これだけだと簡単に動作したように思われてしまうかも知れませんが、1時間くらい上手く動きませんでした。 ダイアグラム上はPluto SDR Sourceのコネクションを削除しているはずなのに、 なぜかリアルタイムで受信したものがOFDM処理されてしまい、画面上では接続されているFile Sourceからのデータが処理されませんでした。

どうやら、“Generate the flow graph"ボタンで失敗していたのが原因だったようです。これが失敗する場合でも、なぜか"Execute the flow graph"は押せることがあるようです。でも、実際にExecuteしても、Generateする前のダイアグラムが実行されるようで、画面ではPluto SDR Sourceは外しているのに、繋がっている状態の(古い)ダイアグラムが実行されていたようです。

思わぬ所でハマってしまいましたが、ひとまず生データの取得はできました。44秒位ですが、約360MB程度のファイルサイズです。1.0*10^6(サンプリング周波数)64/638(複素1サンプルバイト数)*44(秒)=352MBで、おおよそ計算どおりです。

また、家のアンテナ端子だと受信感度がイマイチだったので、DX アンテナ US10WBと、AliExpressにて、SMAオス<-> Fメスの変換アダプタを購入しました。無線機器を繋ごうとすると、こういったアダプタが(種類として)沢山必要になりますね。

データ型について

File Sinkで生成されるファイルは、GNU RadioではComplexとなっています。 Pythonのgr-utilsのgr_plot_iqのコードを見てみると、 numpy.complex64となっているようです。numpyのData Typesのページから、結局、実部、虚部の順に32-bit floatとしてバイナリデータとなっているようです。