지상에서 포물선 운동을 하는 물체의 속도 벡터를 모두 모아 시작을 같게 만들면 속도 벡터의 끝이 그리는 자취(hodograph라 함)는 직선이 된다. 이는 속도의 차이가 가속도이고 지상에서 중력가속도는 크기와 방향이 일정하기 때문이다.
$$\vec {v}(t) =\vec {v}_0 - g\hat {j} t$$
태양계에서 행성은 태양을 한 초점으로 하는 타원 궤도상에서 움직인다. 그러면 행성의 각 지점에서 속도 벡터의 시작점을 한 지점에 모을 때 벡터의 끝점이 그리는 궤적은 무엇일까? 원 궤도를 그리면 속력이 일정하므로 당연히 원이 될 것으로 예상할 수 있지만, 타원 궤도에서는 속력은 에너지 보존을 고려하면 일정할 수 없다. 그런데 타원궤도를 그리는 경우에도 속도의 hodograph는 원궤도에서와 마찬가지로 벡터 공간에서 원으로 표현된다. 왜 그럴까?
각운동량의 크기가 일정하므로 속도의 증가분이 회전각의 증가분과 비례함을 의미하고, 기하학적으로는 속도 벡터를 모아둔 벡터 공간에서 각이 증가하면서 움직이는 벡터의 끝이 반지름 $GMm/L$인 원주 위에 있을 때만 가능한 관계이다. 벡터 공간에서 벡터의 시작점이 각을 재는 원점이 아니다. hodograph의 반지름은 근일점에서 속력 $v_P$와 원일점에서 속력 $v_A$의 절반임으로 표현되는데 $v_R = (v_A + v_P)/2 = \frac {GMm}{L}$임을 근일점, 원일점에서 역학적 에너지 보존과 각운동량 보존식을 이용해서 확인할 수 있다.
태양계의 행성 운동에서 각운동량은 보존이 된다(Kepler의 제2법칙). 이는 태양이 행성에 작용하는 힘인 만유인력이 중심력의 형태를 띠고 있기 때문이다. 식으로는 태양에서 행성까지 위치 벡터를 $\vec {r}$이라면 태양이 행성에 작용하는 만유인력은 $\vec {F} = \frac {GMm}{r^3}\vec {r}$로 쓸 수 있으므로 만유인력이 만드는 토크가 $\vec {\tau} = \vec {r}\times \vec {F}=0$임을 쉽게 알 수 있다. 따라서 각운동량 보존은 자명해진다.
만약 행성의 위치를 재는 원점을 태양이 아니라 다른 지점으로 잡으면 어떻게 될까? 이 경우 만유인력의 방향과 위치 벡터의 방향이 나란하지 않으므로 토크가 0이 안된다. 그럼 각운동량은 원점을 어디로 잡는가에 따라 보존되기도 하고 안되기도 하는 물리량일까? 무엇을 놓치고 있는 것일까?
사람이 걷는 동작은 복잡하지만 몇 가지 가정을 하면 단순한 강체의 운동으로 근사를 할 수 있다. 우선 사람이 걷는 동안 항상 한 발은 땅을 딛고 있다. 그리고 땅을 딛고 있는 발을 기준으로 몸의 질량중심은 원호를 그리면서 거의 일정한 속도로 움직이게 된다. 사람을 질량중심에 모든 질량이 뭉친 점으로 근사를 하면 걷는 동작은 거의 질량이 없는 막대(다리: 길이=$L$)에 매달린 역진자(inverted pendulum) 운동과 비슷하다고 볼 수 있다. 또한 원호를 따라 중심이 이동하는 속력은 거의 일정하다고 근사할 수 있다.
이 경우에 질량중심은 등속 원운동을 한다고 볼 수 있고, 뉴턴의 제2법칙에 의해서 질량중심에서 발 쪽을 향하는 힘 성분이 구심력 역할을 한다. 구심력 역학을 할 수 있는 힘은 중력, 수직항력, 그리고 마찰력이 있다.
사람을 너무 단순화시킨 감은 있지만 왜 다리가 긴 사람이 더 빨리 걸을 수 있는지를 물리 법칙으로 추정할 수 있다는 것을 보여준다. 사람을 좀 더 복잡한 강체로 근사하더라도 $\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}$이므로 합리적인 추정이 된다.
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}
이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 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에서 찾을 수 있다.
$F_n$ 인덱스의 이진비트가 twiddlefactor를 어떻게 곱해야 하는지 알려주고, 비트를 역으로 읽으면 계수가 나타나는 순서를 알 수 있게 한다.
#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;
}
Edge를 보존하는 low-pass filter로 주어진 픽셀 값을 그 위치에서 edge 방향을 따라 구한 평균값으로 대체한다. 방향은 8방향 (45도) 기준이다. Edge의 방향은 그림처럼 중심 픽셀을 기준으로 주변의 8개를 픽셀을 이용해서 45도씩 회전하면서 차이(평균 변화율)를 계산하여 가장 작은 값을 주는 방향으로 잡는다. 중심 픽셀이 line[1][1]일 때
0도 에지: abs(line[1][0] - line[1][2])가 최소
45도 에지: abs(line[0][0] - line[2][2])가 최소
90도 에지: abs(|ine[0][1] - line[2][1])가 최소
135도 에지: abs(line[0][2] - line[2][0])가 최소
Edge 방향의 3 픽셀 평균값으로 중심 픽셀 값을 대체하면 edge는 보존하면서 mean filter를 적용한 효과를 얻을 수 있다.
