スポンサーサイト 

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

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 ▲
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。