FFTW : C/C++ 고속 푸리에 변환 라이브러리

라이브러리 소개

현재 진행중인 연구에서 사용중인 C언어 푸리에 변환 라이브러리인 FFTW에 대해 간략히 포스팅해볼까 합니다. FFTW는 Fastest Fourier Transform in the West 의 약자인데요. 저 같은 토종 아시안은 이해 못하는 개그에서 파생된 이름이라는 카더라가 있습니다만, 중요한 내용은 아니므로 패스하겠어요.

 

FFTW Home Page

 

위에 링크된 FFTW 웹사이트에서 다운로드 가능하며, 리눅스나 유닉스에서는 소스 파일을 빌드하여 libfftw3.a 등의 라이브러리 파일을 만든 다음, 이를 링크하는 방식으로 이용할 수 있습니다. macOS 사용자의 경우 Homebrew를 통해서도 설치를 할 수 있습니다. 터미널 콘솔을 연 다음 brew install fftw 를 입력하면 되고, 설치가 제대로 완료되었다면 brew list 를 입력했을때 목록에서 fftw를 확인할 수 있을 것입니다.

 

반응형

 

윈도우에서 구동할 경우 DLL 라이브러리를 웹사이트에서 다운로드 받아서 사용하거나, vcpkg를 사용해서 설치할 수 있습니다. vcpkg를 통해 라이브러리를 설치하고 비주얼 스튜디오에서 사용하는 방법은 다음 포스팅에 소개되어 있습니다.

 

vcpkg로 비주얼 스튜디오 라이브러리 설치하기

 

사용법 (버전 3.3.9)

FFTW 라이브러리를 사용하는 데 있어서 필요한 변수형으로는 fftw_complex 및 fftw_plan 이 있는데요. fftw_complex 는 복소수 변수형으로, 변환 대상이 되는 함수와 변환된 함수가 fftw_complex 변수형으로 이루어진 배열에 저장되는 방식을 취하고 있습니다. fftw_plan 은 고속 푸리에 변환을 수행하기 위한 정보를 담고 있는 변수형이며, FFTW를 초기화하기 위한 함수 및 푸리에 변환을 실행하는 fftw_execute 함수와 관련이 있습니다. 라이브러리를 사용하기 위해서는 먼저 헤더파일을 C언어나 C++ 소스파일에 포함시켜야 합니다.

 

#include<fftw3.h>

 

FFTW 라이브러리를 사용한 고속 푸리에 변환은 전체적으로 다음과 같은 구조를 가지고 있습니다.

 

// 이산화된 격자의 갯수
int nbin = 100;

// 푸리에 변환 대상이 되는 함수를 저장하는 배열
fftw_complex *array_fft_ini;
// 변환된 함수를 저장하는 배열
fftw_complex *array_fft_fin;

// 동적 메모리 할당
array_fft_ini =
    (fftw_complex *)fftw_malloc(nbin *
        sizeof(fftw_complex));
array_fft_fin =
    (fftw_complex *)fftw_malloc(nbin *
        sizeof(fftw_complex));

/* 푸리에 변환을 위한 fftw_plan형 변수를
 * 선언하고 초기화 */
fftw_plan fft_plan;
fft_plan =
    fftw_plan_dft_1d(nbin,
        array_fft_ini, array_fft_fin,
        FFTW_BACKWARD, FFTW_MEASURE);

// 변환 대상이 되는 함수를 지정
for (int iq = 0; iq < nbin; iq++) {
    array_fft_ini[iq][0] = [실수부분];
    array_fft_ini[iq][1] = [허수부분];
}

// 고속 푸리에 변환 실행
fftw_execute(fft_plan);

// 변환된 함수
for (int ix = 0; ix < nbin; ix++) {
    [실수부분] = array_fft_fin[ix][0];
    [허수부분] = array_fft_fin[ix][0];
}

// 동적으로 할당된 메모리 해제
fftw_destroy_plan(fft_plan);
fftw_free(array_fft_ini);
fftw_free(array_fft_fin);

 

가장 먼저 동적 메모리 할당을 위한 fftw_malloc 함수를 사용하여, 푸리에 변환의 대상이 되는 함수와 변환된 함수를 저장하기 위한 배열 (array_fft_ini 및 array_fft_fin) 을 구성할 필요가 있습니다.

그 다음으로 fftw_plan 변수형을 가진 변수를 선언하고 이를 초기화해야 하는데요. 일반적인 1차원 이산 푸리에 변환의 경우 fftw_plan_dft_1d 함수를 호출하면 되는데, 이 함수는 격자의 개수 nbin, 배열의 포인터 및 정방향인지 역방향인지 여부 (FFTW_FORWARD / FFTW_BACKWARD) 등을 매개변수로 받고 있습니다. 2차원이나 3차원 푸리에 변환의 경우 각각 fftw_plan_dft_2d 및 fftw_plan_dft_3d 함수들을 호출하여 초기화할 수 있습니다.

그 다음 단계는 변환 대상이 되는 함수의 실수 및 허수부분을 array_fft_ini 배열에 저장하는 것입니다. 그러면 푸리에 변환을 수행하는 fftw_execute 함수를 실행했을 때, 변환된 함수들의 값이 array_fft_fin 배열에 저장되는데요. 필요에 따라 이 배열에 저장된 실수와 허수 값들을 사용하면 되겠습니다.

 

마지막으로 FFTW 라이브러리의 사용이 끝나고 나면, fftw_destroy_plan 및 fftw_free 함수들을 사용해서 동적으로 할당받은 메모리를 해제해줘야 합니다.

 

FFTW 라이브러리를 사용하면서 가장 주의해야 할 점을 꼽으라면, 아마도 변환 대상이 되는 함수들이 특정한 정의역에서 존재한다는 것, 또는 함수들이 0이 아닌 유의미한 값을 갖는 특정한 범위를 가정하고 있다는 게 아닐까 싶습니다. 예를 들어서 FFTW를 이용한 푸리에 변환으로 연결되는 두 함수 X와 Y가 있다고 가정할 때, X와 Y의 정의역 모두 최소값은 0이며 최대값들 사이에는 특정한 관계가 성립됩니다. 좀 더 자세히 말하자면,

두 최대값들의 곱 = 2 * 원주율 * 격자 간격의 갯수

의 등식이 성립하죠. 따라서 올바른 계산을 하기 위해서는, 현재 다루고자 하는 문제에 맞게 보정을 해 줄 필요가 있습니다.

 

 

푸리에 변환으로 연결되어 있는 두 함수 phi와 psi가 있고 이들은 각각 변수 q와 x의 함수들로 존재한다고 할 때, q 또는 x의 최소값이 0이 아닌 경우, 변환대상이 되는 함수 phi를 FFTW에 그대로 집어넣으면 안되고, 위에 언급된 대로 별도의 복소수를 곱해줄 필요가 있습니다. 마찬가지로 변환된 함수 psi역시 FFTW에서 얻은 함수에 또다른 복소수를 곱해줘야 올바른 값을 얻을 수 있지요. q와 x의 격자간격 역시 서로 독립적이지 않은 관계로 둘 중 하나는 임의대로 정할 수 있지만, 그러면 나머지 하나는 자동으로 결정됩니다.

 

여기서 q와 x는 물리적으로 각진동수와 시간 혹은 파동수와 공간에 대응될 수 있는데요. 이러한 개념들이 낮설게 느껴지거나, 파동을 수학적으로 다루기 위한 도구인 푸리에 변환에 대해 더 자세한 사항이 궁금하시다면 다음 포스팅을 읽어보는 것을 추천드립니다.

 

물리학 상식 : 파동의 이해

 

예시 : 가우스 함수

푸리에 변환을 사용하기에 가장 만만한 함수로는 가우스 함수가 있죠. (응?) 그래서 FFTW를 이용해서 푸리에 변환을 해보았습니다. 참고로 변환대상이 되는 함수(phi)가 실수함수일지라도 가우스 함수의 평균이 0이 아니라면, 변환된 함수(psi)는 가우스 함수에 삼각함수를 곱한 형태의 복소수 함수가 되는 것이 특징입니다. 수치해석으로 구한 결과를 해석적인 계산과 비교해 볼 텐데, 이는 새로운 도구나 방법론의 유효성을 검증하는 절차 중의 하나입니다. 제가 잘 알고 있는 문제에 대한 답을 재현할 수 있어야, 미지의 문제에 대해 컴퓨터가 내놓은 답을 신뢰할 수 있기 때문이죠.

 

 

 

FFTW를 이용해 얻은 함수가 해석적인 계산으로 얻은 것과 잘 일치한다는 걸 볼 수 있습니다. 이 포스팅에서는 1차원 푸리에 변환에 대해 다루었지만, 다차원 푸리에 변환 역시도 가능하기 때문에 상당히 유용한 라이브러리라 할 수 있겠습니다. 마지막으로 위에 나온 예시에 사용된 C++ 코드를 접은글로 첨부합니다. FFTW를 설치하고 소스코드를 컴파일 한 뒤에, -lfftw3 플래그를 이용해서 라이브러리를 링크하여 실행파일을 얻을 수 있습니다.

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<fftw3.h>

/* comment out the following line,
 * if you do not want to write a file */
#define WRITE_FILE

static const double PI = 3.141592653589793;

double func_gaussian(double x, double mean_x, double sigma_x);

int main(int argc, char *argv[]) {
    // mean and width of a Gaussian function in the q-space
    double mean_q = 1.;
    double sigma_q = 0.5;

    // lattice spacing in the q-space
    double delta_q = 0.02;

    double q_min = mean_q - 10. * sigma_q;
    double q_max = mean_q + 10. * sigma_q;
    int nbin_q = (int)floor((q_max - q_min) / delta_q + 1e-8);
    if (nbin_q % 2 == 1) {
        nbin_q += 1;
    }
    delta_q = (q_max - q_min) / (double)nbin_q;
    double *qbin;
    double **func_q;
    qbin = new double [nbin_q + 1];
    func_q = new double *[nbin_q + 1];
    for (int iq = 0; iq <= nbin_q; iq++) {
        qbin[iq] = q_min + delta_q * (double)iq;
        func_q[iq] = new double [2];
        // We are doing Fourier transform of a real Gaussian function.
        // real part
        func_q[iq][0] = func_gaussian(qbin[iq], mean_q, sigma_q);
        // imaginary part
        func_q[iq][1] = 0.;
    }

    int nbin_x = nbin_q;
    /* lattice spacing in the x-space,
     * which is depending on one in the q-space */
    double delta_x = 2. * PI / (delta_q * (double)nbin_q);
    double x_min = -0.5 * delta_x * (double)nbin_x;
    double x_max = 0.5 * delta_x * (double)nbin_x;
    double *xbin;
    double **func_x;
    xbin = new double [nbin_x + 1];
    func_x = new double *[nbin_x + 1];
    for (int ix = 0; ix <= nbin_x; ix++) {
        xbin[ix] = x_min + delta_x * (double)ix;
        func_x[ix] = new double [2];
    }

    // the input array containing a function to be tranformed
    fftw_complex *array_fft_ini;
    // the output array containing the outcome
    fftw_complex *array_fft_fin;
    array_fft_ini =
        (fftw_complex *)fftw_malloc(nbin_q * sizeof(fftw_complex));
    array_fft_fin =
        (fftw_complex *)fftw_malloc(nbin_x * sizeof(fftw_complex));

    // declare and initialize FFT transformer
    fftw_plan fft_plan;
    fft_plan =
        fftw_plan_dft_1d(nbin_q, array_fft_ini, array_fft_fin,
                         FFTW_BACKWARD, FFTW_MEASURE);

    // define the input array
    for (int iq = 0; iq < nbin_q; iq++) {
        double factor_in[2];
        factor_in[0] = cos(x_min * delta_q * (double)iq);
        factor_in[1] = sin(x_min * delta_q * (double)iq);
        array_fft_ini[iq][0] =
            func_q[iq][0] * factor_in[0] -
            func_q[iq][1] * factor_in[1];  // real part
        array_fft_ini[iq][1] =
            func_q[iq][1] * factor_in[0] +
            func_q[iq][0] * factor_in[1];  // imaginary part
    }

    /* perform FFTW to take the input array
     * and compute the output array */
    fftw_execute(fft_plan);

    // collect the result from FFTW arrays
    for (int ix = 0; ix < nbin_x; ix++) {
        double factor_out[2];
        factor_out[0] = cos(q_min * xbin[ix]) * delta_q / (2. * PI);
        factor_out[1] = sin(q_min * xbin[ix]) * delta_q / (2. * PI);
        func_x[ix][0] =
            array_fft_fin[ix][0] * factor_out[0] -
            array_fft_fin[ix][1] * factor_out[1];  // real part
        func_x[ix][1] =
            array_fft_fin[ix][1] * factor_out[0] +
            array_fft_fin[ix][0] * factor_out[1];  // imaginary part
    }
    func_x[nbin_x][0] = func_x[0][0];
    func_x[nbin_x][1] = func_x[0][1];

    // write an output file
    #ifdef WRITE_FILE
    FILE *fout;
    fout = fopen("func_gaussian_fft.dat", "w");
    fprintf(fout,
        "#  q  phi(q).re  phi(q).im  x  psi(x).re  psi(x).im\n");
    for (int iq = 0; iq <= nbin_q; iq++) {
        int ix = iq;
        fprintf(fout, "%e  %e  %e  %e  %e  %e\n",
            qbin[iq], func_q[iq][0], func_q[iq][1],
            xbin[ix], func_x[ix][0], func_x[ix][1]);
    }
    fclose(fout);
    #endif

    // deallocate the FFTW stuffs
    fftw_destroy_plan(fft_plan);
    fftw_free(array_fft_ini);
    fftw_free(array_fft_fin);

    for (int iq = 0; iq <= nbin_q; iq++) {
        delete [] func_q[iq];
    }
    delete [] func_q;
    delete [] qbin;

    for (int ix = 0; ix <= nbin_x; ix++) {
        delete [] func_x[ix];
    }
    delete [] func_x;
    delete [] xbin;
}

double func_gaussian(double x, double mean_x, double sigma_x) {
    double dx_scaled = (x - mean_x) / sigma_x;
    double func =
        exp(-0.5 * fabs(dx_scaled * dx_scaled)) / (sigma_x * sqrt(2. * PI));
    return func;
}

 

푸리에 변환은 과학 및 공학 분야에서 광범위하게 사용되는데요. 소리나 전파 등의 시그널에 대한 정보를 이미지화 한 다음 패턴 분석에 사용할 수 있습니다. 이를 스펙트로그램이라 부르는데, 자세한 사항이 궁금하시다면 다음 포스팅을 참고하면 좋습니다.

 

푸리에 변환과 C++ 스펙트로그램

 

앞에서 예시로 언급한 가우스 함수는 정규 분포 (normal distribution)라는 확률 분포를 따르는 연속적인 변수의 확률밀도함수이기도 합니다. 가우스 함수의 정의와 특징에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

수학 상식 : 정규 분포와 가우스 함수

 

---

 

푸리에 변환을 위한 새로운 클래스를 정의하고, 변환대상이 되는 함수들을 매개변수로 넘겨주는 방식으로 코드를 재구성한다면 더 편리하게 사용할 수 있습니다.

 

C/C++ 에서 함수를 매개변수로 사용하기

 

메인 함수의 명령행 인자에 대해서는 다음 포스팅을 참고하면 좋습니다.

 

Command-Line Arguments (C/C++ 공통)

 

---

 

위에 나온 그래프를 그리기 위해서 GNUPLOT (그누플롯)이라는 프로그램을 사용했는데요. 다양한 기능을 가지고 있으면서도 무료로 사용가능하기 때문에, 기본적인 사용 방법을 알아두면 유용합니다.

 

GNUPLOT : 다용도 그래프 유틸리티