JP3709291B2 - 高速複素フーリエ変換方法及び装置 - Google Patents

高速複素フーリエ変換方法及び装置 Download PDF

Info

Publication number
JP3709291B2
JP3709291B2 JP29173498A JP29173498A JP3709291B2 JP 3709291 B2 JP3709291 B2 JP 3709291B2 JP 29173498 A JP29173498 A JP 29173498A JP 29173498 A JP29173498 A JP 29173498A JP 3709291 B2 JP3709291 B2 JP 3709291B2
Authority
JP
Japan
Prior art keywords
data
butterfly
point
same
register
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP29173498A
Other languages
English (en)
Other versions
JP2000122999A (ja
Inventor
恒平 撫原
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NEC Corp
Original Assignee
NEC Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by NEC Corp filed Critical NEC Corp
Priority to JP29173498A priority Critical patent/JP3709291B2/ja
Publication of JP2000122999A publication Critical patent/JP2000122999A/ja
Application granted granted Critical
Publication of JP3709291B2 publication Critical patent/JP3709291B2/ja
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Landscapes

  • Complex Calculations (AREA)

Description

【0001】
【発明の属する技術分野】
本発明は、高速フーリエ変換に関し、特にマイクロプロセッサを用いた高速フーリエ変換の高速演算方法及び装置に関する。
【0002】
【従来の技術】
フーリエ変換 (Fourier Transform) は、時間領域の信号を周波数領域へ変換する手法である。ディジタル信号処理では、離散時間信号に対するフーリエ変換である 離散フーリエ変換 (Discrete Fourier Transform、DFT) が、信号の周波数解析、音声認識、ブロック・フィルタの高速演算などに広く応用されている。しかし、DFT は原理的に演算量が多いため、実際のシステムでは、高速演算手法である高速フーリエ変換 (Fast Fourier Transform、FFT) が用いられる。
【0003】
FFT では、DFT を小さい DFT へ繰り返し分割することで演算量を軽減する。分割の仕方により Radix-2 FFT、Radix-4 FFT など種々の手法がある。 Radix-2 FFT では、入力データ点数を N とすると、N 点 DFT を 2 個の N/2 点 DFT に分割する。次に各 N/2 点 DFT を 2 個の N/4 点 DFT へ分割する。これを繰り返し、最終的に N/2 個の 2 点 DFT まで分割し、この 2 点 DFT を直接計算する。Radix-4 では、N 点 DFT を順次 1/4 に分割し、最終的に 4 点 DFT に分割する。 2 点 DFT、 4 点 DFT は、そのシグナル・フローの形状から「バタフライ演算」とも呼ばれる。
【0004】
FFT では乗算を多用するため、高速な FFT が必要なシステムでは、乗算器を内蔵したディジタル信号処理プロセッサ (Digital Signal Processor、 DSP) が採用されてきた。 最近では、パーソナル・コンピュータや情報家電機器などのマイクロプロセッサ応用製品で、ソフトウェアによるマルチメディア信号処理を行いたいという要求があり、汎用マイクロプロセッサでも乗算器を内蔵するものがある。
【0005】
さらに、SIMD (Single Instruction stream-Multiple Data streams) 並列演算命令を導入し、さらに高いマルチメディア信号処理性能を実現したマイクロプロセッサもある。このようなマイクロプロセッサの拡張例としては、日経エレクトロニクス、No. 661、 1996 年 5 月 6 日、 105 から 119 ページに 「86 系 マイクロプロセサのマルチメディア命令 MMX を解剖」と題して掲載 された Intel の MMX 命令セット、 Kazumasa Suzuki 他により IEEE MICRO Magazine、 Vol。 18、 No. 2、 1998 年 4 月、 36 から 47 ページに "V830R/AV: Embedded multimedia RISC processor" と題して掲載された NEC の V830R 命令セットがあげられる。
【0006】
これらのマルチメディア信号処理向け命令セットでは、 例えば 64 ビット・レジスタ上に 16 ビットの値 4 個を格納し、 1 命令でこれら 4 個の値に対し 4 並列演算を行える。このような並列演算命令を以下 SIMD 命令と呼ぶ。 典型的な SIMD 命令の例を図 13 に示す。 この例では、64 ビット長のレジスタ r1 (1301)、 r2 (1302) 上に ハーフワード長のデータ 4 個 (a0、 b0、 a1、 b1、 a2、 b2、 a3、 b3) が詰め込まれており、r1 (1301)、 r2 (1302) 上の同じハーフワード位置間で 独立した 4 個の加算 (a0+b0、 a1+b1、 a2+b2、 a3+b3) を行い、 レジスタ r3 (1303) 上に 4 個のハーフワード長の結果を得ている。 ここでハーフワードは 16 ビット・データ長を意味するとする。
【0007】
SIMD 命令は、異なるレジスタ上の同じハーフワード位置間で、同一の演算を複数組のデータに対し並列に行うというその特徴から、 まったく同じシグナル・フローに従って処理する複数組の独立したデータを並列に処理する場合に最も適している。 逆に、同じレジスタ上の異なるハーフワード位置間で演算を行う場合、 SIMD 演算命令のほかにフォーマット変換命令を併用する必要が生じ、オーバーヘッドが発生する。
【0008】
これら SIMD 命令を用い FFT を高速化する試みがすでに発表されている。Intel はアプリケーション・ノート AP-555、 発注番号 243040-001、 1996 年 3 月、 "Using MMX Instructions to Perform Complex 16-bit FFT" において、 SIMD 命令を用いて Radix-2 FFT を高速化した例を示している。 しかし、このアプリケーション・ノートに示された例では、次に詳しく示すように、 演算量に比べてデータ変換のオーバヘッドが大きく、高速な FFT 演算が行えない。Intel のアプリケーション・ノートに示された SIMD 命令による FFT 演算手法 (以下、従来法と呼ぶ) を詳しく説明する。
【0009】
図 11 に、Radix-2 FFT の構成要素である 2 点バタフライ演算の シグナル・フローを示す。 2 点 FFT の入力を X0、 X1、出力を Y0、 Y1 とする。 このように 図11 では入出力を区別するために、入力と出力に別の変数名を つけているが、実際には Y0 を X0 と、Y1 を X1 と同じ場所に格納し、次段のバタフライ演算を行うことが多い。 X0、 X1、 Y0、 Y1 は複素数であり、実部と虚部を区別する必要があるときは、 例えば X0 の実部を xr0、 虚部を xi0、虚数単位を j とおき X0 = xr0 + j xi0 と書くことにする。 シグナル・フロー中のひとつの枝は回転因子と呼ばれる係数 W を持つ。 c は係数 W の実部、s は係数 W の虚部とする。
【0010】
一般的に実際に計算するためには、実部と虚部を別に取り扱う必要がある。 そこで図 11 の 2 点 FFT のシグナル・フローを実部と虚部に分けて 書き下した式を数 1 に示す。
【0011】
【数1】
Figure 0003709291
【0012】
数式 1 の演算を SIMD 命令を用い並列化した例を図 12 に示す。 前提として、2 点バタフライ演算の入力データ xr0、 xi0、 xr1、 xi0 は メイン・メモリ 1201 上のデータ配列 X (1202) に、 回転因子 c、 s は係数テーブル W (1203) に格納されているものとする。 2 点バタフライ演算の結果はふたたびデータ配列 X (1202) に格納する。
【0013】
まず、データ配列 X (1202) から xr1、 xi1 を 64 ビット・レジスタ 2 個 (1204、 1205) の下位 2 ハーフワードに ロードし、フォーマット変換命令 PUNPACKLdq により、1 個の 64 ビット・ レジスタ 1206 に 4 個のハーフワード・データを詰め込む。
【0014】
つぎに 4 個のハーフワード・データを格納したレジスタ 1206 と、 4 個のハーフワード長の係数を係数テーブル 1203 からロードした レジスタ 1207 の間で積和演算 PMADDwd を行う。 積和演算の結果 1208 は 32 ビット長で得られるが、これを PSRAd 命令で右シフトし、 ハーフワードに切り詰める (1209)。
【0015】
ハーフワード長に切り詰めた積和演算結果 1215 は、 零 (1214) から減算することで -1 倍する (1216)。 これともとのハーフワード長の積和演算結果 1209 の間で PACKSSdw 命令を用い フォーマット変換を行い、64 ビット・レジスタ 1211 上に 4 個の積和演算結果を 格納する。
【0016】
さらに、 別の 64 ビット・レジスタ (1217、 1218) の下位 2 ハーフワードに xr0、 xi0 をロードし、 フォーマット変換命令 PUNPACKLdq により、1 個の 64 ビット・ レジスタ 1219 に 4 個のハーフワード・データを詰め込む。 これをレジスタ 1211 と SIMD 加算命令で4 並列加算すると、 最終結果を得る (1213)。 最終結果は、次段の計算のため、あるいは最終結果を残すため、 メイン・メモリ上のデータ配列に格納する (1220)。 yr0 は xr0、 yi0 は xi0、 yr1 は xr1、 yi1 は xi1 が収められていた 位置へ格納する。
【0017】
このように従来法では、1 個の 2 点バタフライ演算から並列性を抽出し、 4 並列 SIMD 命令を適用し、高速化している。
【0018】
【発明が解決しようとする課題】
図 12 に示した従来法は、以下に述べる問題があり、SIMD 命令 に適していない。
【0019】
1 個の 2 点バタフライ演算に入力するデータはわずか 2 組の複素数、 実部、虚部を別に扱ってもデータ点数として 4 点にすぎない。 この少ないデータ点数の中から、典型的な SIMD 演算命令を無駄なく適用する ために 4 並列のデータ並列性を抽出するのはかなり困難である。
【0020】
図 12 に示した例では、あるデータの実部と虚部という異質のデータを 同一レジスタ上で混在させるデータ・フォーマットを採用することにより (1206)、 4 並列のデータ並列性を確保している。 このデータ・フォーマットを実現するには、メイン・メモリ上からロードした データに対しデータ・フォーマット変換命令 PUNPKLdq を適用する必要がある という欠点がある。 このデータ・フォーマット変換命令は、xr0、 xi0 の組に対してと xr1、 xi1 の組に対しての計 2 回適用されている。
【0021】
また、最下位ハーフワードの乗算結果と最下位から 2 番目のハーフワード の乗算結果を加算し 32 ビットの結果を得ているため、 並列度が 2 に落ちているところがある (1208)。 さらに 32 ビットの結果を 16 ビット長にフォーマット変換するための 右シフト命令が余計に必要になるという欠点がある。
【0022】
このように、従来法が採用している並列性の抽出法は、 1 個の 2 点バタフライ演算という狭い範囲から並列性を抽出し、 実部と虚部を同一レジスタ上に混在するデータ・フォーマットを採用し 並列演算をしているため、 データ・フォーマット変換のオーバーヘッドが大きい、 また並列度が途中で 2 に制限されるという欠点があり、 SIMD 命令に適していない。
【0023】
本発明の目的は、SIMD 命令に適した高速な FFT 演算手法を提案する ことにある。 具体的には、Radix-4 FFT アルゴリズムを、 4 並列の SIMD 演算命令セット上で常に 4 並列で行え、 データのフォーマット変換のオーバーヘッドを削減した 高速な演算手法を提案する。
【0024】
【課題を解決するための手段】
本発明は、Radix-4 FFT の全段において、4 個の 4 点バタフライ演算を SIMD命令 を用い並列に処理する。 より具体的には、独立した入力データの実部を保持する配列 (図7の 701) と 入力データの虚部を保持する配列 (図 7 の 702)、 それぞれの配列から連続した 4 要素を取り出し保持する手段 (図 2 )、 必要な係数を生成し保持する手段 (図 3 )、および図 13 に示すような SIMD 並列演算命令を実行する機構から構成される。 最終段においては、係数の生成および保持手段は不要であるが、 4 行 4 列の行列の転置を高速に行う機構 ( 図 4 、 405、 406、 408、 409) が 必要である。 Radix-4 FFT では、隣接する 4 要素は独立したバタフライ演算に属するため、 隣接する 4 要素を同一レジスタに入れ、 通常のスカラー演算と同様の 4 点バタフライ演算のシグナル・フローに従い、 SIMD 演算命令を用いて処理を行うことで、4 個の 4 点バタフライ演算を並列処理 できる。 最終段においては、隣接する 4 要素間でバタフライ演算を行うため、 前段までとは異なり、 隣接する 4 個のバタフライ演算への入力データを 1 個のレジスタにロードする。 すなわち 4 行 4 列の行列の転置と同等の操作により 3 要素おきに 4 個のデータを レジスタにロードし、 通常のスカラー演算と同様の 4 点バタフライ演算のシグナル・フローに従い、 SIMD 演算命令を用いて処理を行うことで、4 個の 4 点バタフライ演算を並列処理 する。
【0025】
【発明の実施の形態】
本発明の実施例を図1から図9、数2、 3、 4 を用いて説明する。
【0026】
【数2】
Figure 0003709291
【0027】
【数3】
Figure 0003709291
【0028】
【数4】
Figure 0003709291
【0029】
まず Radix-4 FFT について説明する。 Radix-4 FFT は、従来例で用いていた Radix-2 FFT より演算量が少なく 高速である。 また、 R. Meyer とK. Schwarz が ICASSP (International Conference on Acoustics、 Speech and Signal Processing) 1990 年、 予稿集 1503 から 1506 ページで ``FFT Implementation on DSP Chips - Theory and Practice'' と題して 述べているように、 Radix-4 FFT は SR (Split Radix) 型 FFT に比べ、 乗算の回数は多少多いもののデータのアドレス計算が単純で加減算の回数 が少ないため、 乗算も加減算も同じ 1クロックで実行できる最近のマイクロプロセッサ上 で高速に実行できる。
【0030】
N 点の DFT を Radix-4 FFT アルゴリズムで求める時の手順を説明する。 ここで N は 4 のべき乗とする。 N 点の DFT をまず 4 個の N/4 点 DFT に分割する。 つぎにそれぞれの N/4 点 DFT を 4 個の N/16 点 DFT に分割する。 この 4 分割を繰り返し、最終的に N/4 個の 4 点 DFT まで分割し、 この 4 点 DFT を直接計算する。
【0031】
N 点 FFT は log4N 段の 4 点バタフライ演算から構成される。 また各段は N/4 個の 4 点バタフライ演算から構成される。 従って、N 点 FFT 全体は N/4 log4N 個の 4 点バタフライ演算を含む。
【0032】
4 点 DFT (バタフライ演算) のシグナル・フローを 図 5 に示す。 これを実部、虚部に分け書き下した数式を数 2 に示す。 図 5 中、係数 W は数 3 に示すとおりである。 また数 2 では簡単化のため、 数 4 に示す中間項 r1、 s1、 r2、 s2、 r3、 s3 を使っている。 以下、4 点 DFT の内部演算ではなく、4 点 DFT の入出力に注目するため、 図 5 に示した 4 点 DFT を 図 6 のように簡略に表記する。
【0033】
各段のバタフライ演算への入力データ、および係数を説明する。 このとき 図7に示すように、 N 点 FFT への入力データである N 個の複素数は、 実部が配列 xr[0 〜 (N-1)] (701)、虚部が配列 xi[0 〜 (N-1)] (702) に別々に 収められているとする。 また、係数 W は、虚部と実部が交互に収められているとする (703)。
【0034】
各段のバタフライ演算への入力データ、および係数の関係を 図 8 に示す。 第 i 段 (1 ≦ i ≦ log4N) を構成するバタフライ演算へ入力する データ ( 図 5中 X0、 X1、 X2、 X3) の、配列 xr (701)、 xi (702) 上での 添字の間隔 Sd は N/(4i) である。 また、係数の間隔 ( 図 7 中 m) の間隔 Sc は 4i-1 である。
【0035】
64 点 Radix-4 FFT の例で説明する。 64 点 Radix-4 FFT のシグナル・フローを 図 9 に示す。 64 点 Radix-4 FFT は、4 点バタフライ演算 3 段から構成され、各段は 16 個の バタフライ演算を含む。 第 1 段の 16 個のバタフライ演算のデータ間隔 Sd は 16、 係数間隔 Sc は 1 である。
【0036】
すなわち、例えば第 1 のバタフライ演算 901 における X0 の実部は xr[0]、虚部は xi[0]、 X1 の実部は xr[16]、 虚部は xi[16]、 X2 の実部は xr[32]、虚部は xi[32]、 X3 の実部は xr[48]、 虚部は xi[48] である。 このバタフライ演算の係数 W は数 3 で n=0 とおいた場合に相当し、 Wn = W2n = W3n = W0 であり、 実部 が cos(0)、 虚部が sin(0) の値を 図 7 に示した係数テーブル W (703) から取得する。
【0037】
次に、第 2 のバタフライ演算 902 における X0 の実部は xr[1]、虚部は xi[1]、 X1 の実部は xr[17]、 虚部は xi[17]、 X2 の実部は xr[33]、 虚部は xi[33]、 X3 の実部は xr[49]、 虚部は xi[49] である。 このバタフライ演算の係数 W は数 3 で n を第 1 段における 係数間隔 1 とおいた場合に相当し、 W1 の実部は W[2]、 虚部は W[3]、 W2 の実部は W[4]、 虚部は W[5]、 W3 の実部は W[6]、 虚部は W[7] の値を 図 7 に示した係数テーブル W (703) から取得する。
【0038】
64 点 FFT の第 2 段の 16 個のバタフライ演算演算のデータ間隔 Sd は 4、 係数間隔 Sc も 4 である。 すなわち、例えば第 3 のバタフライ演算 903 における X0 の実部は xr[0]、虚部は xi[0]、 X1 の実部は xr[4]、 虚部は xi[4]、 X2 の実部は xr[8]、 虚部は xi[8]、 X3 の実部は xr[12]、 虚部は xi[13]である。 このバタフライ演算の係数 W は数式 3 で n=0 とおいた場合に相当し、 Wn = W2n = W3n = W0 であり、 実部 が cos(0)、 虚部が sin(0) の値を 図 7 に示した係数テーブル W (703) から取得する。
【0039】
次に第 4 のバタフライ演算 904 における X0 の実部は xr[1]、虚部は xi[1]、 X1 の実部は xr[5]、 虚部は xi[5]、 X2 の実部は xr[9]、 虚部は xi[9]、 X3 の実部は xr[13]、 虚部は xi[13] である。 このバタフライ演算の係数 W は数 3 で n を第 2 段における 係数間隔 4 とおいた場合に相当し、 W4 の実部は W[8]、 虚部は W[9]、 W8 の実部は W[16]、 虚部は W[17]、 W12 の実部は W[24]、 虚部は W[25] の値を 図 7 に示した係数テーブル W (703) から取得する。
【0040】
第 2 段では、ブロック 0 (907)、 ブロック 1 (908)、 ブロック 2 (909)、 ブロック 3 (910)で係数を共通に使える。 数 3 に示した W の定義と、sin、 cos の周期性から、 第 2 段では 4 バタフライ演算ごとに 同じ係数を用いたバタフライ演算が出現するからである。
【0041】
64 点 Radix-4 FFT の第 3 段の 16 個のバタフライ演算のデータ間隔 Sd は 1、 係数間隔 Sc は 16 である。 すなわち、64 点 Radix-4 FFT の最終段である第 3 段では、隣接する 4 要素間 で 4 点バタフライ演算を行う。例えば第 5のバタフライ演算 905 における X0 の実部は xr[0]、虚部は xi[0]、 X1 の実部は xr[1]、 虚部は xi[1]、 X2 の実部は xr[2]、 虚部は xi[2]、 X3 の実部は xr[3]、 虚部は xi[3] である。 また、第 6 のバタフライ演算 906 における X0 の実部は xr[4]、虚部は xi[4]、 X1 の実部は xr[5]、 虚部は xi[5]、 X2 の実部は xr[6]、 虚部は xi[6]、 X3 の実部は xr[7]、 虚部は xi[7] である。 係数は sin、 cos の周期性から単純になる。
【0042】
以上、Radix-4 FFT のシグナルフローと 64 点 Radix-4 FFT における例を示した。 次にこの Radix-4 FFT を SIMD 命令を用いて 4 並列化する本発明の手法を示す。
【0043】
図 1 に、本発明の Radix-4 FFT の制御の流れをフローチャートで示す。 N を FFT 点数とする。実行開始 (101) 後、 データ距離 Sd の初期値を N に、 係数距離 Sc の初期値を 1 に設定する (102)。
【0044】
Sd と Sc の初期設定後、 変数 k についてのループ (120) を実行する。このループは最終段以外の各段の バタフライ演算に対応しており、N 点 Radix-4 FFT では (log4N) - 1 回 繰り返される。 これは k を N に初期化し (103)、ループ末尾で条件 114 が成立している限り、 k を更新しながらループ内を実行することで行う。
【0045】
変数 k についてのループ (120) 内部では、 図 8 に示したとおり データ距離 Sd をバタフライ演算 1 段ごとに 1/4 にする。 さらに係数の配列 W をアクセスする時の添字となる変数 a を 0 に初期化する (104)。 その後、変数 j についてのループ (121) を実行する。
【0046】
変数 j についてのループ (121) は係数 W のアクセスに対応する。 N 点 FFTの、例えば第 1 段では、N/4 個のバタフライ演算を実行するため、 後述するように 4 個のバタフライ演算を並列に処理する本発明では、 変数 j のループはN/16 回 繰り返される。 第 2 段では、第 1 段と同様 N/4 個のバタフライ演算を実行するが、 4 個のバタフライ演算間で係数の組が共有できるため、 4 個のバタフライ演算を並列に処理する本発明では、 変数 j のループは第 1 段の 1/4 の N/64 回繰り返される。 これは j を 0 に初期化し (105)、ループ末尾で条件 113 が成立している限り j を更新しながらループ内の係数アクセスを実行することで行う。 係数アクセスの詳細については、 図 3 を参照しながら後述する。 係数アクセス後には係数アクセス時の添字 a を係数距離 Sc だけ増加させる。
【0047】
変数 i0 についてのループ (122) は、第 1 の 4 並列バタフライ演算 (110) および演算に必要なデータのアクセス (109)、 演算結果の書き込み (111) に対応する。 N 点 FFT の、例えば第 1 段では、1 組の係数は 1 組のバタフライ演算でのみ 用いられるため、j についてのループ 1 回につき、i0 についてのループも 1 回実行する。第 2 段では、1 組の係数を 4 組のバタフライ演算で再利用 できるため、j についてのループ 1 回につき、i0 についてのループを 4 回実行する。 これは、i0 を j に初期化し (108)、ループ末尾で条件 112 が成立している限り i0 を更新しながらデータ・アクセス (109)、第 1 の 4 並列バタフライ演算 (110)、 結果書き込み (111) を繰り返すことで行う。
【0048】
データ・アクセス (109) の詳細について 図 2 を参照しながら説明する。 本発明の 並列 Radix-4 FFT では、 配列 xr (201)中、添字 i0 から始まる 連続した 4 要素 (202) を収めたレジスタ xr0 (206)、 配列 xi (218)中、添字 i0 から始まる 連続した 4 要素 (211)を収めたレジスタ xi0 (215)、 配列 xr (201)中、添字 i0+Sd から始まる 連続した 4 要素 (203) を収めたレジスタ xr1 (207)、 配列 xi (208)中、添字 i0+Sd から始まる 連続した 4 要素 (212) を収めたレジスタ xi1 (216)、 配列 xr (201)中、添字 i0+2*Sd から始まる 連続した 4 要素 (204) を収めたレジスタ xr2 (208)、 配列 xi (208)中、添字 i0+2*Sd から始まる 連続した 4 要素 (213)を収めたレジスタ xi2 (217)、 配列 xr (201)中、添字 i0+3*Sd から始まる 連続した 4 要素 (205) を収めたレジスタ xr3 (209)、 配列 xi (208)中、添字 i0+3*Sd から始まる 連続した 4 要素 (214) を収めたレジスタ xi3 (218) を入力とする、独立した 4 組のバタフライ演算を SIMD 命令を用いて 並列に行う。 配列 xr、 xi 中の連続した 4 要素は、64 ビット長のロード命令 1 命令で レジスタに格納できる。
【0049】
xr0、 xi0、 xr1、 xi1、 xr2、 xi2、 xr3、 xi3 に対し、 次に述べる係数アクセス (106) において取得した s1、 c1、 s2、 c2、 s3、 c3 とともに、 数 2 に従って 4 点バタフライ演算を行う。 4 並列 SIMD 演算に用いるデータを格納する 64 ビット・レジスタの 各ハーフワード位置は互いに独立であり、 同一のハーフワード位置間でのみ乗算、加減算などの演算が行われるため、 一般のスカラー型のデータに対する 4 点バタフライ演算の定義式である 数 2 をそのまま SIMD データ型に適用すると、4 組のバタフライ演算 が同時に求められる。
【0050】
64 点 Radix-4 FFT を例にとり、データ・アクセスを具体的に説明する。 図 10は、64 点 Radix-4 FFT 第 1 段の、第 1、第 2 回めの i0 についての ループ (122) 内部の実行において、並列処理されるバタフライ演算を示している。 第 1 回めの実行では、入力データとして xr[0〜3]、 xi[0〜3]、 xr[16〜19]、 xi[16〜19]、 xr[32〜35]、 xi[32〜35]、 xr[48〜51]、 xi[48〜51] をレジスタ xr0 (1001)、 xi0 (1002)、 xr1 (1003)、 xi1 (1004)、 xr2 (1005)、 xi2 (1006)、 xr3 (1007) xi3 (1008) に格納し、 数 2 のシグナル・フローに従って 4 点バタフライ演算を実行する。 これにより、 図 9 に示した 64 点 FFT の第 1 段のバタフライ演算 16 個 のうち、 図 10 中実線で示したバタフライ演算 4 個 (1009、 1010、 1011、 1012) が 並列に計算できる。 次の繰り返しでは、入力データとして xr[4〜7]、 xi[4〜7]、 xr[20〜23]、 xi[20〜23]、 xr[36〜39]、 xi[36〜39]、 xr[52〜55]、 xi[52〜55] をレジスタに格納し、 64点 FFT の第 1 段のバタフライ演算のうち、 図 10中破線で示したバタフライ演算 4 個が並列に計算できる。 以下同様に 64 点 Radix-4 FFT では、第 1 段のバタフライ演算、 第 2 段のバタフライ演算各 16 個を 4 個ずつ計算する。
【0051】
次に係数アクセスの詳細について、図 3 を参照しながら説明する。 1 組のバタフライ演算では、 図 5 に示したように、3 個の複素数の 係数 W が必要である。 本発明では、4 組のバタフライ演算を並列に行うので、 係数アクセス (106) では、4 組のバタフライ演算に必要な 12 個の複素数 の係数を、係数テーブル W (703) からレジスタに読み出す。
【0052】
あるバタフライ演算が数 3 に示した Wn として 係数テーブル中の a 番目の複素数の要素を必要としているとする。 このときこの複素数の虚部は係数テーブル (703) 中 W[a*2]、 実部は W[a*2+1] に格納されている。 このバタフライ演算は同時に W2n として、虚部が W[a*4]、 実部が W[a*4+1] である複素数を、 W3n として、虚部が W[a*6]、 実部が W[a*6+1] である複素数を必要とする。
【0053】
このバタフライ演算と同時に演算される 2 番目のバタフライ演算は、係数として 係数テーブル W (703) 上で係数距離 Sc だけ離れた位置に格納されている W[(a+Sc)*2]、 W[(a+Sc)*2+1]、 W[(a+Sc)*4]、 W[(a+Sc)*4+1]、 W[(a+Sc)*6]、 W[(a+Sc)*6+1] を必要とする。 同様に、3 番目のバタフライ演算は、係数として 係数テーブル W (703) 上でさらに係数距離 Sc だけ離れた位置に格納されている W[(a+2*Sc)*2]、 W[(a+2*Sc)*2+1]、 W[(a+2*Sc)*4]、 W[(a+2*Sc)*4+1]、 W[(a+2*Sc)*6]、 W[(a+2*Sc)*6+1] を必要とする。 同様に、4 番目のバタフライ演算は、係数として 係数テーブル W (703) 上でさらに係数距離 Sc だけ離れた位置に格納されている W[(a+3*Sc)*2]、 W[(a+3*Sc)*2+1]、 W[(a+3*Sc)*4]、 W[(a+3*Sc)*4+1]、 W[(a+3*Sc)*6]、 W[(a+3*Sc)*6+1] を必要とする。
【0054】
図 3 に示すように、 係数テーブル上の必要な値は、いったん 2 個ずつレジスタに格納され、 6 個のレジスタに詰め込まれる。 各レジスタ上には、各係数の実部または虚部に対応する値が、 4 バタフライ演算分格納される。 このフォーマット変換は、例えば NEC V830R プロセッサの場合 2 段階の フォーマット変換命令 (314、 315、 316、 317) で行える。 4 バタフライ演算分のW1に対応する実部の値を格納したレジスタを c1 (309)、 虚部の値を格納したレジスタを s1 (308)、 W2に対応する実部の値を格納したレジスタを c2 (311)、 虚部の値を格納したレジスタを s2 (310)、 W3に対応する実部の値を格納したレジスタを c3 (313)、 虚部の値を格納したレジスタを s3 (312)とおく。
【0055】
最終段の演算について、図 4 を参照しながら説明する。 図 9 で示したように、最終段より前では、 配列 xr (701)、 xi (712) の連続した 4 要素は独立したバタフライ演算に属しているため、 連続 4 要素をレジスタに格納し、これら 4 個の要素に対し 同じシグナル・フローに従って演算することにより、 4 個のバタフライ演算を並列に演算できた。 最終段でも、前段までと同様、4 個のバタフライ演算を同一のシグナル・フローに 従い並列に処理する。 しかし、最終段では 図 9に示した 64 点 FFT の例でもわかるように、 データ配列 xr (701)、 xi (712) の連続した 4 要素は 同じバタフライ演算に属しており、依存関係があるため、 異なった並列化手法をとる必要がある。
【0056】
そこで最終段では、データのフォーマット変換を行い、 隣接する 4 個のバタフライ演算を SIMD 命令を用いて並列に行う。 最終段の処理の流れを 図 4 に示す。 最初に j を 0 に初期化し (401)、データ配列 xr (701)、 xi (702) の、 添字 j から始まる各 16 個の要素の実部を レジスタ xr0、 xr1、 xr2、 xr3 にロードし (403)、 虚部をレジスタ xi0、 xi1、 xi2、 xi3 に ロードすると (404)、 連続する 4 要素が同一レジスタにロードされる。 次に実部、虚部各 4 本のレジスタ上に置いた 16 個の要素を、4 行 4 列の行列 とみなし、行列の転置を行う (405、 406)。 このデータ・フォーマット変換により ひとつのバタフライ演算に必要な、連続 4 要素は別のレジスタに格納され、 独立した 4個のバタフライ演算に必要な 4 要素は同一レジスタの異なる ハーフワード位置に格納できる。 このフォーマット変換のオーバーヘッドは小さい。 例えば NEC V830R プロセッサでは、フォーマット変換命令 8 命令 で 4 個のレジスタ上の 16 要素の転置が行える。
【0057】
転置後のレジスタ xr0、 xi0、 xr1、 xi1、 xr2、 xi2、 xr3、 xi3 に対し、 図 5 に示したシグナル・フローに従って 4 点バタフライ演算を行うと、 隣接する 4 個のバタフライ演算に対し並列演算が行える (407)。 バタフライ演算の結果は再び転置し (408、 409)、配列 xr (701)、 xi (702) に最終結果として書き込む (410)。
【0058】
本発明の実施例は、 64 ビット・レジスタ上に 16 ビット・データを 4 個格納し、 4 並列演算を行うという、現在一般的な SIMD 命令セットを用い説明したが、 本発明はレジスタのビット長、各データの長さ、並列度が異なるプロセッサに 対しても容易に適用できる。
【0059】
【発明の効果】
本発明の Radix-4 FFT 演算手法によれば、4 並列 SIMD 命令を持つプロセッサ上 で常に 4 個のバタフライ演算を並列に行うことができ、高い並列性を達成できる。 また、同一レジスタに実部のみもしくは虚部のみを収めるという規則的な データ構造を採用しているため、データのフォーマット変換のオーバーヘッドを 最小化できる。 結果として、高速な Radix-4 FFT の演算が可能である。 例えば、 NEC V830R プロセッサ上で本発明の手法により 256 点 FFT を実行した場合、 所要クロック数は4,093 クロックである。 これは V830R と類似した SIMD 命令セットを持つ Intel PentiumII プロセッサ上で 従来の手法により 256 点 FFT を実行した場合の所要クロック数5,522 クロック に比べ、35% も高速である。
【図面の簡単な説明】
【図 1】本発明の FFT 方法を示すフローチャートである。
【図 2】データアクセス方法を示す図である。
【図 3】係数アクセス方法を示す図である。
【図 4】本発明の FFT の最終段の処理を示すフローチャートである。
【図 5】4 点バタフライ演算のシグナル・フロー図である。
【図 6】の簡略化した 4 点バタフライ演算のシグナル・フロー図である。
【図 7】本発明の FFT におけるデータ構造を示す図である。
【図 8】バタフライ演算で使用するデータの間隔と係数の間隔を示す表である。
【図 9】64 点 Radix-4 FFT のシグナル・フロー図である。
【図 10】本発明の方法で 64 点 FFT を 4 並列化した場合のシグナル・フロー図である。
【図 11】2 点バタフライ演算のシグナル・フロー図である。
【図 12】従来の並列 FFT 演算手法を示す図である。
【図 13】SIMD 演算命令の動作例を説明する図である。
【符号の説明】
201 入力データ X の実部を収めた配列 xr
206、207、208、209 xr の 4 要素を収めた 64 ビット・レジスタ
214、215、216、217 64 ビット・レジスタ
218 入力データ X の虚部を収めた配列 xi
301 係数テーブル W
302、303、304、305 係数テーブルのデータを 2 個ロードした 64 ビット・レジスタ
306、307 フォーマット変換の中間結果を保持する 64 ビット・レジスタ
308、309 4 組の係数を保持する 64 ビット・レジスタ
310、311 4 組の係数を保持する 64 ビット・レジスタ
312、313 4 組の係数を保持する 64 ビット・レジスタ
701 入力データ X の実部を収めた配列 xr
702 入力データ X の実部を収めた配列 xi
703 係数テーブル W
901 第 1 のバタフライ演算
902 第 2 のバタフライ演算
903 第 3 のバタフライ演算
904 第 4 のバタフライ演算
905 第 5 のバタフライ演算
906 第 6 のバタフライ演算
907 第 0 ブロック
908 第 1 ブロック
909 第 2 ブロック
910 第 3 ブロック
1001 最初の繰り返しで使用する xr0 の値を保持する 64 ビット・レジスタ
1002 最初の繰り返しで使用する xi0 の値を保持する 64 ビット・レジスタ
1003 最初の繰り返しで使用する xr1 の値を保持する 64 ビット・レジスタ
1004 最初の繰り返しで使用する xi1 の値を保持する 64 ビット・レジスタ
1005 最初の繰り返しで使用する xr2 の値を保持する 64 ビット・レジスタ
1006 最初の繰り返しで使用する xi2 の値を保持する 64 ビット・レジスタ
1007 最初の繰り返しで使用する xr3 の値を保持する 64 ビット・レジスタ
1008 最初の繰り返しで使用する xi3 の値を保持する 64 ビット・レジスタ
1009 第 7 の 4 点バタフライ演算
1010 第 8 の 4 点バタフライ演算
1011 第 9 の 4 点バタフライ演算
1012 第 10 の 4 点バタフライ演算
1201 メイン・メモリ
1202 データ配列 X
1203 係数テーブル W
1204、1205、1206、1207、1208、1209、1210、1211、1212、1213、1214、1215、1216、 1217、1218、1219 64 ビット・レジスタ
1301 64 ビット・レジスタ
1302 64 ビット・レジスタ
1303 64 ビット・レジスタ

Claims (2)

  1. 複数のデータを同一のレジスタに格納し,それぞれのデータに対し同一の演算を並列に行う,いわゆる SIMD 型演算命令を備えるマイクロプロセッサを用いた高速複素フーリエ変換方法において,
    入力データである複素数の実部と虚部をそれぞれ独立したデータ格納場所に格納するステップと,
    最終段以外のバタフライ演算では隣接する 4 組のバタフライ演算への入力データである 16 個の複素数の実部と虚部計 32 個のデータを,バタフライ演算への入力データ位置が同一かつデータ格納場所で隣接する 4 個の実部データ 4 組と,バタフライ演算への入力データ位置が同一かつデータ格納場所で隣接する 4個の虚データ 4 組に分割し,これら計 8 組の実部,虚部のデータをそれぞれ異なるレジスタにロードして4 点バタフライ演算を行い,ふたたび前記格納場所に格納するステップと,
    最終段のバタフライ演算では最終段の入力データとして前記データ格納場所からバタフライ演算への 4 個の入力データ位置それぞれについて,実部,虚部とも前記データ格納場所から 3 要素おきに 4 個の要素を 1 組として1 個のレジスタにロードすることを繰り返し,これら計 8 組のデータをそれぞれ異なるレジスタにロードし,4 点バタフライ演算を行い,ふたたび前記格納場所に格納するステップを備え,
    全段で独立した 4 個の 4 点バタフライ演算を並列に処理することを特徴とする高速複素フーリエ変換方法。
  2. 複数のデータを同一のレジスタに格納し,それぞれのデータに対し同一の演算を並列に行う,いわゆる
    SIMD 型演算命令を備えるマイクロプロセッサを用いた高速複素フーリエ変換装置において,
    入力データである複素数の実部と虚部をそれぞれ独立したデータ格納場所に格納する手段と,
    最終段以外のバタフライ演算では隣接する 4 組のバタフライ演算への入力データである 16 個の複素数の実部と虚部計 32 個のデータを,バタフライ演算への入力データ位置が同一かつデータ格納場所で隣接する 4 個の実部データ 4 組と,バタフライ演算への入力データ位置が同一かつデータ格納場所で隣接する 4個の虚部データ 4 組に分割し,これら計 8 組の実部,虚部のデータをそれぞれ異なるレジスタにロードして4 点バタフライ演算を行い,ふたたび前記格納場所に格納する手段と,
    最終段のバタフライ演算では最終段の入力データとして前記データ格納場所からバタフライ演算への 4 個の入力データ位置それぞれについて,実部,虚部とも前記データ格納場所から 3 要素おきに 4 個の要素を 1 組として1 個のレジスタにロードすることを繰り返し,これら計 8 組のデータをそれぞれ異なるレジスタにロードし,4 点バタフライ演算を行い,ふたたび前記格納場所に格納する手段を備え,
    全段で独立した 4 個の 4 点バタフライ演算を並列に処理することを特徴とする高速複素フーリエ変換装置.
JP29173498A 1998-10-14 1998-10-14 高速複素フーリエ変換方法及び装置 Expired - Fee Related JP3709291B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP29173498A JP3709291B2 (ja) 1998-10-14 1998-10-14 高速複素フーリエ変換方法及び装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP29173498A JP3709291B2 (ja) 1998-10-14 1998-10-14 高速複素フーリエ変換方法及び装置

Publications (2)

Publication Number Publication Date
JP2000122999A JP2000122999A (ja) 2000-04-28
JP3709291B2 true JP3709291B2 (ja) 2005-10-26

Family

ID=17772718

Family Applications (1)

Application Number Title Priority Date Filing Date
JP29173498A Expired - Fee Related JP3709291B2 (ja) 1998-10-14 1998-10-14 高速複素フーリエ変換方法及び装置

Country Status (1)

Country Link
JP (1) JP3709291B2 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100692997B1 (ko) * 2001-04-17 2007-03-12 삼성전자주식회사 패스트퓨리에변환 장치
US8229014B2 (en) 2005-03-11 2012-07-24 Qualcomm Incorporated Fast fourier transform processing in an OFDM system
US8266196B2 (en) * 2005-03-11 2012-09-11 Qualcomm Incorporated Fast Fourier transform twiddle multiplication
KR20120077164A (ko) 2010-12-30 2012-07-10 삼성전자주식회사 Simd 구조를 사용하는 복소수 연산을 위한 사용하는 장치 및 방법

Also Published As

Publication number Publication date
JP2000122999A (ja) 2000-04-28

Similar Documents

Publication Publication Date Title
US6304887B1 (en) FFT-based parallel system for array processing with low latency
JP2000285105A (ja) 行列ベクトル乗算命令を用いて高速フーリエ変換を実行するための方法及びシステム
US20010032227A1 (en) Butterfly-processing element for efficient fast fourier transform method and apparatus
KR100538605B1 (ko) 데이터 처리 장치, 행렬 변환 방법 및 컴퓨터 판독가능한 매체
US7761495B2 (en) Fourier transform processor
AU579621B2 (en) Computer and method for discrete transforms
US6993547B2 (en) Address generator for fast fourier transform processor
Jiang et al. Twiddle-factor-based FFT algorithm with reduced memory access
JP3709291B2 (ja) 高速複素フーリエ変換方法及び装置
US20060075010A1 (en) Fast fourier transform method and apparatus
Cui-xiang et al. Some new parallel fast Fourier transform algorithms
US7774397B2 (en) FFT/IFFT processor
Nadehara et al. Radix-4 FFT implementation using SIMD multimedia instructions
EP1447752A2 (en) Method and system for multi-processor FFT/IFFT with minimum inter-processor data communication
EP1426872A2 (en) Linear scalable FFT/IFFT computation in a multi-processor system
Takala et al. Butterfly unit supporting radix-4 and radix-2 FFT
Sorensen et al. Efficient FFT algorithms for DSP processors using tensor product decompositions
CN104572578B (zh) 用于显著改进微控制器中fft性能的新颖方法
WO1999053417A1 (en) Device for converting series of data elements
US20180373676A1 (en) Apparatus and Methods of Providing an Efficient Radix-R Fast Fourier Transform
JP2000231552A (ja) 高速フーリエ変換方法
US20040003017A1 (en) Method for performing complex number multiplication and fast fourier
EP0988605A2 (en) Device for converting series of data elements
Pei et al. Efficient bit and digital reversal algorithm using vector calculation
Yeh Parallel implementation of the fast Fourier transform on two TMS320C25 digital signal processors

Legal Events

Date Code Title Description
A02 Decision of refusal

Free format text: JAPANESE INTERMEDIATE CODE: A02

Effective date: 20020205

RD01 Notification of change of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7421

Effective date: 20050310

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20050808

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20080812

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20090812

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20090812

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100812

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110812

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110812

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120812

Year of fee payment: 7

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130812

Year of fee payment: 8

LAPS Cancellation because of no payment of annual fees