[Bluestein's FFT] 2^n의 길이가 아닌 데이터의 Fast Fourier Transform

Algorithm

[Bluestein's FFT] 2^n의 길이가 아닌 데이터의 Fast Fourier Transform

빠릿베짱이 2015. 6. 12. 08:31

출처 : http://en.wikipedia.org/wiki/Chirp_Z-transform

FFT를 수행하고 싶은데, 2^n승이 아닐 경우 어떻게 처리해야 하는가?

----> FFT data size is not power of two.

1. Convolution Theorem : [link]

- $\mathcal{F}\{f*g\} = \mathcal{F}\{f\} \cdot \mathcal{F}\{g\}$

- $f*g= \mathcal{F}^{-1}\big\{\mathcal{F}\{f\}\cdot\mathcal{F}\{g\}\big\}$

- f와 g의 컨볼루션의 결과의 푸리에 변환과 f와 g의 각각 푸리에 변환의 곱과 같다.

- 위의 두 가지 성질을 이용해서 2^n 승이 아닌 데이터 길이에 대해 FFT를 계산할 수 있다.

chirpZscheme

<출처 : http://www.katjaas.nl/chirpZ/chirpZ.html>

Bluestein's algorithm[4] expresses the CZT as a convolution and implements it efficiently using FFT/IFFT. As the DFT is a special case of the CZT, this allows the efficient calculation of discrete Fourier transform (DFT) of arbitrary sizes, including prime sizes. (The other algorithm for FFTs of prime sizes, Rader's algorithm, also works by rewriting the DFT as a convolution.) It was conceived in 1968 by Leo Bluestein.[5] Bluestein's algorithm can be used to compute more general transforms than the DFT, based on the (unilateral) z-transform (Rabiner et al., 1969).

블루 스틴 알고리즘은 컨볼루션으로 CZT를 표현하며, FFT와 IFFT를 사용하여 효율적으로 구현할 수 있다. 이산 푸리에 변환은 CZT의 특수한 경우이다, CZT는 임의의 크기의 이산 푸리에 변환을 매우 효율적으로 계산이 가능하다. (소수를 포함하여) - 소수의 FFT 를 위한 다른 알고리즘으로 Rader's algorithm 또한 컨볼루션으로 DFT를 수행한다. 이것은 Leo Bluestein에 의해 1968년에 고안되었다. 블루스테인 알고리즘은 이산 푸리에 변환보다 z 변환에 기반하여 좀 더 일반화된 변환으로 사용할 수 있다.

Recall that the DFT is defined by the formula - DTF는 다음과 같이 정의된다.

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk } \qquad k = 0,\dots,N-1.

If we replace the product nk in the exponent by the identity - exponential 지수의 nk를 다른 표현으로 변경하면.

n k = \frac{-(k-n)^2}{2} + \frac{n^2}{2} + \frac{k^2}{2}

we thus obtain: - 이를 DFT 수식에 적용하여 정리하면,

X_k = e^{-\frac{\pi i}{N} k^2 } \sum_{n=0}^{N-1} \left( x_n e^{-\frac{\pi i}{N} n^2 } \right) e^{\frac{\pi i}{N} (k-n)^2 } \qquad k = 0,\dots,N-1.

This summation is precisely a convolution of the two sequences a_n and b_n defined by:

- 위의 수식에서 summation은 a_n 과 b_n 의 두 시계열의 컨볼루션이다.

a_n = x_n e^{-\frac{\pi i}{N} n^2 }

b_n = e^{\frac{\pi i}{N} n^2 },

with the output of the convolution multiplied by N phase factors b_k^*. That is:

- 결국, 컨볼루션의 출력은 N 위상 b_k^* 와 곱해진다. 따라서,

X_k = b_k^* \sum_{n=0}^{N-1} a_n b_{k-n} \qquad k = 0,\dots,N-1.

This convolution, in turn, can be performed with a pair of FFTs (plus the pre-computed FFT of complex chirp bn) via the convolution theorem. The key point is that these FFTs are not of the same length N: such a convolution can be computed exactly from FFTs only by zero-padding it to a length greater than or equal to 2N–1. In particular, one can pad to a power of two or some other highly composite size, for which the FFT can be efficiently performed by e.g. the Cooley–Tukey algorithm in O(N log N) time. Thus, Bluestein's algorithm provides an O(N log N) way to compute prime-size DFTs, albeit several times slower than the Cooley–Tukey algorithm for composite sizes.

- 이 컨볼루션은 컨볼루션 이론을 통해 한 쌍의 FFT로 수행할 수 있다. 여기서 중요한 점은 이러한 FFT는 같은 길이 N이 아니다. 그 컨볼루션은 2N-1과 같거나 더 큰 길이의 제로 패딩을 수행한 데이터에 의해 FFT를 계산할 수 있다. 특히, 효율적인 FFT(Cooley-Tukey 알고리즘)를 수행하기 위해, 2의 승수 또는 어떤 HCN(Highly composite number)의 크기로 패딩할 수 있다. 따라서 블루스테인 알고리즘은 Composite 크기를 위한 Cooley-Tukey 알고리즘보다 몇 배 느리긴 하지만, Prime-size의 DFT를 계산하는데 O(NlongN)의 방법을 제공한다.

The use of zero-padding for the convolution in Bluestein's algorithm deserves some additional comment. Suppose we zero-pad to a length M ≥ 2N–1. This means that an is extended to an array An of length M, where An = an for 0 ≤ n < N and An = 0 otherwise—the usual meaning of "zero-padding". However, because of the bk–n term in the convolution, both positive and negative values of n are required for bn (noting that b–n = bn). The periodic boundaries implied by the DFT of the zero-padded array mean that –n is equivalent to M–n. Thus, bn is extended to an array Bn of length M, where B0 = b0, Bn = BM–n = bn for 0 < n < N, and Bn = 0 otherwise. A and B are then FFTed, multiplied pointwise, and inverse FFTed to obtain the convolution of a and b, according to the usual convolution theorem.

- 블루스테인 알고리즘에서 컨볼루션을 위한 제로 패딩의 사용에 대해서는 몇가지 추가적으로 설명할 가치가 있다. M>2N-1 에 대해 제로 패딩을 가정해보자. a_n은 길이 M의 배열 A_n으로 확장된다. 여기서 0 < n < N 에 대해서는 A_n= a_n 이며, 반대의 경우에는 A_n은 0이다.( N<= n < M 인 경우)

그러나, 컨볼루션에서 b_(k-n)의 항으로 인해, b_n에 대해 양수와 음수 모두 필요하게 된다. 제로패팅된 배열의 DFT에 의한 주기적인 경계는 -n이 M-n과 같다는 것을 의미한다. 따라서 b_n은 M 길이의 B_n으로 확장된다. 여기서 B_0= b_0이며, 0< b_n < N에 대해서는 B_n = B_m-n = b_n이다, 또한 나머지 B_n은 0이다. 따라서 A와 B를 FFT 수행한 후, 각각 곱하고, 다시 역 FFT를 수행하면, 컨볼루션 이론에 의해 a와 b의 컨볼루션 결과를 얻을 수 있다.

[Code] : http://www.nayuki.io/page/free-small-fft-in-multiple-languages

Matlab code version 1 :

%데이터 열기

fp = fopen('c:\0.ssm');

fread(fp, 44100, 'int16', 0);

fread(fp, 44100+13093, 'int16', 0);

signal = fread(fp, 22050, 'int16', 0);

fclose(fp);

%정규화

signal = signal/32767;

figure(1);

plot(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = complex(signal, 0);

n = [0:22049];

temp = pi * n.*n / 22050 ;

costable = cos(temp)';

sintable = sin(temp)';

areal = zeros(65536,1);

aimag = zeros(65536,1);

areal(1:22050,:) = real(data) .* costable + imag(data) .* sintable;

aimag(1:22050,:) =-real(data) .* sintable + imag(data) .* costable;

breal = zeros(65536,1);

bimag = zeros(65536,1);

breal(1,1) = costable(1,1);

bimag(1,1) = sintable(1,1);

for i=2:22050

breal(i,1) = breal(65536-i+2) = costable(i);

bimag(i,1) = bimag(65536-i+2) = sintable(i);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

afft = fft(complex(areal, aimag));

bfft = fft(complex(breal, bimag));

conv = afft.*bfft;

conv = ifft(conv);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

outreal = imag(conv);

outimag = real(conv);

outreal = outreal(1:22050) .* costable + outimag(1:22050) .* sintable;

outimag = -outreal(1:22050) .* sintable + outimag(1:22050) .* costable;

mfft2 = 10*log10(outreal.*outreal*2 + outimag.*outimag*2 );

figure(2)

plot(mfft2(1700:2300), 'r');

bluestein_ver2.m

%데이터 열기

fp = fopen('c:\0.ssm');

fread(fp, 44100, 'int16', 0);

fread(fp, 44100+13093, 'int16', 0);

signal = fread(fp, 22050, 'int16', 0);

size(signal)

fclose(fp);

%정규화

signal = signal/32767;

figure(1);

plot(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = complex(signal, 0);

n = [0:22049];

temp = pi * n.*n / 22050 ;

W2= exp(j*temp)';

a = zeros(65536,1);

a(1:22050) = data.*W2;

b = zeros(65536,1);

b(1,1) = 1/W2(1,1);

for i=2:22050

b(i,1) = b(65536-i+2) = 1/W2(i);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

afft = fft(a);

bfft = fft(b);

conv = afft.*bfft;

conv = ifft(conv);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

outreal = imag(conv);

outimag = real(conv);

outreal = outreal(1:22050) .* costable + outimag(1:22050) .* sintable;

outimag = -outreal(1:22050) .* sintable + outimag(1:22050) .* costable;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mfft = 10*log10(outreal.*outreal*2 + outimag.*outimag*2 );

figure(2)

plot(mfft(1700:2300), 'r');

Matlab code version 2 :

Result

Signal Frequency

저작자표시 비영리 변경금지

'Algorithm' 카테고리의 다른 글

MSER ( Maximally stable extremal regions ) (0)	2015.08.04
마이크를 이용해서 소리의 위치 검출 (0)	2015.06.12
Computer Vision Resource (0)	2015.05.15
두 분포 간 거리 (distance between distribution) (0)	2013.09.10
Best Fit Plane ( 평면 방정식 찾기 ) (0)	2013.07.12

현재글[Bluestein's FFT] 2^n의 길이가 아닌 데이터의 Fast Fourier Transform

Live, like today is the last day to live.

추적, tracking, Face Detection, Kinect, 수유 맛집, Scala, 라즈베리파이, Deep Learning, 제스처 인식, RANSAC, 셀인, 제주도맛집, AdaBoost, 셀프 인테리어, Spark, Face, OpenNI, Face Alignment, 협제맛집, ISO26262,

Today :
Yesterday :

Jooo's life story