Go to Qiita Advent Calendar 2024 Top
LoginSignup
12
17

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@YSRKEN

大浦版FFTをC++とJavaに移植してみた

Last updated at Posted at 2015-11-02

#概要
　題名通りです。元がC言語なのでそこまで困難でもなかったのですが、修正する際に気をつけたことをメモしておきます。移植したコードを入手したい人は、次のGitHubのリンクから入手して下さい。

　　Ooura-FFT-Library-by-Other-Language

※今回移植したのはfft4g.c(FFT パッケージ in C - 高速版 (基数 4, 2))だけです。他も同じような方法で移植可能です。

#元々のソースの使い方
　端的に書くとこんな感じです。若干トリッキーですが、その性能は折り紙付きです。

sample1.cpp
#include <cmath>
#include <iostream>

#include "fft4g.c"

const int n = 4;		//データ数
const int n_sqrt = 2;	//√(データ数)

int main(){
	// 変換用バッファの初期化
	int ip[2 + n_sqrt];
	double w[n / 2];
	// 入力データを用意
	//・1 + 2i, 3 + 4i, 5 + 6i, 7 + 8i
	//・仕様上、データ数の2倍だけ用意する必要がある
	double a[] = {1, 2, 3, 4, 5, 6, 7, 8};
	if(sizeof(a) / sizeof(a[0]) != n * 2) return -1;
	// 変換開始
	ip[0] = 0;
	cdft(n * 2, 1, a, ip, w);
	// 結果を出力
	// ・4 + 5i, 0 - 2i, -1 - 1i, -2
	// ・生データだと正規化されてないので、データ数で割る必要がある
	// ・Mathematicaだと出力を√(データ数)で割らないと正規化出来ない
	for(int i = 0; i < n; i++){
		std::cout << "[" << i << "] " << (a[i * 2] / n) << " " << (a[i * 2 + 1] / n) << "i" << std::endl;
	}
}

　ここで誰もが思うのが、ipとwが邪魔だということです。「要素数nが決まればipとwに必要な領域は自動で決まる」とか「最初に使う時だけはip[0] = 0にする」とかを毎回考えるのは面倒ですよね？

#クラスとしてまとめる
　というわけで、オブジェクト指向らしく全体をクラスで包んでしまいましょう。フーリエ変換をする際は、その都度インスタンスを生成すればいいというわけです。また、オブジェクト指向言語の性質を利用して次のような利点も生まれます。

  • 初期化・破棄が楽になる(コンストラクタ・デストラクタのご利益)
  • 使用するメソッド以外は隠蔽できる(カプセル化)
  • FFTの中身が変わっても、インターフェース部分は統一できる(継承)

　まず最初に移植対象としたのはC++ですが、やることは意外に多かったです。具体的にはこんな感じ。

  • 全体をclassで包む
  • publicにするメンバ関数と、privateにするメンバ変数・関数を振り分ける
  • 各関数冒頭に付いてくるおびただしいプロトタイプ宣言を削除
  • 変数スコープを考えつつ、宣言位置を関数ブロック先頭から移動……させようとしたがパフォーマンス低下を心配して撤回
  • 配列w(cosとsinのテーブル)を用意する際、要求されるサイズがcdft・rdftとddct・ddstとdfct・dfstで異なるので、あえて大目に取っておく

　後は、メンバ関数と関数引数は言語側で自動判別してくれるので特に弄る必要はないでしょう。実際、この程度の条件で、目的のcppファイルが完成しました。

#ポインタが使えない場合に備える
　C→C++と来てその次はJava……だったのですが、その場合、修正個所が次のように追加されます。

  • 各種数学関数が、例えばsin関数だとMath.sinとなるように書き直される
  • private指定とpublic指定が各メソッドごとに掛けないとダメなのでその辺を書き直す
  • Javaだとポインタが使えないので、例えば配列aに対して「a + 2」といった表記ができない(C++だとaのアドレスを2つ分ずらした場所になる)ので、+演算子でズラす分を引数で指定するといった作戦を取る。具体的にはこんな感じ
sample2.java
/* 変更前 */
makect(nc, ip, w + nw);
private void makect(int nc, int[] ip, double[] c){
    c[nch] = 0.5 * c[0];
}
/* 変更後 */
makect(nc, ip, w, nw);
private void makect(int nc, int[] ip, double[] c, int nw){
    c[nch + nw] = 0.5 * c[0 + nw];
}

#テストコードで検証
　一通り移植が終わった後、テストコードで検証しました。元々付属していたテストコードであるtestxg.cをC++とJavaに移植してから検証します。

testxg.cpp
#include "fft4g.hpp"
#include <cstdio>
#define MAX(x,y) ((x) > (y) ? (x) : (y))

/* random number generator, 0 <= RND < 1 */
#define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200.0))

#ifndef NMAX
#define NMAX 8192
#define NMAXSQRT 64
#endif

void putdata(int nini, int nend, double *a);
double errorcheck(int nini, int nend, double scale, double *a);

int main()
{
    int n, ip[NMAXSQRT + 2];
    double a[NMAX + 1], t[NMAX / 2 + 1], err;

    printf("data length n=? (must be 2^m)\n");
    scanf("%d", &n);
    ip[0] = 0;

	fft4g fft(n);

    /* check of CDFT */
    putdata(0, n - 1, a);
    fft.cdft(1, a);
    fft.cdft(-1, a);
    err = errorcheck(0, n - 1, 2.0 / n, a);
    printf("cdft err= %g \n", err);

    /* check of RDFT */
    putdata(0, n - 1, a);
    fft.rdft(1, a);
    fft.rdft(-1, a);
    err = errorcheck(0, n - 1, 2.0 / n, a);
    printf("rdft err= %g \n", err);

    /* check of DDCT */
    putdata(0, n - 1, a);
    fft.ddct(1, a);
    fft.ddct(-1, a);
    a[0] *= 0.5;
    err = errorcheck(0, n - 1, 2.0 / n, a);
    printf("ddct err= %g \n", err);

    /* check of DDST */
    putdata(0, n - 1, a);
    fft.ddst(1, a);
    fft.ddst(-1, a);
    a[0] *= 0.5;
    err = errorcheck(0, n - 1, 2.0 / n, a);
    printf("ddst err= %g \n", err);

    /* check of DFCT */
    putdata(0, n, a);
    a[0] *= 0.5;
    a[n] *= 0.5;
    fft.dfct(a, t);
    a[0] *= 0.5;
    a[n] *= 0.5;
    fft.dfct(a, t);
    err = errorcheck(0, n, 2.0 / n, a);
    printf("dfct err= %g \n", err);

    /* check of DFST */
    putdata(1, n - 1, a);
    fft.dfst(a, t);
    fft.dfst(a, t);
    err = errorcheck(1, n - 1, 2.0 / n, a);
    printf("dfst err= %g \n", err);

    return 0;
}


void putdata(int nini, int nend, double *a)
{
    int j, seed = 0;

    for (j = nini; j <= nend; j++) {
        a[j] = RND(&seed);
    }
}


double errorcheck(int nini, int nend, double scale, double *a)
{
    int j, seed = 0;
    double err = 0, e;

    for (j = nini; j <= nend; j++) {
        e = RND(&seed) - a[j] * scale;
        err = MAX(err, fabs(e));
    }
    return err;
}
testxg.java
import java.util.Scanner;

public class testxg{
	static int seed;
	public static void main(String args[]){
		// 入力
		System.out.println("data length n=? (must be 2^m)");
		Scanner scan = new Scanner(System.in);
		int n = Integer.parseInt(scan.next());
		// 処理
		fft4g fft = new fft4g(n);
		double[] a = new double[n + 1], t = new double[n / 2 + 1];
		double err;

		/* check of CDFT */
		putdata(0, n - 1, a);
		fft.cdft(1, a);
		fft.cdft(-1, a);
		err = errorcheck(0, n - 1, 2.0 / n, a);
		System.out.println("cdft err= " + err);

		/* check of RDFT */
		putdata(0, n - 1, a);
		fft.rdft(1, a);
		fft.rdft(-1, a);
		err = errorcheck(0, n - 1, 2.0 / n, a);
		System.out.println("rdft err= " + err);

		/* check of DDCT */
		putdata(0, n - 1, a);
		fft.ddct(1, a);
		fft.ddct(-1, a);
		a[0] *= 0.5;
		err = errorcheck(0, n - 1, 2.0 / n, a);
		System.out.println("ddct err= " + err);

		/* check of DDST */
		putdata(0, n - 1, a);
		fft.ddst(1, a);
		fft.ddst(-1, a);
		a[0] *= 0.5;
		err = errorcheck(0, n - 1, 2.0 / n, a);
		System.out.println("ddst err= " + err);

		/* check of DFCT */
		putdata(0, n, a);
		a[0] *= 0.5;
		a[n] *= 0.5;
		fft.dfct(a, t);
		a[0] *= 0.5;
		a[n] *= 0.5;
		fft.dfct(a, t);
		err = errorcheck(0, n, 2.0 / n, a);
		System.out.println("dfct err= " + err);

		/* check of DFST */
		putdata(1, n - 1, a);
		fft.dfst(a, t);
		fft.dfst(a, t);
		err = errorcheck(1, n - 1, 2.0 / n, a);
		System.out.println("dfst err= " + err);
	}
	private static void putdata(int nini, int nend, double[] a){
		seed = 0;
		for(int j = nini; j <= nend; j++){
			a[j] = rnd();
		}
	}
	private static double errorcheck(int nini, int nend, double scale, double[] a){
		seed = 0;
		double err = 0;
		for(int j = nini; j <= nend; j++) {
			double e = rnd() - a[j] * scale;
			err = Math.max(err, Math.abs(e));
		}
		return err;
	}
	private static double rnd(){
		seed = (seed * 7141 + 54773) % 259200;
		return 1.0 * seed / 259200;
	}
}

#テスト結果
　n=1024とした場合、どちらも良好な精度を得ました。移植は成功のようです。

// C++
cdft err= 5.55112e-016
rdft err= 6.66134e-016
ddct err= 5.55112e-016
ddst err= 6.66134e-016
dfct err= 4.44089e-016
dfst err= 6.66134e-016
// Java
cdft err= 4.440892098500626E-16
rdft err= 4.440892098500626E-16
ddct err= 6.661338147750939E-16
ddst err= 4.718447854656915E-16
dfct err= 5.551115123125783E-16
dfst err= 6.661338147750939E-16

#大浦版FFTのC++・Java移植版
　GitHubに移行しました。

12
17
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
12
17

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?