Geometry & Recognition

사람이 걷는 동작은 복잡하지만 몇 가지 가정을 하면 단순한 강체의 운동으로 근사를 할 수 있다. 우선 사람이 걷는 동안 항상 한 발은 땅을 딛고 있다. 그리고 땅을 딛고 있는 발을 기준으로 몸의 질량중심은 원호를 그리면서 거의 일정한 속도로 움직이게 된다. 사람을 질량중심에 모든 질량이 뭉친 점으로 근사를 하면 걷는 동작은 거의 질량이 없는 막대(다리: 길이=$L$)에 매달린 역진자(inverted pendulum) 운동과 비슷하다고 볼 수 있다. 또한 원호를 따라 중심이 이동하는 속력은 거의 일정하다고 근사할 수 있다.

이 경우에 질량중심은 등속 원운동을 한다고 볼 수 있고, 뉴턴의 제2법칙에 의해서 질량중심에서 발 쪽을 향하는 힘 성분이 구심력 역할을 한다. 구심력 역학을 할 수 있는 힘은 중력, 수직항력, 그리고 마찰력이 있다.

수직에서 $\theta$만큼 기울어졌을 때 등속원운동의 가속도 벡터는

$$ \vec{a} =\frac{V^2}{L} (-\sin \theta \hat{x} -\cos \theta \hat{y})$$

이고,  운동 방정식의 수식 성분을 보면 중력과 수직항력이 기여하므로

$$ \sum F_y = N - mg = m a_y = -\frac{m V^2}{L} \cos \theta $$

움직이는 속력이 너무 빠르면 달리는 동작이 되고, 이 경우 양발이 땅에서 떨어지게 된다. 발이 땅에서 떨어지지 않으려면  수직항력 $N\ge0$인 조건을 만족하도록 걷는 속력을 조절해야 한다. 

$$N \ge  0 \quad \rightarrow  \quad V \le  \sqrt{\frac{gL }{\cos \theta}}$$

원호를 그리며 움직이는 동안 이 관계가 성립해야 하므로 걷기가 가능한 최대속력은:

$$V_\text{max} = \sqrt{gL}$$

사람을 너무 단순화시킨 감은 있지만 왜 다리가 긴 사람이 더 빨리 걸을 수 있는지를 물리 법칙으로 추정할 수 있다는 것을 보여준다. 사람을 좀 더 복잡한 강체로 근사하더라도 $\sqrt{gL}$ 앞의 factor만 바뀔 것이다.  사람이 걷는 동안 땅을 딛지 않는 발은 약간 구부러진 상태로 엉덩이를 축으로 일종의 물리진자처럼 행동하는데 이를 이용해도 역시 같은 추정치를 얻을 수 있다.

다리 길이가 $L=1\text{m}$이면,  $V_\text{max}= \sqrt{(1 \text{m}) \times (9.8\text{m/s}^2)} = 3.13\text{m/s}=11.3\text{km/h}$. 경보 세계기록이 대략 $15\text{km/h}$이므로 합리적인 추정이 된다.

youtu.be/HypJY1XWkWY

FFT는 입력 데이터의 개수$(N)$가 2의 지수승으로 주어질 때  $O(N\log N)$의 연산만으로 빠르게 DFT을 수행하는 알고리즘이다. 주어진 $N$개의 data $\{F_0, F_1, ..., F_{N-1}\}$의  DFT $\{a_0, a_1,..., a_{N-1}\}$는

$$ a_k = \sum_{n = 0}^{N-1} F_n \exp\Big( -i \frac{2\pi nk }{N} \Big)=\sum_{n = 0}^{N-1} F_n (W^n )^k , \quad W = \exp\Big( -i \frac{2\pi}{N}\Big)$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:

\begin{align} a_k &= \sum _{n=0}^{N/2-1} F_{2n} (W^{2n})^k+W^{k} \sum_{n=0}^{N/2-1} F_{2n+1}(W^{2n})^k\\ &=\text{(N/2개 even 항에 대한 DFT)} + W^k \times \text{(N/2개 odd 항에 대한 DFT)}\\ &= e_k + W^k o_k \end{align}

더불어 주기성 때문에 $a_k$의 계산도 절반만 하면 된다: 

$$a_{k+ N/2} = e_k - W^k o_k, \quad k = 0,1,..., N/2-1.$$

이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 divide and conquer 기법을 적용할 수 있는 구조이므로 매우 빠르게 연산을 수행할 수 있다(Cooley-Tukey algorithm). 역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im}  \rightarrow  - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22에서 찾을 수 있다.

N=8일 때 전개:

\begin{align}&a_k \\ &= \sum_{n=0}^7 F_n (W^{n})^k  \\ &= \sum_{n=0}^3 F_{2n} (W^{2n})^k +  W^k \sum_{n=0}^3 F_{2n+1}( W^{2n})^k   \\ &= \left[ \sum _{n=0}^1 F_{4n} (W^{4n})^k  + W^{2k} \sum _{n=0}^1 F_{4n +2} (W^{4n})^k\right]+W^k \left[\sum _{n=0}^1 F_{4n+1} (W^{4n})^k + W^{2k} \sum_{n=0}^1 F_{4n+3} (W^{4n})^k\right]\\ &= \left[ (F_0 + W^{4k} F_4) + W^{2k} (F_2 + W^{4k} F_6) \right]  + W^k \left[(F_1 + W^{4k} F_5 )+W^{2k} (F_3 + W^{4k} F_7) \right]\\ &= \left[ (F_{000} + W^{4k} F_{100}) + W^{2k} (F_{010} + W^{4k} F_{110}) \right] \\ &\qquad\qquad\qquad\qquad+ W^k \left[(F_{001} + W^{4k} F_{101} )+W^{2k} (F_{011} + W^{4k} F_{111}) \right]\end{align}

$F_n$ 인덱스의 이진비트가 twiddle factor를 어떻게 곱해야 하는지 알려주고, 비트를 역으로 읽으면 계수가 나타나는 순서를 알 수 있게 한다.

#define TWOPI 6.28318530717958647693
void split(int N, double* data) {
    if (N == 8) {
        double t1 = data[4]; data[4] = data[1];
        double t2 = data[2]; data[2] = t1; data[1] = t2;
        t1 = data[5]; data[5] = data[3];
        t2 = data[6]; data[6] = t1; data[3] = t2; 
        return;
    }
    double *tmp = new double [N / 2];
    for (int i = 0; i < N / 2; i++)   //copy odd elements to tmp buffer;
        tmp[i] = data[2 * i + 1];
    for (int i = 0; i < N / 2; i++)   // move even elements to lower half;
        data[i] = data[2 * i];            
    for (int i = 0; i < N / 2; i++)   // move odd elements to upper half;
        data[i + N / 2] = tmp[i];     
    delete [] tmp;
}
void fft4(double *re, double *im) {
    double tr1 = re[0]-re[2]+im[1]-im[3];
    double ti1 = im[0]-im[2]-re[1]+re[3];
    double tr2 = re[0]+re[2]-re[1]-re[3];
    double ti2 = im[0]+im[2]-im[1]-im[3];
    double ti3 = im[0]-im[2]+re[1]-re[3];
    double tr3 = re[0]-re[2]-im[1]+im[3];
    re[0] = re[0]+re[2]+re[1]+re[3];
    im[0] = im[0]+im[2]+im[1]+im[3];
    re[1] = tr1; im[1] = ti1;
    re[2] = tr2; im[2] = ti2;
    re[3] = tr3; im[3] = ti3;
}
int fft ( int N, double* re, double* im ) {
    if (N <= 0 || N & (N-1)) return 0; // N is not a power of 2;
    if ( N < 2 ) return 1;
    else if (N == 2) {
        double tr1 = re[0]-re[1], ti1 = im[0]-im[1];
        re[0] = re[0]+re[1]; im[0] = im[0]+im[1]; 
        re[1] = tr1; im[1] = ti1; 
        return 1;
    }  else if (N == 4) {
        fft4(re, im);
        return 1;
    }
        
    split ( N, re);
    split ( N, im);
    fft ( N / 2, &re[    0], &im[    0] );
    fft ( N / 2, &re[N / 2], &im[N / 2] );
    for ( int k = 0; k < N / 2; k++ ) {
        double Ere = re[k];
        double Eim = im[k];
        double Ore  = re[k + N / 2];
        double Oim  = im[k + N / 2];
        double Wre =  cos ( TWOPI * k / N );  //real part of twiddle factor:W
        double Wim =  -sin ( TWOPI * k / N ); //imag part of twiddlw factor:W
        double WOre = Wre * Ore - Wim * Oim;
        double WOim = Wre * Oim + Wim * Ore;
        re[k] = Ere + WOre;
        im[k] = Eim + WOim;
        re[k + N / 2] = Ere - WOre;
        im[k + N / 2] = Eim - WOim;
    }
    return 1;
};

6개의 주파수 $\{f_i = 2, 5, 9, 11, 21, 29\text{Hz}\}$가 섞인 신호 $s(t)$를 $[0,1]$ 초 구간에서 일정한 간격으로 sampling 해서 64개 data를 얻고(Nyquist frequency = 32Hz), 이에 대해 FFT를 수행한다.

$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$

int fft_test_main() {
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for (int i = 0; i < samples; i++) {
        double signal = 0;
        for (int k = 0; k < nfreqs; k++) 
            signal += (k + 1) * cos(TWOPI * freq[k] * i / samples);
        re[i] = signal;
        im[i] = 0.;
    }
    fft(samples, &re[0], &im[0]);
    for (int i = 0; i < samples; i++) TRACE("%d, %f, %f\n", i, re[i], im[i]);
    return 0;
}

FFT 결과:

FFT2D