山口屋~活動日誌~

私生活で主な出来事をピックアップ

フーリエ変換 FFT 基数 Rader Winograd Bluestein chirp z-transform

2023-11-23 06:10:26 | ソフトウェア開発
FFTの基数は2でなければならないわけではなく3,4,5など他の数であってもよい。
○Ooura's Mathematical Software Packages:FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法
FFT内でよく登場する複素数乗算について、通常の実数乗算4回と実数加算2回で行うのではなく、変形して実数乗算3回と実数加算3回で行う方法も知られている。
○Ooura's Mathematical Software Packages:FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法

基数3,4,5のアルゴリズム例は下記の文献に記載がある。
○高橋大介、金田康正:分散メモリ型並列計算機による2,3,5基底一次元FFTの実現と評価、情報処理学会論文誌、Vol.39、No.3、pp.519-528、1998→情報処理学会(PDF)
基数8のアルゴリズム例でわかりやすい文献はなかなか見つからず、未確認だが載っていると思われる文献も含める。
○高橋大介、金田康正:積和演算命令に向いた8基底FFTカーネルの提案、情報処理学会論文誌、Vol.41、No.7、pp.2019-2026、2000→情報処理学会(PDF)
○ S. Goedecker : Fast Radix 2, 3, 4, and 5 Kernels for Fast Fourier Transformations on Computers with Overlapping Multiply - Add Instructions, SIAM Journal on Scientific Computing, Vol.18, No.6, pp.1605-1611, 1997 →非公式?(PDF)
○Charles Van Loan : Computational Frameworks for the Fast Fourier Transform → Society for Industrial and Applied Mathematics
○Daisuke Takahashi : Fast Fourier Transform Algorithms for Parallel Computers → Springer

発展の過程にはいろいろ歴史があるようだが、まだまだ調査中。
●基数2,4の組み合わせ
○W. M. Gentleman, G. Sande : Fast Fourier Transforms for fun and profit, Proceedings of 1966 Fall Joint Computer Conference, American Federation of Information Processing Societies Conference Proceedings, Vol.29, pp.563-578, 1966
●基数8
○G. D. Bergland : A fast Fourier transform algorithm using base 8 iterations, Mathematics of Computation, Vol.22, No.102, pp.275-279, 1968→American Mathematical Society(PDF)
●PFA:Prime Factor Algorithm
○I. J. Good : The Interaction Algorithm and Practical Fourier Analysis, Journal of the Royal Statistical Society Series B, Vol. 20, pp. 361-372, 1958
○I. J. Good : The Interaction Algorithm and Practical Fourier Analysis, Journal of the Royal Statistical Society Series B, Vol. 22, pp. 372-375, 1960
○S. Winograd : On computing the Discrete Fourier transform, Mathematics of Computation, Vol.32, No.141, pp.175-199, 1978→American Mathematical Society(PDF)
→Prime Factor Algorithm について、ルーチンを工夫して乗算回数を減らすものらしい。
○中国剰余定理(CRT:Chinese Remainder Theorem)
→Prime Factor Algorithm について、分解後の回転因子を上手く求めるものっぽい?
○大野豊、磯田和男:数値計算ハンドブック、pp.929-945、1990、オーム社
○Qiita:Prime-Factor FFT
→Prime Factor Algorithm について、工夫を終えたあとのもの。
○鳥居達生、杉浦洋:基数2のFFTに基づく任意項数の離散型Fourier変換、情報処理学会論文誌、Vol.25、No.1、pp.30-36、1984→情報処理学会(PDF)
○Qiita:任意要素数の高速フーリエ変換
→1次元フーリエ変換をゼロ埋めで2のべき乗のフーリエ変換に書き換える?配列の生成を伴うため、要素数がある程度大きい場合に適用するとのこと。
○鳥居達生:離散形フーリエ変換の高速算法の動向、計測と制御、Vol.25、No.12、pp.1067-1073、1986→J-STAGE(PDF)

Numeric Theory Translation:数論変換(整数環高速フーリエ変換)、Fast Molulo Transformation:高速剰余変換、とかも仲間なのかな?
○Senの競技プログラミング備忘録:NTT(数論変換)のやさしい解説
○後 保範(Ushiro Yasunori)のホームページ:FMT(高速剰余変換)の原理

任意の要素数のFFTとしては、すでに chirp z-transform あるいは Bluestein アルゴリズム と呼ばれているものがあったようだ。
○Wikipedia(英語版):Chirp Z-transform
○L. R. Rabiner, R. W. Schafer, and C. M. Rader : The chirp z-transform algorithm and its application, Bell System Technical Journal, Vol.48, No.5, pp.1249-1292, 1969, also published in IEEE Tr. Audio & Electroacoustics, vol. 17, no. 2, pp. 86-92, 1969
○L. Bluestein : A linear filtering approach to the computation of discrete Fourier transform, IEEE Transactions on Audio and Electroacoustics, Vol.18, No.4, pp.451-455, 1970
○T. T. Wang : The segmented chirp z-transform and its application in spectrum analysis, IEEE Transactions on Instrumentation and Measurement, Vol.39, No.2, pp.318–323, 1990
○腰も砕けよ 膝も折れよ:Chirp z変換によるFFT
○記:FFT (Chirp z-変換アルゴリズム)

なお、z-transform の逆変換は50年近く謎のままだったが最近になって解明されたらしい。
○Vladimir Sukhoy & Alexander Stoytchev : Generalizing the inverse FFT off the unit circle, Scientific Reports, 2019→nature.com(PDF)
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

三角関数 n倍角 Chebyshev多項式 漸化式 FFT

2023-11-23 05:57:15 | ソフトウェア開発
Chebyshev多項式を用いると、三角関数のn倍角を漸化式で表すことができる。三角関数のn倍角を順に求めていく計算では、三角関数を計算するルーチンの呼出は最初の1回だけで、残りは漸化式で逐次的に数値計算することができる。
一松信:初等関数の数値計算、pp.129-132、1974、教育出版
上記文献によれば元ネタは下記文献らしい。
C. W. Clenshaw : A note on the summation of Chebyshev series, Mathematical Tables and other Aids to Computation, Vol.9, No.51, pp.118-120, 1955

Chebyshev多項式と漸化式は、例えば下記を参照。
潮田康夫:チェビシェフの多項式とn倍角の公式,数研通信,No.69,pp.8-11,2011→数研出版(PDF)
cosのn倍角はcosの多項式で表されるが、sinのn倍角はcosの多項式とsinの積で表される点に注意。

Chebyshev多項式を用いた漸化式で、三角関数のn倍角を逐次求めるC言語のソースコードは、FFTの一部として下記に記載がある。
奥村晴彦:C言語による標準アルゴリズム事典、pp.346-348、2018、技術評論社
Qiita:技術評論社 Software Technology シリーズ一覧
ソースコードに解説が無く難解だったため、当該部分のみ引用して補足させてもらう。

(1) 三角関数を計算してFFTの回転因子を求める等、下準備をする。
t = sin(PI / n);
dc = 2 * t * t; // = 1 - cos (2 * PI / n) ※0の近傍での桁落ちを避けるため、2 * sin(PI / n) * sin(PI / n) で計算している。
ds = sqrt(dc * (2 - dc)); // = sin (2 * PI / n)
t = 2 * dc;

(2) k倍角(c[k], s[k])を漸化式で求め、次の(k+1)倍角を求めるための下準備をする。
c[k] = c[k-1] - dc; dc += t * c[k]; // 更新前の dc は -c[k] + c[k-1] の意味合い
s[k] = s[k-1] + ds; ds -= t * s[k]; // 更新前の ds は s[k] - s[k-1] の意味合い

漸化式を用いた計算では桁落ちに注意しなければならない場合もある。
吉田年雄:数値解析の基礎・基本、理工系数学の基礎・基本7、pp.125-127、2005、牧野書店(廃業)
特性方程式の解が初期値となる cos の計算は桁落ちが起きていないか心配であったが、問題ないと解説されている。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする