실습 B. DFT와 FFT

DFT는 입력을 N개의 복소 지수와 하나씩 내적하는 계산이다. 이 내적을 정의대로 모두 하면 곱셈이 $N^2$번 필요하고, FFT는 같은 값을 $\tfrac N2\log_2 N$번 정도로 얻는다. 이 실습에서는 두 방법을 직접 코드로 쓰고, 값이 같다는 것과 시간이 다르다는 것을 확인한다.

이 실습의 자리: 시뮬레이션 — 실습 A에서 신호를 배열로 만들고 스펙트럼을 그려 보았고, 이 실습은 그 스펙트럼을 내는 DFT를 직접 계산하고 빠르게 계산하는 방법을 쓰며, 다음 실습 C에서 표본화와 창이 스펙트럼에 하는 일을 본다.

준비

활동은 셋이다. 활동 1과 3은 함수를 직접 쓰고, 활동 2는 측정 결과를 보고 설명한다. 함수를 쓰는 활동에는 함수의 이름, 입력, 출력, 계산 순서가 있다. 그 순서대로 쓰고, 바로 아래 확인 셀로 맞는지 본다.

이론과 기호는 2장 DFT와 2장 FFT의 계산량을 따른다.

X[k]=\sum_{n=0}^{N-1}x[n]\,e^{-j2\pi kn/N},\qquad k=0,\ldots,N-1

그림과 표 옆의 기록에 넣기 단추를 누르면 그 결과가 기록에 들어가고, 글 상자에 쓴 설명도 함께 저장된다. 다 마치면 도구 막대의 저장을 누른다. 이 실습은 Pluto를 쓰지 않는다.

활동 1. DFT를 행렬 곱으로 만들기

$X[k]$ 하나는 입력 $x$ 와 복소 지수 $e^{-j2\pi kn/N}$ ( $n=0,\ldots,N-1$ )의 내적이다. 이 지수 N개를 $k$ 번째 행에 늘어놓은 N×N 행렬을 $W$ 라 하면, N개의 내적을 한 번에 쓴 것이 행렬 곱 $\mathbf X=W\mathbf x$ 이다. 즉 $W[k,n]=e^{-j2\pi kn/N}$ 이다.

쓸 함수 dft_matrix(N)

입력: 정수 N
출력: N×N 복소 배열 W, 원소는 $W[k,n]=e^{-j2\pi kn/N}$

계산 순서

번호 배열 n = np.arange(N)을 만든다. 행 번호 $k$ 와 열 번호 $n$ 은 같은 0…N−1이다.
모든 $(k, n)$ 짝의 곱 $kn$ 을 담은 N×N 배열을 만든다. np.outer(n, n)이 그 배열이다([k, n] 자리에 k*n).
그 배열을 지수에 넣는다: np.exp(-2j * np.pi * kn / N). 반복문은 필요 없다.

def dft_matrix(N):
    """N×N DFT 행렬, W[k, n] = e^{-j2πkn/N}"""
    # TODO: 위의 계산 순서대로 dft_matrix 함수를 구현하세요.
    raise NotImplementedError('dft_matrix 함수를 구현하세요.')

아래 셀은 $N=8$ 로 확인한다. 첫째, 행 1을 복소 평면에 찍는다. DFT는 음의 지수를 쓰므로 행 1의 점은 시계 방향으로 한 바퀴 돌아야 한다. 둘째, 2장의 입력 $x=[1,2,0,-1,0,0,0,0]$ 으로 W @ x와 np.fft.fft(x)를 비교한다. 두 값의 차가 $10^{-12}$ 보다 작으면 맞게 쓴 것이다.

W = dft_matrix(8)
show.iq(W[1], title='W의 행 1: n = 0…7')
x = np.array([1, 2, 0, -1, 0, 0, 0, 0])
X = W @ x
max_difference = np.max(np.abs(X - np.fft.fft(x)))
print('W @ x =', np.round(X, 3))
print(
    f'최대 차 |W@x − fft(x)| = {max_difference:.1e}',
    '→ 맞음' if max_difference < 1e-12 else '→ 다시 본다',
)

차가 0이 아니라 아주 작은 값으로 남는 것은 부동소수점 반올림 때문이다. FFT는 DFT와 같은 합을 계산 순서만 바꿔 구한다. 그러니 값이 같아야 하는 것은 당연하고, 남은 물음은 '얼마나 빠른가'이다.

활동 2. N이 두 배가 되면 시간은 몇 배가 되는가

행렬 곱 W @ x는 출력 N칸 각각에 입력 N개를 곱해 더하므로 곱셈이 $N^2$ 번이다. FFT는 같은 결과를 대략 $\tfrac N2\log_2 N$ 번의 곱셈으로 얻는다(활동 3에서 그 이유를 코드로 본다). 그러면 N을 두 배로 할 때마다 행렬 곱의 시간은 약 4배, FFT의 시간은 2배를 조금 넘게 늘어나야 한다.

아래 셀은 N = 64…4096에서 두 방법의 시간을 잰다. 측정은 dsp.seconds_per_call이 한다. 이 함수는 합이 0.05초를 넘을 때까지 같은 계산을 되풀이한 뒤 한 번의 평균 시간을 돌려준다. 행렬을 만드는 시간은 재지 않는다. 몇 초 걸리는 동안 페이지가 멈춘 것처럼 보인다.

sizes = [64, 128, 256, 512, 1024, 2048, 4096]
generator = np.random.default_rng(7)
dft_times, fft_times = [], []
for N in sizes:
    dft_weights = dft_matrix(N)
    signal = generator.standard_normal(N) + 1j * generator.standard_normal(N)
    dft_times.append(dsp.seconds_per_call(lambda signal: dft_weights @ signal, signal))
    fft_times.append(dsp.seconds_per_call(np.fft.fft, signal))
    del dft_weights                                   # N = 4096의 W는 256 MiB라 바로 지운다
show.line(dft_times, fft_times, x=sizes, labels=['W @ x', 'np.fft.fft'],
          points=True, logx=True, logy=True,
          title='계산 시간', xlabel='N', ylabel='한 번의 시간', yunit='s')
show.table([{'N': N, 'W@x (µs)': round(dft_time * 1e6, 1),
             'fft (µs)': round(fft_time * 1e6, 1), '시간의 비': round(dft_time / fft_time)}
            for N, dft_time, fft_time in zip(sizes, dft_times, fft_times)], title='이 브라우저에서 잰 시간')

로그–로그 그림에서는 $N^2$ 이 기울기 2인 직선, $N\log N$ 은 기울기가 1보다 조금 큰 직선으로 보인다. 큰 N 쪽에서 두 선의 기울기를 비교한다. 작은 N에서는 함수를 부르고 배열을 만드는 고정 시간이 계산 시간보다 커서 두 선이 거의 평평하다.

큰 N에서 두 선의 기울기는 각각 얼마쯤인가? 작은 N에서는 왜 그 기울기가 나오지 않는가?

활동 3. FFT를 직접 쓰기

FFT는 한 가지 관찰에서 나온다. DFT의 합을 짝수 번째 표본과 홀수 번째 표본으로 나누면, 각각은 길이 N/2인 DFT가 된다. 짝수 표본 $x[0], x[2], \ldots$ 의 N/2점 DFT를 $E[k]$ , 홀수 표본 $x[1], x[3], \ldots$ 의 N/2점 DFT를 $O[k]$ 라 하자. 홀수 표본은 한 칸 늦게 시작하므로 그 몫에는 회전 인자 $w_N^k=e^{-j2\pi k/N}$ 이 곱해진다. 그래서 $k=0,\ldots,N/2-1$ 에서

X[k]=E[k]+w_N^k\,O[k],\qquad X[k+N/2]=E[k]-w_N^k\,O[k]

이다. 아래 식은 $w_N^{k+N/2}=-w_N^k$ 이기 때문에 위 식과 부호만 다르다. 길이 N의 DFT 하나가 길이 N/2의 DFT 둘과 곱셈 N/2번으로 바뀐 것이다. $E$ 와 $O$ 도 같은 방법으로 다시 나눌 수 있고, 길이가 1이 될 때까지 나누면 길이 1인 수열의 DFT는 그 표본 자신이다. 이렇게 자기 자신을 부르는 함수가 radix-2 FFT이다.

쓸 함수 fft_rec(x)

입력: 복소 배열 x, 길이 N은 2의 거듭제곱
출력: 길이 N의 복소 배열 X, np.fft.fft(x)와 같은 값

계산 순서

N = 1이면 x를 복소수 배열로 바꿔 그대로 돌려준다(x.astype(complex)).
E = fft_rec(x[0::2]), O = fft_rec(x[1::2])로 짝수·홀수 표본의 DFT를 구한다.
회전 인자 tw = np.exp(-2j * np.pi * np.arange(N // 2) / N)을 만든다.
앞 절반 E + tw * O와 뒤 절반 E - tw * O를 np.concatenate로 이어 돌려준다.

def fft_rec(x):
    """radix-2 FFT. len(x)는 2의 거듭제곱"""
    # TODO: 위의 계산 순서대로 fft_rec 함수를 구현하세요.
    raise NotImplementedError('fft_rec 함수를 구현하세요.')

아래 셀은 길이 1024의 합성 복소 수열로 np.fft.fft와 값을 비교하고, 세 구현의 시간을 잰다.

test_signal = generator.standard_normal(1024) + 1j * generator.standard_normal(1024)
max_difference = np.max(np.abs(fft_rec(test_signal) - np.fft.fft(test_signal)))
print(
    f'N = 1024 최대 차 |fft_rec − np.fft.fft| = {max_difference:.1e}',
    '→ 맞음' if max_difference < 1e-9 else '→ 다시 본다',
)
show.table([{'N': N, 'W@x (ms)': round(dft_times[sizes.index(N)] * 1e3, 2),
             'fft_rec (ms)': round(
                 dsp.seconds_per_call(fft_rec, generator.standard_normal(N) + 0j) * 1e3, 2),
             'np.fft.fft (ms)': round(fft_times[sizes.index(N)] * 1e3, 3)}
            for N in (256, 1024, 4096)],
           title='세 구현의 시간')

fft_rec는 계산량이 $N\log N$ 인데도 작은 N에서는 행렬 곱보다 느릴 수 있다. Python에서 함수를 한 번 부르는 비용이 곱셈 한 번보다 훨씬 크고, fft_rec는 길이 N에서 자기 자신을 약 2N번 부르기 때문이다. N이 커질수록 계산량의 차가 그 비용을 이긴다. np.fft.fft는 같은 알고리즘을 컴파일된 코드로 돌려 그 비용이 없다.

N = 4096에서 fft_rec와 W @ x 중 무엇이 빨랐는가? 계산량과 함수 호출 비용으로 설명한다.