고속 푸리에 변환

이 단원에서는 고속 푸리에 변환(FFTs)을 위한 함수들을 기술합니다. 라이브러리에서는 radix-2 방법과 mixed-radix 방법을 제공합니다. 효율성을 위해서 각 기능들은 복소수용 함수와 실수용 함수로 나뉘어져 있습니다. Mixed-radix 함수들은 FFTPACK 라이브러리를 재구현한 형태입니다. FFTPACK의 포트란 코드는 Netlib에서 확인할 수 있습니다(FFTPACK에는 sine, cosine 변환 기능도 제공하지만 GSL 에서 더이상 제공하지 않습니다.). 본 문서에서 기술되는 알고리즘들에 대한 더 자세한 정보와 파생사항들은 “GSL FFT Algorithms” 문서를 참고할 수 있습니다. (참고 문헌과 추가 자료 참고)

수학적 정의

고속 푸리에 변환은 이산 푸리에 변환 (DFT)를 효율적으로 계산하기 위한 알고리즘들을 의미합니다.

\[x_j = \sum_{k=0}^{n-1} z_k \exp(-2 \pi i j k / n)\]

DFT는 연속 푸리에 변환을 근사할 필요가 있을 때 사용됩니다. 일반적으로 공간과 시간에 대해 이산 간격으로 샘플들이 주어졌을 때 사용합니다. 기초적인 이산 푸리에 변환은 행렬-벡터곱 \(W \vec{z}\) 으로 계산됩니다. 일반적인 행렬-곱의 계산 복잡도는 \(n\) 개의 데이터에 대해 \(O(n^2)\) 의 복잡도를 가집니다. 고속 푸리에 변환은 분할 정복 전략을 사용해 행렬 \(W\) 를 여러개의 작은 하위-행렬들과 그에 상응하는 길이 \(n\) 의 정수들로 분할해 연산을 수행합니다. \(n\) 이 여러개의 정수들의 곱 \(f_1 f_2 \dots fm\) 으로 표현될 수 있다면 DFT는 \(O(n \sum f_i)\) 의 시간 복잡도를 가집니다. radix-2 FFT 방법은 \(O(n\log_2 n)\) 의 복잡도를 가집니다.

모든 FFT 함수들은 3가지 형태의 변환을 제공합니다. 변환, 역변환, 그리고 조정 계수가 없는 역변환입니다. 이들은 같은 수학적 정의에 기반해 있습니다.

고속 푸리에 변환은 \(x = \text{FFT(z)}\) 로 다음과 같이 정의됩니다.

\[x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi ijk /n)\]

역변환은 \(x= \text{IFFT}(z)\) 로 다음과 같이 정의됩니다.

\[z_j = \sum_{k=0}^{n-1} x_k \exp(2\pi ijk /n)\]

계수 \(1/n\) 로 완전한 역변환을 얻을 수 있습니다.

예를 들어서, gsl_fft_complex_forward() 와 gsl_fft_complex_inverse() 를 합성하면 완전히 같은 값을 얻을 수 있습니다(수치적 오류가 같이 나올 수 있습니다).

일반적으로, 변환-역변환 과정에서 지수 함수의 지수 부호를 \(\pm\) 중 자유롭게 선택할 수 있습니다. GSL 에서는 FFTPACK과 같은 규약을 따릅니다. 이 규약은 변환에서 음의 부호를 역변환에서는 양의 부호를 사용합니다. 이러한 부호 규칙의 이점은 역변환으로 본래의 함수를 다시 만들 수 있다는 점입니다. 반면, Numerical Recipes 부호 규약이 반대입니다. 변환에 양의 부호를 역변환에 음의 부호를 사용합니다.

조정 계수 없는 역 변환은 backwards FFT 로 불립니다. 이는 말 그대로 역 변환에서 크기를 조정하는 계수를 생략한 변환입니다.

\[z_j^{backwards} = \sum_{k=0}^{n-1} x_k \exp(2\pi ijk /n)\]

결과물의 구체적인 크기가 중요하지 않을 때 조정계수가 없는 변환은 역변환에서 불필요한 나눗셈을 줄일 수 있습니다.

복소수 FFTs 개요

복소수 FFT 함수들은 부동 소수점 배열을 입출력에 사용합니다. 이러한 배열은 실수-허수부가 번갈아 할당된 형태로 사용됩니다. 예를 들어서, 다음은 길이 \(6\) 의 복소 FFT 함수 입력용 배열입니다.

double x[3*2];
gsl_compleX_packed_array data = x;

이 배열은 길이 \(3\) 의 복소수 배열 z[3] 와 같이 사용할 수 있습니다.

data[0] = Re(z[0])
data[1] = Re(z[0])
data[2] = Re(z[1])
data[3] = Re(z[1])
data[4] = Re(z[2])
data[5] = Re(z[2])

배열의 인덱스는 DFT의 수학 정의에 사용된 인덱스와 동일한 순서를 가집니다. 별도의 인덱스 변환을 할 필요가 없습니다.

strid 인자는 사용자 함수를 사용할 때, z[i] 대신 z[strode*i] 형태로 변환을 할 수 있게 해줍니다. \(1\) 보다 큰 이 값은 행렬의 열들에 대해 FFT를 적용할 수 있게 해줍니다. \(1\) 값을 가지면, 배열 사이에 건너뛰는 값 없이 변환을 할 수 있습니다.

벡터 인자들에 대해 FFT를 적용할 수 있습니다. 예를 들어 gsl_vector_complex * v 사용하면 다음과 같은 형태로 이 단원의 함수들을 사용할 수 있습니다.

gsl_complex_packed_array data = v -> data;
size_t *stride* = v->stride;
size_t *n* = v->size;

물리학등의 응용에서 기억해 두어야 할 점은, DFT의 인덱스들이 물리적 주파수에 그대로 대응되지는 않는다는 점입니다. 만약 DFT의 시간 간격이 \(\Delta\) 라면, 주파수 공간은 \(-1/(2\Delta)\) 에서 \(1/(2\Delta)\) 로 주파수의 부호가 \(\pm\) 을 모두 가집니다. \(+\) 주파수는 배열의 시작 지점에서 중간지점까지, \(-\) 의 주파수는 배열의 끝에서 중간까지 공간에 저장됩니다.

다음은 배열 data 자세한 정보와 대응되는 시간-공간의 \(z\) , 주파수-공간의 \(x\) 를 나타낸 표입니다.

index    z               x = FFT(z)

0        z(t = 0)        x(f = 0)
1        z(t = 1)        x(f = 1/(n Delta))
2        z(t = 2)        x(f = 2/(n Delta))
.        ........        ..................
n/2      z(t = n/2)      x(f = +1/(2 Delta),
                                 -1/(2 Delta))
.        ........        ..................
n-3      z(t = n-3)      x(f = -3/(n Delta))
n-2      z(t = n-2)      x(f = -2/(n Delta))
n-1      z(t = n-1)      x(f = -1/(n Delta))

\(n\) 이 짝수일 때, \(n/2\) 지점은 양수와 음수 부호의 주파수 \((+1/(2\Delta), -1/(2\Delta))\) 값을 모두 가집니다. 이 둘은 같습니다. \(n\) 이 홀수라면, 위의 표와 같은 구조를 가집니다. 하지만 \(n/2\) 는 존재하지 않습니다.

복소수 Radix-2 FFTs

이 단원에서 서술하는 radix-2 알고리즘은 간단하고 작지만, 가장 효율적인 방법은 아닙니다. 이 방법은 Cooley-Tukey 알고리즘을 이용해서 길이가 \(2^m, m \in \mathbb{N}\) 인 복소수 데이터의 FFTs를 계산합니다. 계산 과정에서 해당 배열 외에 다른 저장소가 필요 없습니다.

대응되는 자기 정렬 mixed-radix 방법은 계산 과정에서 별도의 공간을 사용해 더 좋은 효율을 보여줍니다.

모든 함수들은 헤더파일 gsl_fft_complex.h 에 기술되어 있습니다.

int gsl_fft_complex_radix2_forward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_transform(gsl_complex_packed_array data, size_t stride, size_t n, gsl_fft_direction sign)

int gsl_fft_complex_radix2_backward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_inverse(gsl_complex_packed_array data, size_t stride, size_t n)

변환, 역변환, 그리고 조정계수 없는 역변환을 주어진 길이 n 과 stride 이용해 복소수 배열 data 에 적용합니다. 이 함수들은 실시간 선택 알고리즘인 radix-2를 사용하며 결과값은 data 그대로 저장됩니다. 배열의 길이 n 은 반드시 \(2^m, m \in \mathbb{N}\) 여야 합니다. transform 이 붙은 함수에 관해, sign 인자는 forward (\(+1\)) 나 backward (\(-1\))가 될 수 있습니다.

오류 없이 변환이 정상적으로 완료되면 각 함수들은 GSL_SUCESS 값을 반환합니다. 만약 n \(2\) 의 배수가 아닐 때에는 GSL_ECOM 값을 반환합니다.

int gsl_fft_complex_radix2_dif_forward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_dif_transform(gsl_complex_packed_array data, size_t stride, size_t n, gsl_fft_direction sign)

int gsl_fft_complex_radix2_dif_backward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_dif_inverse(gsl_complex_packed_array data, size_t stride, size_t n)

Radix-2 변환들의 decimation-in-frequency 형태 입니다.

다음은 짧은 펄스파의 128 길이 샘플 데이터를 고속 푸리에 변환하는 예제 프로그램입니다. 푸리에 변환 결과를 실수로 얻기 위해 펄스는 시간 공간에서 원점을 중심으로 대칭인 범주 (\(-10 \dots 10\))에 대해 정의되어 있습니다. 여기서 음수인 시간은 배열의 끝에서 반올림 됩니다.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
{
  int i; double data[2*128];

  for (i = 0; i < 128; i++)
    {
       REAL(data,i) = 0.0; IMAG(data,i) = 0.0;
    }

  REAL(data,0) = 1.0;

  for (i = 1; i <= 10; i++)
    {
       REAL(data,i) = REAL(data,128-i) = 1.0;
    }

  for (i = 0; i < 128; i++)
    {
      printf ("%d %e %e\n", i,
              REAL(data,i), IMAG(data,i));
    }
  printf ("\n\n");

  gsl_fft_complex_radix2_forward (data, 1, 128);

  for (i = 0; i < 128; i++)
    {
      printf ("%d %e %e\n", i,
              REAL(data,i)/sqrt(128),
              IMAG(data,i)/sqrt(128));
    }

  return 0;
}

프로그램에서 암묵적으로 기본 오류 관리자(오류를 발견하면 프로그램을 정지합니다)를 사용함을 가정함을 알고 있어야합니다. 만약 안전한 오류 관리자를 사용하지 않는다면, gsl_fft_complex_radix2_forward() 의 반환값을 검사해야합니다.

변환된 값들은 \(1/\sqrt n\) 로 크기 조정이 이루어집니다. 따라서 변환된 결과 값은 입력된 초기 값과 잘 맞아 떨어집니다. 입력 값을 적절히 조정해 허수부를 \(0\) 으로 만들었기 때문에 실수 값만을 볼 수 있습니다. \(t=128\) 단위로 변환이 이루어져, DFT는 연속 푸리에 변환을 근사해 변형된 sine 함수를 보여줍니다.

\[\int_{-a}^{+a} e^{-2 \pi i k x} dx = {\sin(2\pi k a) \over\pi k}\]

예제 프로그램의 결과는 다음과 같이 그래프로 그릴 수 있습니다 그림 1.

그림 1 펄스 신호와 그 신호의 이산 푸리에 변환, 예제 프로그램의 결과

복소수 mixed-radix FFT

이 단원은 복소수 자료에 대한 mixed-radix FFT 알고리즘에 대해 서술합니다. Mixed-radix 함수들을 이용해 임의의 길이를 가지는 복소수 자료들의 고속 푸리에 변환을 얻을 수 있습니다. Paul Swartzrauber가 작성한 포트란 라이브러리 FFTPACK의 재구현체를 포함하고 있습니다. 이 방법에 관한 이론은 Clive Temperton 이 작성한 리뷰 논문 “Self-sorting Mixed-radix FFTs”에서 설명하고 있습니다. 함수들은 FFTPACK과 인덱스 표기와 기초 알고리즘을 공유합니다.

Mixed-radix 알고리즘은 부분 변환 모듈에 기반해 있습니다. 잘 최적화된, 소규모 자료를 대상으로한 고속 푸리에 변환 모듈을 사용해서 더 큰 규모의 고속 푸리에 변환을 만들어내는 방법입니다. 2, 3, 4, 5, 6 그리고 7 단계에 대한 최적화된 연산 모듈을 포함하고 있습니다. 4와 6 규모의 고속 푸리에 변환 모듈은 \(2*2\) 과 \(2*3\) 으로 2, 3 단계의 변환 모듈을 합성한 것보다 빠른 속도를 보여줍니다.

모듈 내에 구현되지 않은 길이는 일반적인 길이 \(n\) 구현체로 계산해야 합니다. 해당 구현체는 Singleton 방법을 사용해 이산 푸리에 변환 효율적으로 계산합니다. 해당 모듈은 \(O(n^2)\) 의 길이를 가지는 일반적인 자료를 사용하더라도 낮은 단계로 분해 할 수 있습니다. 예를 들어서 143 길이의 변환을 하고자 한다면 이는 \(11*13\) 의 두 단계로 분해할 수 있습니다. 큰 소수를 인수로 가질 때가 가장 피해야하는 경우입니다. 예로, \(2*3*99991\) 등이 있습니다. 이러한 경우 \(O(n^2)\) 복잡도를 가지는 알고리즘의 연산에 대부분의 구동 시간을 소모하게 됩니다. 이러한 상황에 빠졌다면 GSL 베포 파일에 들어있는 “GSL FFT Algorithms” 문서를 참고하십시오.

Mixed-radix의 초기화 함수 gsl_fft_complex_wavetable_alloc 는 주어진 길이 \(n\) 의 크기를 가지게 됩니다. \(f_i\) 는 \(n\) 의 인자들입니다. 사용자가 작성한 프로그램에서 이러한 분해가 잘 이루어지지 않았을 때, 변환의 속도 저하를 경고하고 싶을 수도 있습니다. 라이브러리 내부에 있는 작은 규모의 소인수 분해 모듈에서 분해할 수 없는 길이의 모듈을 자주 사용한다면, “GSL FFT Algorithms” 문서를 참고해서 다른 인자들에 대한 지원을 추가할 수 있습니다.

이 단원에 기술된 함수들은 gsl_fft_complex.h 에 기술되어 있습니다.

gsl_fft_complex_wavetable *gsl_fft_complex_wavetable_alloc(size_t n)

길이 n 의 복소수 고속 푸리에 변환을 위한 삼각 함수 참고표를 생성합니다. 오류가 없다면 새로 할당된 gsl_fft_complex_wavetable 의 포인터를 반환합니다. 오류가 발생하면 Null 포인터를 반환합니다. 길이 n 은 여러 부분 변환들로 분해되고, 파동 참고표에 대응되는 삼각 함수 계수들이 저장됩니다. 삼각 함수 계수들은 정확도를 위해 sin 와 cos 를 사용해 직접 계산됩니다. 빠른 계산을 위해 재귀적 방법을 사용해 참고표를 작성할 수도 있습니다. 하지만, 프로그램에서 동일한 길이의 고속 푸리에 변환을 반복해서 수행한다면, 이 계산은 한번만 수행해도 최종 결과에 영향을 미치지 않습니다.

파동 참고표 구조체는 같은 길이를 사용하는 모든 변환에 반복적으로 사용할 수 있습니다. 이 표는 다른 고속 푸리에 변환 함수에 사용된다고 해서 값들이 변하지 않습니다. 같은 파동 참고표를 같은 길이를 가지는 값들의 변환, 역변환, 그리고 조정 계수가 없는 역변환 에 모두 사용할 수 있습니다.

void gsl_fft_complex_wavetable_free(gsl_fft_complex_wavetable *wavetable)

파동 참고표 wavetable 가 할당된 메모리를 해제합니다. 이 파동 참고표는 해당하는 길이의 고속 푸리에 변환이 더 이상 필요 없을 때 해제될 수 있습니다.

이 함수들의 연산은 gsl_fft_complex_wavetable 구조체에 적용됩니다. 이 구조체는 고속 푸리에 변환을 위한 내부 계수들을 가지고 있습니다. 해당 구조체의 내부 변수들을 직접 설정하는 것이 필수적이지는 않으나 해당 값들을 알고 있으면 몇몇 상황에서 유용할 수 있습니다. 예를 들어서 고속 푸리에 변환 길이에 대한 소인수 분해가 주어져 있으면, 실행 시간의 추정치나 수치적 오차를 계산하는 데 사용할 수 있습니다.

파동 참고표 구조체는 헤더 파일 gsl_fft_complex.h 에 정의되어 있습니다.

type gsl_fft_complex_wavetable

이 구조체는 mixed-radix fft 알고리즘을 위한 인수 분해와 삼각 함수 참고표를 담고 있습니다. 다음의 구성으로 이루어져 있습니다.

size_t n	복소수 지점들의 갯수
size_t nf	길이:code:n 이 분해된 인수들의 갯수
size_t factor[64]	인수들의 배열. 첫번째 nf 원소만이 사용됩니다.
gsl_complex * trig	길이 n 변환에 대해 사전 할당된 삼각 함수 참고표를 가르키는 포인터
gsl_complex * twiddle[64]	This is an array of pointers into trig , giving the twiddle factors for each pass.

type gsl_fft_complex_workspace

mixed-radix 알고리즘은 변환 중간 단계를 유지하기 위한 별도의 작업 공간이 필요합니다.

gsl_fft_complex_workspace *gsl_fft_complex_workspace_alloc(size_t n)

길이 n 의 복소수 값들의 변환을 위한 작업 공간을 반환합니다.

void gsl_fft_complex_workspace_free(gsl_fft_complex_workspace *workspace)

작업 공간 workspace 가 할당된 메모리를 해제합니다. 이 작업 공간은 동일한 길이의 변환이 더 이상 필요 없을 때 해제할 수 있습니다.

다음 함수들은 실제 변환을 계산합니다.

int gsl_fft_complex_forward(gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable *wavetable, gsl_fft_complex_workspace *work)

int gsl_fft_complex_transform(gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable *wavetable, gsl_fft_complex_workspace *work, gsl_fft_direction sign)

int gsl_fft_complex_backward(gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable *wavetable, gsl_fft_complex_workspace *work)

int gsl_fft_complex_inverse(gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable *wavetable, gsl_fft_complex_workspace *work)

이 함수들은 각각 고속 푸리에 변환, 역변환 그리고 조정 계수 없는 역변환을 주어진 길이 n 과 stride 걸음 크기를 가지는 복소수 값들 data 에 대해 계산합니다. mixed-radix decimation-in-frequency 알고리즘을 사용합니다. 길이 n 에 대한 제약은 존재하지 않습니다. 길이 2, 3, 4, 5, 6 그리고 7 에 대해 최적화된 부분 변환 모듈이 있습니다. 나머지 길이들은 \(O(n^2)\) 의 시간 복잡도를 가지는 상대적으로 느린 방법으로 계산됩니다. 호출 과정에서 변환하는 값들의 길이에 맞는 파동 참고표 wavetable 와 작업 공간 work 를 인자로 제공해야합니다. wavetable 에는 참고용 삼각 함수 계수들을 포함하고 있습니다. transform 이 붙은 함수에 관해, sign 인자는 forward ( \(+1\) ) 나 backward ( \(-1\) )가 될 수 있습니다.

오류가 없으면 0 값을 반환합니다. 다음은 이 함수를 위해 정의된 gsl_errno 의 오류 조건 입니다.

GSL_EDOM	자료의 길이 n 이 양의 정수가 아닐 때, (예 n=0 ).
GSL_EINVAL	자료의 길이 n 와 계산을 위해 사용하는 파동 참고표 wavetable 가 맞지 않을 때.

다음 프로그램은 Mixed-radix 알고리즘을 이용해 짧은 펄스에 대해, 길이 630 ( \(1=2*3*3*5*7\) )의 고속 푸리에 변환을 계산하는 예제를 보여줍니다.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
{
  int i;
  const int n = 630;
  double data[2*n];

  gsl_fft_complex_wavetable * wavetable;
  gsl_fft_complex_workspace * workspace;

  for (i = 0; i < n; i++)
    {
      REAL(data,i) = 0.0;
      IMAG(data,i) = 0.0;
    }

  data[0] = 1.0;

  for (i = 1; i <= 10; i++)
    {
      REAL(data,i) = REAL(data,n-i) = 1.0;
    }

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e %e\n", i, REAL(data,i),
                                IMAG(data,i));
    }
  printf ("\n");

  wavetable = gsl_fft_complex_wavetable_alloc (n);
  workspace = gsl_fft_complex_workspace_alloc (n);

  for (i = 0; i < (int) wavetable->nf; i++)
    {
       printf ("# factor %d: %zu\n", i,
               wavetable->factor[i]);
    }

  gsl_fft_complex_forward (data, 1, n,
                           wavetable, workspace);

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e %e\n", i, REAL(data,i),
                                IMAG(data,i));
    }

  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);
  return 0;
}

프로그램에서 암묵적으로 기본 오류 관리자(오류를 발견하면 프로그램을 정지합니다)를 사용함을 가정함을 알고 있어야합니다. 만약, 안전한 오류 관리자를 사용하지 않는다면 gsl 속 모든 함수들의 반환값을 검사해야합니다.

실수 FFTs 개요

실수 값들에 적용하는 고속 푸리에 변환은 복소수인 경우와 비슷하지만, 변환과 역변환에서 중요한 차이가 있습니다. 바로 실수 값들의 푸리에 변환은 실수 값이 아니라는 점입니다. 변환 된 값들은 다음과 같은 대칭성을 가지는 복소 수열로 나타나게 됩니다.

\[z_k = z_{n-k}^*\]

이러한 대칭성을 가지는 수열을 컬례 복소 나 반 복소성을 가지고 있다고 부릅니다. 이러한 구조적 차이점은 변환(실수->반 복소), 역변환(반 복소 -> 실수)에 대해 이전의 복소 푸리에 변환과 다른 설계를 필요로 합니다. 때문에, 실수 푸리에 변환의 함수들은 두 개의 범주로 나뉘어집니다. 실수 배열을 처리하는 gsl_fft_real 함수들과 반 복소 배열을 처리하는 gsl_fft_halfcomplex 함수들입니다.

gsl_fft_real 함수들은 실수 배열의 주파수 계수를 계산합니다. 실수 배열 \(x\) 에 대한, 반 복소 계수 \(c\) 는 푸리에 분석에 따라 다음과 같이 주어집니다.

\[c_k = \sum_{j=0}^{n-1} x_j \exp(-2 \pi i j k /n)\]

gsl_fft_halfcomplex 함수들은 역변환과 조정 계수 없는 역변환을 계산합니다. 이 함수들은 반 복소 주파수 계수 \(c\) 로 부터 실수 배열을 재구성합니다.

\[x_j = {1 \over n} \sum_{k=0}^{n-1} c_k \exp(2 \pi i j k /n)\]

반 복소 배열의 대칭성으로 인해 배열 전체의 절반만 저장해도 됩니다. 나머지 절반은 반 복소 대칭성으로 계산 과정에서 재구성 됩니다. 이 방법은 짝수, 홀수를 가리지 않고 모든 크기에 대해 사용할 수 있습니다. 짝수 크기의 배열에 대해, 중간값 \(k/2\) 도 실수가 됩니다. 따라서 반 복소 배열의 계산에 필요한 것은 n 크기의 실수 저장 공간뿐입니다. 계산 결과인 실수 배열도 입력 배열과 같은 크기의 배열에 저장됩니다.

정확한 저장 공간의 관리는 알고리즘에 따라 다릅니다. radix-2와 mixed-radix에 따라 차이가 있습니다. radix-2 방법은 실수와 허수부가 저장되는 위치를 가능한 한 멀리 떨어지게 저장하도록 강제합니다. 반면, mixed-radix 방법은 이러한 제약이 없으며, 실수와 허수부가 인접 위치에 저장합니다. 이는 메모리 접근에서 인접성 향상에 좋습니다.

실수 Radix-2 FFTs

이 단원은 radix-2 고속 푸리에 변환 알고리즘을 실수 값들에 적용시키는 함수들에 대해 다룹니다. 이들은 Cooley-Tukey 알고리즘을 사용해 \(2\) 의 거듭 제곱 크기의 길이를 가지는 값들에 대해 계산합니다.

이 radix-2 고속 푸리에 변환 함수들은 헤더 파일 gsl_fft_real.h 에 저장되어 있습니다.

int gsl_fft_real_radix2_transform(double data[], size_t stride, size_t n)

radix-2 고속 푸리에 변환을 주어진 stride 간격과 길이 n 를 가지는 실수 배열 data 대해 계산합니다. 계산 결과는 반 복소 배열로 각 실수-허수 값들은 정해진 위치 규약에 따라 저장됩니다. 길이 \(n\) 의 배열에 대해, \(k < n/2\) 대해, \(k\) 번째 복소수의 실수 값은 \(k\) 번째 배열에 저장되고, 대응 되는 허수 값은 \(n-k\) 번째 배열에 저장됩니다. \(k > n/2\) 인 항들은 대칭성 \(z_k = z_{n-k}^*\) 에 의해 재구성 됩니다. \(k=0\) 이나 \(k= n/2\) 인 경우는 항상 실수가 되며, 따로 처리됩니다. 실수 값은 \(0\) 과 \(n/2\) 위치에 저장되며, 허수 값은 \(0\) 이므로 따로 저장되지 않습니다.

다음 표는 어느 실수 배열의 계산 값 data 허수부가 \(0\) 인 동일한 복소수 배열을 이용해 얻은 결과를 비교한 표입니다. ( stride \(= 1\) 로 가정합니다.)

complex[0].real    =    data[0]
complex[0].imag    =    0
complex[1].real    =    data[1]
complex[1].imag    =    data[n-1]
...............         ................
complex[k].real    =    data[k]
complex[k].imag    =    data[n-k]
...............         ................
complex[n/2].real  =    data[n/2]
complex[n/2].imag  =    0
...............         ................
complex[k'].real   =    data[k]       k' = n - k
complex[k'].imag   =   -data[n-k]
...............         ................
complex[n-1].real  =    data[1]
complex[n-1].imag  =   -data[n-1]

계산 결과는 나중에 기술될 함수 gsl_fft_halfcomplex_radix2_unpack() 로 완전한 복소수 배열로 바꿀 수 있습니다.

반 복소 배열을 위한 radix-2 고속 푸리에 변환 함수들은 헤더 파일 gsl_fft_halfcomplex.h 에 기술되어 있습니다.

int gsl_fft_halfcomplex_radix2_inverse(double data[], size_t stride, size_t n)

int gsl_fft_halfcomplex_radix2_backward(double data[], size_t stride, size_t n)

radix-2 고속 푸리에 변환의 역변환과 조정 계수 없는 역변환을 주어진 stride 간격과 길이 n 가지는 반 복소 배열 data 대해 계산합니다. 이때, data gsl_fft_real_radix2() 서 사용된 결과 배열과 동일한 저장 규칙을 가지고 있습니다. 결과로 나오는 실수 배열은 기본 배열 순서를 따릅니다.

int gsl_fft_halfcomplex_radix2_unpack(const double halfcomplex_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)

gsl_fft_real_radix2_transform() 함수의 계산 결과로 나오는 halfcomplex_coefficient 배열을 일반 복소수 complex_coefficient 배열로 변환합니다. 이 함수는 복소수 배열을 대칭성 \(z_k = z_{n-k}^*\) 를 이용해, 중복되는 요소들을 재구성합니다.

알고리즘은 다음과 같습니다.

complex_coefficient[0].real = halfcomplex_coefficient[0];
complex_coefficient[0].imag = 0.0;

for (i = 1; i < n - i; i++)
  {
    double hc_real = halfcomplex_coefficient[i*stride];
    double hc_imag = halfcomplex_coefficient[(n-i)*stride];
    complex_coefficient[i*stride].real = hc_real;
    complex_coefficient[i*stride].imag = hc_imag;
    complex_coefficient[(n - i)*stride].real = hc_real;
    complex_coefficient[(n - i)*stride].imag = -hc_imag;
  }

if (i == n - i)
  {
    complex_coefficient[i*stride].real = halfcomplex_coefficient[(n - 1)*stride];
    complex_coefficient[i*stride].imag = 0.0;
  }

실수 mixed-radix FFTs

이 단원에서는 실수 값들을 위한 mixed-radix 고속 푸리에 알고리즘을 기술합니다. 고속 푸리에 변환을 위한 mixed-radix 함수들은 모든 크기의 값에 대해 적용할 수 있습니다. 이 함수들은 Paul Swarztrauber가 포트란 FFRPACK 라이브러리에 작성한 실수 고속 푸리에 변환 기능들을 재구현한 형태입니다. 이 알고리즘에 사용된 이론들은 Clive Temperton이 작성한 “Fast Mixed-Radix Real Fourier Transforms”을 참고할 수 있습니다. 이 단원에 작성된 함수들은 FFTPACK에 있는 기반 알고리즘과 같은 인덱스 규약을 가집니다.

이 단원의 함수들은 FFTPACK의 반 복소 배열 저장 규약을 사용합니다. 이 규약에서 실수 배열의 반 복소 변환은 \(0\) 에서 증가하는 순서로, 실수-허수부가 번갈아 이웃해 가며 저장됩니다. 실수로 확인된 값에 대해서는 허수부가 저장되지 않습니다. \(0\) 의 주파수를 가지는 부분에는 허수부가 저장되지 않습니다(). 크기가 짝수인 입력 배열에서는, \(n/2\) 의 허수부가 저장되지 않습니다. 이는 \(z_k = z_{n-k}^*\) 로 인해 자명하게 실수이기 때문입니다.

저장 규약은 다음 예시로 보는 것이 가장 정확합니다. 아래의 표는 \(n=5\) 의 홀수 길이 배열에 대한 계산 결과를 나타냅니다. 각 두 열은 대응되는 반-복소열 값 \(5\) 개를 나타내고 있습니다. 이 값들은 각각 gsl_fft_real_transform() 반환값 halfcomplex[] 나태납니다. complex[] 동일한 실수열이 gsl_fft_complex_backward() 복소수열로 전달 되었을 때 반환됩니다. 이 실수열은 허수부가 0 복소수열로 취급됩니다.

complex[0].real  =  halfcomplex[0]
complex[0].imag  =  0
complex[1].real  =  halfcomplex[1]
complex[1].imag  =  halfcomplex[2]
complex[2].real  =  halfcomplex[3]
complex[2].imag  =  halfcomplex[4]
complex[3].real  =  halfcomplex[3]
complex[3].imag  = -halfcomplex[4]
complex[4].real  =  halfcomplex[1]
complex[4].imag  = -halfcomplex[2]

complex 배열의 뒷 부분 , complex[3] complex[4] 는 대칭 조건으로 채워집니다. 대칭성에 의해 주파수 성분이 0 부분의 허수부 complex[0].imag \(0\) 으로 정해집니다.

다음 표는 길이가 짝수인 배열의 결과를 나타내줍니다. \(n=6\) 인 경우입니다. 짝수 길이를 가지는 경우 완전히 실수인 성분이 두 개 생기게 됩니다.

complex[0].real  =  halfcomplex[0]
complex[0].imag  =  0
complex[1].real  =  halfcomplex[1]
complex[1].imag  =  halfcomplex[2]
complex[2].real  =  halfcomplex[3]
complex[2].imag  =  halfcomplex[4]
complex[3].real  =  halfcomplex[5]
complex[3].imag  =  0
complex[4].real  =  halfcomplex[3]
complex[4].imag  = -halfcomplex[4]
complex[5].real  =  halfcomplex[1]
complex[5].imag  = -halfcomplex[2]

complex 배열의 뒷부분, complex[4] complex[5] 대칭 조건에 의해 결정됩니다. complex[0].imag complex[3].imag \(0\) 으로 정해집니다.

이 함수들은 헤더 파일 gsl_fft_real.h gsl_fft_halfcomplex.h 정의되어 있습니다.

type gsl_fft_real_wavetable

type gsl_fft_halfcomplex_wavetable

FFT를 위한 고정된 크기의 파동 자료표를 위한 구조체입니다.

gsl_fft_real_wavetable *gsl_fft_real_wavetable_alloc(size_t n)

gsl_fft_halfcomplex_wavetable *gsl_fft_halfcomplex_wavetable_alloc(size_t n)

math:n 길이의 실수 원소들에 대한 FFT에 쓰이는 고정된 크기의 삼각함수 파동 자료표를 제공합니다. 이 함수들은 새로 할당된 구조체를 가르키는 포인터를 반환하고 오류가 발견되면 Null 인터를 반환합니다. 길이 n 부분 변환들의 곱으로 분해되며, 인자와 삼각함수 계수들은 파동 자료표에 저장됩니다. 삼각함수 계수들은 정확도를 위해 sin cos 호출해 직접 계산됩니다. 속도를 위해 재귀적 관계를 사용해 계산할 수도 있습니다. 하지만, 응용 프로그램에서 많은 횟수의 FFT를 같은 길이의 자료들에 대해 계산한다면 이 파동 자료표는 최종 결과에 영향을 미치지 않습니다.

파동 자료표 구조체는 같은 크기의 자료를 변환하는데 반복적으로 사용될 수 있습니다. 이 표는 다른 FFT 함수들의 호출에서 변하지 않습니다. 순방향의 실수, 반복소 변환에 대해 적절한 유형의 파동 자료표를 사용해야 합니다.

void gsl_fft_real_wavetable_free(gsl_fft_real_wavetable *wavetable)

void gsl_fft_halfcomplex_wavetable_free(gsl_fft_halfcomplex_wavetable *wavetable)

파동 자료표 wavetable 에 할당된 메모리를 해제합니다. 파동 자료표의 해제는 동일한 크기의 FFT가 더 이상 필요 없을 때 해제될 수 있습니다.

Mixed-radix 알고리즘은 변환의 중간 단계를 유지하기 위해 추가적인 작업 공간을 필요로 합니다.

type gsl_fft_real_workspace

실수 FFT를 계산하기 위한 계수들을 가지고 있는 작업 공간입니다.

gsl_fft_real_workspace *gsl_fft_real_workspace_alloc(size_t n)

길이 n 의 실수 변환을 위한 작업 공간을 할당합니다. 같은 작업공간이 순방향 실수나 역방향 반 복소 변환에 대해 사용될 수 있습니다.

void gsl_fft_real_workspace_free(gsl_fft_real_workspace *workspace)

작업 공간 workspace 에 할당된 메모리를 해제합니다. 작업 공간은 동일한 크기의 FFT가 더이상 필요 없을 때 해제될 수 있습니다.

다음 함수들은 실수, 반복소 자료들에 대한 변환을 계산합니다.

int gsl_fft_real_transform(double data[], size_t stride, size_t n, const gsl_fft_real_wavetable *wavetable, gsl_fft_real_workspace *work)

int gsl_fft_halfcomplex_transform(double data[], size_t stride, size_t n, const gsl_fft_halfcomplex_wavetable *wavetable, gsl_fft_real_workspace *work)

주어진 길이 n 의 실수나 반복소 배열 값 data 대한 FFT를 계산합니다. mixed-radix decimation-in-frequency 알고리즘을 사용합니다. gsl_fft_real_transform() 에 대해, data 시간 순으로 배열된 실수 자료로 취급됩니다. gsl_fft_halfcomplex_transform() 에 대해, data 는 위에 서술한 반-복소 순서로 배열된 푸리에 계수를 가지고 있습니다. n 대한 제약은 없습니다. 효율을 위해 모듈들은 길이 \(2,3,4\) 와 \(5\) 에 대한 부분 변환을 제공합니다. 일반적인 \(n\) 길이 모듈에 대해 나머지 값들은 \(O(n^2)\) 의 속도로 느리게 계산됩니다. 호출시 삼각함수 파동 자료표를 포함하는 wavetable 작업 공간 work 와 같이 제공해야합니다.

int gsl_fft_real_unpack(const double real_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)

실수 배열 real_coefficient 을 동일한 복소수 배열 complex_coefficient 로 변환합니다. 이 복소수 배열은 허수부가 모두 \(0\) 이고 실수부가 실수 배열과 같은 배열이며, gsl_fft_complex 함수들에 대해 사용할 수 있습니다.

다음과 같이 간단하게 구현되어 있습니다.

for (i = 0; i < n; i++)
{
  complex_coefficient[i*stride].real = real_coefficient[i*stride];
  complex_coefficient[i*stride].imag = 0.0;
}

int gsl_fft_halfcomplex_unpack(const double halfcomplex_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)

gsl_fft_real_transform() 서 반환된 반복소 계수 배열 halfcompelex_coefficient 와 일반적인 복소수 배열 complex_coefficient 을 변환합니다. 이 함수는 복소수 배열을 \(z_k = z_{n-k}^*\) 대칭성을 사용해 중복된 요소들을 재 구성합니다. 이 변환에 사용된 알고리즘은 다음과 같습니다.

complex_coefficient[0].real = halfcomplex_coefficient[0];
complex_coefficient[0].imag = 0.0;

for (i = 1; i < n - i; i++)
  {
    double hc_real = halfcomplex_coefficient[(2 * i - 1)*stride];
    double hc_imag = halfcomplex_coefficient[(2 * i)*stride];
    complex_coefficient[i*stride].real = hc_real;
    complex_coefficient[i*stride].imag = hc_imag;
    complex_coefficient[(n - i)*stride].real = hc_real;
    complex_coefficient[(n - i)*stride].imag = -hc_imag;
  }

if (i == n - i)
  {
    complex_coefficient[i*stride].real = halfcomplex_coefficient[(n - 1)*stride];
    complex_coefficient[i*stride].imag = 0.0;
  }

gsl_fft_real_transform() 과 gsl_fft_halfcomplex_inverse() 의 사용을 보여주는 예제 프로그램이 있습니다. 이 프로그램은 사각형 모양의 실수 신호를 생성합니다. 신호는 푸리에 변환을 통해 주파수 공간으로 변환돼고, 가장 낮은 \(10\) 개의 주파수는 푸리에 계수 배열에서 제거되어 gsl_fft_real_transform() 에 의해 반환됩니다.

나머지 푸리에 계수들은 다시 역변환되어 시간 공간으로 되돌아옵니다. 이는 필터를 거친 사각 신호를 나타냅니다. 푸리에 계수는 반 복소 대칭성을 사용해 저장되므로 양의 신호와 음의 신호는 모두 제거되고 최종 필터링된 신호도 실수가 됩니다.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

int
main (void)
{
  int i, n = 100;
  double data[n];

  gsl_fft_real_wavetable * real;
  gsl_fft_halfcomplex_wavetable * hc;
  gsl_fft_real_workspace * work;

  for (i = 0; i < n; i++)
    {
      data[i] = 0.0;
    }

  for (i = n / 3; i < 2 * n / 3; i++)
    {
      data[i] = 1.0;
    }

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e\n", i, data[i]);
    }
  printf ("\n");

  work = gsl_fft_real_workspace_alloc (n);
  real = gsl_fft_real_wavetable_alloc (n);

  gsl_fft_real_transform (data, 1, n,
                          real, work);

  gsl_fft_real_wavetable_free (real);

  for (i = 11; i < n; i++)
    {
      data[i] = 0;
    }

  hc = gsl_fft_halfcomplex_wavetable_alloc (n);

  gsl_fft_halfcomplex_inverse (data, 1, n,
                               hc, work);
  gsl_fft_halfcomplex_wavetable_free (hc);

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e\n", i, data[i]);
    }

  gsl_fft_real_workspace_free (work);
  return 0;
}

이 프로그램의 결과는 다음 그림 그림 2 과 같습니다.

그림 2 실수 펄스 신호의 저주파 통과 필터의 적용, 예제 프로그램의 결과.

참고 문헌과 추가 자료

FFT 입문에 다음의 리뷰 논문을 추천합니다.

• P. Duhamel and M. Vetterli. Fast Fourier transforms: A tutorial review and a state of the art. Signal Processing, 19:259-299, 1990.

GSL 구현체에 쓰인 알고리즘을 참고하기 위해서 “GSL FFT Algorithms” 문서를 참고할 수 있습니다 1 . 이 문서는 FFT에 관한 일반적인 정보들과 각각의 구현체와 그 알고리즘에 관한 여러 파생 정보들을 기술하고 있습니다. 이에 더해, 관련 참고 문헌들도 함께 제공하고 있습니다. 편의를 위해 중요한 참고문헌들을 이곳에 함께 제공합니다.

FFT를 예시 프로그램과 함께 제공하는 몇몇 입문 서적들이 있습니다. Brigham 저 “The Fast Fourier Transform”과 Burrus and Parks 저 “DFT/FFT and Convolution Algorithms”이 있습니다.

• 5. Oran Brigham. “The Fast Fourier Transform”. Prentice Hall, 1974.

• 3. 19. Burrus and T. W. Parks. “DFT/FFT and Convolution Algorithms”, Wiley, 1984.

이 입문 서적들은 radix-2 FFT 방법을 자세히 서술하고 있습니다. FFTPACK 구현체의 핵심 알고리즘인 Mixed-radix 알고리즘은 Clive Temperton의 논문에서 잘 다루어져 있습니다.

• Clive Temperton, Self-sorting mixed-radix fast Fourier transforms, Journal of Computational Physics, 52(1):1-23, 1983.

실수 데이터에 관한 FFTs의 파생 정보들은 다음의 두 논문에서 다루고 있습니다.

• Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney Burrus. Real-valued fast Fourier transform algorithms. “IEEE Transactions on Acoustics, Speech, and Signal Processing”, ASSP-35(6):849-863, 1987.

• Clive Temperton. Fast mixed-radix real Fourier transforms. “Journal of Computational Physics”, 52:340-350, 1983.

1979년도에, IEEE에서는 포트란 FFT 프로그램을 중점적으로 리뷰한 개괄 문서을 출판했습니다. 해당 문서는 “Programs for Digital Sigal Processinng”이란 제목으로 여러 다른 FFT 알고리즘들의 구현에 좋은 참고 문헌이 되어줍니다.

• Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal Processing Committee, editors. Programs for Digital Signal Processing. IEEE Press, 1979.

대규모 데이터에 FFT를 적용할 때, 이러한 작업 전용으로 만들어진 FFTW 라이브러리를 쓰는 것을 추천합니다. 이 라이브러리는 Frigo와 Johnson이 작성했습니다. FFTW 라이브러리는 구동 하드웨어에 맞추어 최고의 효율을 내도록 스스로를 최적화 합니다. 이 라이브러리는 GNU GPL하에서 사용가능합니다.

• FFTW Website, http://www.fftw.org/

FFTPACK의 소스 코드는 http://www.netlib.org/fftpack/ 에서 확인할 수 있습니다.

1

이 문서는 GSL 라이브러리에 포함되어 있습니다. doc/fftalgorithms.tex 경로에 tex 문서로 되어있으며, pdf는 이곳 에서 내려 받을 수 있습니다.