PythonからCで書かれたFFTW3を使いたいアホの子

誰がアホやねん！
まあ、それはおいといて……
今回はCの高速FFTライブラリ、FFTW3のラッパーを自分で作ってみたいと思ったりしました。
もうすでにあるじゃん！ってツッコミは無効化です。だから、アホの子だって言ったんです。
じゃ、メモ始めます。
※Windowsは今回、はごってます。我が校にいじめはありません。
※今回、コメントは英語になってます。日本語コメントはめんどいです。
※コメント英語の用法ミスは、ご自由に脳内変換してください。

動作環境：

Ubuntu10.04(64bit)、FreeBSD 8.1(64bit)
FFTW3がCから使えてれば良い

使用イメージ：

まず、どんな風に使うラッパーかと言うと、

#!/usr/bin/env python
#coding: utf-8
import fftw3py
import math

pi=math.atan(1)*4

def make_test_data(num):
    sig = []
    for i in range(num):
        sig.insert(len(sig),0.25*math.sin(pi*i/16))
    return sig

signal=make_test_data(64)

fftw3py.init(len(signal))
fftw3py.inject(signal,"with-hanning")
freq=fftw3py.exe()
fftw3py.destory()

print signal
print freq

こんなイメージ。大雑把にフローを書くと、

make_test_data()で入力信号作る。(別に外部からでもいいけど割愛)
fftw3py(ラッパーモジュール名).init()で初期化。
fftw3py.inject()でCにおけるデータ型に変換。(with-hanningオプションつけるとハニング窓を適用)
fftw3py.exe()で実行。fftw3py.destory()でメモリクリア。(実際の計算はCで行っているので、メモリクリアは必須。)

ちなみに、destoryはdestroyの最上級……ではなく、ただのスペルミスです。
細かいことは気にせずに、つまり、こういう使い方ができるラッパーを書きます。

ラッパー本体：

さて、ラッパー本体の実装には2つのファイルを作ります。今回、モジュール名をfftw3pyとしたので、fftw3py.cとfftw3py.hを作ることになります。

Pythonではimportすると、モジュール名のpyファイルをパスから探します。つまり、ファイル名はモジュール名と同一である必要があることに注意してください。なお、素人の言うことなので、違った場合は嘲笑でもしといてください。

ではでは、それらのファイルの中身はというと、
fftw3py.h

extern PyObject* FFTW3PY_init(PyObject* self, PyObject* args);
extern PyObject* FFTW3PY_inject(PyObject* self, PyObject* args);
extern PyObject* FFTW3PY_exe(PyObject* self, PyObject* args);
extern PyObject* FFTW3PY_destory(PyObject* self, PyObject* args);

extern void initfftw3py();

fftw3py.c

#include <math.h>
#include "Python.h"
#include "fftw3.h"
#include "fftw3py.h"

/***** global variables *****/
static double *in;
static fftw_complex *out;
static fftw_plan plan;
static int sample_points;

/***** Initialization function *****/
PyObject* FFTW3PY_init(PyObject* self, PyObject* args){
    /***** translate python variable to c variable *****/
    if(!PyArg_ParseTuple(args, "i", &sample_points))
        return NULL; // return NULL mean Error
    in = (double *)malloc(sizeof(double)*sample_points);
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*sample_points);
    plan = fftw_plan_dft_r2c_1d(sample_points, in, out, FFTW_ESTIMATE);

    // If you don't have any return value, you must return null character.
    return Py_BuildValue("");
}

/***** set FFT variables *****/
PyObject* FFTW3PY_inject(PyObject* self, PyObject* args){
    int i;
    double h;
    PyObject* plist;
    const char* string = NULL;

    // Parsing args.
    if(!PyArg_ParseTuple(args,"Os",&plist,&string))
        return NULL;

    // set input signal's array from python's list object.
    for(i=0;i<sample_points;i++){
        if(!strcmp(string,"with-hanning"))
            h = (1 - cos(M_PI*i/sample_points))/2;
        else
            h=1;
        in[i] = h*PyFloat_AsDouble(PyList_GetItem(plist,i));
    }

    return Py_BuildValue("");
}

PyObject* FFTW3PY_exe(PyObject* self, PyObject* args){
    int i;
    PyObject* result = PyList_New(sample_points);

    fftw_execute(plan);
    for(i=0;i<(sample_points/2);i++){
        PyList_SetItem(result,i,PyFloat_FromDouble(out[i][0])); //real
        PyList_SetItem(result,i+sample_points/2,PyFloat_FromDouble(out[i][1])); //image
    }

    return result;
}

PyObject* FFTW3PY_destory(PyObject* self, PyObject* args){
    free(in);
    fftw_free(out);
    fftw_destroy_plan(plan);

    return Py_BuildValue("");
}

/***** Define relations for each methods *****/
static PyMethodDef methods[] = {
    {"init", FFTW3PY_init, METH_VARARGS},
    {"inject", FFTW3PY_inject, METH_VARARGS},
    {"exe",FFTW3PY_exe, METH_VARARGS},
    {"destory", FFTW3PY_destory, METH_VARARGS},
    {NULL},
};

/***** Constructor *****/
void initfftw3py(){
    Py_InitModule("fftw3py", methods);
}

です。細かい説明は抜きます。たぶん、コードから読み取っていただく方が素人の説明よりよほど分かりやすいでしょうから。
ただ、要するにPythonとCで異なるデータの扱い方に関して、双方の橋渡し的な部分が今回の肝ということです。データの扱いさえその言語のものに変換してしまえば、あとはその言語でプログラミングするのとほとんど変わりません。
続いて、コンパイルです。

コンパイル、実行：

飽きてきたので、Makefileのせます。
Makefile

.SUFFIXES: .c .so

DIST_FILES= *~ \#*\#*
TARGETS= *.so
FLAGS=`pkg-config --cflags fftw3`
LIBS=`pkg-config --libs fftw3`
PYHEAD=`../SetPyHeadEnv`

.c.so:
        @gcc -shared -fpic ${FLAGS} ${LIBS} -I${PYHEAD} $< -o `echo $@|sed s/.so//`module.so

fftw3py.so:

test:
        @python test.py

clean:
        rm -rf ${DIST_FILES}
distclean: clean
        rm -rf ${TARGETS}

ここで、突然出てきたSetPyHeadEnvとは、FreeBSDとUbuntuの環境差を吸収させるためのクッションみたいなもんです。いや、自分も忘れてたんですが……
SetPyHeadEnv(shell script)

#!/bin/sh
PYVERSION=2.6

case `uname` in
    FreeBSD)
        echo "/usr/local/include/python$PYVERSION"
        #sed "s/PYHEAD=[a-z1-9/.]*/PYHEAD=\/usr\/local\/include\/python$PYVERSION/" Makefile > Makefile.x
        ;;
    Linux)
        echo "/usr/include/python$PYVERSION"
        #sed "s/PYHEAD=[a-z1-9/.]*/PYHEAD=\/usr\/include\/python$PYVERSION/" Makefile > Makefile.x
        ;;
    *)
        ;;
esac

以上。

で、makeするとfftw3pymodule.soが生成されます。これが、モジュールになります。あと、色々コンパイラに怒られますが、無視してます。ツンデレですから。

実行は、たとえばこんなイメージで使いたいんだってところで書いたPythonコードをtest.pyとかで保存したら、

% python test.py

でオーケーです。どんな結果が表示されるかって言うと、

[0.0, 0.048772580504032062, 0.095670858091272445, 0.13889255825490054, 0.20786740307563631,..., -0.095670858091272654, -0.04877258050403202]
[5.7211887261098457e-18, -1.8421947292037881e-16,..., 1.4136105444788528e-16, 0.0, 1.9294112024003265e-17]

みたいな。最初のリストが入力信号、次のリストがFFTW3が返してくれた結果となります。

こんな具合で十分かな？これを後で見て、自分が分かるだろうか……
まあ、分かれ、自分！