(57)【特許請求の範囲】
【請求項１】
ベクトル計算機を用いて要素数ＮのデータＸ（Ｎ）の複素高速フーリエ変換を行う処理方法において，
前記変換対象となる要素数Ｎを，Ｎ＝ｐ・ｑ・ｐの形の３つの因子に分解する処理を行う過程と，
Ｗp ×Ｉpq〔ただし，×はクロネッカー積，Ｗp は要素数ｐのフーリエ変換行列，Ｉpqは次数がｐｑの単位行列。〕による第１の変換であって，該第１の変換は，前記長さＮ＝ｐ・ｑ・ｐのデータＸ（Ｎ）をｐｑ・ｐなる２次元データＸ′（０：ｐｑ−１，０：ｐ−１）とみなして，該２次元データに対する２次元目に関して長さｐのフーリエ変換をｐ・ｑ本多重に行う変換である処理と，
（Ｄq p ×Ｉp ），Ｄpqp またはＤpqp ，（Ｄq p ×Ｉp ）〔ただし，Ｄq p はｄｉａｇ（ωj sm），ｊ＝ｓｑ＋ｍなる対角行列。ｊは０からｐ＊ｑ−１まで動き，ｍは０からｑ−１まで動き，ｓは０からｐ−１まで動く。ω＝ｅｘｐ（SN２πｉ／（ｐ＊ｑ）），SNは＋１か−１，Ｉp は次数がｐの単位行列。〕による第１の回転であって，該第１の回転は，前記第１の変換による処理結果の長さＮ＝ｐ・ｑ・ｐのデータＹ（Ｎ）を，ｐ・ｑ・ｐなる３次元データＹ′（０：ｐ−１，０：ｑ−１，０：ｐ−１）とみなして，該３次元データの２次元と３次元からなるｐ個のｑ・ｐの２次元データに対して回転因子を掛ける演算を，多重度ｐの１次元目に対して行うものであり，前記回転因子を掛ける演算は，前記１次元目のｐ個に多重に，Ｙ′（ｉ，ｓ，ｍ）にωsmを各ｓ，ｍに対して掛けるものである処理と，
Ｉp ×Ｗq ×Ｉp による第２の変換と途中結果の配列の転置の組み合わせであって，該第２の変換と途中結果の配列の転置の組み合わせは，前記第１の回転の処理結果の長さＮ＝ｐ・ｑ・ｐのデータＺ（Ｎ）を，ｐ・ｑ・ｐなる３次元データＺ′（０：ｐ−１，０：ｑ−１，０：ｐ−１）とみなして，２次元目の長さｑのデータに対するフーリエ変換を，１次元目と３次元目のｐ・ｐ組多重に行うものであり，この結果をメモリに格納するときに，１次元目と３次元目を転置して入れ換えて格納するものである処理と，
（Ｄq p ×Ｉp ），Ｄpqp またはＤpqp ，（Ｄq p ×Ｉp ）による第２の回転であって，該第２の回転は，前記第２の変換と途中結果の配列の転置の組み合わせの結果に対し，前記第１の回転と同様の処理を行うものである処理と，
Ｗp ×Ｉpqによる第３の変換であって，該第３の変換は，前記第２の回転による処理結果に対し，前記第１の変換と同様の処理を行うものである処理とを，
前記分解した要素に対してベクトル処理により実行する過程とを含み，
要素数を再帰的に分解して計算することを特徴とするベクトル計算機による複素高速フーリエ変換処理方法。
【請求項２】
ベクトル計算機を用いて要素数ＮのデータＸ（Ｎ）の複素高速フーリエ変換を行う処理方法において，
前記変換対象となる要素数Ｎを，Ｎ＝ｐ・ｑ・ｐの形の３つの因子に分解する処理を行う過程と，
Ｗp ×Ｉpq〔ただし，×はクロネッカー積，Ｗp は要素数ｐのフーリエ変換行列，Ｉpqは次数がｐｑの単位行列。〕による変換の逆変換であって，該変換の逆変換は，前記要素数ＮのデータＸ（Ｎ）に対する請求項１記載の第３の変換の逆変換である処理と，
（Ｄq p ×Ｉp ），Ｄpqp またはＤpqp ，（Ｄq p ×Ｉp ）〔ただし，Ｄq p はｄｉａｇ（ωj sm），ｊ＝ｓｑ＋ｍなる対角行列。ｊは０からｐ＊ｑ−１まで動き，ｍは０からｑ−１まで動き，ｓは０からｐ−１まで動く。ω＝ｅｘｐ（−SN２πｉ／（ｐ＊ｑ）），SNは＋１か−１，Ｉp は次数がｐの単位行列。〕による回転の逆変換であって，該回転の逆変換は，前記変換の逆変換の処理結果に対する請求項１記載の第２の回転のωをω-1に置き換えた逆変換である処理と，
Ｉp ×Ｗq ×Ｉp による第２の変換の逆変換と途中結果の配列の転置の逆変換の組み合わせであって，該第２の変換の逆変換と途中結果の配列の転置の逆変換の組み合わせは，前記回転の逆変換の処理結果に対する請求項１記載の第２の変換と途中結果の配列の転置の組み合わせにおける第２の変換の逆変換と転置の逆変換である処理と，
（Ｄq p ×Ｉp ），Ｄpqp またはＤpqp ，（Ｄq p ×Ｉp ）による回転の逆変換であって，該回転の逆変換は，前記第２の変換の逆変換と途中結果の配列の転置の逆変換の組み合わせの処理結果に対する請求項１記載の第１の回転のωをω-1に置き換えた逆変換である処理と，
Ｗp ×Ｉpqによる変換の逆変換であって，該変換の逆変換は，前記回転の逆変換の処理結果に対する請求項１記載の第１の変換の逆変換である処理とを，
前記分解した要素に対してベクトル処理により実行する過程とを含み，
要素数を再帰的に分解して計算することを特徴とするベクトル計算機による複素高速フーリエ変換処理方法。
【請求項３】
請求項１または請求項２記載のベクトル計算機による複素高速フーリエ変換処理方法において，
Ｎ＝ｐ・ｑ・ｐの形で要素数を再帰的に分解して計算を行うときに，中間ノードをひとまとめにして，各ｑとｐを決定することを特徴とするベクトル計算機による複素高速フーリエ変換処理方法。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は，ベクトル計算機を用いて複素高速フーリエ変換（ＦＦＴ：Fast Fourier Transform) を行う処理方法に係り，特に，高速にかつ作業領域をほとんど必要としないで処理できるようにしたベクトル計算機による複素高速フーリエ変換処理方法に関する。
【０００２】
主に，科学技術計算を高速に行う計算機として，ベクトル計算機が多く用いられるようになってきており，各種の計算処理において，そのベクトル処理を活かした高速な演算を可能とする技術が望まれている。このような科学技術計算の分野において，特に，複素ＦＦＴは演算量が多く，一般に作業領域も膨大な作業領域を必要とすることが多いので，処理の高速化を図るとともに，膨大な作業領域を必要としない処理技術が必要とされる。
【０００３】
【従来の技術】
ベクトル計算機により複素高速フーリエ変換を行う処理手段は，例えばベクトル計算機向けの科学技術計算用サブールチンライブラリとして提供され，ＦＯＲＴＲＡＮプログラムなどからの呼び出しによってフーリエ変換の処理を行うようになっている。従来からこの計算を高速化するための技術が種々考えられており，例えばＰｅａｓｅのアルゴリズム（ｉｓｏｇｅｏｍｅｔｒｉｃなアルゴリズム），Ｓｔｏｃｋｈａｍのアルゴリズムなどが知られている。これらの参考文献としては，例えば以下のものがある。
〔参考文献〕島崎著：計算機科学／ソフトウェア技術講座９「スーパーコンピュータとプログラミング」，第８章“高速フーリエ変換とベクトル化”共立出版（株），１９８９。
【０００４】
これらの方法では，高速化のためにベクトル長を確保して演算を行うことを目的としているが，Ｓｔｏｃｋｈａｍのアルゴリズムでは，アクセスは一定間隔のデータのアクセスとなる。また，Ｐｅａｓｅのアルゴリズムでは，インデックスをリストとして利用するリストベクトル処理となる。
【０００５】
図７に，それらのアルゴリズムの中核になるループの手続きの例を示す。説明をわかりやすくするため，図ではＦＯＲＴＲＡＮの記述形式で表している。
図７（Ａ）は，Ｓｔｏｃｋｈａｍのアルゴリズムであって，同図からもわかるように，配列Ｃのアクセスは，奇数と偶数の要素をそれぞれアクセスするため，連続域をアクセスすることはできない。また，この中核の部分を実現するため，ω2**(L+1)sg(r) をテーブルに格納しておく必要がある。
【０００６】
ただし，ω2**(L+1)＝ｅｘｐ（２πｉ／２L+1 ）であり，ω2**(L+1)は，１の２L+1 乗根である。ｓは符号で＋１か−１である。
ｇ（ｒ）は，ｉｎｔｅｇｅｒ〔ｒ／Ｎ・２＊＊−（Ｌ＋１）〕である。
【０００７】
図７（Ｂ）は，Ｐｅａｓｅのアルゴリズムの場合の例を示している。Ｐｅａｓｅのアルゴリズムでは，二重ループに関するｊ1 ，ｊ3 に関するリストを実行時に作成し，またはデータとして実行前に保持しておいて，リストベクトルによるベクトル処理を行うようにしていた。このアルゴリズムは，間接指標付きのベクトル演算が高速であるベクトル計算機では高速である。しかし，一般には指標付きのベクトルロードよりも連続域のベクトルロードのほうが高速である。また，インデックスリストを作成する領域も必要となる。
【０００８】
【発明が解決しようとする課題】
本発明は上記問題点の解決を図り，いわゆるインプレース（ｉｎ  ｐｌａｃｅ）かつセルフ・ソーティング（ｓｅｌｆ−ｓｏｒｔｉｎｇ）で，作業領域をほとんど必要とせず，しかも変換のサンプル数が大きくなるにつれて従来の方式よりも高速な複素高速フーリエ変換処理方法を提供することを目的としている。
【０００９】
【課題を解決するための手段】
ＳｔｏｃｋｈａｍのアルゴリズムやＰｅａｓｅのアルゴリズムなどの従来の方法のように，一定間隔でデータをアクセスする方法や，作業領域が多く必要となる間接指標のベクトルアクセスに依存する方法でなく，連続アクセス中心の方法のほうが高速となり得る。
【００１０】
そこで，ベクトル長ＮがＮ＝ｔ・ｕと２つの因子に分解されるとき，ＦＦＴの分解定理を利用することにより，連続アクセス中心の方法に変換することを考える。この変換は，２次元のＦＦＴと似た形をしており，真ん中のローテーションと行列の転置の操作の前後に，ＦＦＴを２次元目の変数に対して行う変換である。この方法は，一般にＮはｔとｕという２つの因子に分解される場合のアルゴリズムで，行列の転置を行うとき作業領域を必要とし，結果としてこの作業領域の大きさは，入力データと同じ大きさを必要とする。したがって，ｉｎ  ｐｌａｃｅなアルゴリズム（余分な作業配列領域なしで演算を実行できるアルゴリズム）ではない。
【００１１】
本発明は，ｉｎ  ｐｌａｃｅにするため，この方法を２回連続して行うことを主要な特徴とする。すなわち，Ｎ＝ｒ・ｓと分解し，さらにｒ＝ｔ・ｕと分解する。ここで，ｔ＝ｓとなるようにｕを選べるとき，２回の転置を連続して行うことを，一回のｉｎ  ｐｌａｃｅな置換で行うことができる。したがって，ｉｎ
ｐｌａｃｅなアルゴリズムを実現することができる。
【００１２】
これは，作業領域が少なく，連続アクセス中心の１次元のＦＦＴであり，特に連続アクセスを高速に実行する能力を持つベクトル計算機に適した処理方法を実現できる。すなわち，この方法は，作業領域は少なく，ｉｎ  ｐｌａｃｅかつｓｅｌｆ−ｓｏｒｔｉｎｇという，よい性質を持っている。また，異なるＦＦＴのアルゴリズムを結合する方式も可能にする。
【００１３】
本発明による配列の処理方法は，以下のとおりである。
図１（Ａ）に示すように，要素数Ｎ＝ｐ2 ｑの１次元配列Ｘ（Ｎ）は，ｐ×ｑ×ｐの３次元配列Ｘ（ｐ，ｑ，ｐ）とみなすことができる。また，ｐｑ×ｐの２次元配列Ｘ（ｐｑ，ｐ）ともみなせる。これに着目して，以下のように処理する。
【００１４】
次のフーリエ変換を行うものとする。
αk ＝Σａj ｅｘｐ（SN２πｉｊｋ／Ｎ）
ただし，Σはｊ＝０からＮ−１までの総和，SNは符号であり１または−１である。ｉは虚数単位である。α＝ＷN ・ａとも記述できる。
【００１５】
ＷN はＮ＝ｐ2 ｑとなるとき，図１（Ｂ）に示すような変換によって，以下の手順で計算する。
▲１▼  変換１（Ｗp ×Ｉpq）→回転１→変換２（Ｉp ×Ｗq ×Ｉp ）→転置→回転２→変換３（Ｗp ×Ｉpq）
▲２▼  変換１→回転１→転置→変換２→回転２→変換３
変換１および変換３は，Ｗp ×Ｉpqのクロネッカー積で表される変換である。変換２は，Ｉp ×Ｗq ×Ｉp のクロネッカー積で表される変換である。転置は，１次元目と３次元目を転置するものである。
【００１６】
回転因子１，回転因子２の組は，行列（Ｄq p ×Ｉp ），Ｄpqp またはＤpqp ，（Ｄq ×Ｉp ）のいずれかの組である。または，各部品の逆変換を逆の順序で作用したものが考えられる。すなわち，次の手順である。
【００１７】
▲３▼  変換３の逆変換→回転２の逆変換→転置の逆変換→変換２の逆変換→回転１の逆変換→変換１の逆変換
▲４▼  変換３の逆変換→回転２の逆変換→変換２の逆変換→転置の逆変換→回転１の逆変換→変換１の逆変換
各々の変換，回転の逆変換は，ω＝ω-1と置き換えたものである。
【００１８】
以上の変換の中の転置は，Ａ（ｐ，ｑ，ｐ）の１次元目と３次元目を転置するもので，各次元長が同じであり，作業領域を必要とせずに行うことができる。
行列Ｗp ×Ｉpqは，長さｐｑの連続ベクトルごとを要素とする要素数ｐのＦＦＴとなり，連続アクセス中心のＦＦＴとなる。
【００１９】
これでは，ベクトル長は，ほぼ〜（Ｎ）1/2 程度のベクトル長にとどまる。さらにベクトル長を長くするために，変換Ｗp ×ＩpqのＷp に同じことを行う。
このことで, ｐ＝ｒ2 ＊ｓと分解できたとき，変換１および変換３のＷp は，同様に変換１と変換３を生じる。この変換は，Ｗr ×Ｉrsとなる。つまりＷr ×Ｉrs×Ｉpq＝Ｗr ×Ｉrspqとなりベクトル長が延びる。このように，ベクトル長ｐｑの演算が，ベクトル長ｒｓｐｑの演算に変換され，かつ連続アクセスとなる。
【００２０】
変換２は，Ｉr ×Ｗs ×Ｉr ×Ｉpq＝Ｉr ×Ｗs ×Ｉrpq となり，長さｒｐｑのベクトルを一要素とみた要素数ｓにＦＦＴをｒ組み行うこととなりベクトル長は長くなる。
【００２１】
また転置は，Ａ（ｐ，ｑ，ｐ）がＡ（ｐ，ｑ，ｒ，ｓ，ｒ）となり，３次元目と５次元目の転置となりベクトル長ｐ＊ｑの転置となる。他に回転因子はｐｑのベクトル長のものと，ｒｐｑベクトル長に分かれる。
【００２２】
変換１，変換３のベクトル長は，
Ｎ＝ｐ1 2 ｑ1 ，ｐ1 ＝ｐ2 2 ｑ2 ，……，ｐk-1 ＝ｐk 2 ｑk 
のように分解されたとき，末端の変換１，変換３のベクトル長はｖk ＝ｐ1 ｑ1 ｐ2 ｑ2 ……ｐk ｑk となり，その中間ノード（深さｉ）のベクトル長は，ｖi-1 ＊ｐi とｖi-1 の２通りが現れる. すなわち，各ノードの分解の深さｉに対応させて，ｖi とノードに対応するｐi ，ｑi がわかれば，計算することができることになる。
【００２３】
例えば，Ｎ＝２11の場合，
Ｎ＝２11＝（（２2 ）2 ＊２）2 ＊２  または  （（２2 ）2 ＊１）2 ＊８
などと分解できる。
【００２４】
図２にその分解の様子を示す。繰り返し適用すると，計算は木構造になり，変換Ｗr ×Ｉrspqがノード▲１▼，▲３▼，▲５▼，▲７▼に対応する。ベクトル長は，ｒｓｐｑ＝２9 になる。連続アクセスでベクトル長は長くなる。転置のベクトル長＝ｐ＊ｑは，図２（Ａ）に示すＴＹＰＥ１では２6 であるが，図２（Ｂ）に示すＴＹＰＥ２では２7 となり長くなる。これで準備ができたので，処理の順序で実行すればよい。
【００２５】
分解がｋまでできたとき，２k+1 −１個のノードができる。
処理順に番号を１〜２k+1 −１まで振ったとき，この番号が２j で割切れると，ｋ＋１−ｊのノードに対応するので，処理を行う。このとき，深さに対応するテーブルに，深さｉのベクトル長とｐi ，ｑi を格納しておくことにより計算することができる。
【００２６】
【作用】
本発明の動作原理は，以下のとおりである。
１次元のＦＦＴを計算する上で，サンプル数を３つの因子に分解する。ｎ＝ｔ・ｕ・ｔと考えて，３つの因子に分解することを考える。
【００２７】
まず，以下のＤ１〜Ｄ３を定義する。
（定義Ｄ１）
Ａ×Ｂ＝（ａijＢ）で定義されるクロネッカー積を表す。
【００２８】
Ｂが（ｑ，ｑ）の行列，Ａは（ｐ，ｐ）の行列のとき，
〔Ａ×Ｂ〕（ｔｑ＋ｕ，ｒｑ＋ｓ）＝〔Ａ〕（ｔ，ｒ）・〔Ｂ〕（ｕ，ｓ）
（定義Ｄ２）
0003647887.tif 000002
これは，ｒｐ＋ｓとｓｑ＋ｒを入れ換える
（定義Ｄ３）
0003647887.tif 000003
以下の関係Ａ１〜Ａ３（補助定理）を利用する。
（Ａ１）    Ｐq p Ｄq p ＝Ｄp q Ｐq p 
ＡP ，Ｂq はそれぞれ（ｐ，ｐ），（ｑ，ｑ）の行列を表す。
（Ａ２）    Ｐq p ・（Ａp ×Ｂq ）＝（Ｂq ×Ａp ）・Ｐq p 
（Ａ３）    （Ａp ・Ｂp ）×Ｉq ＝（Ａp ×Ｉq ）・（Ｂp ×Ｉq ）
ＦＦＴに関する分解定理Ｔ１およびＴ２がある。
〔定理〕
0003647887.tif 000004
〔定理終わり〕
これをＮ＝ｐｑｒと分解される場合に（Ｔ１）および（Ｔ２）を２回連続に作用させて適用し, 補助定理（Ａ１），（Ａ２），（Ａ３）を適用して，
0003647887.tif 000005
上記１番目の式において，『（Ｐr q ×Ｉp ）Ｐqrp 』の部分は，
▲１▼  ｓ＋ｔ×ｐをｔ＋ｓ×ｑｒへ変換したあと，
▲２▼  ｔ＝ｕ＋ｖ×ｑをｖ＋ｕ×ｒへ変換することである。
すなわち，ｓ＋ｕ×ｐ＋ｖ×ｐｑ  →  ｖ＋ｕ×ｒ＋ｓ×ｑｒの変換を行うことになり，ｉｎ  ｐｌａｃｅな置換となる。
【００２９】
同様に，２番目の式において，『Ｐr pq（Ｐq p ×Ｉr ）』の部分は，
▲１▼  ｖ＋（ｕ＋ｓ×ｑ）×ｒをｖ＋（ｓ＋ｕ×ｐ）×ｒへ変換したあと，
▲２▼  ｖ＋（ｓ＋ｕ×ｐ）×ｒをｓ＋ｕ×ｐ＋ｖ×ｐｑへ変換することである。すなわち，ｖ＋ｕ×ｒ＋ｓ×ｑｒ  →  ｓ＋ｕ×ｐ＋ｖ×ｐｑの変換を行うことである。
【００３０】
これらはｐ＝ｒの場合，ともに入力データの領域だけを使ってｉｎ  ｐｌａｃｅに行うことができる。
ｐ＝ｒの場合を展開すると，以下の式のようになる。
0003647887.tif 000006
これにより，ｉｎ  ｐｌａｃｅでｓｅｌｆ−ｓｏｒｔｉｎｇなアルゴリズムとなる。
【００３１】
上記の（Ｃ）の式を別の表現で記述すると，以下のようになる。
ｙk ＝Σｅｘｐ（−２πｉｊｋ／ｎ）×ｘj 
（Σはｊ＝０からｎ−１までの総和，ｋ＝０，…，ｎ−１，ｉ：虚数単位）
ｎ＝ｎ0 ｎ1 ｎ2 と分解できるとき，
0003647887.tif 000007
これを次のように５段階のステップに分解して行うことができる。
【００３２】
特に，（Ｓ４）の式は，上記式（Ｄ）の  Ｄqpp Ｐp pq（Ｐq p ×Ｉp ）の部分を，対角行列の積と，添字の転置を一緒に行うことで，べクトルのロード／ストアを節約することができ，効率よく処理することができる。
ステップ（Ｓ１）
0003647887.tif 000008
ステップ（Ｓ２）
0003647887.tif 000009
ステップ（Ｓ３）
0003647887.tif 000010
ステップ（Ｓ４）
0003647887.tif 000011
ステップ（Ｓ５）
0003647887.tif 000012
【００３３】
【実施例】
図３は，本発明の実施例に用いるベクトル計算機の一例を示す。
図中，２０は主記憶装置（ＭＳＵ），３０はベクトル計算を行うベクトルユニット（ＶＵ），４０はスカラ計算を行うとともに，ベクトルユニット３０に対しベクトル命令を送るスカラユニットを表す。
【００３４】
ベクトルユニット３０において，３１−１，３１−２はロードストアパイプラインであり，ベクトル命令で指定された主記憶装置２０のデータのアクセスを処理するもの，３２はアドレス変換レジスタ，３３はマスク演算のためのマスクデータを格納するマスクレジスタ，３４は処理されるベクトルデータが格納されるベクトルレジスタ，３５−１，３５−２はマスクパイプラインであり，マスクデータを対象として演算を実行するためのもの，３６−１，３６−２はユニバーサルパイプラインであり，乗算命令，加算命令，乗算・加算命令および論理演算命令を実行する汎用の演算パイプライン，３７は除算命令を実行するための除算パイプライン，３８はスカラユニット４０からの指示によりベクトル命令の実行制御を行うベクトル命令制御部を表す。
【００３５】
また，スカラユニット４０において，４１は主記憶装置２０からの命令およびデータを蓄えるバッファストレジ（キャッシュメモリ），４２は仮想アドレスと実アドレスの変換を行うアドレス変換部，４３はレジスタファイル，４５はアクセスするアドレスを生成するアドレス生成部，４５は命令を実行する演算部，４６はスカラ命令とベクトル命令とを切り分けて命令の制御を行う命令制御／プロセス制御部を表す。
【００３６】
ベクトルレジスタ３４は，例えば０から２５５までの番号を持つ最大２５６個で構成され，大量のベクトルデータを保持する。ベクトルロード／ベクトルストア命令の実行と，ベクトル演算命令の実行とを独立に行うことができる。ベクトルレジスタ３４は，プログラムの指定でベクトルデータを格納できるようになっており，種々のベクトル演算命令に対して最適なベクトルレジスタ３４の構成を選択できるようになっている。この種のベクトル計算機は広く知られているので，これ以上の詳しい説明は省略するが，一般に一定の間隔でデータをアクセスするよりも，連続したデータをアクセスして処理するほうが高速に処理することができる。
【００３７】
図４および図５は，本発明の一実施例であって，次の複素高速フーリエ変換を行う処理フローを示している。
αk ＝Σａj ｅｘｐ（SN２πｉｊｋ／Ｎ）
ただし，Σはｊ＝０からＮ−１までの総和，SNは符号であり１または−１である。ｉは虚数単位である。α＝ＷN ・ａとも記述できる。
【００３８】
図４は，処理７，処理９の変換の符号を
ＷN （ｋ，ｊ）＝ωjk    ｊ，ｋ＝０，・・・Ｎ−１
ω＝ｅｘｐ（SNｓ２πｉ／Ｎ）でｓ＝１
としたときの処理フローである。
【００３９】
また，図５は，処理７，処理９の変換の符号を
ＷN （ｋ，ｊ）＝ωjk    ｊ，ｋ＝０，・・・Ｎ−１
ω＝ｅｘｐ（SNｓ２πｉ／Ｎ）でｓ＝−１
としたときの処理フローである。
【００４０】
以下，図４および図５に示す処理１〜処理１１に従って，本発明の実施例を詳細に説明する。
〔処理１〕
Ｎをｐi ，ｑi （ｉ＝１，…，ｋ）に分解する。すなわち，
Ｎ＝ｐ1 2 ｑ1 ，ｐ1 ＝ｐ2 2 ｑ2 ，……，ｐk-1 ＝ｐk 2 ｑk   と分解する。
【００４１】
〔処理２〕
ｖ1 ＝ｐ1 ｑ1 ，ｖ2 ＝ｖ1 ＊ｐ2 ｑ2 ，……，ｖk ＝ｖk-1 ＊ｐk ｑk 
を計算する。
【００４２】
Ｎ＝ｖk ＊ｐk 
ｐk が充分小さくなるまで分解する。
〔処理３〕
ｊ＝１に初期化し，ｊ＞２K+1 −１になるまで以下の処理を繰り返す。
【００４３】
〔処理４〕
ｍｏｄ（ｊ，２i ）＝０なる最大のｉを見つける。
〔処理５〕
ｉが０なら処理９へ進み，０でないなら処理６〜処理８を実行する。なお，処理６を行ったときは処理８を実行し，処理６Ａを行ったときは処理８Ａを実行する。
【００４４】
〔処理６〕
ｈ＝ｋ＋１−ｉ番目の処理１での分解に対応した回転１と転置／変換２の処理を以下のように行う。
【００４５】
ｈが１でないときのベクトル長は，ｖh-1 となり，ｖh-1 ×ｐh で，ｐh-1 ＝ｐh 2 ｑh となる。
回転を表す行列  Ｄqhph×Ｉvh-1*ph を次式（６．１）で計算する。
【００４６】
Ｄqhph＝ｄｉａｇ（ωj sm）    ｊ＝ｓｑh ＋ｍ    式（６．１）  なる対角行列である。
ｊは０からｐh ＊ｑh −１まで動く。ｍは０からｑh ─１，ｓは０からｐh ─１まで動く。
【００４７】
これらの対角要素を１次元ベクトルに格納して，図６（Ａ）に示すようなルーチンでベクトル処理する。
ベクトル長は，ｖh-1 ＊ｐh となる。
【００４８】
ω＝ｅｘｐ（SN２πｉ／（ｐh ＊ｑh ））として，ＲR ，Ｒi に式（６．１）から計算した，実部，虚部を１次元配列に格納する。出力は入力と同じ領域に作成される。
【００４９】
図６（Ａ）に示すルーチンの呼び出し文は，以下のとおりである。
ＣＡＬＬ  ＲＯＴ（Ａr ，Ａi ，ｖh-1 ＊ｐh ，ｐh ＊ｑh ，ＲR ，Ｒi ）
〔処理６Ａ〕
Ｄph*qh ph×ＩVh-1を同様にルーチンＲＯＴを呼んで行う。
【００５０】
ω＝ｅｘｐ（SN２πｉ／（ｐh ＊ｑh ＊ｐh ））
式（６．１）より，Ｄph*qh phを計算して，ＲR ，Ｒi に実部，虚部を１次元配列として格納する。
【００５１】
このときの図６（Ａ）に示す文の呼び出し文は，以下のとおりである。
ＣＡＬＬ  ＲＯＴ（Ａr ，Ａi ，ｖh-1 ，ｐh ＊ｑh ＊ｐh ，ＲR ，Ｒi ）
〔処理７〕
転置と変換について説明する。この順序はどちらが先でも同じである。このため，この転置は，変換２または回転１，回転２の処理の中に結合することができる。この結果，ベクトルとしてのロード／ストア命令を減じることができる。
【００５２】
ｈが１でないときのベクトル長は，ｖh-1 となり，この部分での転置を行う。Ａ（ｖh-1 ，ｐh ，ｑh ，ｐh ）の配列で２次元目と４次元目を入れ換える。このとき図６（Ｂ）に示すようなサブプログラムとなる。最内ループはベクトル処理される。出力は入力と同じ領域に作成される。
【００５３】
同時に，Ａ（ｖh-1 ，ｐh ，ｑh ，ｐh ）の変換はｑh について行う。これは，
αijkm＝Σａijnmωkn      ｎ，ｋ＝０，……，ｑh −１
（ただし，Σはｎ＝０からｑh −１までの総和）
ここで，ω＝ｅｘｐ（SNｓ２πｉ／ｑh ），符号ｓはアルゴリズムにより決められたもので計算する。
【００５４】
このとき，ｑh が１なら，変換は不要となる。
ｑh が１以外のときは，図６（Ｂ）に示すルーチンＴＲＰＳで，ｉに関するループを展開したものを各ｑh に対して作る。
【００５５】
〔処理８〕
最後に，Ｄph*qh ph×ＩVh-1を同様に図６（Ａ）に示すルーチンＲＯＴを呼んで行う。
【００５６】
ω＝ｅｘｐ（SN２πｉ／（ｐh ＊ｑh ＊ｐh ））
式（６．１）より，Ｄph*qh phを計算して，ＲR ，Ｒi に実部，虚部を１次元配列として格納する。
【００５７】
このときの図６（Ａ）に示すルーチンＲＯＴの呼び出し文は，次のようになる。
ＣＡＬＬ  ＲＯＴ（Ａr ，Ａi ，ｖh-1 ，ｐh ＊ｑh ＊ｐh ，ＲR ，Ｒi ）
〔処理８Ａ〕
Ｄqhph×Ｉvh-1*ph の計算を同様に図６（Ａ）に示すルーチンＲＯＴを呼んで行う。
【００５８】
ω＝ｅｘｐ（SN２πｉ／（ｐh ＊ｑh ））
式（６．１）より，Ｄphphを計算して，ＲR ，Ｒi に実部，虚部を１次元配列として格納する。
【００５９】
図６（Ａ）に示すルーチンＲＯＴの呼び出し文は，以下のとおりである。
ＣＡＬＬ  ＲＯＴ（Ａr ，Ａi ，ｖh-1 ＊ｐh ，ｐh ＊ｑh ，ＲR ，Ｒi ）
〔処理９〕
Ｗpk×ＩVkの中に現れる×は行列のクロネッカー積を表す。
【００６０】
クロネッカー積は，Ａが（ｐ，ｐ）の行列でＢが（ｑ，ｑ）の行列であり，添字は０，・・・，ｐ−１と０，・・・，ｑ−１の範囲である。
Ａ×Ｂは（ａijＢ）を表す新たな行列で，各要素が行列Ｂにａijを掛けた行列を要素に持つ。すなわち，大きさが（ｐｑ，ｐｑ）の行列である。
【００６１】
Ａ×Ｂtq+u,rq+s ＝Ａt,r ・Ｂu,s となる。
Ａ×Ｂの行列の添字は，０，１，・・・，ｐｑ−１である。
長さＮのベクトルを長さｖk の連続ベクトルに，ｐk 等分して，長さｐk のフーリエ変換をｖk 組行うことに相当する。
【００６２】
αik＝Σａij・ωkj  （ただし，Σはｊ＝０からｐｋ−１までの総和）
ω＝ｅｘｐ（SNｓ２πｉ／ｐk ），ｓは符号で各アルゴリズムで決まった値，
ｉ＝０，・・・，ｖk −１。
この部分を連続アクセスし，ベクトル処理して上記計算を行う。
【００６３】
この処理部はｐk の大きさに応じたものを持たせる。
ｐk が２のべき数の場合，つまり２，４，８等の大きさに応じて計算する。その処理ルーチンの例を図６（Ｃ）に示す。２に対して説明すると，Ｎ＝２のＦＦＴを１次域を利用し，ｉｎ  ｐｌａｃｅな計算部を作成して  各要素をベクトルとみて行う装置で計算する。４，８となった場合も同様に行う。
【００６４】
出力は入力と同じ領域に作成される。
上記処理７において，連続ベクトルでベクトル長はｖk の処理となる。従来の方法が等間隔のデータのアクセスやアクセスする要素の位置をインデックスのリストとして持ちベクトル処理するリストベクトル処理に比べ連続アクセス中心になる結果，従来の方法に比べ約２倍高速である。
【００６５】
上記処理における回転１，回転２，変換１〜変換３，転置などの処理の意味をさらに詳しく説明すると以下のとおりである。
［回転１］
回転を表す対角行列Ｄqhph×Ｉvh-1*ph （×はクロネッカー積）による回転を行うことは，ベクトルデータｂm （ｍ＝０，…，Ｎ−１）を長さｖh-1 ＊ｐh の連続ベクトルデータに，ｐh ＊ｑh 等分して，次の長さｐh ＊ｑh の回転をｖh-1*ph組行うことに相当する。ａij＝ｂu （ただし，ｕ＝ｊ＊ｖh-1*ph＋ｉ，ｉ＝１，…，ｖh-1*ph−１，ｊ＝０，…，ｐh ＊ｑh −１）である。結果は同じ長さのベクトルデータβu となる。
【００６６】
αij＝ωj sm×ａij    ただし，ｊ＝ｓｑh ＊ｐh ＋ｍである。
ｊは０からｐh ＊ｑh −１まで動き，ｓは０からｑh −１，ｍは０からｐh −１まで動く。ω＝ｅｘｐ（SNｓ２πｉ／（ｐh ＊ｑh ））である。ｅｘｐ（）中のｓは符号で各アルゴリズムで決まった値であり，これを＋１か−１に変えることで正変換か逆変換にすることができる。Ｎ＝ｖh-1 ＊ｐh ＊ｑh ＊ｐh である。
【００６７】
［回転２］
回転を表す対角行列Ｄph*qh ph×Ｉvh-1（×はクロネッカー積）による回転を行うことは，ベクトルデータｂm （ｍ＝０，…，Ｎ−１）を長さｖh-1 の連続ベクトルデータに，ｐh ＊ｑh ＊ｐh 等分して，次の長さｐh ＊ｑh ＊ｐh の回転をｖh-1 組行うことに相当する。ａij＝ｂu （ただし，ｕ＝ｊ＊ｖh-1 ＋ｉ，ｉ＝１，…，ｖh-1 −１，ｊ＝０，…，ｐh ＊ｑh ＊ｐh −１）である。結果は同じ長さのベクトルデータβu となる。
【００６８】
αij＝ωj sm×ａij    ただし，ｊ＝ｓｑh ＋ｍである。
ｊは０からｐh ＊ｑh ＊ｐh −１まで動き，ｓは０からｐh −１，ｍは０からｐh ＊ｑh −１まで動く。ω＝ｅｘｐ（SNｓ２πｉ／（ｐh ＊ｑh ＊ｐh ））である。ｓは符号で各アルゴリズムで決まった値で，これを＋１か−１に変えることで正変換か逆変換にすることができる。Ｎ＝ｖh-1 ＊ｐh ＊ｑh ＊ｐh である。
【００６９】
［変換２］
この変換は，行列Ｉph×Ｗqh×Ｉvh-1*ph （×はクロネッカー積）で表される。フーリエ変換は，変換するべきベクトルデータｂm （ｍ＝０，…，Ｎ−１）を長さｖh-1 ＊ｐh の連続ベクトルデータに，ｐh ＊ｑh 等分して，ｐh ＊ｑh 組をさらにｐh 等分し，次の長さｑh のフーリエ変換をｖh-1 ＊ｐh ＊ｐh 組行うことに相当する。ａijm ＝ｂu （ただし，ｕ＝ｍ＊ｖh-1 ＊ｐh ＊ｑh ＋ｊ＊ｖh-1 ＊ｐh ＋１，ｉ＝１，…，ｖh-1 ＊ｐh −１，ｊ＝０，…，ｑk −１，ｍ＝０，…，ｐk −１）である。結果は同じ長さのベクトルデータβt となる（ｔ＝０，…，Ｎ−１）。
【００７０】
αikm ＝Σａijm ・ωkj
ただし，Σはｊ＝０からｑk −１までの総和，ｋ＝０，…，ｑk −１。
ω＝ｅｘｐ（SNｓ２πｉ／ｐk ）
ｓは符号で各アルゴリズムで決まった値で，これを＋１か−１に変えることで正変換か逆変換にすることができる。Ｎ＝ｖh-1 ＊ｐh ＊ｑh ＊ｐh である。
【００７１】
［転置］
転置は，ベクトルデータｂm （ｍ＝０，…，Ｎ−１）を長さｖh-1 の連続ベクトルデータに，ｐh ＊ｑh ＊ｐh 等分して，さらにｐh ＊ｑh 等分し，これをｑh 等分した組に対して，
ａijkm＝ｂu （ただし，ｕ＝ｍ＊ｖh-1 ＊ｐh ＊ｑh ＋ｋ＊ｖh-1 ＊ｐh ＋ｊ＊ｖh-1 ＋１，ｉ＝０，…，ｖh-1 −１，ｊ＝０，…，ｐh −１，ｋ＝０，…，ｑh −１，ｍ＝０，…，ｐh −１）としたとき，αimkj＝ａijkmとすることをいう。Ｎ＝ｖh-1 ＊ｐh ＊ｑh ＊ｐh である。
【００７２】
［変換１および変換３］
この変換は，行列Ｗpk×Ｉvk（×はクロネッカー積）で表される。フーリエ変換は，変換するべきベクトルデータｂm （ｍ＝０，…，Ｎ−１）を長さｖk の連続ベクトルデータに，ｐk 等分して，次の長さｐk のフーリエ変換をｖk 組行うことに相当する。ａij＝ｂu （ただし，ｕ＝ｊ＊ｖk ＋１，ｉ＝１，…，ｖk −１，ｊ＝０，…，ｐk −１）である。結果は同じ長さのベクトルデータβt となる（ｔ＝０，…，Ｎ−１）。
【００７３】
αik＝Σａij・ωkj
ただし，Σはｊ＝０からｐk −１までの総和，ｉ＝０，…，ｖk −１。
ω＝ｅｘｐ（SNｓ２πｉ／ｐk ），ｓは符号で各アルゴリズムで決まった値で，これを＋１か−１に変えることで正変換か逆変換にすることができる。Ｎ＝ｖk ＊ｐh である。
【００７４】
【発明の効果】
以上説明したように，本発明によれば，連続アクセス中心のベクトル処理が可能になり，連続アクセスが高速なベクトル計算機の処理性能を向上させることができる。また，インプレースかつセルフ・ソーティングで作業領域をほとんど必要としない処理が可能になる。
【図面の簡単な説明】
【図１】本発明の原理説明図である。
【図２】本発明による分解例説明図である。
【図３】本発明の実施に用いるベクトル計算機の例を示す図である。
【図４】本発明の実施例処理フローチャートである。
【図５】本発明の実施例処理フローチャートである。
【図６】本発明の実施例で用いる部品化したプログラムの例を示す図である。
【図７】従来技術の例を説明するための図である。
【符号の説明】
１０    ベクトル計算機
１１    メモリ
【図１】
0003647887.tif 000013
【図２】
0003647887.tif 000014
【図３】
0003647887.tif 000015
【図４】
0003647887.tif 000016
【図５】
0003647887.tif 000017
【図６】
0003647887.tif 000018
【図７】
0003647887.tif 000019