Signal Processing and It’s Application

April. 2011 written by 이주형 (Lee Ju Hyung)

신호 처리(Signal Processing)는 자연계나 물성계에 존재하는 신호자체의 분석을 통해 신호 속에 내재되어 있는 다양한 정보를 추출할 뿐만 아니라 원하는 형태로 대상 신호를 변형하는 것을 말한다. 신호는 아날로그 또는 디지털로 표현되며, 주로 음성, 이미지, 영상 처리와 같은 분야에서 쓰이고 있지만, 활용하기에 따라 더 많은 분야에 적용시킬 수 있을 것이다.

이 글에서는 중요한 수학적 개념들을 먼저 살펴 보고, 마지막에 데모를 통해 활용 가능성을 짚어보기로 한다.

1. n 차원 벡터 (N-Dimensional Vector)

서로 독립된 n 개의 값들의 나열, ex) 1차원, 2차원, 3차원, … n차원

각 차원의 값은 서로 영향을 주지 않음.

3D 그래픽스에서 주로 사용하는 것은 3차원 벡터

2. 벡터 공간 연산 (Vector Space Operation)

3 차원에서의 성질들을 n 차원으로 확장 !

2.1 Length

벡터 자기 자신의 길이를 간단히 알 수 있음

    \[\text{length}(\mathbf{a}) = \sqrt{\mathbf{a} \cdot \mathbf{a}}\]

2.2. Normal

자기 자신의 길이로 모든 성분을 나누어 단위벡터(normal vector) = (길이가 1 인 벡터) 로 만들 수 있음

2.3 Projection

두개의 벡터는 서로 투영 (projection) 한 벡터를 구할 수 있음

    \[\text{proj}_{\mathbf{a}}(\mathbf{b}) = \frac{(\mathbf{a} \cdot \mathbf{b}) \mathbf{a}}{ |\mathbf{a}|^2 }\]

    \[\text{proj}_{\mathbf{b}}(\mathbf{a}) = \frac{(\mathbf{b} \cdot \mathbf{a}) \mathbf{b}}{ |\mathbf{b}|^2 }\]

2.4 Orthogonal 판단

두 벡터가 직교(orthogonal) 하는지 여부는 투영했을 때 길이가 0 인지 확인하여 알수 있음

간단히 \mathbf{a} \cdot \mathbf{b} = 0 을 확인

2.5 Orthonormal

n 차원의 단위직교(orthonormal) 벡터들은 무수히 많지만 n 개 만큼만 존재

2.6 투영하여 성분값 알아내기

임의의 벡터를 단위직교(orthonormal) 벡터들과 dot product 하면 각각의 단위 성분을 얻을 수 있음

3. Sine, Cosine 함수의 직교성

  • \sin(x + \alpha)[0, 2\pi] 구간안에서 적분값이 0 이다 (\alpha 는 임의의 상수)
  • \cos(x + \alpha)[0, 2\pi] 구간안에서 적분값이 0 이다
  • \sin, \cos 함수는 주기를 정수배하더라도 마찬가지로 적분값이 0 이다.
  • 함수들을 곱해도 즉, \sin(k_1 x + \alpha) \cos(k_2 x + \beta)[0, 2\pi] 구간안에서 적분값이 0 이다 (k_1, k_2 는 임의의 정수)
  • 두 함수의 곱의 적분은 두 함수 벡터의 내적(dot product) 이다

따라서

    \[\begin{aligned} & ...,\; \sin(-3x),\; \sin(-2x),\; 0,\; \sin(x),\; \sin(2x),\; \sin(3x),\; ... \\ & ...,\; \cos(-3x),\; \cos(-2x),\; \cos(-x),\; 1,\; \cos(x),\; \cos(2x),\; \cos(3x),\; ...  \end{aligned}\]

위 함수들은 [0, 2\pi] 구간안에서 서로 orthogonal 하다. (\sin, \cos 함수의 성질을 따라 구간을 [-\pi, \pi] 로 잡아도 상관없다)

4. 푸리에 급수 (Fourier Series)

  • 주기(periodical) 함수의 주기를 정수배해서 더해도 함수의 기본 주기는 불변
  • [0, N] 구간안에서 k 주기를 갖는 \sin(\frac{2 k \pi x}{N} + \alpha_k) 함수들을 합성(combination) 하면 임의의 주기 함수를 만들어 낼 수 있음 (k0 보다 큰 서로 다른 정수)
  • 임의의 signal 은 구간을 정해서 주기함수처럼 다룰 수 있음
  • 따라서 임의의 복잡한 signal 을 정수 frequency 를 갖는 여러 \sin 주기 함수(sinusoid) 로 분해가 가능 (반대로 합성도 가능)

    \[f(x) = a_0 \sin \left( \frac{0\pi x}{N} + \alpha_0 \right) + a_1 \sin \left( \frac{2\pi x}{N} + \alpha_1 \right) + a_2 \sin \left( \frac{4\pi x}{N} + \alpha_2 \right) + a_3 \sin \left( \frac{6\pi x}{N} + \alpha_3 \right) + ...\]

이것을 푸리에 급수라고 하고 일반적인 형태로 sin 과 cos 을 사용하여 아래와 같이 쓴다.

(1)   \[f(x) = a_0 + \sum_{k=1}^\infty \left\{ a_k \cos \left( \frac{2 k \pi x}{N} \right) + b_k \sin \left( \frac{2 k \pi x}{N} \right) \right\} \]

4.1 푸리에 급수의 계수(coefficient) 구하기

합성하는 방법을 알았으니 이제 분해하는 방법을 알아보자

푸리에 급수에서 각각의 주기 함수 성분은 서로 직교(orthogonal) 한다는 성질을 이용해서 각 계수들을 알아낼 수 있다.

먼저, 구간 [0, N] 에서 orthonormal 한 형태의 주기함수들을 구해보자 (3D 벡터를 normalize 하는 것과 계산 원리는 똑같다)

normalized 된 \cos 주기함수:

    \[\begin{aligned} g_k(x) &= \frac{ \cos \left( \frac{ 2 k \pi x }{ N } \right) }{ \sqrt{ \int_0^N \cos^2 \left( \frac{ 2 k \pi x }{ N } \right) dx } } \\ \Downarrow \\ g_0(x) &= \frac{1}{\sqrt{N}} \\ g_1(x) &= \sqrt{\frac{2}{N}} \cos\left(\frac{2\pi x}{N}\right) \\ g_2(x) &= \sqrt{\frac{2}{N}} \cos\left(\frac{4\pi x}{N}\right) \\ g_3(x) &= \sqrt{\frac{2}{N}} \cos\left(\frac{6\pi x}{N}\right) \\ \vdots \\ \end{aligned}\]

normalized 된 \sin 주기함수:

    \[\begin{aligned} h_k(x) &= \frac{ \sin \left( \frac{ 2 k \pi x }{ N } \right) }{ \sqrt{ \int_0^N \sin^2 \left( \frac{ 2k \pi x }{ N } \right) dx } } \\ \Downarrow \\ h_1(x) &= \sqrt{\frac{2}{N}} \sin\left(\frac{2\pi x}{N}\right) \\ h_2(x) &= \sqrt{\frac{2}{N}} \sin\left(\frac{4\pi x}{N}\right) \\ h_3(x) &= \sqrt{\frac{2}{N}} \sin\left(\frac{6\pi x}{N}\right) \\ \vdots \\ \end{aligned}\]

이제 임의의 signal 함수 f(x) 를 방금 구한 orthonormal 함수와 dot product 하여 각 계수들을 구할 수 있다

    \[\begin{aligned} c_0 &= \int_0^N f(x) \frac{1}{\sqrt{N}} dx \\ c_1 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \cos\left(\frac{2\pi x}{N}\right) dx \\ c_2 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \cos\left(\frac{4\pi x}{N}\right) dx \\ c_3 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \cos\left(\frac{6\pi x}{N}\right) dx \\ \vdots \\ s_1 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \sin\left(\frac{2\pi x}{N}\right) dx \\ s_2 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \sin\left(\frac{4\pi x}{N}\right) dx \\ s_3 &= \int_0^N f(x) \sqrt{\frac{2}{N}} \sin\left(\frac{6\pi x}{N}\right) dx \\ \vdots \\ \end{aligned}\]

c_k,\; s_k 는 orthonormal 함수들의 계수들이고, 푸리에 급수를 처음 설명할 때 나왔던 식 (1)a_k,\; b_k 계수들을 구해보자.

    \[f(x) = a_0 + \sum_{k=1}^\infty \left\{ a_k \cos \left( \frac{2 k \pi x}{N} \right) + b_k \sin \left( \frac{2 k \pi x}{N} \right) \right\}\]

    \[\begin{aligned} a_0 &= \frac{c_0 g_0(x)}{1} = \frac{1}{N} \int_0^N f(x) dx \\ a_1 &= \frac{c_1 g_1(x)}{\cos\left(\frac{2\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \cos\left(\frac{2\pi x}{N}\right) dx \\ a_2 &= \frac{c_2 g_2(x)}{\cos\left(\frac{4\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \cos\left(\frac{4\pi x}{N}\right) dx \\ a_3 &= \frac{c_3 g_3(x)}{\cos\left(\frac{6\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \cos\left(\frac{6\pi x}{N}\right) dx \\ \vdots \\ b_1 &= \frac{s_1 h_1(x)}{\sin\left(\frac{2\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \sin\left(\frac{2\pi x}{N}\right) dx \\ b_2 &= \frac{s_2 h_2(x)}{\sin\left(\frac{4\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \sin\left(\frac{4\pi x}{N}\right) dx \\ b_3 &= \frac{s_3 h_3(x)}{\sin\left(\frac{6\pi x}{N}\right)} = \frac{2}{N} \int_0^N f(x) \sin\left(\frac{6\pi x}{N}\right) dx \\ \vdots \\ \end{aligned}\]

최종적으로 정리하면..

    \[f(x) = \frac{a_0}{2} + \sum_{k=1}^\infty \left\{ a_k \cos \left( \frac{2 k \pi x}{N} \right) + b_k \sin \left( \frac{2 k \pi x}{N} \right) \right\}\]

a_0 항을 a_k 항에 포함시키기 위해 원래의 a_0 값에 2 를 곱하고 위의 식에서는 다시 2 로 나눴다.

(3)   \[a_k = \frac{2}{N} \int_0^N f(x) \cos\left(\frac{2 k \pi x}{N}\right) dx \]

(4)   \[b_k = \frac{2}{N} \int_0^N f(x) \sin\left(\frac{2 k \pi x}{N}\right) dx \]

4.2 오일러 공식 (Euler’s Formula)

맥클로린 급수를 이용하면 어떠한 함수라도 급수의 형태로 전개할 수 있다.

    \[\begin{aligned} f(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 + ... \\ &= \sum_{n=0}^\infty \frac{1}{n!} f^{(n)} (0) x^n \end{aligned}\]

이를 이용해서 \cos(x), \sin(x), e^x 함수를 전개해 보자.

  • \cos(x)n 번 미분하여 x0 을 대입하면, 1, 0, -1, 0, 1, 0, -1, 0, ... (즉, 4 번 미분하면 0 번 미분한 것과 같다)
  • \sin(x)n 번 미분하여 x0 을 대입하면, 0, 1, 0, -1, 0, 1, 0, -1, ... (즉, 4 번 미분하면 0 번 미분한 것과 같다)
  • 따라서 \sin(x), \cos(x) 함수는 맥클로린 급수에서 각 항이 겹치지 않는다.
  • e^x 는 미분해도 e^x, x0 을 대입하면 항상 1

    \[\begin{aligned} \cos(x) &= 1 - \frac{1}{2!}x^2 + \frac{1}{4!}x^4 - \frac{1}{6!}x^6 + ... \\ \sin(x) &= x - \frac{1}{3!}x^3 + \frac{1}{5!}x^5 - \frac{1}{7!}x^7 + ... \\ e^x &= 1 + \frac{1}{1!}x + \frac{1}{2!}x^2 + \frac{1}{3!}x^3 + \frac{1}{4!}x^4 + ... \\ \end{aligned}\]

이제 e^{i x} 를 전개해 보면

e^{i x}\cos(x), \sin(x) 함수와 같이 4번 미분하면 0 번 미분한 것과 같다, x0 을 대입하면, 1, i, -1, -i, 1, i, -1, -i,... (반복)

따라서

    \[e^{i x} = \cos(x) + i \sin(x)\]

인 것을 알 수 있으며, 이게 오일러 공식이다. (x\pi 를 대입하면 그 유명한 오일러 등식(e^{i\pi} - 1 = 0 이 된다)

오일러 공식으로 부터 \sin(x), \cos(x) 같은 2\pi 주기 함수와 자연로그 상수 e, 그리고 허수 i 의 관계를 알 수 있다

    \[\begin{aligned} e^{i x} &= \cos(x) + i \sin(x) \\ e^{-i x} &= \cos(x) - i \sin(x) \\ \cos(x) &= \frac{e^{i x} + e^{-i x}}{2} \\ \sin(x) &= \frac{e^{i x} - e^{-i x}}{2i} \end{aligned}\]

euler_graph

4.3 푸리에 급수의 복소 표현

푸리에 급수식을 오일러 식을 이용해 고쳐써 보자. 먼저 푸리에 급수의 계수 공식 부터..

(5)   \[\begin{aligned} a_k &= \frac{2}{N} \int_0^N f(x) \cos\left(\frac{2 k \pi x}{N}\right) dx \\ &= \frac{2}{N} \int_0^N f(x) \frac{1}{2} \left( e^{i \frac{2 k \pi x}{N}} + e^{-i \frac{2 k \pi x}{N}} \right) dx \\ &= \frac{1}{N} \int_0^N f(x) \left( e^{i \frac{2 k \pi x}{N}} + e^{-i \frac{2 k \pi x}{N}} \right) dx  \end{aligned}\]

(6)   \[\begin{aligned} b_k &= \frac{2}{N} \int_0^N f(x) \sin\left(\frac{2 k \pi x}{N}\right) dx \\ &= \frac{2}{N} \int_0^N f(x) \frac{1}{2 i} \left( e^{i \frac{2 k \pi x}{N}} - e^{-i \frac{2 k \pi x}{N}} \right) dx \\ &= \frac{1}{N i} \int_0^N f(x) \left( e^{i \frac{2 k \pi x}{N}} - e^{-i \frac{2 k \pi x}{N}} \right) dx  \end{aligned}\]

푸리에 급수 공식도 바꿔보자

    \[\begin{aligned} f(x) &= \frac{a_0}{2} + \sum_{k=1}^\infty \left\{ a_k \cos \left( \frac{2 k \pi x}{N} \right) + b_k \sin \left( \frac{2 k \pi x}{N} \right) \right\} \\ &= \frac{a_0}{2} + \sum_{k=1}^\infty \left\{ \frac{a_k}{2} \left( e^{\frac{2 k \pi i x}{N} } + e^{-\frac{2 k \pi i x}{N} } \right) + \frac{b_k}{2 i} \left( e^{\frac{2 k \pi i x}{N} } - e^{-\frac{2 k \pi i x}{N} } \right) \right\} \\ &= \frac{a_0}{2} + \sum_{k=1}^\infty \left\{ \frac{1}{2} \left( a_k - i b_k \right) e^{\frac{2 k \pi i x}{N}} + \frac{1}{2} \left( a_k + i b_k \right) e^{-\frac{2 k \pi i x}{N} } \right\} \\ \end{aligned}\]

안쪽의 항을 식 (5)(6) 을 이용하여 고쳐쓰면

(7)   \[\frac{1}{2} \left( a_k - i b_k \right) = \frac{1}{N} \int_{0}^{N} f(x) e^{\frac{-2 k \pi i x}{N}} dx \]

(8)   \[\frac{1}{2} \left( a_k + i b_k \right) = \frac{1}{N} \int_{0}^{N} f(x) e^{\frac{2 k \pi i x}{N}} dx \]

(7)C_k 로 놓으면 식 (8)C_{-k} 이다. (둘은 복소평면에서 실수축으로 서로 대칭인 관계에 있다 – 켤레 복소수)

이제 다시 f(x) 를 정리하면..

    \[\begin{aligned} f(x) &= \frac{a_0}{2} + \sum_{k=1}^\infty \left\{ C_k e^{\frac{2 k \pi i x}{N}} + C_{-k} e^{-\frac{2 k \pi i x}{N} } \right\} \\ &= C_0 + \sum_{k=1}^\infty C_k e^{\frac{2 k \pi i x}{N}} + \sum_{k=-1}^{-\infty} C_{k} e^{\frac{2 k \pi i x}{N}} \\ &= \sum_{k=-\infty}^\infty C_k e^{\frac{2 k \pi i x}{N}} \\ \end{aligned}\]

이로써 푸리에 급수의 복소 표현을 얻었다.

(9)   \[f(x) = \sum_{k=-\infty}^\infty C_k e^{\frac{2 k \pi i x}{N}} \]

(10)   \[C_k = \frac{1}{N} \int_{0}^{N} f(x) e^{\frac{-2 k \pi i x}{N}} dx \]

5. 푸리에 변환 (Fourier Transform)

지금까지는 함수 f(x) 를 구간 [0, N] 에서 정의된 함수로 보았지만, 이제 구간을 [-\infty, \infty] 으로 확장시켜 보자.

먼저 기존의 식에서 구간을 [-\frac{N}{2}, \frac{N}{2}] 로 바꾸고, 식 (10) 를 식 (9) 에 대입해보자.

    \[f(x) = \sum_{k=-\infty}^\infty \left( \frac{1}{N} \int_{-\frac{N}{2}}^{\frac{N}{2}} f(x) e^{-2 \pi i t x} dx \right) e^{2 \pi i t x},\quad t = \frac{k}{N}\]

이제 N 을 무한대로 보내보자.

    \[\begin{aligned} f(x) &= \lim_{N\to\infty} \sum_{k=-\infty}^\infty \left( \frac{1}{N} \int_{-\frac{N}{2}}^{\frac{N}{2}} f(x) e^{-2 \pi i t x} dx \right) e^{2 \pi i t x} \\ &= \lim_{N\to\infty} \frac{1}{N} \sum_{k=-\infty}^\infty \left( \int_{-\frac{N}{2}}^{\frac{N}{2}} f(x) e^{-2 \pi i t x} dx \right) e^{2 \pi i t x} \\ &= \int_{-\infty}^\infty \left( \int_{-\infty}^{\infty} f(x) e^{-2 \pi i t x} dx \right) e^{2 \pi i t x} dt \\ \end{aligned}\]

괄호 안쪽의 식을 t 에 대한 함수 g(t) 로 놓으면

    \[g(t) = \int_{-\infty}^{\infty} f(x) e^{-2 \pi i t x} dx\]

    \[f(x) = \int_{-\infty}^\infty g(t) e^{2 \pi i t x} dt\]

가 되고, g(t) 를 함수 f(x) 에 있어서의 푸리에 변환 (Fourier Transform), f(x) 를 푸리에 역변환 (Inverse Fourier Transform) 이라고 한다.

5.1 DFT (Discrete Fourier Transform)

연속적인(continuous) 푸리에 변환 공식을 구간 [0, N-1] 의 N 개의 막대 그래프 형태로 생각하여 재구성

    \[\begin{aligned} C_k &= \frac{1}{N} \sum_{x=0}^{N-1} f(x) e^{\frac{-2 k \pi i x}{N}} \\ &= \frac{1}{N} \sum_{x=0}^{N-1} f(x) \left\{ \cos\left( \frac{2 k \pi x}{N}\right) - i \sin\left( \frac{2 k \pi x}{N}\right) \right\}  \end{aligned}\]

    \[\begin{aligned} f(x) &= \sum_{k=0}^{N-1} C_k e^{\frac{2 k \pi i x}{N}} \\ &= \sum_{k=0}^{N-1} C_k \left\{ \cos\left( \frac{2 k \pi x}{N} \right) + i \sin\left( \frac{2 k \pi x}{N} \right) \right\}  \end{aligned}\]

여기서 C_k 를 구할때 필요한 e^{\frac{-2 k \pi i}{N}} 를 주목해보자. e^{\frac{-2 k \pi i}{N}} 는 복소평면의 단위원을 N 등분한 것을 k 가 증가함에 따라 시계방향으로 도는 그림이다. 아래는 N = 8 일때를 나타낸 그림이다.

rotation-operator-w8

k \ge \frac{N}{2} 부터 시계 반대방향으로 도는 켤레복소수라고 생각해도 된다. 만약 f(x) 가 실수함수라면 C_kk = \frac{N}{2} 부터 실수부는 같고 허수부만 부호가 다른 대칭형태가 된다는 것을 알 수 있다.

또한 e^{\frac{-2 k \pi i}{N}} 또는 e^{\frac{2 k \pi i}{N}} 는 복소평면을 회전하는 값이므로 계산과정에서 중복되는 부분이 생기게 된다. 이를 이용해 연산을 줄이는 알고리즘이 FFT 다.

5.2 FFT (Fast Fourier Transform)

  • N 개의 샘플값에 DFT 연산을 하면 N^2 번의 복소수 곱셈이 필요 – 실시간으로 사용하기엔 부담스러움
  • N 값을 특정 기수(radix) 에 맞추고, 공식을 분할 후 중복요소를 찾아낸다.
  • FFT 알고리즘은 그 변종이 굉장히 많지만, 여기서는 가장 일반적인 Radix-2 FFT 와, Radix-4 FFT 만 다룬다.

5.2.1 Radix-2 FFT

  • N = 2^n 일때 사용할 수 있다.
  • N^2 보다 훨씬 적은 \frac{N}{2} log_2 N 번의 복소수 곱셈이 필요하다.

DFT 식을 좀더 간단하게 표현하기 위해 먼저 아래의 복소수를 정의하자.

    \[W_N^{k} = e^{\frac{-2 k \pi i}{N}}\]

W_N^{k} 는 복소평면의 단위원을 N등분 한것이고, e 의 지수가 음수이니 k 가 증가함에 따라 시계방향으로 단위원을 도는 그림을 상상하면 된다. 그럼 다음의 성질을 갖는다는 것을 쉽게 알 수 있다.

(12)   \[W_N^{ N } = W_N^{ -N } = 1 \]

(13)   \[W_N^{ \frac{N}{2} } = W_N^{ -\frac{N}{2} }  = -1 \]

(14)   \[W_N^{ \frac{N}{4} } = -i \]

(15)   \[W_N^{ -\frac{N}{4} } = i \]

이제 푸리에 급수의 계수 공식을 x 가 홀수일 때 와 짝수일 때 별로 나눠본다. (\frac{1}{N} 은 일단 생략)

    \[\begin{aligned} C_k &= \sum_{x=0}^{N-1}f(x) W_N^{k x} \\ &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_N^{2 k x} + \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_N^{k (2x + 1)} \\ &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_N^{2 k x} + W_N^{k} \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_N^{2 k x} \\ \end{aligned}\]

W_N^2 = W_{\frac{N}{2}} 이므로,

    \[C_k = \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_{\frac{N}{2}}^{k x} + W_N^{k} \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_{\frac{N}{2}}^{k x}\]

위 식에서

    \[\begin{aligned} Y_k &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_{\frac{N}{2}}^{k x} \\ Z_k &= \sum_{x=0}^{\frac{N}{2}-1} f(2x + 1) W_{\frac{N}{2}}^{k x}  \\ \end{aligned}\]

이라고 놓으면 C_{k} 는,

(16)   \[C_k = Y_k + W_N^{k} Z_k \]

이제 C_{k+\frac{N}{2}} 를 구해보면..

(17)   \[\begin{aligned} C_{k+\frac{N}{2}} &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_{\frac{N}{2}}^{(k + \frac{N}{2}) x} + W_N^{(k + \frac{N}{2})} \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_{\frac{N}{2}}^{(k + \frac{N}{2}) x} \\ &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_{\frac{N}{2}}^{k x} W_{\frac{N}{2}}^{\frac{N}{2}x} + W_N^{k} W_N^{\frac{N}{2}} \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_{\frac{N}{2}}^{k x} W_{\frac{N}{2}}^{\frac{N}{2}x}\\ &= \sum_{x=0}^{\frac{N}{2}-1} f(2x) W_{\frac{N}{2}}^{k x} - W_N^{k} \sum_{x=0}^{\frac{N}{2}-1} f(2x+1) W_{\frac{N}{2}}^{k x} \\ &= Y_k - W_N^{k} Z_k  \\ \end{aligned}\]

위의 식 (16), (17)Cooley-Tukey butterfly operation 이라고 하고 다이어그램으로 그리면 아래처럼 나비모양으로 나타낼 수 있다.

Cooley-Tukey-butterfly

Y_k, Z_k 는 재귀적으로 계속 쪼갤 수 있다. log_2 N 번 만큼 쪼개면 Radix-2 FFT 알고리즘이 된다.

ex) N = 8 일때 Radix-2 FFT 연산과정 (8^2 = 64번의 복소수 곱셈 연산이 \frac{8}{2}log_2{8} = 12번의 복소수 곱셈 연산으로 줄었다)

FFT-butterfly

5.1.2 Radix-4 FFT

  • N = 4^n 일때 사용할 수 있다.
  • Radix-2 FFT 와 방식은 비슷하나, 항을 4개로 분할한다.
  • Radix-2 FFT 보다 약간 적은 \frac{3N}{4} log_4 N 번의 복소수 곱셈이 필요하다.

    \[\begin{aligned} C_k &= \sum_{x=0}^{N-1}f(x) W_N^{k x} \\ &= \sum_{x=0}^{\frac{N}{4}-1} f(4x) W_N^{4 k x} + \sum_{x=0}^{\frac{N}{4}-1} f(4x+1) W_N^{k (4x + 1)} + \sum_{x=0}^{\frac{N}{4}-1} f(4x+2) W_N^{k (4x + 2)} + \sum_{x=0}^{\frac{N}{4}-1} f(4x+3) W_N^{k (4x + 3)} \\ &= \sum_{x=0}^{\frac{N}{4}-1} f(4x) W_N^{4 k x} + W_N^{k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+1) W_N^{4 k x} + W_N^{2k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+2) W_N^{4 k x} + W_N^{3k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+3) W_N^{4 k x} \\ &= \sum_{x=0}^{\frac{N}{4}-1} f(4x) W_{\frac{N}{4}}^{k x} + W_N^{k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+1) W_{\frac{N}{4}}^{k x} + W_N^{2k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+2) W_{\frac{N}{4}}^{k x} + W_N^{3k} \sum_{x=0}^{\frac{N}{4}-1} f(4x+3) W_{\frac{N}{4}}^{k x} \end{aligned}\]

위 식에서

    \[\begin{aligned} Y_k &= \sum_{x=0}^{\frac{N}{4}-1} f(4x) W_{\frac{N}{4}}^{k x} \\ Z_k &= \sum_{x=0}^{\frac{N}{4}-1} f(4x + 1) W_{\frac{N}{4}}^{k x}  \\ U_k &= \sum_{x=0}^{\frac{N}{4}-1} f(4x + 2) W_{\frac{N}{4}}^{k x}  \\ V_k &= \sum_{x=0}^{\frac{N}{4}-1} f(4x + 3) W_{\frac{N}{4}}^{k x}  \\ \end{aligned}\]

이라고 놓고 C_{k}, C_{k+\frac{N}{4}}, C_{k+\frac{2N}{4}}, C_{k+\frac{3N}{4}} 을 구해보면..

    \[\begin{aligned} C_{k} &= Y_k + W_N^{k} Z_k + W_N^{2k} U_k + W_N^{3k} V_k \\ C_{k+\frac{N}{4}} &= Y_k + W_N^{k+\frac{N}{4}} Z_k + W_N^{2(k+\frac{N}{4})} U_k + W_N^{3(k+\frac{N}{4})} V_k \\ C_{k+\frac{2N}{4}} &= Y_k + W_N^{k+\frac{2N}{4}} Z_k + W_N^{2(k+\frac{2N}{4})} U_k + W_N^{3(k+\frac{2N}{4})} V_k \\ C_{k+\frac{3N}{4}} &= Y_k + W_N^{k+\frac{3N}{4}} Z_k + W_N^{2(k+\frac{3N}{4})} U_k + W_N^{3(k+\frac{3N}{4})} V_k \\ \end{aligned}\]

(12), (13), (14), (15) 에 의해서

    \[\begin{aligned} C_{k} &= Y_k + W_N^{k} Z_k + W_N^{2k} U_k + W_N^{3k} V_k \\ C_{k+\frac{N}{4}} &= Y_k - i W_N^{k} Z_k - W_N^{2k} U_k + i W_N^{3k} V_k \\ C_{k+\frac{2N}{4}} &= Y_k - W_N^{k} Z_k + W_N^{2k} U_k - W_N^{3k} V_k \\ C_{k+\frac{3N}{4}} &= Y_k + i W_N^{k} Z_k - W_N^{2k} U_k - i W_N^{3k} V_k \end{aligned}\]

다시 정리하면

    \[\begin{aligned} C_{k} &= (Y_k + W_N^{2k} U_k) + (W_N^{k} Z_k + W_N^{3k} V_k) \\ C_{k+\frac{N}{4}} &= (Y_k  - W_N^{2k} U_k) - i (W_N^{k} Z_k - W_N^{3k} V_k) \\ C_{k+\frac{2N}{4}} &= (Y_k + W_N^{2k} U_k) - (W_N^{k} Z_k + W_N^{3k} V_k) \\ C_{k+\frac{3N}{4}} &= (Y_k - W_N^{2k} U_k) + i (W_N^{k} Z_k - W_N^{3k} V_k) \end{aligned}\]

Radix-4-butterfly

Y_k, Z_k, U_k, V_k 를 재귀적으로 계속 쪼갤 수 있다. log_4 N 번 만큼 쪼개면 Radix-4 FFT 알고리즘이 된다.

ex) N = 16 일때 Radix-4 FFT 연산과정 (16^2 = 256번의 복소수 곱셈 연산이 \frac{16 \times 3}{4}log_4{16} = 24번의 복소수 곱셈 연산으로 줄었다. Radix-2 FFT 였다면 \frac{16}{2}log_2{16} = 32)

FFT-radix4

6. DCT (Discrete Cosine Transform)

  • DCT 는 DFT 의 cosine 성분만을 사용 (sine 성분이 없으므로 2 \pi 가 아닌 \pi 를 기본 주기로 해도 된다)
  • DFT 와는 달리 허수 성분이 없으므로 계산이 간단해서 더 널리 쓰임 (JPEG 이미지 압축, MPEG, Theora, H.264 비디오 압축)

0 부터 N-1 까지 N 개의 데이터에 대하여 정수 x 는 막대의 중간점인 x + \frac{1}{2} 이라고 놓고 다시 써보면..

    \[C_k = u_k \sum_{x=0}^{N-1} f(x) \cos\left(\frac{k \pi (2x + 1)}{2N}\right)\]

    \[f(x) = \sum_{k=0}^{N-1} u_k C_k \cos\left(\frac{k \pi (2x + 1)}{2N}\right)\]

    \[u_k = \begin{cases}  \sqrt{\frac{1}{N}} & \quad \text{if}\;k = 0 \\  \sqrt{\frac{2}{N}} & \quad \text{if}\;k \neq 0 \\  \end{cases}\]

상수 C_k 성분을 구하는 과정을 FDCT (Forward DCT) 라고 하고, 반대로 C_k 에서 f(x) 를 구하는 과정을 IDCT (Inverse DCT) 라고 한다

(※ u_k 를 양쪽에서 같은 값으로 만들어 FDCT, IDCT 가 대칭이 되게 만들었다)

1 차원 DCT 를 알았으니 2 차원에서 적용시켜 보자

    \[\begin{aligned} C_{kl} &= u_k v_l \sum_{y=0}^{M-1} \sum_{x=0}^{N-1} f(x, y) \cos\left(\frac{k \pi (2x + 1)}{2N}\right) \cos\left(\frac{l \pi (2y + 1)}{2M}\right) \\ &= v_l \sum_{y=0}^{M-1} \left\{ u_k \sum_{x=0}^{N-1} f(x, y) \cos\left(\frac{k \pi (2x + 1)}{2N}\right) \right\} \cos\left(\frac{l \pi (2y + 1)}{2M}\right) \end{aligned}\]

    \[\begin{aligned} f(x, y) &= \sum_{l=0}^{M-1} \sum_{k=0}^{N-1} u_k v_l C_{kl} \cos\left(\frac{k \pi (2x + 1)}{2N}\right) \cos\left(\frac{l \pi (2y + 1)}{2M}\right) \\ &= \sum_{l=0}^{M-1} v_l \left\{ \sum_{k=0}^{N-1} u_k C_{kl} \cos\left(\frac{k \pi (2x + 1)}{2N}\right) \right\} \cos\left(\frac{l \pi (2y + 1)}{2M}\right) \end{aligned}\]

    \[u_k = \begin{cases}  \sqrt{\frac{1}{N}} & \quad \text{if}\;k = 0 \\  \sqrt{\frac{2}{N}} & \quad \text{if}\;k \neq 0 \\  \end{cases}\]

    \[v_l = \begin{cases}  \sqrt{\frac{1}{M}} & \quad \text{if}\;l = 0 \\  \sqrt{\frac{2}{M}} & \quad \text{if}\;l \neq 0 \\  \end{cases}\]

6.1 빠른 DCT 알고리즘 #1 – FFT 이용

임의의 함수 f(x) 가 있다고 할때 x = 0 을 기준으로 대칭인 \hat{f}(x) 를 정의한다.

    \[\hat{f}(x) = \begin{cases}  f(x) & \quad (0 \le x \le N-1) \\ f(-x - 1) & \quad (-N \le x \le -1) \end{cases}\]

구간 [-N, N] 에서 \hat{f}(x) 의 DFT 공식은

    \[C_k = \frac{1}{2N} \sum_{x=-N}^{N-1} \hat{f}(x) e^{-\frac{\pi k i \hat{x}}{N}}\]

\hat{x} = x + 0.5 라고 놓고 위의 식을 \hat{x} 로 다시 쓰면

    \[\begin{aligned} C_k &= \frac{1}{2N} \sum_{\hat{x}=-N+0.5}^{N-0.5} \hat{f}(\hat{x} - 0.5) e^{-\frac{\pi k i (\hat{x} - 0.5)}{N}} \\ &= \frac{1}{2N} \sum_{\hat{x}=-N+0.5}^{N-0.5} \hat{f}(\hat{x} - 0.5) e^{-\frac{\pi k i \hat{x}}{N}} e^{-\frac{\pi k i}{2N}} \\ &= \frac{1}{2N} \sum_{\hat{x}=-N+0.5}^{N-0.5} \hat{f}(\hat{x} - 0.5) \left( \cos \left( \frac{\pi k \hat{x}}{N} \right) - i \sin \left( \frac{\pi k \hat{x}}{N} \right) \right) e^{-\frac{\pi k i}{2N}} \\ &= \frac{1}{N} \sum_{\hat{x}=0.5}^{N-0.5} \hat{f}(\hat{x} - 0.5) \cos \left( \frac{\pi k \hat{x}}{N} \right) e^{-\frac{\pi k i}{2N}} \\ &= \frac{1}{N} \sum_{x=0}^{N-1} \hat{f}(x) \cos \left( \frac{\pi k (x + 0.5)}{N} \right) e^{-\frac{\pi k i}{2N}} \\ &= DCT ( f(x) ) e^{-\frac{\pi k i}{2N}} \end{aligned}\]

따라서 DCT 는 DFT 를 이용해서 구할 수 있다.

    \[DCT ( f(x) ) = DFT ( \hat{f}(x) [-N, N] ) e^{\frac{\pi k i}{2N}}\]

6.2 빠른 DCT 알고리즘 #2 – DCT 행렬 분해 (Matrix Decomposition)

  • cosine 의 대칭성을 이용
  • 8 개의 샘플에 22 번의 곱셈과 28 번의 덧셈 – 계산 의존도가 없으므로 병렬로 실행 가능 (GPU 연산에 적합한 형태)

    \[\begin{bmatrix} C_0 \\ C_1 \\ C_2 \\ C_3 \\ C_4 \\ C_5 \\ C_6 \\ C_7 \end{bmatrix} = \frac{1}{ \sqrt{8} } \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ a & c & d & f & -f & -d & -c & -a \\ b & e & -e & -b & -b & -e & e & b \\ c & -f & -a & -d & d & a & f & -c \\ 1 & -1 & -1 & 1 & 1 & -1 & -1 & 1 \\ d & -a & f & c & -c & -f & a & d \\ e & -b & b & -e & -e & b & -b & e \\ f & -d & c & -a & a & -c & d & -f \end{bmatrix} \begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ f_3 \\ f_4 \\ f_5 \\ f_6 \\ f_7 \end{bmatrix}\]

    \[a = \sqrt{2} \cos \left( \frac{\pi}{16} \right),\;  b = \sqrt{2} \cos \left( \frac{2 \pi}{16} \right),\; c = \sqrt{2} \cos \left( \frac{3 \pi}{16} \right),\;  d = \sqrt{2} \cos \left( \frac{5 \pi}{16} \right),\; e = \sqrt{2} \cos \left( \frac{6 \pi}{16} \right),\; f = \sqrt{2} \cos \left( \frac{7 \pi}{16} \right)\]

위의 8-DCT 행렬의 대칭성을 이용해서 아래처럼 분해가 가능하다

    \[\begin{bmatrix} C_0 \\ C_2 \\ C_4 \\ C_6 \end{bmatrix} = \frac{1}{ \sqrt{8} } \begin{bmatrix} 1 & 1 & 1 & 1 \\ b & e & -e & -b \\ 1 & -1 & -1 & 1 \\ e & -b & b & -e \end{bmatrix} \begin{bmatrix} f_0 + f_7 \\ f_1 + f_6 \\ f_2 + f_5 \\ f_3 + f_4 \end{bmatrix}\]

    \[\begin{bmatrix} C_1 \\ C_3 \\ C_5 \\ C_7 \end{bmatrix} = \frac{1}{ \sqrt{8} } \begin{bmatrix} a & -c & d & -f \\ c & f & -a & -d \\ d & a & f & -c \\ f & d & c & a \end{bmatrix} \begin{bmatrix} f_0 - f_7 \\ f_6 - f_1 \\ f_2 - f_5 \\ f_4 - f_3 \end{bmatrix}\]

위의 첫번째 식은 같은 방식으로 더 분해가 가능하다.

여기서

    \[\begin{aligned} f_{07} = f_0 + f_7 \\ f_{16} = f_1 + f_6 \\ f_{25} = f_2 + f_5 \\ f_{34} = f_3 + f_4 \\ \\ f_{0734} = f_{07} + f_{34} \\ f_{1625} = f_{16} + f_{25}  \end{aligned}\]

라고 하면

    \[\begin{aligned} \begin{bmatrix} C_0 \\ C_4 \end{bmatrix} &= \frac{1}{ \sqrt{8} } \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} f_{07} + f_{34} \\ f_{16} + f_{25} \end{bmatrix} \\ &= \frac{1}{ \sqrt{8} }  \begin{bmatrix} f_{0734} + f_{1625} \\ f_{0734} - f_{1625} \end{bmatrix} \end{aligned}\]

    \[\begin{bmatrix} C_2 \\ C_6 \end{bmatrix} = \frac{1}{ \sqrt{8} } \begin{bmatrix} b & e \\ e & -b \end{bmatrix} \begin{bmatrix} f_{07} - f_{34} \\ f_{16} - f_{25} \end{bmatrix}\]

7. 데모 (Demo)

8. 참고 (Reference)

Leave a Reply

Your email address will not be published. Required fields are marked *