FIT2012（第 11 回情報科学技術フォーラム）社団法人

電子情報通信学会

信学技報
TECHNICAL REPORT OF IEICE.

THE INSTITUTE OF ELECTRONICS,
RC-005
INFORMATION AND COMMUNICATION ENGINEERS

ポリフェーズ・フィルタ・バンクを用いた基数 2k FFT に関して
電波望遠鏡用分光器への適用
中原

啓貴†

中西

裕之††

笹尾

勤†††

† 鹿児島大学 大学院 理工学研究科 電気電子工学専攻 〒 890–0065 鹿児島県鹿児島市郡元 1–21–40†† 鹿児島大学 大学院 理工学研究科 物理宇宙専攻 〒 890–0065 鹿児島県鹿児島市郡元 1–21–40††† 九州工業大学 大学院 情報工学研究科 情報創成工学専攻 〒 820–8502 福岡県飯塚市大字川津 680–4
あらまし 電波望遠鏡における分光器 (spectrometer) は宇宙から受信した電波をスペクトルに変換する重要な箇所であり, 広帯域かつ高速なフーリエ変換が必要とされている. 本論文ではポリフェーズ・フィルタ・バンクを用いた基数

2k FFT 回路について述べる. 既存の基数 22 FFT 回路の解析を行い, メモリが実装のボトルネックになることを示す.
提案回路は, ポリフェーズ・フィルタ・バンクを用いて荒い帯域分割を行った後, 基数 2k FFT 回路を用いてスペクトルに変換する. 4 並列化した提案回路を Altera 社 FPGA ボードに実装し, サンプル点 227 の FFT を 0.048 秒で行えることを示す. 提案回路は, サンプル点 227 FFT に関して, SETI spectrometer よりも 20.81 倍高速であった.

On a Radix 2k FFT Using a Polyphase Filter Bank
Application to a Spectrometer for a Radio Telescope
Hiroki NAKAHARA† , Hiroyuki NAKANISHI†† , and Tsutomu SASAO†††† Faculty of Engineering, Kagoshima University, 1-21-40, Korimoto, Kagoshima 890–0065, Japan†† Faculty of Science, Kagoshima University, 1-21-40, Korimoto, Kagoshima 890–0065, Japan††† Department of Computer Science and Electronics, Kyushu Institute of Technology
680–4, Kawazu, Iizuka, Fukuoka, 820–8502 Japan
Abstract In radio telescope, the spectrometer analyzes the radio frequency (RF) received from celestial objects
at the frequency domain by performing the fast Fourier transform (FFT). The FFT can be realized by the radix 2k
FFT (R2k FFT). In Radio Astronomy, the number of points for the FFT is larger than that for the general purpose
one. Thus, the size of memory using a conventional method is too large to implement. We use the polyphase filterbank to reduce the size of the R2k FFT. We implemented the 227 point FFT by the R2k F F T with the polyphasefilter bank. Compared with the SETI spectrometer for 227 -FFT, the four parallelized proposed ones for 227 -FFT is
20.81 times faster.

1. は じ め に
1. 1 電波望遠鏡における分光器 [1], [9], [18]
電波天文学 (Radio astronomy) は, 天体から放射される電波を電波望遠鏡で観測し, 研究を行う天文学の一分野である. 1931
年, ベル研究所のカール・ジャンスキーによる銀河系からの電波 (波長 14.6m) の初観測以降, 電波望遠鏡による宇宙の観測が行われており, 初期宇宙に存在する生まれたての銀河, 星の誕
生, 太陽系のような惑星系の誕生, 宇宙にある物質の進化の解明が進んでいる.
図 1 に電波望遠鏡の構成を示す. 天体からの電波 (RF: Radio
Frequency) を受信アンテナ上のいくつかの鏡によって反射し,受信機へ伝送する. 次に, フィードホーン (Feedhone) によって電磁波を空間モードから導波管モードに変換し, 受信機へと転送する. 受信機は, 増幅器とミクサを用いて扱いやすい
信号 (IF: Intermediate Frequency) に変換する. 最後に, 分光器 (Spectrometer) を用いて受信電圧からパワースペクトラ
ムに変換し, 測定を行う. 従って, 分光器の周波数帯域 (f ) と分解能 (∆f ) が重要なパラメータである. 初期の分光器は各周波数を直接解析し, 同期回路を用いて集計するフィルタ・バ
ンク方式を採用していた. フィルタ・バンク方式では分解能が低いため, レーザー光の回折効果を利用した音響光学偏向素子
による周波数解析を行う音響光学型スペクトル分析器 (AOS:
Acousto-Optical Spectrometer) が提案された. AOS は圧電素

49
( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）Feedhone

Voltage(Analog)
A/Dᄌ឵ེ
Voltage(Digital)

ฃାᯏ

䊐䊷䊥䉣
ᄌ឵ེ

ฃା䉝䊮䊁䊅

ಽశེ
(Spectrometer)

表 1 FFT 解析点数 N の比較.
アプリケーション
一般
用途 [28]ベースバンド3GPP
OFDM
28
CT スキャナ 210 -213
高音質サウンド
29

電波
望遠鏡

OH 輝線
ゼーマン効果SETISpectrometer

䉴䊕䉪䊃䊦
ⶄ⚛⛘ኻ୯
䊡䊆䉾䊃
䉴䊕䉪䊃䊦
ᒝᐲ

N
210 -211䉴䊕䉪䊃䊦ᒝᐲ

220
227

図 1 電波望遠鏡とデジタル分光器.

子等のアナログ素子による非線形特性の調整が困難である. 現行の分光器は, 性能のばらつきが殆どなく調整も不要なデジタ
ル分光器が用いられている. 図 1 にデジタル分光器を示す. デジタル分光器では, まず, 受信機からの中間周波数信号が A/D変換器でサンプルされデジタル信号に変換される. 次に, フーリエ変換 (FFT) を行い, スペクトルを得る. 最後に, 複素絶対値 (Magnitude of a complex number) を求め, スペクトル強度を得る. 近年, 高速 A/D 変換器が開発されており, 高周波の信号を扱えるようになった. 自己相関演算はフーリエ変換と比較して軽い処理であることから, デジタル分光器における処理の重い演算はフーリエ変換である.

1. 2 電波望遠鏡における要求
広帯域解析: 近年, 広帯域周波数解析の有用性が報告されている [19]. 第一の有用性は, 広帯域のスペクトルを得ることにより基本物理パラメータの測定範囲が広がることである. 例えば,輝線観測では星間ガスに存在する色々な種類の原子, 分子, 電離
ガスからの輝線を一度に測定することで, 新しい星間ガス成分の探査, 温度測定が可能となる. 第二の有用性は, 望遠鏡の感度が向上することである. 観測データのノイズは帯域を f とする√と 1/ f に比例する. 従って, f を大きくするほど感度は向上
し, より暗い天体を検出することが可能となる.
表 1 に電波望遠鏡用分光器と一般用途の FFT 解析点数 N の比較を示す. フーリエ変換器を FPGA で実現した場合, 既存の手法では O(log2 N ) の乗算器, LUT, 及び O(N ) の組込みメモリが必要である [30]. 表 1 より N が非常に大きい電波望遠鏡では, 一般用途の帯域にチューニングした FFT では組込みメモリ
がボトルネックとなり, 実現できない. 従って, 組込みメモリをコンパクトにする手法が要求される.
高速処理: 電波望遠鏡は開口合成法を用いて天体から受信した RF から電波図を作成する. 開口合成法は地球の自転を利用
してアンテナの位置を移動させ, 高詳細な図を得る手法であり,スキャナがヘッダをずらしながら図面を得るのと同じ原理である. 従って, 1 データを解析する時間を短くするほど高詳細な図を得ることができる. よって, 電波望遠鏡では高速な FFT が求められる.
省面積: 天体が放射する電波は大気中の水分子や酸素分子によって吸収されるため, 空気が薄く乾燥した標高 4,000m〜
5,000m のハワイ島マウナケア山頂, チリ・アンデス山脈のアタ
カマ砂漠（注 1）等に電波望遠鏡が設置されている. また, 短波より

（注 1）：現在, 国際共同プロジェクトで計画中の ALMA (アタカマ大型ミリ波サブミリ波干渉計) 計画では, チリのアタカマ砂漠にパラボナアンテナ 80 台を設置し, 2012 年から本格的な観測が始まる予定である.

波長が長い電波は電離層で反射されるため, 宇宙空間上の衛星で受信する必要があり, 1997 年に電波天文衛星 HALCA が打ち上げられた [8]. 従って, 電波望遠鏡は高原や宇宙空間に設置するため, 広い敷地を必要とするスーパーコンピュータは使えない. アンテナのみを設置して保存したデータを転送するという手段 (VLBI) も考えられるが, ごく一部のデータしか観測できない.
低消費電力: 高原や宇宙空間では電源や電力送信ケーブルの敷設が限られるため, 低消費電力であることが必須である. 例
えば, NASA の MARVEL (Mars Volcanic Emission and Life)
ミッションでは消費電力を 10 W 以下にすることが要求されている. GPU は価格対性能比が FPGA よりも優れているが [11],消費電力が数 100 W と大きい.
従って, これらの条件を満たす FPGA ボード上に分光器を実装する.

1. 3 関 連 研 究
現在, 国際連携プログラム SKA (Square Kilometer Array)が集光面積 1 平方キロメートルのアンテナを使った電波望遠鏡を開発中である [23]. 地球外生命体探査 (SETI: Search for
Extra-Terrestrial Intelligence) のプロジェクトの 1 つである
SERENDIP プロジェクト [24] では地球外生命体からの信号を解析するため, SETI Spectrometer は f = 200 MHz 帯域の電波を ∆f = 2 Hz の分解能で 1 秒間に解析できるフーリエ変換器を用いている [20].
Cooley と Tukey はフーリエ変換を高速に行う高速フーリエ変換 (FFT: Fast Fourier Transform) を提案した [5].
Duhamel ら は Split-Radix に基 く FFT を 提 案 した [6]. 積
和演算に適した FFT 演算に関して, Linzer らは基数 2,4,8,
及 び Split-Radix の FFT を [16], Karner ら は 基 数 2,3,5 のFFT を [13], Goedecker は基数 2,3,4,5 の FFT を提案した [7].
Kondo らは SETI Spectrometer と同等の分光器を NVIDIA 社
Tesla S1070 (GPU4 台) [26] 上に実装した [14]. 専用回路実装に関しては, He らが基数 22 に基く FFT 回路を提案した [10].
また, FFT の誤差に関する解析や [15], [21], FPGA に実装した
場合のハードウェア量の解析 [30], 及び, FFT の実装でボトル
ネックとなる回転因子のメモリを CORDIC [2], [27] や区分線形近似回路 [12], [22], [25] で実現する手法が提案されている.
1. 4 提案手法と本論文の貢献点
提案手法は実時間で行う広帯域フーリエ変換を FPGA ボー
ド上に実現する. 既存の FFT 回路はメモリ量が実装のボトル
ネックとなるため, ポリフェーズ・フィルタ・バンクを用いて,荒い帯域分割を行った後に各帯域に対して FFT を行う. 提案
手法は分割した帯域に対して FFT を逐次行うため, FFT のサイズを大幅に削減可能である. 本論文の貢献点は以下の通り.
従来 FFT 回路の面積を解析: 基数 2k FFT 回路の面積を解析
し, メモリ量がボトルネックになることを示した.
ポリフェーズ・フィルタ・バンクと基数 2k FFT を組合せ, コ
ンパクトかつ高速な広帯域 FFT を実現: 227 点 FFT を 1 台のFPGA ボード上に実装し, 性能を明らかにした. 筆者の知る限りでは, 1 台の FPGA ボード上に 227 点もの FFT を実装した報告はない. 商用の FFT ライブラリでは高々216 点 FFT までしか実現できない [3]. ポリフェーズ・フィルタ・バンクを用いた FFT が提案されているが, 商用の FFT (R2k FFT) を使っており, FFT の基数を 2k に拡張しポリフェーズ・フィルタ・バ
ンクと組合せたのは本報告が初めてである.
既存手法と比較を行い, 高速かつ低消費電力であることを実証:

50( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）既存の電波望遠鏡用 SETI spectrometer と比較を行い, 提案手法は同一帯域で 20.81 倍高速であることを示した. GPU を 4台用いた GPU spectrometer と比較を行い, 4.37 倍高速かつ約
1 桁低消費電力であることを示した.
第 2 章では基数 2k の高速フーリエ変換器について述べ, 面積の解析を行う. 第 3 章ではポリフェーズ・フィルタ・バンクについて述べる. 第 4 章では実験結果を示し, 第 5 章で本論文のまとめを行う.

Butterfly Operator(BF)

TF

R

R

R

R

BF

-

2. 基数 2k 高速フーリエ変換器 (R2k F F T )

R

R

R

R

R

R

R

R

R

R

R

∑

N −1
nm
, (m = 0, 1, . . . , N − 1)
x(n)WN(1)n=0nmで定義される. WN
は回転因子 (Twiddle Factor) といい
2πnm= e−j N nm = cos(
WN

2π
2π
nm) + j sin( nm)
N
N

(2)

で定義される. j は虚数単位を表す. 本論文では, 入力信号 x(n)は複素数のストリーム入力, N は 2 のべき乗であると仮定する.

2. 2 高速フーリエ変換
Cooley と Tukey らは乗算器を O(N log2 N ) 個に削減する
N 点 高 速 フ ー リ エ 変 換 (Fast Fourier Transform: N FFT) を提案した. 本論文では回路実現が容易な時間間引き (decimation-in-time)（注 2）FFT を用いる.
N 個の入力を偶数番目と奇数番目に分けると, 式 (1) は

∑

N/2−1

X(m) =

∑

N/2−1
2nm
x(2n)WN+

n=0(2n+1)mx(2n + 1)WN

n=0

R

となる.
mn
WN/2
を用いると

∑

N/2−1

n=0

∑

N/2−1
nm
m
x(2n)WN/2
+ WN

00
01
10
11

Twiddle Factor Memory (TF)

BF

TF

図 6 基数 23 バタフライ演算回路 (R23 Butterf ly).

を用いると乗算を共有化できる (図 3 (b)).
入 力 側 か ら 出 力 側 方 向 に ス テ ー ジ の イ ン デック ス を
1, 2, . . . , log2 N とする. このとき, ステージ 1 では隣接した2 点をバタフライ演算器に入力し, 隣接した 2 点の出力を得る (図 2 (a)). ステージ 2 では 2 点離れた 2 点をバタフライ演
算器に入力し, 2 点離れた 2 点の出力を得る (図 2 (b)). ステー
ジ 3 では 4 点離れた 2 点をバタフライ演算器に入力し, 4 点離れた 2 点の出力を得る (図 2 (c)). すなわち, 各ステージ間でタイミングを調整する調整回路が必要となる. 調整回路をシフトレジスタとセレクタで構成する手法が提案されており, バタフ
ライ演算器と組合わせて N -FFT を実現できる [10]. これを基数 2k バタフライ演算回路 (R2k Butterf ly) という.
［例 2.1］ 図 6 に図 2 に示した 8-FFT を実現する基数 23 バタ
フライ演算回路を示す. ただし, バタフライ演算は図 3 (b) を用いている.
基数 2k バタフライ演算回路のハードウェア量を見積もる. 乗
算器数 M ulR2k（注 3）はM ulR2k = 4k = 4log2 N : O(log2 N ),2mn
ここで WN= e−j(2π/N )2mn = e−j(2π/(N/2))mn =

=

0
1

00
01
10
11

Timing Adjuster

2. 1 離散フーリエ変換
x(n), (n = 0, 1, . . . , N − 1) を入力信号とすると N 点離散
フーリエ変換 (以降, N -DFT と表記) はX(m) =

TF
01(2n+1)mx(2n + 1)WN/2

である. 加算器数 AddR2k は各バタフライ演算器の 2 個の複素加算器に 4 個, 1 個の複素乗算器に 2 個必要であるから

AddR2k = (4 + 2)k = 6log2 N : O(log2 N ),である. レジスタ数 RegR2k はステージ i 毎に 2i + 2i − 1 個必要である. ただし, ステージ 1 は不要であるため

n=0m= DF Teven (m) + WN
DF Todd (m),

− 1 である. 式 (1) は N -DFTである. ここで m = 0, 1, . . . , N
2-DFTに分解できることを示している. N -DFT にが 2 個の N
2対して上式を 1-DFT になるまで s = log2 N 回再帰的に適用することで N -FFT が得られる. このとき, s をステージ数と定義する. N -FFT では 1 ステージ毎に N個の乗算器が必要であ2
る. N -FFT は s ステージで構成されるため, O(N log2 N ) 個の乗算器が必要である. N -DFT では乗算器数が O(N 2 ) 個であったので, N -FFT を用いることでフーリエ変換をコンパクトに実現できる.
2. 3 基数 2k バタフライ演算器 (R2k Butterf ly )N -FFT ではバタフライ演算が基本演算単位となる. 図 3 にバタフライ演算を示す. バタフライ演算器は複素加算回路 (図 4),及び複素乗算回路 (図 5) から構成される. また, 回転因子の値は回転因子メモリで保持する. 図 3 (a) に示すバタフライ演算では,
m+N/2
m
= −WN
乗算が 2 回必要であるが, 回転因子の対称性 WNRegR2k =k∑(2i + 2i − 1) : O(2k ),i=2である. セレクタに用いる 2 マルチプレクサ数 M uxR2k はス
テージ i 毎に 2i 個の信号を選択するセレクタが 2 個必要であるから

M uxR2k = 2 ×

k∑(2i − 1) : O(2k ),i=2である. 回転因子メモリ量 T widdleM emR2k は

T widdleM emR2k = N s : O(N log2 N )である. 従って, N -FFT を基数 2k バタフライ演算器で実現した場合, O(2k ) 個必要なレジスタと 2 マルチプレクサが実装のボトルネックとなる.
（注 3）：本論文では FPGA の乗算器ブロック 1 個で乗算を実現するため, 乗算（注 2）：周波数間引き (decimation-in-frequency)FFT も存在する.

回数と乗算器数を同一と見なして見積る.

51( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）Stage 1

Stage 2

Stage 3

x
X(0)

x(0)
x(4)
x(2)
x(6)
x(1)
x(5)

-1

-1

W84

-1

X(2)

-1

X(3)
-1

W

W80
W81

-1

-1 W82-1W84-1 W83
(b)

-1

-1
0
8

x(3)
x(7)

X(1)

W80

(a)

X = x + yWNm

WNm

(a)

y
W

X(4)

m+ N / 2N

Y = x + yW

= −W

m
N

X = x + yWNm

x

X(5)
-1
X(6)

(b)

X(7)

m+ N / 2
N

WNm

y

-1

(c)

a

a
b

+

a+c

b

㽢

c
d

+

b+d c

㽢

(a+bj)+(c+dj)

Y = x − yWNm

図4

複素加算器.

Transpose-1d+ (bc+ad)j㽢(a+bj)㽢(c+dj)

図5

複素乗算器.

N/q points

W80

x(2)

-1

-1

W

-1

N/q FFT

X(1)N-points
input signal

X(2)
-1

-1

2
8

x(1)

1
8

X(3)

Filter 1

W

N/q FFT

X(4)Filter 2
N points

X(5)

x(5)
-1

W82

x(3)
x(7)

ac-bd

X(0)

x(0)

x(6)

-

図 3 バタフライ演算器.

図 2 8-FFT のデータフロー.

x(4)

㽢

-1

-1

W82

-1

図7

-1
X(6)
-1

X(7)

W83
転 置 の 例.

(a)Transpose Memory
Radix 2k
Butterfly
Unit

N/q FFT

Filter q

Radix 2k
Butterfly
Unit

(b)

(c)

(d)

図 9 フィルタ・バンクを用いた R2k FFT.

Radix 2k
Butterfly
Unit

セレクタに用いる 2 マルチプレクサ数 M uxR2k F F T は

図 8 基数 2k FFT.

M uxR2k F F T = qk∑(2i − 1) : O(2k ),

(4)

i=2

基数 2k バタフライ演算回路のレジスタとセレクタが増加する原因は, ステージ毎に異なる距離の 2 点を計算するからである. 従って, ステージ間に各点のインデックスを置き換える回
路を挿入することで, 隣接した 2 点のバタフライ演算で実現でき, 調整回路を削減することができる. インデックスを置き換えることを転置 (Transpose) といい, 転置メモリで実現する.
［例 2.2］ 図 7 に図 2 に示したデータフローのステージ 2 とス
テージ 3 の間に転置メモリを挿入した例を示す. 転置を行うことで, ステージ 3 の演算は隣接した 2 点のバタフライ演算で計算でき, タイミング調整回路は不要である.
s = dlog2 N e をステージ数とする. 図 8 に, 基数 2k バタフラ
イ演算回路間に転置メモリを挿入した基数 2k FFT (R2k F F T )を示す. 基数 2k FFT 回路では, k ステージ毎に転置メモリを配置するため, ks 個の転置メモリが必要である. 従って, k を大きくすれば転置メモリ数を削減できる.
ステージ数 s = log2 N の N -FFT は基数 2k バタフライ演算器を q = ks 個用いて実現できる. 基数 2k 転置 FFT 回路の乗算器数 M ulR2k F F T は

M ulR2k F F T = 4s = 4log2 N : O(log2 N ),
加算器数 AddR2k F F T は

AddR2k F F T = 6s = 6log2 N : O(log2 N ),レジスタ数 RegR2k F F T はRegR2k F F T = qk∑(2i + 2i − 1) : O(2k ),

i=2

(3)

であり, 転置メモリ量 T ransM emR2k F F T は

T ranM emR2k F F T = N q : O(N )(5)である. N -FFT のステージ数 s は一定であるから, 転置メモリ
数を削減するために k を大きくすると, 式 (3), (4) が示すようにレジスタと 2 マルチプレクサが増加してしまう. 従って, 与えられたハードウェアの制約を満たす適切な k を決めなければならない. 本論文では, k を実験的に求める. 回転因子メモリ量
T widdleM emR2k は

T widdleM emR2k F F T = N s : O(N log2 N )

(6)

である.
以上の議論より R2k F F T において, 実装のボトルネックとなるのは N  k より回転因子メモリ量である. 従って, 本論文
ではポリフェーズ・フィルタ・バンクを用いて N そのものを削減することを提案する.

3. ポ リ フェー ズ・フィル タ・バ ン ク を 用 い た
R2k F F T
3. 1 フィルタ・バンクを用いた R2k F F T の削減
R2k F F T の解析結果より, 面積は解析点数 N に比例して増加する. 従って, N を削減することで面積を削減できる. 図 9にフィルタ・バンクを用いた R2k F F T を示す. 周波数解析を行う帯域を ∆f とし, 解析点数を N とする. まず, q 個のフィルタを用いて荒い帯域分割を行う (図 9 (a)). 各フィルタにより, 帯
域は ∆f /q に分割される. 次に, ダウン・サンプリングを行い,

52( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）
t=8

x

Z-1

h0

h1

䉺䉾䊒ᢙ:

Z-1

Z-1

Z-1

Z-1

h2

h3

h4

h5

Z-1
h6

Z-1
Z-3
h7

㰿

Z-1

Z-3
h5

h1
㸣M

Z-2

y

Z-3

Z-3

3. 2 ポリフェーズ・フィルタ [17]ここではフィルタ・バンクを構成するポリフェーズ・フィル
タについて述べる. 間引き回路はフィルタとダウン・サンプリ
ングで構成される. デジタルフィルタでは FIR (有限インパルス応答) フィルタと IIR (無限インパルス応答) フィルタの 2 種
類に大別される. IIR フィルタはフィードバックループを用いるため, インパルス応答は再帰的となり無限に持続する. 一方,
FIR フィルタはフィードバックループを用いないため, インパ
ルス応答は有限となる. 従って, 同数のタップ数に関しては IIR
フィルタの方が FIR フィルタよりも高精度なフィルタを構成できる. しかし, FIR フィルタは全ての周波数範囲において群
遅延が完全に一定であることや安定性が高いことで, IIR フィ
ルタよりも優れている. 電波望遠鏡では長期間の測定が行われ,安定した出力が求められることから本論文では FIR フィルタを採用する. 図 10 (a) に FIR フィルタを示す. 時刻 t における入力を x(t) とすると, FIR フィルタでは遅延要素の z 変換に相当するレジスタ z −1 を通過し, タップ係数を乗算した値全てを加算する. 係数の値 hi とタップ数 t を変えることにより, 任意の周波数応答特性を実現することができる.
入力信号に対して直接ダウン・サンプリングを行った場合,
高周波成分による影響 (エイリアス) を受け, 余分な成分が含まれてしまうことがある. 従って, フィルタリングによる高周波成分の除去が必要となる. 図 10 に FIR フィルタを用いた間引き回路を示す. FIR フィルタに M 点ダウン・サンプリング回路 (図 10 (b)) を付加することで間引き回路を実現できる. M
点ダウン・サンプリング回路は M 点サンプル毎に M − 1 点サ
ンプルを捨てる.
図 10 に示した間引き回路は図 11 に示すサブフィルタ回路に分割できる. 図 11 に示すサブフィルタも FIR フィルタである. ダウン・サンプリングを行ってからフィルタリングを行っても出力は変わらないので, M 点ダウン・サンプリング回路を入力側に移動する (図 12). 図 12 において, 各サブフィルタへ入力を分配している回路 (図 12 灰色部分) を回転スイッチで置き換えることで図 13 に示すポリフェーズ・フィルタが得られる. 元の FIR フィルタのタップ数を t とすると, ポリフェーズ・
フィルタは M 分割したタップ数 t/M のサブフィルタを時分割実行していると言える. 従って, 図 14 に示す各サブフィルタを逐次模擬する時分割方式ポリフェーズ・フィルタ (以降, 単に䉰䊑䊐䉞䊦䉺ᢙ:
M

㸣M

Z-3
h7

h3

N 点のデータを N/q 点に削減する (図 9 (b)). フィルタリングとダウン・サンプリングを行うことを間引き (Decimation)という. 次に, N/q-FFT を用いてフィルタリングした帯域のみ
周波数変換を行う (図 9 (c)). FFT 演算を行う前にダウン・サ
ンプリングを行っているので周波数領域に変換された N/q 点のデータを周波数シフトし元の帯域に戻す (図 9 (d)). これらの操作を各帯域毎に行って N 点スペクトラムを構成する.

㰿

h6

h2

図 10 FIR フィルタを用いた間引き回路.

䉰䊑䊐䉞䊦䉺

h4h0(a) FIR䊐䉞䊦䉺(b) Mὐ䉻䉡䊮䊶䉰䊮䊒䊥䊮䉫

tM

図 11 サブフィルタ分割.
Z-3

㸣M

h4

h0
㸣M

Z-1

Z-3
h5

h1
㸣M

Z-2

Z-3

㸣M

Z-3

㰿

h6

h2
Z-3

h7

h3

図 12 ダウン・サンプリング回路の移動.
Z-3
h4

h0

࿁ォ䉴䉟䉾䉼
(Rotating Switch)

Z-3
h5

h1
Z-3

㰿

h6

h2
Z-3
h3

図 13

h7

ポリフェーズ・フィルタ.

౉ജ䊋䉾䊐䉜
䊜䊝䊥

x

㰿

図 14

h0

h4

h1

h5

h2

h6

h3

h7

ଥᢙ䊜䊝䊥

䊧䉳䉴䉺

y

reset

時分割方式ポリフェーズ・フィルタ.

ポリフェーズ・フィルタと表記) で実現できる. この回路では,入力バッファメモリで遅延要素 z −1 による過去の入力を保持す
る. 係数メモリに各時刻毎の係数を格納する. そして, 乗算器と加算器を用いて各時刻毎のサブフィルタの出力を求める. レジ
スタに計算結果を格納する.
ポリフェーズ・フィルタのハードウェア量を見積もる. 元のFIR フィルタのタップ数を t とし M 点ダウン・サンプリングを行うとする. 乗算器数 M ulP oly は53( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）䊘䊥䊐䉢䊷䉵䊶
䊐䉞䊦䉺䊶䊋䊮䉪

䊋䉾䊐䉜䊜䊝䊥

Polyphase Filter 1

y1(0) y1(1)Polyphase Filter 2

y2(0) y2(1)Polyphase Filter qyq(0) yq(1)䊥䊒䊦
0.8dBએਅ

N/q point
R2kFFTᏪၞᄖ
ᷫ⴮㊂32dBએਅ๟ᵄᢙ䉲䊐䊃

図 15

ポリフェーズ・フィルタ・バンクを用いた R2k FFT.

M ulP oly =

t
M

ㅢㆊၞ

ㆫ⒖ၞ
(7)

図 16

ポリフェーズ・フィルタの振幅特性.

ㆫ⒖ၞ
Tであり, 加算器 AddP oly は

AddP oly = 1

㒖ᱛၞ

ㅢㆊၞ
P

(8)

であり, レジスタ数 RegP oly は

RegP oly = 1

3dBᷫ⴮
䈚䈩੤Ꮕ

(9)

であり, 入力バッファメモリ M emBufP oly は

M emBufP oly = t

ᏪၞT+PᲤ䈮FFT

(10)
図 17

であり, 係数メモリ M emCoefP oly は

M emCoefP oly = t

(11)

である.

3. 3 ポリフェーズ・フィルタ・バンクを用いた R2k FFT
図 15 にポリフェーズ・フィルタ・バンクを用いた R2k F F Tを示す. 入力信号は q 個のポリフェーズ・フィルタからなるポリフェーズ・フィルタ・バンクにより q 個の帯域に分割される.
転置メモリに各フィルタからの出力 yi (t), i = 1, 2, . . . q が各ア
ドレス (列) に格納される. R2k F F T は転置メモリの行に格納された各フィルタの出力を逐次読み出し, N/q-FFT を行う. 最後に, 元の帯域に戻すために周波数シフトを行う. 周波数シフト
回路は FFT の出力をシフトするだけなので, R2k F F T やポリ
フェーズ・フィルタ・バンクと比較して無視できるほど小さい.
ポリフェーズ・フィルタ・バンクのハードウェア量を見積もる. 元の FIR フィルタのタップ数を t とし M 点ダウン・サン
Nとプリングを行うとする. ポリフェーズ・フィルタ数を q = M
すると乗算器数 M ulP F B は式 (7) より
M ulP F B =t
q : O(q)M(12)

であり, 加算器 AddP F B は式 (8) より

AddP F B = q : O(q)(13)

であり, レジスタ数 RegP F B は式 (9) より

RegP F B = q : O(q)(14)

であり, 入力バッファメモリ M emBufP F B は式 (10) より

M emBufP F B = qt : O(q)(15)

であり, 係数メモリ M emCoefP F B は式 (11) より

M emCoefP F B = qt : O(q)(16)

である. 従って, ポリフェーズ・フィルタ数 q に対してハード

ポリフェーズ・フィルタ・バンクの振幅特性.

ウェア量は線形的に増加する. フィルタの性能を示すタップ数
t はなるべく小さくしたいがフィルタの精度に影響するので, 本
論文では実験的に求める.
各フィルタの出力は独立なので, 複数のフィルタ毎に転置メ
モリと R2k F F T を用意することで容易に並列化可能である.
従って, FPGA の空きリソースを利用して高速な分光器を構成できる.

4. 実 験 結 果
4. 1 実 装 環 境
基数 2k 転置 FFT 回路を Altera 社 DE4 Development and
Education Board に実装した. 搭載 FPGA は Altera 社 Stratix
IV GX(ALUT 数: 424,960 個, レジスタ数: 424,960 個, M9k(9k
ビット組込みメモリ) 数: 1,280 個, 18 × 18DSP 数: 1,024 個)であり, DDR2SO-DIMM コネクタを 2 個, SSRAM を 1 個,
PCIExpress (x8) コネクタを 1 個, 及び A/D 変換器と接続するための HSMC コネクタを 2 個持つ. FPGA 合成ツールはAltera 社 Quartus II v.11.0 を用いた. また, DDR2SDRAM
コントローラは Altera 社の MegaCore Function を用いて生成した.
式 (3), (4) より, レジスタ数と MUX 数は O(2k ) で増加する.
各ステージの基数 2k がなるべく等しくなるように, k を求めた. 実装に用いた FPGA ボードは DDR2SO-DIMM を 2 個実装可能なので, 3 つの基数 2k バタフライ演算回路間に転置メモ
リとして用いた. 最終段の転置メモリは実装していない. 従って, k = d log32 N e とした.
4. 2 ポリフェーズ・フィルタ・バンクの設計
ポリフェーズ・フィルタの元となる FIR フィルタの係数とタップ数を実験的に求めた. 図 16 にポリフェーズ・フィルタの振幅特性を示す. 衛星搭載用フィルタ・バンクの設計事例 [29] を引用し, 通過域のリプルを 0.8 dB 以下, 帯域外減衰量を 32 dB
以下とした. また, 設計の自由度を考慮して窓関数は Kaiser 窓54( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）㪉㪌㪇㪇㪉㪅㪎㪋㪇㪇

㪇㪅㪎㪍㪉㪊

㪇㪅㪌㪇㪇㪇

㪌㪇㪇㪇㪅㪈㪇㪌㪉
㪇㪅㪇㪇㪇㪇

図 18

䊘䊥䊐䉢䊷䉵䊶䊐䉞䊦䉺䊶䊋䊮䉪㪩2k㪝㪝㪫

㪊㪇㪇㪇㪇

㪈㪋㪇㪇
㪈㪉㪇㪇
㪛㪪㪧䊑䊨䉾䉪ᢙ

㪉㪌㪇㪇㪇
㪉㪇㪇㪇㪇
㪈㪌㪇㪇㪇
㪈㪇㪇㪇㪇
㪌㪇㪇㪇

数.

㪊㪉
㪍㪋
㪈㪉㪏
㪉㪌㪍

䉺䉾䊒ᢙ

㪈㪇㪇㪇
㪏㪇㪇
㪍㪇㪇
㪋㪇㪇
㪉㪇㪇

㪈㪇㪉㪋

図 21 DSPを選択した. 図 17 にポリフェーズ・フィルタ・バンクの振幅特性を示す. この特性も衛星搭載用フィルタ・バンクの設計事例を引用した. タップ数 t を固定し, 各フィルタの遷移域が 3 dB減衰する点が交差するように各フィルタバンクの係数を求めた.
従って, 提案回路は帯域 T + P で Nq -FFT を行う.
ポリフェーズ・フィルタ・バンクの解析より, タップ数 t はなるべく小さい値が望ましいがタップ数を小さくするとフィルタの遷移域が広がり, 他の帯域の成分が混入するため FFT の誤差が増えてしまう. 図 18 に直接 FFT を行った結果に対して,
ポリフェーズ・フィルタ・バンクを行い FFT を行った結果との平均誤差を求めた結果を示す. ただし, この結果は Altera 社
DSP ブロック (丸め誤差回路を含む) を用いた結果であり, 実
装に依存することに注意されたい. 図 18 より, タップ数 t を増やすことで平均誤差を削減できる.

4. 3 ポリフェーズ・フィルタ数 q とハードウェア量の関係
ポリフェーズ・フィルタ数 q とタップ数 t を変化させ, FPGAのハードウェア量を測定した. 図 19 に ALUT 数を示す. 元の解析点数を N とする. ポリフェーズ・フィルタ数 q に対して提案回路の FFT の解析点数は N/q となるので, R2k で消費する
ALUT 数は線形的に減少する. しかし, ポリフェーズ・フィル
タで加算器を用いるため, ポリフェーズ・フィルタ・バンク数 qが指数関数的に増えると加算器が消費する ALUT 数も指数関数的に増加する. 従って, ALUT 数に関してトレード・オフがあり図 19 より, q =64〜256 でバランスが取れていることがわ
かる. 図 20 に組込みメモリ (M9K) 数を示す. タップ数 t を増やすとポリフェーズ・フィルタで間引きできるサンプル数 M が

表2

㪈㪇㪉㪋

䊘䊥䊐䉢䊷䉵䊶䊐䉞䊦䉺ᢙ q

図 19 ALUT 数.

㪌㪈㪉

㪉㪌㪍

㪈㪉㪏

㪍㪋

㪏

㪉

㪇
㪊㪉

䊘䊥䊐䉢䊷䉵䊶䊐䉞䊦䉺ᢙ q

㪌㪈㪉

㪉㪌㪍

㪈㪉㪏

㪍㪋

㪊㪉

㪈㪍

㪏

㪋

㪉

㪇

㪋

㪘㪣㪬㪫ᢙ

㪈㪇㪉㪋

䊘䊥䊐䉢䊷䉵䊶䊐䉞䊦䉺ᢙ q

提案手法の誤差解析.

図 20 M9K

㪊㪌㪇㪇㪇

㪌㪈㪉

䉺䉾䊒ᢙ t

㪉㪌㪍

㪇

㪉㪌㪍

㪈㪉㪏

㪈㪉㪏

㪍㪋

㪍㪋

㪊㪉

㪊㪉

㪉

㪈㪍

㪈㪇㪇㪇

㪈㪍

㪈㪅㪇㪇㪇㪇

㪊㪉
㪍㪋
㪈㪉㪏
㪉㪌㪍

㪈㪌㪇㪇

㪏

㪈㪅㪊㪍㪌㪉

㪈㪅㪌㪇㪇㪇

㪋㪉㪅㪇㪇㪇㪇䉺䉾䊒ᢙ㪉㪇㪇㪇㪈㪅㪐㪎㪋㪇

㪤㪐㪢ᢙᐔဋ⺋Ꮕ 㪲㩼㪴

㪉㪅㪌㪇㪇㪇

㪈㪍㪊㪅㪇㪇㪇㪇

数.

実 装 結 果.

ALUT 数 M9k 数

DSP
ブロック数

ポリフェーズ・フィルタ・バンク

2048

256

128

2k FFT(4 台)

80,9360

320

DDR2 SDRAM

6,036

6

0

89,020262

448

コントローラ (2 台)合計増えるため, M9K を削減できる. しかし, タップ数 t 以上は間引きできないのでこれ以上は M9K 数が増える. 従って, M9K
に関してはフィルタ数 q を増やしたとしても常に削減できるとは限らない. 図 21 に DSP ブロック (乗算器+加算器) 数を示
す. M9K と同様に, 間引きできるサンプル数の限界によりタッ
プ数 t までしか削減できない. ただし, M ９ K, DSP ブロックの実験結果は FPGA の構成要素の変更により変わる値であることに注意されたい.

4. 4 設計した分光器
SETI spectrometer と同じ 227 FFT を行う分光器を設計した.
図 19 より, FPGA 実装のボトルネックとなる ALUT 数のバラ
ンスが取れているタップ数 t = 128, ポリフェーズ・フィルタ数
q = 128 とし, R2k F F T を 4 並列に実装し, 高速化を行った. 実
装結果を表 2 に示す. また, FPGA 以外に DDR2-SODIMM を2 個, SSRAM を 1 個を使用した. 最大動作周波数は 375.2 MHzであったが, DDR2SO-DIMM を 667 MHz で動作させるので,
PLL を用いてシステムクロックを 333 MHz とした. 提案回路55
( 第 1 分冊 )Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.

FIT2012（第 11 回情報科学技術フォーラム）は 1 クロックで 2 個のデータを処理できるので, 333 MHz 動作
1
させる場合は 333×2
20 秒で 2 個のデータを処理できる. 提案回
路は 4 並列化しているので, N = 227 サンプルを FFT する時
227
1
1
間は 333×2
20 × 2 × 4 = 0.048 [sec] である.
FFT 演算では N 点のデータに対して N 回加算が行われる
ため, 最悪 log2 N ビット桁上げが生じてしまう. 従って入力信号の精度を r ビット (符号なし) とすると, 最悪の場合, FFT の
出力は r + log2 N + 1 (符号あり) ビットとなる [7], [21]. 提案
回路は A/D 変換器からの入力を 8 ビットとし, 512 個のポリフェーズ・フィルタでフィルタリングした 218 点のデータに対して FFT 演算を逐次行い, 227 点のスペクトラムを得る. よって,最悪の場合の出力は 8 + log2 218 + 1 = 27 ビットとなる. 実装
した回路は, 内部で 18 ビットの実数と 18 ビットの虚数に拡張
し FFT を行っている. R2k F F T は 3 ステージ構成で実装したので, 各ステージの出力で３ビットシフトしスケーリング調整を行っている. また, 誤差を減らすため, 乗算器を実装した DSP
ブロック出力に丸め回路を用いている. SETI Spectrometer もFFT 内部で 36 ビットの複素数 (実部 18 ビット, 虚部 18 ビット) を用いており, スケーリングを行っているので実用上問題はない.

4. 5 他の手法との比較
ポリフェーズ・フィルタ・バンクと R22 F F T を用いた SETI
spectrometer が提案されている [20]. SETI spectrometer は
227 FFT 演算を 1.00 秒で行うため, 提案手法は 20.81 倍高速である. また, SETI spectrometer は R22 F F T を用いるため, 5台の FPGA を必要とする. 一方, 提案回路は R2k F F T を用いるため FFT 回路部をコンパクトに実現でき, 空きリソースを並列化に適用できる.

[5]

[6]
[7]

[8]
[9]
[10]
[11]

[12]

[13]
[14]

[15]

[16]

[17]

5. ま と め[18]

本論文ではポリフェーズ・フィルタ・バンクを用いた R2k FFT
回路について述べた. 既存の基数 22 FFT 回路の解析を行い, メ
モリが実装のボトルネックになることを示した. ポリフェーズ・
フィルタ・バンクを用いて荒い帯域分割を行った後, R2k FFT回路を用いてスペクトルに変換する. 4 並列化した提案回路をAltera 社 FPGA ボードに実装し, サンプル点 227 の FFT 演
算を 0.048 秒で行えることを示した. 提案回路は, サンプル点
227 FFT 演算に関して, SETI spectrometer よりも 20.81 倍高速であった.
本論文で実装した回路のボトルネックは外付けメモリの帯域である. メモリ数を増加し, 帯域を増やすことで更なる高速化が可能である. 今後の課題は, 電波望遠鏡に適したメモリ帯域が広い FPGA ボードを開発し, その性能を明らかにすることである. また, 本研究の成果による FFT ライブラリは電波望遠
鏡国際開発コミュニティである CASPER [4] で公開する予定で
ある.

6. 謝

[19]

[20]

[21]
[22]

[23]
[24]
[25]

[26]
[27]

辞

本研究は, 一部, 日本学術振興会・科学研究費補助金, 及び日本学術振興会・頭脳循環プロジェクト補助金による.
文
献

[28]
[29]

赤羽 賢司, 海部 宣男, 田原 博人, “宇宙電波天文学,” 共立出版株式会社,
1988.
R. Andraka, “A survey of CORDIC algorithms for FPGA based
computers,” FPGA1998, pp. 191-200, 1998.
Altera Corp., “FFT MegaCore Function v11.0”, 2011.
CASPER: Collaboration for Astronomy Signal Processing and

[30]

[1]
[2]
[3]
[4]

56
( 第 1 分冊 )Electronics Research,https://casper.berkeley.edu/
J. W. Cooley and J. W. Tukey, “An algorithm for the machine
computation of complex fourier series,” Mathematics of Computation, Vol. 19, pp. 297-301, 1965.
P. Duhamel and H. Hollmann, “Split-radix FFT algorithm,”
Electron. Lett., Vol. 20, pp. 14-16, 1984.
S. Goedecker, “Fast radix 2,3,4 and 5 kernels for fast fouriertransformations on computers with overlapping multiply-add
instructions,” SIAM Journal Sci. Comput., Vol. 18, pp. 16051611, 1997.
HALCA: Highly Advanced Laboratory for Communications andAstronomy, http://www.vsop.isas.jaxa.jp/top.html
半田 利弘, “基礎からわかる天文学: 太陽系から銀河, 観測技術や宇宙論まで,” 誠文堂新光社, 2011.
S. He and M. Torkelson, “A new approach to pipeline FFT processor,” IPPS1996, pp. 766-770, 1996.
D. H. Jones, A. Powell, C. S. Bouganis, P. Y. K. Cheung, “GPUVersus FPGA for High Productivity Computing,” FPL2010,
pp. 119-124, 2010.
V. K. Jain, S. A. Wadekar, and L. Lin, “A universal nonlinear
component and its application to WSI,” IEEE Trans. Components, Hybrids, and Manufacturing Technology, Vol. 16, No. 7,
pp. 656-664, Nov., 1993.
H. Karner, et al., “Multiply-add optimized FFT kernels,” Math.
Models and Methods in Appl. Sci., Vol. 11, pp. 105-117, 2001.
H. Kondo, E. Heien, M. Okita, D. Werthimer, K. Hagihara, “Amulti-GPU spectrometer system for real-time wide bandwidth
radio signal analysis,” ISPA2010, pp. 594-604, 2010.
W. R. Knight and R. Kaiser, “A Simple Fixed-Point Error
Bound for the Fast Fourier Transform,” IEEE Trans. Acoustics, Speech and Signal Proc., Vol. 27, No. 6, pp. 615-620, 1979.
E. N. Linzer and E. Feig, “Implementation of eﬃcient FFT algorithms on fused multiply-add architectures,” IEEE Trans.
Signal Processing, Vol. 41, pp. 93-107, 1993.
R. G. Lyons, “Understanding Digital Signal Processing,” Pearson Education, 1996.
中井直正, 坪井昌人, 福井康雄, “宇宙の観測 II(16): 電波天文学”, 日本
評論者, 2009.
A. Hirota, N. Kuno, N. Sato, H. Nakanishi, T. Tosaki, andK. Sorai, “Variation of molecular gas properties across the spiralarms in IC 342: Large-scale 13CO (1-0) emission”, PASJ2010,
No. 62, pp. 1261-1275, 2010.
A. Parsons et.al, “PetaOp/Second FPGA signal processing for
SETI and radio astronomy,” Proc. of the Asilomar Conferenceon Signals, Systems, and Computers, 2006.
L. R. Rabiner and B. Gold, “Theory and Application of Digital
Signal Processing,” Prentice-Hall Inc., 1975.
T. Sasao, S. Nagayama and J. T. Butler, “Numerical function
generators using LUT cascades,” IEEE Trans. on Comput.,
Vol.56, No.6, June 2007, pp.826-838.
SKA: Square Kilometre Array,http://www.skatelescope.org/
SERENDIP: The Search for extra terrestrial intelligence at UC
Berkeley, http://seti.berkeley.edu/SERENDIP/
M. J. Schulte and J. E. Stine, “Approximating elementary functions with symmetric bipartite tables,” IEEE Trans. Comput.,Vol. 48, No. 8, pp. 842-847, Aug., 1999.
NVIDIA Corp., “TESLA S1070 GPU COMPUTING SYSTEM,” http://www.nvidia.com/
J. E. Volder, “The CORDIC trigonometric computing technique,” IRE Trans. on Electronic Comput., pp. 330-334, 1959.
Xilinx Inc., “LogiCORE IP fast fourier transform v7.1”, 2011.
山下 史洋, 風間 宏志, 中須賀 好典, “衛星搭載用帯域幅可変 FFT フィ
ルタバンクの提案と基本動作特性,” 電子情報通信学会論文誌 B, 通信
J85-B(12), pp. 2290-2299, 2002.
B. Zhou, Y. Peng, and D. Hwang, “Pipeline FFT architectures
optimized for FPGAs,” Int’l Journal of Reconfigurable Computing, Vol. 2009, 2009.

Copyright © 2012 by
The Instiute of Electronics, Information and Communication Engineers andInformation Processing Society of Japan All rights reserved.