complex.hを使ったFFT関数 

用途が激しく限られる、C言語で書かれたFFT(高速フーリエ変換)関数です。Tomy's Homeに掲載されていたものを、C99で追加された標準ライブラリcomplex.hのcomplex型に対応させてみました（※C++の標準ライブラリのものとは違います）。
complex型とは所謂 複素数型のことで、実部と虚部を単一変数で扱えます。これを使えば、コードが大分すっきりするかも。


ただ、未確認ですが、GCCでしか動かないと思います。
用途は…卒論で使うぐらいかな…。
[fft.c]

使い方は、渡す配列がcomplex型という以外ココと同じです。
flagが0のときは通常のFFTで、1のときは逆FFTになります。
fft1のiterは、例えばnが512だったら9というように2の乗数の指数部を指定して下さい。


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

void fft1(complex c[], int n, int iter, int flag)
{
    int i, it, j, j1, j2, k, p, q, xp, xp2;
    double arg, sign, w;
    complex t, ww;

    // 指定サイズチェック
    if (n < 2) {
        fprintf(stderr, "Error : n < 2  in fft1()\n");
        return;
    }
    if (iter <= 0) {
        iter = 0;
        i = n;
        while((i /= 2) != 0) iter++;
    }
    j = 1;
    for (i = 0; i < iter; i++) j *= 2;
    if (n != j) {
        fprintf(stderr, "Error : n != 2 ^ k  in fft1()\n");
        return;
    }

    sign = (flag == 1)? 1.: -1.;
    xp2 = n;
    for (it = 0; it < iter; it++) {
        xp = xp2;
        xp2 = xp / 2;
        w = M_PI / xp2;
        for (k = 0; k < xp2; k++) {
            arg = k * w;
            ww = cos(arg) + I * sign * sin(arg);
            i = k - xp;
            for (j = xp, p = i+xp, q = i+xp+xp2;
                 j <= n; j += xp, p += xp, q += xp) {
                t = c[p] - c[q];
                c[p] = c[p] + c[q];
                c[q] = t * ww;
            }
        }
    }
    j1 = n / 2;
    j2 = n - 1;
    j = 1;
    for (i = 1, p = 0; i <= j2; i++, p++) {
        if (i < j) {
            t = c[p];
            c[p] = c[j-1];
            c[j-1] = t;
        }
        k = j1;
        while (k < j) {
            j -= k;
            k /= 2;
        }
        j += k;
    }
    if (flag == 0) return;
    w = n;
    for (i = 0, p = 0; i < n; i++, p++) c[p] = c[p] / w;
    return;
}


void fft2(complex c[], int m, int n, int flag)
{
    int i, iter, j, k, p;
    complex *w;

    // サイズチェック
    if (m < 2) {
        fprintf(stderr, "Error : Illegal parameter in fft2()\n");
        return;
    }
    iter = 0;
    i = m;
    while ((i /= 2) != 0) iter++;
    j = 1;
    for (i = 1; i <= iter; i++) j *= 2;
    if (m != j) {
        fprintf(stderr, "Error : m != 2 ^ k  in fft2()\n");
        return;
    }

    w = (complex *)malloc(m * sizeof(complex));
    if (!w) {
        fprintf(stderr, "Error : Out of memory  in fft2()\n");
        return;
    }

    for (j = 0; j < n; j++, k += n) {
        for (i = 0, p = j; i < m; i++, p += n) w[i] = c[p];
        fft1(w, n, iter, flag);
        for (i = 0, p = j; i < m; i++, p += n) c[p] = w[i];
    }
    for (i = k = 0; i < m; i++, k += n) {
        for (j = 0, p = k; j < m; j++) w[j] = c[p++];
        fft1(w, n, iter, flag);
        for (j = 0, p = k; j < m; j++) c[p++] = w[j];
    }

    free(w);
}

コメント

C言語どこで勉強したんですか？

Cは主に大学の講義でです。
ただ、このコードは参考元のモノを少し弄っただけなので、自分で書いたとは言えないけれど…

コメントの投稿















管理者にだけ表示を許可する

トラックバック

この記事のトラックバックURL
http://yotama.blog84.fc2.com/tb.php/4-367a530b

FFT

FFTFFT* 高速フーリエ変換 (Fast Fourier Transform)* ファイナルファンタジータクティクス (FINAL FANTASY TACTICS)* フロンティア航空の国際民間航空機関|ICAO航空会社コード。 　FFT 関連いろいろブログ紹介

  • [2007/02/07 05:43]
  • URL |
  • そらの部屋 |
  • TOP ▲