計測自動制御学会東北支部 第 213 回研究集会 (2003.12.12)
資料番号 213-11

Octave を用いた並列信号処理アルゴリズム開発環境の構築
Implementation of a Prototyping Environment for Parallel
Signal Processing Algorithms Using Octave○小原武† , 青木孝文† , 樋口龍雄‡○ Takeshi Obara† , Takafumi Aoki† and Tatsuo Higuchi‡†東北大学大学院情報科学研究科, ‡東北工業大学工学部電子工学科†Graduate School of Information Sciences, Tohoku University,‡Department of Electronics, Tohoku Institute of Technology
キーワード : 並列処理 (parallel processing), 信号処理 (signal processing), Octave (Octave), MPI (MPI)
連絡先 : 〒 980-8579 仙台市青葉区荒巻字青葉 05 東北大学大学院情報科学研究科 青木研究室
小原武, Tel.: (022)217-7169， Fax.: (022)263-9308， E-mail: obara@aoki.ecei.tohoku.ac.jp

1.

して並列 Octave を構築してきた2) ．GNU Octave1)

はじめには，ディジタル信号処理の分野で広く用いられる
近年，科学技術計算の分野において要求される演算能力は増加の一途をたどっている．ディジMatlab と互換性を持つフリーソフトウェアである．オリジナルの Octave には並列計算を行う機

タル信号処理の分野においても，膨大な計算量を能が存在しないので，並列信号処理アルゴリズム必要とする信号処理アルゴリズムの実装や大規模開発環境を構築するために，並列計算ライブラリな信号処理システムの迅速なプロトタイピングななどを用いた拡張が必要である．今回作成した並どを行うことが求められており，これらを高速に列 Octave には，最もよく使われている並列計算短時間で処理する手段として，並列処理技術の有効性が注目されている．一方，近年，PC などのライブラリの 1 つである MPI (Message PassingInterface) 3) ライブラリを利用して並列処理用の汎用品を利用した安価な並列処理技術が注目され，関数群を実装した．並列 Octave は，PC クラスタ

COTS (Commercial Oﬀ-the-Shelf) クラスタなどのなどに手軽に導入できるため，並列信号処理アル名称で普及が進んでいる．ディジタル信号処理のゴリズムのプロトタイピングや実証研究などに有分野においても，これらの安価な並列処理技術を用である．
積極的に活用するための高水準なシステム開発環
本稿では，Octave を用いた並列信号処理アル境が求められている．
ゴリズム開発環境の作成と，分散メモリ型並列計
これに対し，筆者らは，並列計算機上で手軽に算機と共有メモリ型並列計算機上で作成した並列利用できる高水準並列ディジタル信号処理環境とOctave を評価した結果について述べる．

–1–

2.

並列 Octave

利用した．

2.1

GNU Octave

2.2

並列計算環境GNU Octave は，もともと， J. B. Rawlings

一般に並列計算を実行するプラットフォーム (並(The Univ. of Wisconsin-Madison) と J. G. Ekerdt列コンピュータ) は，大きく共有メモリ型と分散メ(The Univ. of Texas) らによる化学反応器の設計にモリ型に分けられる．共有メモリ型には，主に 1関する教科書の付属ツールとして作成された．現在

台のマシンに複数の CPU を持つ SMP (Symmet-は J. W. Eaton (The Univ. of Wisconsin-Madison)ric Multi-Processor) マシンの形で提供される．各らによって開発が継続されている．Octave はディ

CPU は，全てのメモリに対して同じようにアクセジタル信号処理の分野で広く利用される Matlab

スすることができるが，CPU の数が多くなると高と互換性を持つ，いわゆる Matlab クローンの価になる．分散メモリ型は，主にネットワークで一種である．Matlab に慣れたユーザならば，各結合された複数のマシンという形で構成される．

種信号処理アルゴリズムの実装やシステムのプロ各マシンは，一般に安価に作ることができるが，トタイピングを手軽に行うことができる．また，他のマシンのメモリ内のデータを参照する際に，Octave は General Public License (GPL) に基づく通信してデータを送受信しなくてはいけないためフリーソフトウェアであるため，ソースコードのオーバーヘッドが生じる．

改変及び再配布を自由に行うことができる．

並列プログラミングをする際に，MPI (Mes-Octave は，C++を基本言語として，Fortran お

sage Passing Interface) ，PVM (Parallel Virtual

よび各種の数値計算ライブラリ，さらには ﬂex やMachine)4) ，HPF (High Performance Fortran)5) ，bison といった字句・構文解析系生成ツールなどをOpenMP6) といった並列計算用のライブラリやプ用いて構築されている．そのソースコードは，オログラミング言語がよく使われている．MPI は，ブジェクト指向言語の抽象化能力を活用した現代あるプロセスから他のプロセスへデータを明示的的な設計であり，注意深く検討されたクラス構造に送信する通信方式であるメッセージパッシングにに基づいて改変・拡張を系統的に行うことができ基づいた並列計算を行うための規格である．PVMるように配慮されている．また，約 340 ページのは，メッセージパッシングを利用して並列計算をマニュアルをはじめとして，インターネットから行うためのソフトウェアである．プロセス間の通得られるさまざまなドキュメント類が存在する．

信以外にも，フォールトトレラント機能などを備Octave は，組み込み関数および実行時に動的えている．HPF は，Fortran を元とした拡張言語にリンクされる関数を定義するためのマクロ (DE-である．データ分散を行う指示文などを基にして，FUN DLD) が用意されている．このマクロを記述コンパイラが並列システムに応じた計算処理のアした C++ソースコードは，oct ファイルと呼ばれロケーションを行う．OpenMP は，主に共有メモるオブジェクトコードにコンパイルされる．作成

リ型並列計算機で用いられ，C や Fortran などでした oct ファイルを Octave が読み込むことで C++書かれた逐次プログラムに対して指示文などを用ソースコードで記述した機能を実行することがでいて並列プログラミングを行うための規格である．

きる．このようにマクロを利用して，新たな関数今回は，MPI ライブラリを利用して Octave にを定義し，システムへ導入することが可能である．

並列インタフェースを実装した．MPI は，さまざ

並列計算用関数の追加に際しては，oct ファイルをまな要求に応えるために，単一ホスト間の同期・–2–非同期や複数ホスト間の集団命令など多数の関数• データ型に関する取り扱いをユーザから隠を定義している．これらの関数群により，統一的蔽するため，通信時に Octave 上のデータ形なインタフェースと柔軟な並列プログラミングを式を char 型のバイトストリームに変換した実現しており，複数のプラットフォーム上でプロ• Octave では引数の参照が行えないため戻りグラムを作成することができる．MPI が提供して値として取り扱ったいる関数群を以下にまとめる．
追加した Octave の並列計算用関数と，元の MPI• 一対一通信関数の一例を示す．

• 集団通信Octave[retval status] = MPI RECV (size, source, tag, comm)• プロセスグループ通信コンテキスト

出力
入力
入力
入力
入力
出力

• プロセストポロジー• Fortran，C からの呼び出し形式• 環境問い合わせ• プロファイリングインタフェース
本稿では，MPI で定義されている一対一通信関数・集団通信関数・環境問い合わせ関数のうち，利用頻度の高い関数群について Octave の並列計算用関数に追加した．

Octave では，整数型・複素数型・文字・2 次元配列などのデータ形式が抽象化されている．データ形式の抽象化により，ユーザは変数の型を意識する必要がないため，迅速なアルゴリズムのプロトタイピングが可能である．Octave に並列計算用関retvalsize
sourcetagcommstatus

受信結果の保存する変数
受信するバイトストリームのサイズ
送信元のランク
メッセージタグ
コミュニケータ
ステータスオブジェクト

MPI
MPI RECV(buf, count, datatype, source, tag, comm, status)
出力
入力
入力
入力
入力
入力
出力

buf
count
datatypesourcetagcommstatus

受信バッファの先頭アドレス
受信バッファ内の要素数
受信バッファの各要素のデータ型
送信元のランク
メッセージタグ
コミュニケータ
ステータスオブジェクト

Octave で並列処理を行うために追加した並列計算用関数を Table 1に示す．これらの関数群により，Octave 上で同期・非同期の通信処理と集団通信処理を行うことができる．

数を追加するにあたって，このデータ形式の抽象化を失わないようにする必要がある．また，並列2.3

Octave から利用できる並列計算用関数群の定義をMPI 関数群の定義となるべく同じようにすることで，ユーザが新たに並列プログラミング技術を習得する労力を軽減することができる．

並列 Octave の構成

Octave は，バッチファイルによる非インタラクティブな処理に加え，インタラクティブな処理も可能である．このインタラクティブな処理の実行によって，スクリプトのデバッグなどの効率が大Octave 上で MPI 関数群を利用するために，並列 Octave から利用できる並列計算用関数群には以下のような変更を加えた．

幅に向上する．このため，並列処理を行う場合もインタラクティブ性を残すことは有用である．イ
ンタラクティブ性を保存するには，ユーザの操作• MPI 関数群のデータ型である datatype (c.f.

している Octave プロセスから他の Octave プロセ

MPI CHAR，MPI DOUBLE など) に関するスを操作する必要がある．そこで，並列 Octave で引数の入力を省略したのプロセス制御には，socket をベースとしたシン–3–

Table 1

clientOctave並列 Octave で利用することができる並列計算用関数一覧
関数名

mpi init

oct file

server1
Octavesocketoct file

IPC ( Named-Pipe)

IPC ( Named-Pipe)server2
Octave

説明

MPI module

MPI module
oct file

並列実行環境の展開

mpi ﬁnalize

並列実行環境の終了

server3
Octavempi abort

並列実行環境の停止

oct file

mpi comm rank

グループ内のランクを返す

mpi comm size

グループのサイズを返す

mpi send

メッセージの同期送信

mpi recv

メッセージの同期受信

mpi irecv

メッセージの非同期受信

mpi test

非同期受信の確認

mpi wait

非同期受信の完了

mpi bcast

グループ内のブロードキャスト

mpi barrier

グループ内のバリア同期

IPC ( Named-Pipe)

MPI module

IPC ( Named-Pipe)

MPI module

Fig. 1

3.

並列 Octave の構成

性能評価
以上のような並列拡張を行った Octave について，PC クラスタシステム (Fig. 2) と PrimePower

600 (Fig. 3) 上で性能評価を行った．評価には，非グルクライアント・マルチプルサーバーモデルを採用した (Fig. 1)．socket によって確立されたセッ線形ダイナミクスの計算と位相限定相関法を用いた指紋照合を使用した．

ションは，サーバ側の Octave を制御するために用いられる．このモデルを採用することにより，ク

3.1

実験環境

ライアントからの指示でサーバに命令を実行させ性能評価を行った実験環境について述べる．実ることができ，インタラクティブ性と非インタラ

験環境には PC クラスタと PrimePower 600 の 2 種

クティブ性を実現することができる．

類の並列計算機を用いた．PC クラスタは分散メモ作成した並列 Octave は，並列計算用関数の呼リ並列計算機であり，サーバー 1 台，ノード 14 台

び出しを行う oct ファイル群とその呼び出しに対しによって構成される．PrimePower 600 は共有メモて通信機能を提供するための実行ファイルであるリ並列計算機であり，8 CPU を搭載している．PCMPI モジュールによって構成される．並列実行環

クラスタシステムの構成と外観をそれぞれ Table

境を展開すると，Octave と MPI モジュールは別々2 と Fig. 2に，PrimePower 600 の構成と外観をそのプロセスとして動作する．これによって，実行れぞれ Table 3と Fig. 3に示す．

環境展開の柔軟性を確保できる．Octave のデータを送受信するために，socket より高速なプロセス

3.2

非線形ダイナミクスの計算の並列化

間通信 (InterProcess Communication: IPC) によ

並列処理時に通信を要する例として非線形ダるセッションを利用する．Octave と MPI モジューイナミクスの計算を取り上げた．本稿では，興奮

ルとの IPC には名前付きパイプ (Named-Pipe) を波を発生する反応拡散モデルである FitzHugh-南用いた．

雲モデル7) を用いた．FitzHugh-南雲モデルは次式–4–

PC クラスタの構成

Table 2ServerNodeCPU

Intel PentiumIII 1GHz DualIntel PentiumIII 1GHz

MemoryPC800 DRDRAM 1GB

PC133 SDRAM 512MB

Network

1000BASE–SX

100BASE–TX

OSLinux 2.4.20

Number of Nodes

114Fig. 3600)
Fig. 2

共有メモリ型並列計算機 (PrimePower

分散メモリ型並列計算機 (PC クラスタ)

Table 3で表される．
τ (∂u/∂t) =∂v/∂t =2 ∇2 u + u(u − α)(1 − u) − v
u − γv(1)

PrimePower600 の構成

CPU

SPARC64 GP 400MHz

MemoryECC SDRAM 2GB

OSSolaris 5.8

Number of CPUs8

ここで，u, v は物質の濃度，α, , γ, τ はパラメータである．FitzHugh-南雲モデルは，平衡点濃度に保たれている空間上のある地点に閾値以上の濃度を与えると，その点から興奮波を発生する (Fig. 4)．
本稿では，興奮波を発生させるためにパラメータとしてα = 10−6 , γ = 10−1 ,  = 10−1 , τ = 10−3 を用いた．また，連続系のの微分方程式をそのままOctave で計算することはできないため．時間と空濃度を参照しなければならない．並列 Octave で計算する場合は，それぞれのステップごとに周囲の濃度を更新してしまうと計算時間全体に占める通信処理の比率が高くなってしまう．そこで，あらかじめ，冗長に計算区間を区切って通信回数を減少させた．

FitzHugh-南雲モデルによる興奮波伝搬を並列間に関して離散化したモデルを使って計算した 8) ．
式 (1) を解くためには，毎回，隣接する地点の化したときの処理手順を Fig. 5に示す．まず，計–5–

(a)Fig. 4

(b)

(c)

(d)

(e)

FitzHugh-南雲モデルによる興奮波伝搬の様子 ((a) から (e) に進む)あった．また，Speed-up factor は，PC クラスタで
1) procedure FitzHugh-南雲モデルによる興奮波の伝搬
2)

14 並列の時に約 7.6 倍，PrimePower 600 で 8 並列

begin3)for ステップ数 := 0 to 1000 doの時に約 6.7 倍であった．PC クラスタで 8 並列の4)begin

時に通信時間が 73 秒であるのに対し，PrimePower

5)

短冊状に分割した 2 次元の計算区間を6)

各プロセスに対して送信する

7)

if 次の計算に濃度情報の更新が必要 then

8)

begin

9)

隣接する区間を持つプロセス間で10)必要な濃度情報を更新する

11)

end

12)

割り当てられた計算区間を計算する

13)
14)
15)

600 では，17 秒である．PrimePower 600 は通信によるオーバーヘッドが小さいため，台数が増加しても，Speed-up factor が飽和せずに台数効果に近い結果が得られていることがわかる．

3.3

end
分割した計算区間を統合する

end.

位相限定相関法を用いた指紋照合の並列化

並列処理時にで通信を要しない例として位相限定相関法 (Phase-Only Correlation: POC) を用

Fig. 5

FitzHugh-南雲モデルによる興奮波伝搬を並列化したときのの処理手順

いた指紋照合を取り上げた．POC は，画像の周波数領域における位相情報を使った高性能画像照合

算区間をプロセス数だけ短冊状に分割し，さらに技術である．

POC を使った指紋照合 9) は，以下のような手冗長に計算する濃度情報を付け加えて各プロセスに渡す．冗長に割り当てた計算区間の計算が終了順で行われる．

すると，各プロセス間で計算に必要な濃度情報を更新するため，隣接する区間を持っているプロセ

1) 登録画像と入力画像を読み込む
2) -40◦ から 40◦ まで登録画像を回転させる

スどうしで通信を行う．最後に各プロセスの計算

3) 回転した登録画像と入力画像を POC で照合結果を結合する．
以上の処理を PC クラスタと PrimePower 600
上で行った結果を Fig. 6と Fig. 7に示す．並列化を行わない場合は，PC クラスタで約 6246 秒，Prime-Power 600 では約 6752 秒であった．処理を並列化し，最もスコアの高かった画像から回転角度を決定する

4) 回転補正した登録画像と入力画像の共通部分を抽出する

した場合は，PC クラスタで 14 並列の時に約 814秒，PrimePower 600 で 8 並列のときに約 1025 秒で–6–5) POC で共通部分を照合する

7000

14Total elapsed time

5000

Communication time121) procedure POC を用いた指紋照合

Speed-up factor

102)

4000

8

3000

6Region size =1024*1024

2000

4

1000

200
12

3

4

5

6

7

89

Speed-up factor

Elapsed time (sec)6000

10 11 12 13 14

Number of processes

begin3)回転角度 := -40 から 40 度まで 1 度刻み;

4)

begin

5)

for 全ての指紋画像ペアに対して do6)

begin

7)

分担する回転角度の範囲を決定する

8)

for 与えられた回転角度に対して do9)

begin

10)

Fig. 6 興奮波伝搬のシミュレーションにおける経過時間と Speed-up factor (PC クラスタ)
7000
Total elapsed time

7Communication time6

Speed-up factor

5000

54000

43000

3Region size =1024*1024
2000

Speed-up factor

Elapsed time (sec)

12)

86000

21000

0
1234

5

6

78

相関スコアを POC で調べる

end

13)

end

14)

各プロセス間で相関スコアを通信し，15)

最も相関が高くなる回転角度を決定する

16)

照合処理を各プロセスに再配分する

17)

for 与えられた処理に対して do

18)

begin

19)

回転した登録画像と入力画像の共通部分を20)抽出した後に POC で照合スコアを算出する

21)

10

回転させた登録画像と入力画像の11)

end22)

end

23)

end.

Number of processes

Fig. 7 興奮波伝搬のシミュレーションにおける経過時間と Speed-up factor (PrimePower 600)

Fig. 8POC を用いた指紋照合を並列化したときの処理手順

POC を使った指紋照合を並列化したときの処

PC クラスタと PrimePower 600 の実験結果を理手順を Fig. 8に示す．POC を用いた指紋照合にそれぞれ Fig. 9 と Fig. 10に示す．並列化を行わな

おいて，以下の部分が実行時間の大部分を占めてい場合は，PC クラスタで約 3691 秒，PrimePowerいる．

600 で約 8373 秒のであったのに対し，並列化した
場合は，PC クラスタで 14 並列の時に約 330 秒，

1) 登録画像の回転

PrimePower 600 で 8 並列の時に約 1150 秒であった．Speed-up factor は，PC クラスタが 14 並列で2) 回転した登録画像と入力画像の照合

11.2 倍，PrimePower 600 が 8 並列で 7.2 倍であっ

3) POC による共通部分の照合た．POC を用いた指紋照合では，それぞれのデー
これらの部分は，お互いのデータに依存せずに処タに対して処理の依存性がないため，PC クラス理が繰り返し行われるため，容易に処理を分割す

タでも PrimePower 600 でも同様にほぼ理想的なることができる．ただし，3) は，2) の計算結果を台数効果が得られた．

利用するため，2) と 3) の間でプロセス全体の同期待ちが加わる．実験は，384 × 384 の 10 枚の指紋画像を用い，そのすべての組み合わせ (45 通り)に対して行った．

–7–

143500

Total elapsed time

3000

Communication time1210

Speed-up factor

2500

82000
61500

Image size =386*386

1000

Speed-up factor

Elapsed time (sec)

4000

42500
001234587000

Communication time76

Speed-up factor

6000

55000

44000

33000

Image size =386*386

2000

Speed-up factor

Elapsed time (sec)Fig.
9
指紋照合処理における経過時間とSpeed-up factor (PC クラスタ)
Total elapsed time

2
11000
00
12345678

Number of processes

Fig.

10

指紋照合処理における経過時間とspeed-up factor (PrimePower 600)

4.

4) PVM : Parallel Virtual Machine
http://www.csm.ornl.gov/pvm/pvm home.html.
5) HPF: The High Performance Fortran Home Page
http://crpc.rice.edu/HPFF/.
6) OpenMP: Simple, Portable, Scalable SMP Programming
http://www.openmp.org/.
7) J. D. Murray, Mathematical Biology,SpringerVerlag, Berlin, 1993.

6 7 8 9 10 11 12 13 14
Number of processes

80003) Message Passing Interface (MPI) Forum
http://www.mpi-forum.org/.

むすび
本稿では，GNU Octave を用いて並列信号処理環境を構築し，その性能評価を行った結果ついて述べた. 開発した並列 Octave は，並列信号処理アルゴリズムを記述・検証・評価することができる高水準な並列プログラミング環境として活用できると考えられる．今後は，並列計算用の関数群の追加や通信性能の向上などを行っていく予定である．
今回 作 成 した 並 列 Octave の ソー ス は以 下 の

URL にて GPL の下に公開している．
http://www.aoki.ecei.tohoku.ac.jp/octave/

参考文献
1) GNU Octave Home Page
http://www.octave.org/.
2) Parallel Octavehttp://www.aoki.ecei.tohoku.ac.jp/octave/.

–8–8) K. Ito, T. Aoki, and T. Higuchi,“Digital reactiondiﬀusion system —A foundation of bio-inspired texture image processing—,” IEICE Transactions on
Fundamentals of Electronics, Communications andComputer Sciences, Vol. E84-A, No. 8, pp. 1909 –1918, August 2001.
9) K. Ito, H. Nakajima, K. Kobayashi, T. Aoki, andT. Higuchi,“A ﬁngerprint matching algorithm using phase-only correlation,” IEICE Transactions on
Fundamentals of Electronics, Communications andComputer Sciences, March 2004,(to be appeared).