SS1の日記： ガロワ体のマトリックスでハミング （…）

今度は、面接授業でハミング符号理論のお勉強をしてきました。平日の授業なので、すべて参加することはできませんでしたが。主にハミング符号のところだけ受講しました。

誤り訂正符号設計入門 - 藤原英二 - 放送大学面接授業　シラバス
http://syl-web1.code.ouj.ac.jp/ouj-f232/dt-8554.html

講師の藤原英二東工大名誉教授は誤り訂正符号理論の専門家で、中でも拡張ハミング符号の改良である「奇数重み符号(M.H. Hsiao 1970」の応用をいろいろ研究してきた方です。

その中で，でた演習問題(同講義レジュメ、p.36，と書いても参照できる人はいないけど・・・・)がHP電卓で解くには手頃だったので、その結果からレポート書いて持っていったら，先生にウケたのでその話とか。

ハミング符号は、Wikipediaにあるように、一般化すると行列演算として検査行列（H）で表現されます。で、授業では放送大学生向けにわかりやすく、ビット列のANDとXORで論理回路を組む時のような解法で教えてて、んでまあ、電卓でもそういうコードを書いてみて、いろいろいじってたんですが、ふと、ベクトル演算でやってみたらどうなる？　と、書き換えてみたら、こっちのほうが実行時間が短かったという・・・　ベクトルのほうが早いって、どこの地球シミュレータですか！

ちなみにHP48の実効速度は50FLOPS（地球シミュレータのCPU一個分の約1.6×10^‐14倍）くらい。とはいえ、組み込みのMATRIX型を計算に利用すると実効速度きっかりの演算速度を出してくれる。こりゃべんりだ。てなわけで、あれこれ試していると、普通の学生よりハミング符号理論の理解が進んだようでした。

行列を素手で演算するととっても大変なんですが、面倒なところを電卓がやってくれるならこれくらい楽なのも無い。というわけで、HP電卓のマトリックス演算でハミング符号をn次の2元数｛0, 1}にして計算してみた。hp48のMATRIX演算には、いろいろ制限があるので、それに合わせた計算手順になってます。なお，以下の話は先生に添削してもらったわけではないので，間違ってるところもあるかと。

ハミング符号をHP電卓で生成してみる

Wikipediaにあるハミング符号(7,4)の検査行列をHP電卓のMATRIX型で表現するとこんな感じ、

H:[[ 1 0 1 1 1 0 0 ][ 1 1 0 1 0 1 0 ][ 0 1 1 1 0 0 1 ]]　ハミング行列

で、データ列は、こんなかんじ。

D:[ 1 0 1 1 ]　データ列

んでまあ、ハミング符号(De)の生成は、

De='D*TRN(H)'

というかんじなので、このまま演算してくれれば・・・　とおもったが、そこまで簡単でもなかった。いろいろやってみるとHP電卓は、ｎ×ｎのマトリックスとｎ桁の行列みたいな計算はできるけどその逆順の演算や、次元の拡張を伴う計算はできないみたいなのね。んで、生成式はこんな感じに、ハミング行列は，WikipediaによるとH = [ A Im ]で，生成符号がG = [ Ik AT ]なので
HP電卓上で：

[[ 1 0 1 1 ] 1 1 0 1 0 1 1 1 ]　'A'　STO
[ 1 0 1 1 ]　'D'　STO

として。で、'A*D'を評価すると・・・

1:'A*D'　[EVAL]
　　↓
1:[ 3 2 2 ]

ってなる。これがECC (Error Control Code)になる・・・　んだけど、パリティは保存されてるけど2元値になってないので、2元値{0,1}に変換するプログラム（parity）を書いて，上記パラメタで実行する。

1:[ 3 2 2 ] parity
　　↓
1:[ 1 0 0 ]

んであとは，[ D*I₄ D*A ] ってして行列にしてくれればいいんだけど，そもそもHP電卓では単位行列をかけても1次の行列は，

2:[[1 0 0 ]
    [ 0 1 0 ]
        0 0 1]]
1:[ 1 2 3 ]
    *
　　↓
1:[ 1 2 3 ]

てかんじで，そのまま結果の行と列が入れ替わらないなのだね。あと，２つの行列をつなぐとかも組み込み演算子ではできないみたい・・・ そんなわけで，次式

De=[ D*I₄ D*A ]

を計算するプログラムをつくった。

@
@ スタック上の演算結果のマトリックスからパリティを取り出し，
@ 2元値{0, 1}に変換する演算
@
《 →ROW → n 《 1 n START 2 MOD n ROLLD NEXT n →ROW 》》 'parity' STO

@
@ ハミング符号のエラーチェックコード（ECC)の演算
@
《 A D * parity 》 'ECC' STO

@
@ エラーチェックコード（ECC)とデータ部（D)を連結する処理
@
《 D →ROW → n 《 ECC →ROW n + ROW→ 》》 'De' STO

行列の結合はスタック上にいったん分解して，それを組み立てる方式。

で，Deを実行すると，

'De' [EVAL]
        ↓
1:[ 1 0 1 1 1 0 0 ]

という生成されたハミング符号が得られる。こうやって説明すると大変だけど，書いたプログラムは前述の3行だけ。リストやバイナリ演算を使うととてもここまで簡単にはならない。なんでも湘南工科大学ではテスト時に電卓持込み可だけどテスト中にプログラム電卓を片っ端からリセットしてまわる先生がいるらしいが，これだけ短ければ安心だね＾＾。

ハミング符号からシンドロームを求める

さてと，こうして作ったハミング符号の復号はざっと説明しよう。電卓に実装したプログラムは基本的にHの各要素のパラメタ（シンドロームという）を細かくいじるのが目的なのだが，Hの中身はAと単位行列を除きいっしょなので，Aからハミング行列H=[ A I₃ ]を求めて，そこから受信データ Yのシンドローム(S)を演算する。

《 A SIZE 1 GET IDN 》 'Im' STO
《 A →COL → n 《 Im →COL n + COL→ 》》 'H' STO

で，Hを実行すると。一番最初に書いた行列Hが得られて，あとは一ビットの誤りのある受信データYから，シンドロームの計算をしてみる。

[ 1 0 0 0 0 0 0 ] 'E' STO @ エラーベクトルの定義
De E + parity 'Y' STO @ 受信データの定義
《 H Y * parity 》 'S' STO @ シンドロームの計算

で，Sを実行すると。

[ 1 1 0 ] というシンドローム列が得られる。[ 1 1 0 ] というシンドロームは，行列Hの一列目と同じであるため，この値からエラーベクトル（E）を生成することができ，そこから，正しいデータを求めることが可能になる。

符号間のハミング距離を計算してみる

さて。このような1ビット訂正を行うには，ハミング符号の最小符号間距離が3以上なければいけないとのこと。最小符号間距離とは，すべての異なるハミング符号の組み合わせで最小のハミング距離のこと。ハミング距離についてはwikipediaにある。

・・・でまあ，この解説を読んでいると脳がフリーズするわけだが，じっさいにやってみるとそんなに難しくはない。具体的に符号(V1)と符号(V2）の距離を測ってみよう。

[ 1 0 0 0 ] 'D' STO De S @ V1を求める
[ 0 1 0 0 ] 'D' STO De S @ V2を求める
        ↓
2:[ 1 0 0 0 1 1 0 ] @ V1
1:[ 0 1 0 0 0 1 1 ] @ V2

- ABS 2 ^ @ |V1-V2|^2
        ↓
1:4 @ V1-V2のハミング距離は4

で，調べたい符号間距離が出せる。簡単には二つの符号間のエラーベクトル（E)を求めて，そのEの1が立ってる数を数えればいいんだけど。|V1-V2|^2でも大丈夫。理由は… だって距離なんだし（…）。ちなみに，シンドロームのベクトルの１の数を重み(Weight）と呼んで，演習中にいちいち数える必要があったんだけど，これも一個ずつかぞえるプログラムを組むより，こっち使ったほうが速かった（…）。

さて。ハミング符号では，この符号間距離が必ず3以上になるように設計されている。そのとき，最小符号間距離が３である。と呼ぶ。そのため，受信データの1ビットのエラーを訂正することができる。というのは，1ビットエラーの受信データのハミング距離はもともとの符号からは１で，それ以外は2以上になるから。このためエラー前の正しい符号を特定して訂正できるわけだ。

んで，先生のおっしゃるには，「最小符号間距離を5以上にした符号をつくれば，2ビット訂正ができるよ」。というわけで，これが受講後の演習課題になってるわけだが，これが面倒なのよ・・・

そんなわけでエラー訂正にからむ情報理論は，まあ，なかなか面白いなぁと思ったのだが，誤り訂正符号の基礎理論はハミングが考案した1950年から，さっきの奇数重み列の1970年までで，あらかた出尽くしてしまっているそうだ。だから数理科学者になるともうすることが無いらしいし，工学の立場でも研究分野がなかなか見当たらないらしい。

そこで東工大時代の藤原英二氏は，JAXAと多ビット訂正可能な宇宙用メモリーの開発をされたとか。ググると，2件くらい関連特許がありますね。

http://aerospacebiz.jaxa.jp/patent/detail.php?did=1911
http://aerospacebiz.jaxa.jp/patent/detail.php?did=1912

今回の授業でよかったのは，この特許の明細が解読可能になったことかな？
ただし，今回のように「HP電卓で計算してみた」が，できるようになるのは，HP電卓のデータ型に多元ガロワ体が実装されてからになると思う。つーか，数式ソフトで多元ガロワ体の行列サポートしてるやつとか，あるのかな？

--

受講レポートでいつもつけてる自分の放送大学籍番号は，たぶん出席日数不足で単位落としているので省略しました。

あと，あたりまえですが，ハミング符号の演算は実用的な長さ（64ビットとか）になると，バイナリ計算でやるほうがベクトル演算よりもHP電卓のインタープリタでも早くなります。ベクトル演算では演算量が，ハミング符号のサイズそのままに増加するので，データ長が長くなると不利になります。というかベクトル演算が早いのはバイナリ演算では演算結果を読みやすい書式に書き直すコストがかかっていて，その差が出てるのかもしれない。です。

なのでPC上で実用的な誤り訂正符号をシミュレーションする場合は，HP電卓みたいな数式ツールを使うより，VHDLシミュレータをぶん回すほうが早いと思います。ただまあ，論理上の整合性を検討したい場合は，VHDLを使うとコードの整合性の検証に手間をとられることを考えれば，素で定義した行列演算のほうがトータルで速いと思います。それとハード向きの並列演算をソフトで実装した時のコストを体感的に理解するなら，こういうプログラムを書いてみる価値はあると思います。

ちなみに，1ビット誤り訂正と2ビットの誤り検出を行う二つの符号の，ハミングSEC-DED符号と奇数重みSEC‐DED符号を自分がHP電卓に実装したアプリで比較すると，ベクトル演算を用いる場合はハミングSEC-DED符号が早く，バイナリ演算では奇数重みSEC‐DED符号が早くなりました。この差は，パリティ演算にかかるコストと行列演算にかかるコストの並列性の違いが原因だと思います。CPU上では，複数のXOR演算でもレジスタ値なら一発で計算できるので。ただこれも，私の書いた実装ではそうなった・・・ というお話。

--

12/18 読み返したら(笑）がウザかったので，（…）に変更。