C言語でDFT(離散フーリエ変換)のプログラムを作ってみる

2010/09/25
(この記事の文字数: 481)

C言語で DFT のプログラムを書いてみました。

友人の助けがあったおかげですぐに作ることができました。できたDFTのプログラムで適当に作った波形をDFTで周波数空間に変換した後、逆DFTで元の波形に戻したデータをグラフ化してみました。ちなみに元の波形は時間軸をtとすると、次のようなものを使いました。


sin(1*t) + 3*sin(2*t) + 5*sin(5*t) + 3*cos(t);

以下がそれぞれの過程での実部と虚部のグラフです。

最後のグラフは正規化してないですが、これをデータサイズで割って正規化すればほぼ同じ波形に戻るでしょう。

アルゴリズムも改めて勉強しなおしてみればシンプルなものでした。
ただ、基本形の波形であるsinとcosとの相関をとって、相関の強さ＝スペクトル強度となるように周波数空間に置き換えれば良いだけでした。

2次元のDFTも作ったのですが、まだあまり数を試していないので、載せないでおきます。
DFTはできたのものの、やはり遅いので、実用性を考えるとFFTのプログラムがほしいです。次はFFTのプログラムを作ろうと思います。

ちなみに以下が私の作ったDFTのソースコードです。Visual Studio 2008では動作確認しています。


#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define FAILURE (-1)
#define SUCCESS (0)



int dft(double *dest_r, double *dest_i, 
    const double *src_r, const double *src_i, 
    int sampleSize, bool inverse)
{
    double sum_re, sum_im, pi, arg, freqSample;
    int sample;
    int freq;

    if(dest_r == NULL || dest_i == NULL ||
       src_r == NULL || src_i == NULL || sampleSize <= 1){
        return FAILURE;
    }

    // initializing a destination array
    for( sample=0; sample&ltsampleSize; sample++)
    {
        dest_r[sample] = 0.0;
        dest_i[sample] = 0.0;
    }
    
    pi = (inverse == false)? M_PI: -M_PI;

    // main loop
    for( freq=0; freq<sampleSize; freq++)
    {
        sum_re = 0.0;
        sum_im = 0.0;
        arg = ((double)freq/(double)sampleSize) * (2*pi);
        for( sample=0; sample<sampleSize; sample++)
        {
            freqSample = sample * arg;

            sum_re += src_r[sample] * cos( freqSample )
                 + src_i[sample] * sin( freqSample );
            sum_im += src_i[sample] * cos( freqSample )
                 - src_r[sample] * sin( freqSample );
        }
     if(inverse){
          dest_r[freq] = sum_re;
          dest_i[freq] = sum_im;
     }else{
          dest_r[freq] = sum_re/(double)sampleSize;
          dest_i[freq] = sum_im/(double)sampleSize;
     }
    }

    return SUCCESS;
}

int saveArrayAsText1(void *array_src, int arraysize,
          char *const filename)
{
    int i;
    FILE *fp = NULL;

    if(array_src == NULL || arraysize==0)
    { return FAILURE; }

    if((fp=fopen(filename, "w")) == NULL)
    { return FAILURE; }

    for(i=0; i<arraysize; i++)
    {
        fprintf(fp, "%f\n", ((double*)array_src)[i]);
    }
    fclose(fp);

    return SUCCESS;
}

int main()
{
    double *im_r, *im_i, *dest_r, *dest_i;
    int status, datasize;
    int i;

    datasize = 256;

    im_r = (double*)malloc(datasize*sizeof(double));
    im_i = (double*)malloc(datasize*sizeof(double));
    dest_r = (double*)malloc(datasize*sizeof(double));
    dest_i = (double*)malloc(datasize*sizeof(double));

    //initialize
    for(i=0; i&ltdatasize; i++)
    {
        double point = (double)i;
        im_r[i] = sin(point)+3*sin(2*point)
                       +5*sin(5*point)+3*cos(point);
        im_i[i] = 0;
    }

    status = saveArrayAsText1(im_r, datasize,
          "src_r.txt");
    if(status == FAILURE)
    { fprintf(stderr,"dft failed.\n");
   return FAILURE;}

    //dft
    status = dft(dest_r, dest_i, im_r, im_i,
        datasize, false);
    if(status == FAILURE)
    { return FAILURE;}


    status = saveArrayAsText1(dest_r, datasize, 
           "dft_r.txt");
    if(status == FAILURE)
    { fprintf(stderr,"saving file failed.\n");
   return FAILURE;}

    status = saveArrayAsText1(dest_i, datasize,
           "dft_i.txt");
    if(status == FAILURE)
    { fprintf(stderr,"saving file failed.\n");
   return FAILURE;}

    //inverse dft
    status = dft(im_r, im_i, dest_r, dest_i,
        datasize, true);
    if(status == FAILURE)
    { fprintf(stderr,"inverse dft failed.\n");
   return FAILURE;}

    status = saveArrayAsText1(im_r, datasize, 
           "idft_r.txt");
    if(status == FAILURE)
    { fprintf(stderr,"saving file failed.\n");
   return FAILURE;}

    free(im_r);
    free(im_i);
    free(dest_r);
    free(dest_i);

    printf("done!\n");

    return SUCCESS;
}