RustでFFT (1)基本のFFT

フーリエ変換とは

周期関数を三角関数の合成で表す変換。次の式で変換できる。

$$\omega = 2 \pi f $$
$$F(\omega) = \int_{-\infty}^{\infty} f(t)\ e^{- i \omega t},dt$$

また、逆変換も次の式で変換できる。
$$f(t) = {\frac 1 {2 \pi}} \int_{-\infty}^{\infty} F(\omega)\ e^{i \omega t},d\omega$$

なお、正規化係数（上記のフーリエ変換は$1$、逆変換は$\frac {1}{2 \pi}$）はその積が$\frac 1 {2 \pi}$になれば良く、その係数の取り方で別の定義も存在する。

下のページなどに詳しく載っているので参照。
http://www.yukisako.xyz/entry/fourier-transform
http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/fs.html
https://ja.wikipedia.org/wiki/フーリエ変換

離散フーリエ変換

現実世界の時間は連続しているため一定の時間で区切る必要がある。また、値の大きさも一定の間隔で区切る必要がある。

上のフーリエ変換を離散化した式は次の通りとなる。

変換
$$\omega = \frac {2 \pi t}{N} $$
$$F(t)=\sum _{x=0}^{N-1}f(x)e^{-i \omega x}$$

逆変換
$$f(x)={\frac 1 N}\sum _{t=0}^{N-1}F(t)e^{i \omega x}$$

離散フーリエ変換も同様に正規化係数の積が$\frac 1 N$になれば良いため、別の定義もある。

使い道

フーリエ変換及びその変形の使い道は多岐に渡る。例を次に示す。

• 画像、音声、動画の圧縮
• 音声解析
• 多倍長積算

実装の前に

複素数演算について

フーリエ変換は複素数演算である。
Rustは標準では複素数をサポートしない。
したがって、次の選択肢を検討する。

• 外部クレートの使用
• 自前実装

外部クレートを使用する場合、次の大きなリスクがある。

• 悪意のあるコードの混入
• 非公開化

いずれも、JavaScript界隈のnpmjs.comで実際に起こっている。
なお、クレート自身に問題がなくとも、そのクレートが使用しているクレートにも同様の問題が発生する可能性がある為、リスク削減の為には、クレートの依存関係全てをチェックする必要がある。
今回使用を検討しているクレートはnumであり、そのAuthorsはThe Rust Project Developersとなっていることから、上記のリスクは少ないと考えられる。
また、仕様変更されるリスクもあるが、Rust等の新興の言語を使う場合においては享受すべきだろう。

トレイト

テスト・パフォーマンス計測を簡便にする為、以下のトレイトを定義した。
コードはこのトレイトに対して実装していく。

extern crate num_traits;
extern crate num;

use num_traits::float::Float;
use num::Complex;

pub trait AlgorithmName {
    fn name(&self) -> &'static str;
}

pub trait DftAlgorithm<T: Float>: AlgorithmName {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>>;
    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>>;
}

定義式通りの実装

離散フーリエ変換の定義式通りに実装してみる。

use num::{Complex, Float, cast, one};
use num_traits::float::FloatConst;
use std::marker::PhantomData;

pub struct BasicDft<T> {
    phantom: PhantomData<T>,
}

impl<T> BasicDft<T> {
    pub fn new() -> Self {
        Self { phantom: PhantomData }
    }
}

impl<T> AlgorithmName for BasicDft<T> {
    fn name(&self) -> &'static str {
        "定義"
    }
}

impl<T: Float + FloatConst> DftAlgorithm<T> for BasicDft<T> {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        return (0..source.len())
            .map(|i| {
                (1..source.len()).fold(source[0], |x, j| {
                    x +
                        source[j] *
                            Complex::<T>::from_polar(
                                &one(),
                                &(-cast::<_, T>(2 * i * j).unwrap() * T::PI() /
                                      cast(source.len()).unwrap()),
                            )
                })
            })
            .collect::<Vec<_>>();
    }

    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        return (0..source.len())
            .map(|i| {
                (1..source.len())
                    .fold(source[0], |x, j| {
                        x +
                            source[j] *
                                Complex::<T>::from_polar(
                                    &one(),
                                    &(cast::<_, T>(2 * i * j).unwrap() * T::PI() /
                                          cast(source.len()).unwrap()),
                                )
                    })
                    .unscale(cast(source.len()).unwrap())
            })
            .collect::<Vec<_>>();
    }
}

$O(n^2)$の為、要素数が多くなった場合、当然の事ながら、かなり時間がかかる。
したがって、普通は巡回畳み込みを利用したFFTが使用される。
ただ、FFTを実装するにあたり、定義式通りに実装されたこのコードと比較する事により、
確実なテストが可能である。

基本のFFT(2基底)

FFTの概略については以下のページ等を参照
http://zakii.la.coocan.jp/fourie/0_contents.htm
http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/index.html
FFTを実装するにあたっては、京大の大浦助教のサイトは必読である。

以下では要素数が2のべき乗である場合におけるコードを示す。

Cooley-Tukey 型 FFT

Cooley-Tukey 型の特徴は２つ

• In-place演算である（メリット）
• ビットリバース処理が必要である（デメリット）

周波数引き

use num::{Complex, Float, cast, one};
use num_traits::float::FloatConst;
use std::cmp;

pub struct CooleyTukeyDif<T> {
    level: usize,
    ids: Vec<usize>,
    omega: Vec<Complex<T>>,
}

impl<T: Float + FloatConst> CooleyTukeyDif<T> {
    pub fn new() -> Self {
        Self {
            level: usize::max_value(),
            ids: Vec::new(),
            omega: Vec::new(),
        }
    }

    fn initialize(&mut self, len: usize) {
        let lv = len.trailing_zeros() as usize;
        if len != 1 << lv {
            return;
        }
        if self.level != usize::max_value() {
            self.ids = Vec::new();
            self.omega = Vec::new();
        }
        self.level = lv;

        // ビットリバースの計算
        self.ids.push(0);
        for _ in 0..lv {
            for j in 0..self.ids.len() {
                let id = self.ids[j] << 1;
                self.ids[j] = id;
                self.ids.push(id + 1);
            }
        }

        // ωの事前計算
        self.omega.push(one());
        let q = cmp::max(len >> 2, 1);
        let h = cmp::max(len >> 1, 1);
        // 1/4周分を計算
        for i in 1..q {
            self.omega.push(Complex::from_polar(
                &one(),
                &(cast::<_, T>(-2.0).unwrap() * T::PI() / cast(len).unwrap() *
                      cast(i).unwrap()),
            ));
        }

        // 1/4～1/2周分を計算
        for i in q..h {
            let tmp = self.omega[i - q];
            self.omega.push(Complex::new(tmp.im, -tmp.re));
        }

        // 1/2周目から計算
        for i in h..len {
            let tmp = self.omega[i - h];
            self.omega.push(Complex::new(-tmp.re, -tmp.im));
        }

        // 1周ちょうど
        self.omega.push(one());
    }

    fn convert_inner(&mut self, source: &[Complex<T>], is_back: bool) -> Vec<Complex<T>> {
        if self.level == usize::max_value() || source.len() != 1 << self.level {
            self.initialize(source.len());
        }

        // 入力のコピー
        let mut work = source.to_vec();

        // FFT
        let mut po2b = 1 << self.level;
        for i in 0..self.level {
            let po2m = po2b;
            po2b >>= 1;

            for k in 0..po2b {
                let w = if is_back {
                    self.omega[(1 << i) * k].conj()
                } else {
                    self.omega[(1 << i) * k]
                };
                let mut j = k;
                while j < work.len() {
                    let pos_b = j + po2b;
                    let wfa = work[j];
                    let wfb = work[pos_b];
                    work[j] = wfa + wfb;
                    work[pos_b] = (wfa - wfb) * w;
                    j += po2m;
                }
            }
        }

        // 並び替え
        return self.ids
            .iter()
            .map(|&i| if is_back {
                work[i].unscale(cast(work.len()).unwrap())
            } else {
                work[i]
            })
            .collect::<Vec<_>>();
    }
}

impl<T: Float + FloatConst> DftAlgorithm<T> for CooleyTukeyDif<T> {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, false)
    }

    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, true)
    }
}

impl<T> AlgorithmName for CooleyTukeyDif<T> {
    fn name(&self) -> &'static str {
        "Cooley Tukey - DIF"
    }
}

高速化

いくつか高速化の為にテクニックを利用している。

ωの事前計算

現在のCPUで実行する場合、三角関数計算のみである場合、テーブル引きのメリットはほとんどない。
ただし、FFTの場合、1つの角度につき、少なく見積もっても三角関数・乗算がそれぞれ2回必要となる。
CPUのキャッシュメモリサイズによるが、1万要素程度であれば、テーブル引きの方が速い。
また、三角関数の計算も$0$~$\frac \pi 2$の領域のみで充分であり、残りは符号の変換のみで計算できる。

ビットリバース処理と正規化の同時処理

逆フーリエ変換の場合、正規化係数を掛ける必要がある。
ビットリバース処理と別々に行なった場合より、メモリ転送処理が少なくなる為、可能であれば同時に行なった方が良い。

ビット演算の活用

乗除演算、特に除算は非常に遅い。
2のべき乗での乗除は、ビットシフトに置き換えるべきである。

評価

• 入力を変化させたくない場合に限るが、最初のデータ転送が遅い。ただし、このデメリットはアルゴリズムを変形する事で対応できる。
• このコードは高速化の余地はまだあるが、劇的な改善は見込めず、後に示す4基底が大幅に速い為、改良を行わない。

時間引き

use num::{Complex, Float, cast, one};
use num_traits::float::FloatConst;
use std::cmp;
use num_traits::NumAssign;

pub struct CooleyTukeyDit<T> {
    level: usize,
    ids: Vec<usize>,
    omega: Vec<Complex<T>>,
}

impl<T: Float + FloatConst + NumAssign> CooleyTukeyDit<T> {
    pub fn new() -> Self {
        Self {
            level: usize::max_value(),
            ids: Vec::new(),
            omega: Vec::new(),
        }
    }

    fn initialize(&mut self, len: usize) {
        let lv = len.trailing_zeros() as usize;
        if len != 1 << lv {
            return;
        }
        if self.level != usize::max_value() {
            self.ids = Vec::new();
            self.omega = Vec::new();
        }
        self.level = lv;

        // ビットリバースの計算
        self.ids.push(0);
        for _ in 0..lv {
            for j in 0..self.ids.len() {
                let id = self.ids[j] << 1;
                self.ids[j] = id;
                self.ids.push(id + 1);
            }
        }

        // ωの事前計算
        self.omega.push(one());
        let q = cmp::max(len >> 2, 1);
        let h = cmp::max(len >> 1, 1);
        // 1/4周分を計算
        for i in 1..q {
            self.omega.push(Complex::from_polar(
                &one(),
                &(cast::<_, T>(-2.0).unwrap() * T::PI() / cast(len).unwrap() *
                      cast(i).unwrap()),
            ));
        }

        // 1/4～1/2周分を計算
        for i in q..h {
            let tmp = self.omega[i - q];
            self.omega.push(Complex::new(tmp.im, -tmp.re));
        }

        // 1/2周目から計算
        for i in h..len {
            let tmp = self.omega[i - h];
            self.omega.push(Complex::new(-tmp.re, -tmp.im));
        }

        // 1周ちょうど
        self.omega.push(one());
    }

    fn convert_inner(&mut self, source: &[Complex<T>], is_back: bool) -> Vec<Complex<T>> {
        if self.level == usize::max_value() || source.len() != 1 << self.level {
            self.initialize(source.len());
        }

        // 入力の並び替え
        let mut ret = self.ids.iter().map(|&i| if is_back {
            source[i].unscale(cast(source.len()).unwrap()) // 逆変換の場合はこのタイミングで割り戻しておく
        } else {
            source[i]
        }).collect::<Vec<_>>();
        // FFT
        let mut po2 = 1;
        for i in 1..(self.level + 1) {
            let po2m = po2;
            po2 <<= 1;

            for k in 0..po2m {
                let w = if is_back {
                    self.omega[(1 << (self.level - i)) * k].conj()
                } else {
                    self.omega[(1 << (self.level - i)) * k]
                };
                let mut j = k;
                while j < ret.len() {
                    let pos_b = j + po2m;
                    let wfb = ret[pos_b] * w;
                    ret[pos_b] = ret[j] - wfb;
                    ret[j] += wfb;
                    j += po2;
                }
            }
        }
        return ret;
    }
}

impl<T> AlgorithmName for CooleyTukeyDit<T> {
    fn name(&self) -> &'static str {
        "Cooley Tukey - DIT"
    }
}

impl<T: Float + FloatConst + NumAssign> DftAlgorithm<T> for CooleyTukeyDit<T> {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, false)
    }

    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, true)
    }
}

評価

• このコードは高速化の余地はまだあるが、劇的な改善は見込めず、後に示す4基底が大幅に速い為、改良を行わない。

Stockham 型 FFT

Stockham 型の特徴は２つ。Cooley-Tukeyと真逆である。

• Out-place演算である（デメリット）
• ビットリバース処理の必要はない（メリット）

周波数引き

use num::{Complex, Float, cast, one};
use num_traits::float::FloatConst;
use std::cmp;
use std::mem;

pub struct StockhamDif<T> {
    omega: Vec<Complex<T>>,
}

impl<T: Float + FloatConst> StockhamDif<T> {
    pub fn new() -> Self {
        Self { omega: Vec::new() }
    }

    fn initialize(&mut self, len: usize) {
        if len != len.next_power_of_two() {
            return;
        }

        self.omega = Vec::with_capacity(len + 1);

        // ωの事前計算
        self.omega.push(one());
        let q = cmp::max(len >> 2, 1);
        let h = cmp::max(len >> 1, 1);
        // 1/4周分を計算
        for i in 1..q {
            self.omega.push(Complex::from_polar(
                &one(),
                &(cast::<_, T>(-2.0).unwrap() * T::PI() / cast(len).unwrap() *
                      cast(i).unwrap()),
            ));
        }

        // 1/4～1/2周分を計算
        for i in q..h {
            let tmp = self.omega[i - q];
            self.omega.push(Complex::new(tmp.im, -tmp.re));
        }

        // 1/2周目から計算
        for i in h..len {
            let tmp = self.omega[i - h];
            self.omega.push(Complex::new(-tmp.re, -tmp.im));
        }

        // 1周ちょうど
        self.omega.push(one());
    }

    fn convert_inner(&mut self, source: &[Complex<T>], is_back: bool) -> Vec<Complex<T>> {
        let len = source.len();
        if len + 1 != self.omega.len() {
            self.initialize(len);
        }
        let mut ret = source.to_vec();
        let mut tmp = ret.clone();

        // FFT
        let mut local_len = len;
        let half_len = local_len >> 1;
        let mut i = 1;
        while i < source.len() {
            local_len >>= 1;
            for j in 0..local_len {
                let wpos = i * j;
                let w = if is_back {
                    self.omega[len - wpos]
                } else {
                    self.omega[wpos]
                };
                for k in 0..i {
                    let pos = wpos + k;
                    let out_pos = pos + wpos;
                    let xbuf = ret[pos + half_len];
                    tmp[out_pos] = ret[pos] + xbuf;
                    tmp[out_pos + i] = (ret[pos] - xbuf) * w;
                }
            }
            i <<= 1;
            mem::swap(&mut tmp, &mut ret);
        }

        if is_back {
            for i in 0..ret.len() {
                ret[i] = ret[i].unscale(cast(len).unwrap());
            }
        }
        return ret;
    }
}

impl<T> AlgorithmName for StockhamDif<T> {
    fn name(&self) -> &'static str {
        "Stockham - DIF"
    }
}

impl<T: Float + FloatConst> DftAlgorithm<T> for StockhamDif<T> {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, false)
    }

    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, true)
    }
}

評価

これまでのCooley-Tukey型と比較して、上で挙げたメリットに加え、もう一つ、最内ループにおける配列添え字が連続しているといったメリットが存在している。
結果として、ある程度の要素を持つ値のFFTではCooley-Tukey型より速かった。
ただし、要素数が少ない場合は、メモリ確保のデメリットが上回り、Cooley-Tukey型の方が速い。

時間引き

use num::{Complex, Float, cast, one};
use num_traits::float::FloatConst;
use std::cmp;
use std::mem;

pub struct StockhamDit<T> {
    omega: Vec<Complex<T>>,
}

impl<T: Float + FloatConst> StockhamDit<T> {
    pub fn new() -> Self {
        Self { omega: Vec::new() }
    }

    fn initialize(&mut self, len: usize) {
        if len != len.next_power_of_two() {
            return;
        }

        self.omega = Vec::with_capacity(len + 1);

        // ωの事前計算
        self.omega.push(one());
        let q = cmp::max(len >> 2, 1);
        let h = cmp::max(len >> 1, 1);
        // 1/4周分を計算
        for i in 1..q {
            self.omega.push(Complex::from_polar(
                &one(),
                &(cast::<_, T>(-2.0).unwrap() * T::PI() / cast(len).unwrap() *
                      cast(i).unwrap()),
            ));
        }

        // 1/4～1/2周分を計算
        for i in q..h {
            let tmp = self.omega[i - q];
            self.omega.push(Complex::new(tmp.im, -tmp.re));
        }

        // 1/2周目から計算
        for i in h..len {
            let tmp = self.omega[i - h];
            self.omega.push(Complex::new(-tmp.re, -tmp.im));
        }

        // 1周ちょうど
        self.omega.push(one());
    }

    fn convert_inner(&mut self, source: &[Complex<T>], is_back: bool) -> Vec<Complex<T>> {
        if source.len() + 1 != self.omega.len() {
            self.initialize(source.len());
        }
        let mut ret = source.to_vec();
        let mut tmp = ret.clone();

        // FFT
        let mut local_len = ret.len();
        let half_len = local_len >> 1;
        let mut i = 1;
        while i < source.len() {
            local_len >>= 1;
            for j in 0..local_len {
                for k in 0..i {
                    let w = if is_back {
                        self.omega[local_len * k].conj()
                    } else {
                        self.omega[local_len * k]
                    };
                    let pos = i * j + k;
                    let out_pos = 2 * i * j + k;
                    let xbuf = ret[pos + half_len] * w;
                    tmp[out_pos] = ret[pos] + xbuf;
                    tmp[out_pos + i] = ret[pos] - xbuf;
                }
            }
            i <<= 1;
            mem::swap(&mut tmp, &mut ret);
        }

        if is_back {
            for i in 0..ret.len() {
                ret[i] = ret[i].unscale(cast(ret.len()).unwrap());
            }
        }
        return ret;
    }
}

impl<T> AlgorithmName for StockhamDit<T> {
    fn name(&self) -> &'static str {
        "Stockham - DIT"
    }
}

impl<T: Float + FloatConst> DftAlgorithm<T> for StockhamDit<T> {
    fn convert(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, false)
    }

    fn convert_back(&mut self, source: &[Complex<T>]) -> Vec<Complex<T>> {
        self.convert_inner(source, true)
    }
}

評価

こちらもCooley-Tukey型と比較して、最内ループにおける配列添え字が連続している為か、ある程度の要素を持つ値のFFTではCooley-Tukey型より速かった。
同様に、要素数が少ない場合は、メモリ確保のデメリットが上回り、Cooley-Tukey型の方が速い。

終わりに

本当は、1記事で全て載せたかったが、コードが多い為、複数に別けることとした。
近日中に続きを記載するので待ってほしい。