[Algorithm] 푸리에 변환(Fourier Transform)

.1차원 데이터에 대한 이산 푸리에 변환

.이산 푸리에 변환(Discrete Fourier Transform, DFT)

- 원래 신호처리 분야에서 시간축에 따른 신호의 세기를 분석하기 위해 연구

- 이산 함수에 대한 푸리에 변환

- 푸리에 변환에 의해 생성된 함수는 복소수 공간에서 정의, 오일러 공식(Euler formula)

- 복소지수함수를 삼각함수로 변환할 수 있도록 하는 식

- 실제 푸리의 변환 구현 시 함수 F(u)의 실수부(Re)와 허수부(Im)를 따로 고려하여 계산

- 입력 함수인 F(X)가 복소수 함수일 경우

- 

.이산 푸리에 역변환(Invers Discrete Fourier Transform, IDFT)

- 이산 푸리에 변환에 의해 생성된 함수 F(u)는 다시 역변환 과정을 거쳐서 원래의 함수로 변환 가능

#include "StdAfx.h"
#include "IppFourier.h"
 
#include  <math.h>
 
const double PI = 3.14159265358979323846;
// re : 실수부, im : 허수부, N : 데이터 개수, dir : 푸리에 변환 방향(1 : DFT, -1 : IDFT)
void DFT1d(double* re, double* im, int N, int dir)
{
    double* tr = new double[N];
    double* ti = new double[N];
    // 입력 데이터 복사
    memcpy(tr, re, sizeof(double) * N);
    memcpy(ti, im, sizeof(double) * N);
 
    register int i, x;
    double sum_re, sum_im, temp;
 
    for (i = 0; i < N; i++)
    {
        sum_re = sum_im = 0;
        // 푸리에 변환 공식 적용(입력 함수가 복소수 함수일 경우 공식)
        for (x = 0; x < N; x++)
        {
            temp = 2 * dir * PI * ((double)i * x / N);
 
            sum_re += (tr[x] * cos(temp) + ti[x] * sin(temp));
            sum_im += (ti[x] * cos(temp) - tr[x] * sin(temp));
        }
 
        re[i] = sum_re;
        im[i] = sum_im;
    }
  
    if (dir == -1) // IDFT
    {
        for (i = 0; i < N; i++)
        {
            re[i] /= (double)N;
            im[i] /= (double)N;
        }
    }
 
    delete[] tr;
    delete[] ti;
}

- 영상 신호와 같은 자연계에서 발생하는 대부분의 신호는 모두 실수만으로 구성되어 있기 때문에 보통 함수에 전달되는 im 배열의 값은 모두 0으로 설정

- 또, 푸리에 역변환을 이용하여 입력 신호를 다시 복원할 경우 im 배열의 값은 이용하지 않고, 오직 re 배열에 저장된 값만을 참조

.기저 함수(basis funstion)

- 거듭 제곱 함수 또는 위상이 다른 삼각함수와 같이 조합의 기본이 되는 함수

-

.푸리에 변환(Fourier Transform)

- sin 함수와 cos 함수를 기저 함수로 사용하여 입력 신호를 변환하는 방법

-

.고속 푸리에 변환(Fast Fourier Transform, FFT)

- 이산 푸리에 변환을 획기적으로 빠르게 계산하는 방법

- DFT 계산식에서 반복되는 연산을 최소화함으로써 푸리에 변환 계산 속도를 향상시키는 알고리즘

- 원래 1차원 데이터에 대하여 정의된 알고리즘

- 입력 데이터의 개수가 2의 승수로 이루어질 때에만 적용 가능

- 이산 푸리에 변환식을 짝수 번째 데이터와 홀수 번째 데이터로 분리하여 식을 정리

새로운 함수 형태로 정의

아래 식에서는 입력 데이터의 짝수 번째 데이터와 홀수 번째 데이터만을 이용하여 푸리에 변환을 수행.

위 식을 정리하면 아래 식과 같음

이산 푸리에 변환식 정리

.1차원 데이터에 대한 FFT 구현 함수

- 입력 데이터의 순서를 바꾸는 작업은 전통적으로 비트 반전 방법을 사용 (비트의 앞/뒤 순서를 반전)

#include <algorithm>
// re : 실수부, im : 허수부, N : 데이터 개수, dir : FFT 방향(1 : 순방향, -1 : 역방향)
void FFT1d(double* re, double* im, int N, int dir)
{
    register int i, j, k;
 
    //-------------------------------------------------------------------------
    // 입력 데이터 순서 바꾸기 (비트 반전)
    //-------------------------------------------------------------------------
 
    int n2 = N >> 1;
    int nb = 0;
 
    while (N != (1 << nb))
        nb++;
 
    for (i = 0, j = 0; i < N - 1; i++)
    {
        if (i < j)
        {
            std::swap(re[i], re[j]);
            std::swap(im[i], im[j]);
        }
 
        k = n2;
        while (k <= j)
        {
            j -= k;
            k >>= 1;
        }
 
        j += k;
    }
 
    //-------------------------------------------------------------------------
    // 버터플라이(Butterfly) 알고리즘
    //-------------------------------------------------------------------------
 
    int i1, l, l1, l2;
    double c1, c2, t1, t2, u1, u2, z;
 
    c1 = -1.0;
    c2 = 0.0;
    l2 = 1;
 
    for (l = 0; l < nb; l++)
    {
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0;
        u2 = 0.0;
 
        for (j = 0; j < l1; j++)
        {
            for (i = j; i < N; i += l2)
            {
                i1 = i + l1;
                t1 = u1 * re[i1] - u2 * im[i1];
                t2 = u1 * im[i1] + u2 * re[i1];
                re[i1] = re[i] - t1;
                im[i1] = im[i] - t2;
                re[i] += t1;
                im[i] += t2;
            }
 
            z = u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
 
        c2 = sqrt((1.0 - c1) / 2.0);
 
        if (dir == 1) // Forward
            c2 = -c2;
 
        c1 = sqrt((1.0 + c1) / 2.0);
    }
 
    if (dir == -1) // IDFT
    {
        for (i = 0; i < N; i++)
        {
            re[i] /= static_cast<double>(N);
            im[i] /= static_cast<double>(N);
        }
    }
}

 

--

참고 : Visual C++ 영상 처리 프로그래밍