FFT - 고속 푸리에 변환

📝 상세 정리

#

• 이제 슬슬 FFT를 이해할 때가 된것 같다. 정리해보자.

푸리에 변환

#

• 여러 사인파 / 코사인파가 더해진 함수 $h(t)$가 있다.
  • 일단 먼저 사인파만 있다고 가정해보자.
  • 우리는 이 함수 $h(t)$에서 어떤 주파수를 가진 성분들이 있는지 검사하고싶다.
• 더 쉽게, 먼저 사인파 하나에 대해서 생각해보자.
  • 함수 $f(x) = \sin{ax}$가 있다.
    • 우리는 $a$를 모른다고 가정하자.
    • 어떻게 하면 이 $a$를 알아낼 수 있을까?
  • 임의의 사인파 $g(x) = \sin{bx}$를 곱해보자.
    • 삼각함수의 곱셈정리에 의해
    • $f(x)g(x) = \sin{ax}\sin{bx} = \frac{\cos{((a-b)x)} - \cos{((a+b)x)}}{2}$
      • 위와 같이 나타낼 수 있다.
    • 이를 음의 무한대부터 양의 무한대까지 적분하면, 주기함수의 특성을 이용해서 $a, b$가 같은지 다른지 알아낼 수 있다.
      • $a, b$는 주파수기에 $a, b > 0$이 보장된다.
      • $a = b$일 때만 양의 무한대로 발산하고, 그 외의 경우 진동발산한다.
• 위의 특징은 여러 사인파의 합 $h(t)$에 대해서도 간단하게 성립한다.
• 코사인 파가 섞여있다면 코사인파를 곱해서 같은 방법으로 수행해주면 되겠다.
• 이 때, 오일러 공식 $e^{ix} = \cos{x} + i\sin{x}$을 이용해서 실수부 / 허수부로 나눠서 더 쉽게 계산할 수도 있다.
• 이를 식으로 나타내면 다음과 같다.
  • $$\hat{h}(x) = \int_{-\infty}^{\infty}{h(t)e^{-2\pi itx}dt}$$
  • 이는 신호 -> 주파수로의 변환으로 해석하면 된다.
• 반대로 주파수 -> 신호의 변환도 가능하다.
  • $$h(t) = \int_{-\infty}^{\infty}{\hat{h}(x)e^{2\pi itx}dx}$$
  • 이를 푸리에 역변환이라고 하자.
• 여기서 갖고 놀아보자. 이해할때 클로드가 만들어줬다.
  신호 종류: 유한 사인파 (현실적)두 주파수 합성구형파무한 사인파 (이상적)

  ① 시간 도메인 — h(t)

  실제 데이터 (샘플링된 이산 값들)

  ② 주파수 도메인 — |X[k]| (DFT 결과)

  각 주파수 성분의 크기. 스파이크가 있는 곳 = 해당 주파수 존재

• 복소평면에서의 이해도 가능하다. 주파수가 다르다면 상쇄되어 없어진다.

  푸리에 변환 — 복소평면 시각화

  h(t) = sin(2π·f_signal·t)  |  검사 주파수 x 를 바꿔가며 적분값이 어떻게 달라지는지 확인

  신호 주파수 f = 3 Hz 5 Hz 7 Hz

  검사 주파수 x = 1 Hz 2 Hz 3 Hz ✓ 5 Hz 7 Hz

  속도 0.3× 1× 2×

  복소평면 — h(t)·e^{−2πitx} 궤적 합: —

  t = 0.000s Re=0.000 Im=0.000

  시간 도메인 — h(t) 와 검사파

  누적 Re = 0.000 누적 Im = 0.000 |합| = 0.000

이산 푸리에 변환 (DFT)

#

• 현실에서 얻을 수 있는 데이터는 이산적이다.
  • 물론 라플라스 변환마냥 수학적으로 변환해서 계산하는 테크닉적으로 쓰려면 그냥 쓰면 된다.
  • 사인파를 많이 쓸거같은 음악 데이터로 생각해보면, 결국 소리를 받을때도, 저장할때도 이산적으로 저장할 수 밖에 없다.
  • 따라서 이산 데이터에 대해서도 푸리에 변환을 정의해보자.
• 그러면 적분값 대신 시그마 값으로 표현할 수 있겠다.
• $$X[k] = \sum_{n=0}^{N-1}x[n]\exp\left(-i\frac{2\pi kn}{N}\right) = \sum_{n = 0}^{N - 1}{x[n] \omega^{kn}} \,\,\,\,  (\omega = \exp\left(-i\dfrac{2\pi}{N}\right))$$

• 위와 같이 대응되게 생각하면 된다.
  • $\omega$는 단위근이라고 한다.
• 이때 $\omega^N = 1$이므로, $0 \leq k < N$ 인 $k$에 대해서만 유의미한 결과를 얻어낼 수 있다.
  • 따라서 DFT에서는 최대 $N$개의 주파수에 대해서만 검증할 수 있다.
• 여기서도 물론 역변환을 정의할 수 있다. $$x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]\omega^{-kn}$$

고속 푸리에 변환

#

• 위의 푸리에 변환을 계산하는데, 총 $N$개의 주파수에 대해 $N$개 신호의합을 구하니까 $O(N^2)$의 시간복잡도가 걸린다.
• 여기서 관찰을 조금 해보자.
  • 8개의 데이터 $x_0, x_1, x_2, \cdots, x_7$이 있다고 해보자.
  • 이때 $x_n, x_{n+4}$에 곱해지는 계수를 비교하면 $\exp{(-i\pi k)}$ 만큼의 차이가 난다.
  • $k$가 홀수일때는 $-1$, 짝수일때는 $1$이다.
  • 이를 이용해서 DFT의 식을 짝수항과 홀수항으로 나눠서 다음과 같이 쓸 수 있다.  $$\begin{align}  X[k] &= \sum_{n=0}^{N-1} x[n]\exp\left(-i\frac{2\pi kn}{N}\right) \\  &= \sum_{n=0}^{N/2-1} x[2n]\exp\left(-i\frac{2\pi k(2n)}{N}\right) + \sum_{n=0}^{N/2-1} x[2n+1]\exp\left(-i\frac{2\pi k(2n+1)}{N}\right) \\  &= \sum_{n=0}^{N/2-1} x[2n]\exp\left(-i\frac{2\pi kn}{N/2}\right) + \exp\left(-i\frac{2\pi k}{N}\right)\sum_{n=0}^{N/2-1} x[2n+1]\exp\left(-i\frac{2\pi kn}{N/2}\right) \\  &= E[k] + \exp\left(-i\frac{2\pi k}{N}\right)O[k] \\  &= E[k] + \omega^kO[k]  \end{align}$$  - 관찰  - 위 식의 짝수항 / 홀수항 부분은 $N/2$ 크기의 DFT와 같다.  - 이 식은 $k < N/2$ 에 대해서만 정의되지만, $k \geq N/2$에 대해서 처리하는 방법은 위에서 찾았다. $$\begin{align} X[k+N/2] &= \sum_{n=0}^{N-1} x[n]\exp\left(-i\frac{2\pi (k+N/2)n}{N}\right) \\ &= \sum_{n=0}^{N-1} x[n]\exp\left(-i\frac{2\pi kn}{N}\right)\exp(-i\pi n) \\ &= E[k] - \omega^k O[k] \end{align}$$  - 따라서 분할정복을 이용해서 $O(N\log{N})$ 시간에 연산을 수행할 수 있겠다.  - 이를 고속 푸리에 변환, FFT라 하자.

합성곱

#

• 근데 저 수학 씹덕같은걸 왜배웠냐?
  • DFT를 통해 합성곱을, 합성곱으로 다항식의 곱셈을, 다항식의 곱셈으로 두 수의 곱셈을 할 수 있다고 한다.
• 다음과 같은 함수가 있다고 하자.
  • $f(x) = a_{N-1}x^{N-1} + a_{N-2}x^{N-2} + \cdots + a_1x + a_0$
  • 여기서, $x = \omega^k$를 대입하면 그 결과는 DFT와 같다!
  • 이를 이용해서 두 다항식 $f, g$에 대한 DFT로 새로운 함수 $h$에 대한 DFT 결과값을 구하고, 이를 역변환해서 $h$를 찾아낼 수 있겠다.
• 이걸 코드로 옮기면 다음과 같다.

using cd = complex<double>;
const double PI = acos(-1);
vector<cd> fft(vector<cd> a){
    int N = a.size();
    if(N == 1) return a;
    vector<cd> even(N/2), odd(N/2);
    rep(i, 0, N/2){
        even[i] = a[2*i];
        odd[i] = a[2*i+1];
    }
    even = fft(even);
    odd = fft(odd);
    vector<cd> ret(N);
    rep(k, 0, N/2){
        cd t = polar(1.0, -2*PI*k/N) * odd[k];
        ret[k] = even[k] + t;
        ret[k+N/2] = even[k] - t;
    }
    return ret;
}

vector<cd> ifft(vector<cd> a){
    int N = a.size();
    if(N == 1) return a;
    vector<cd> even(N/2), odd(N/2);
    rep(i, 0, N/2){
        even[i] = a[2*i];
        odd[i] = a[2*i+1];
    }
    even = ifft(even);
    odd = ifft(odd);
    vector<cd> ret(N);
    rep(k, 0, N/2){
        cd t = polar(1.0, 2*PI*k/N) * odd[k];
        ret[k] = even[k] + t;
        ret[k+N/2] = even[k] - t;
    }
    return ret;
}

• ifft를 사용할때 마지막에 a의 길이로 나눠줘야 함에 주의하자.
• 사실상 구조가 똑같아서, bool inv같은걸로 최적화하는게 더 좋은거같다

최적화

#

• 하지만 이렇게 재귀함수로 새로 배열을 만들면서 진행하면, 느리고 메모리도 많이 먹는다.
• 인접한 원소들 끼리만 연산할 수 있게, 배열을 재배치하자.
• $[0, 4, 2, 6, 1, 5, 3, 7]$과 같은 모양이면 좋겠다.
• 이는 비트를 좌우로 뒤집는 것으로 수행 가능하다.
• 이렇게 하면 반복문으로 계산할 수 있고, polar 극좌표점 자체도 미리 계산해둬서 쓰면 $NlogN$번에서 $N/2$번으로 줄일 수 있겠다.

namespace FFT{
    void fft(vector<cd> &a, bool inv = false){
        int N = a.size();
        assert(N > 0 && (N & (N-1)) == 0); // N은 2의 거듭제곱이어야 함

        // 비트 뒤집기
        for(int i = 1, j = 0; i < N; i++){
            int bit = N >> 1;
            for(; j & bit; bit >>= 1) j ^= bit;
            j ^= bit;
            if(i < j) swap(a[i], a[j]);
        }

        // 단위근 미리 계산
        double ang = 2 * acos(-1) / N * (inv ? 1 : -1);
        vector<cd> w(N/2); // 단위근
        rep(i, 0, N/2) w[i] = cd(cos(ang*i), sin(ang*i));

        // Cooley-Tukey FFT
        for(int len = 2; len <= N; len <<= 1){
            int step = N / len;
            for(int i = 0; i < N; i += len){
                rep(j, 0, len>>1){
                    cd even = a[i+j], odd = a[i+j+(len>>1)] * w[j*step];
                    a[i+j] = even + odd;
                    a[i+j+(len>>1)] = even - odd;
                }
            }
        }

        // 역 FFT인 경우 결과를 N으로 나눔
        if(inv) rep(i, 0, N) a[i] /= N;
    }
}

• 1067번 문제 기준 608ms -> 92ms로, 6.6배 가까이 줄어들었다!

❔질문 사항

#

🔗 참고 자료

#

• DFT, FFT 이해하기