주어진 함수 $y(x)$를 구간 $0\le x < L$에서 일정한 간격으로 샘플링해서 얻은 $N$개의 data $\{y_n = y(nL/N)|n=0, 1,...,N-1\} $의 discrete Fourier transform(DFT) $Y_k$는
$$ Y_ k = \frac {1}{N} \sum _{n = 0}^{N-1} y_n e^{ -\frac {2\pi i}{N} n k}, $$
로 정의된다. 역변환(IDFT)은 $Y_k$에서 $y_n$을 얻는 과정으로 다음 식으로 계산된다:
$$ y_n = \sum_{k=0}^{N-1} Y_k e^{+\frac {2\pi i}{N} k n}. $$
이산적인 data $y_n$만을 가지고 $x_n = nL/N$에서 미분계수를 구하려면 어떻게 해야 할까? 샘플링 점 사이 구간에서 값을 추정할 수 있도록 보간 함수를 만들어 사용해야 한다. 간단하게는 IDFT 식에서 $n$번째 샘플링 위치인 $x_n = n(L/N)$을 실수 $x$로 대체하여 연속함수를 얻는 방법을 생각할 수 있다. 그리고 IDFT 식에서 $ Y_k e^{+\frac {2\pi i}{N} k n}$ 대신에 추가적으로 위상 요소를 곱해서 $Y_k e^{+\frac {2\pi i}{N} (k + m_k N) n}, ~m_k =\text {integer}$ 로 바꾸더라도 동일한 원 데이터 $y_n$을 얻는다는 사실을 이용하도록 하자. 이 두 가지 관찰을 이용해서 만든 보간 함수는
$$ y(x) = \sum _{k=0}^{N-1} Y_k e^{+\frac{2\pi i}{L} (k +m_k N ) x} $$
이다. 추가 위상이 $y_n$은 바꾸지 않지만 어떤 선택을 하는지에 따라 샘플링 사이 구간에서 다른 진동을 보이므로 미분계수가 영향을 받게 된다. 그러면 어떤 정수 $m_k$을 선택해야 하는가? 여러 가지 선택이 있을 수 있지만 전구간에서 도함수의 제곱 평균을 최소로 만드는 $m_k$을 선택하도록 하자:
$$\underset {m_k} {\text {argmin} } ~ \frac {1}{L} \int_0^L | y'(x)|^2 dx $$
그런데 $$\frac{1}{L} \int_0^L | y'(x)|^2 dx = \Big(\frac {2\pi}{L} \Big)^2 ~\sum_{k=0}^{N-1} (k+m_k N)^2|Y_k|^2 $$ 이므로 각 $k$에 대해서 $|k +m_k N|^2$이 최소가 되도록 선택해야 한다. $k<N/2$일 때는 $m_k=0$을 선택하면 되고, $k>N/2$일 때는 $m_k=-1$을 선택하면 된다. $N=\text{even}$일 때는 $m_{N/2}=0 \text{ or }1$이 같은 결과를 주는데, 두 기여를 균등하게 되도록 선택한다:
$$ Y_{N/2} e^{\frac{2\pi i}{L} (N/2 + m_{N/2} N)x} \rightarrow \frac{1}{2} Y_{N/2}\left( e^{\frac{\pi i}{L} Nx} + e^{-\frac{\pi i}{L} Nx} \right) = Y_{N/2} \cos \Big( \frac{\pi}{L} Nx\Big) $$ 따라서 $m_k$ 선택은 정리하면 다음과 같다:
$$ m_k= \left\{\begin{array}{ll} 0, & 0\le k < N/2 \\ -1, & N/2 <k <N \end{array} \right.$$
따라서 주어진 샘플링 데이터 $\{y_n\}$의 DFT 계수 $\{ Y_n\}$을 구하면 샘플링 데이터를 보간하는 함수(Fourier interpolation)는
$$ y(x) = Y_0 + \sum_{0 <k <N/2} \left( Y_k e^{+\frac {2\pi i}{L} k x} + Y_{N-k} e^{ -\frac {2\pi i}{L} k x} \right) + Y_{N/2} \cos\Big( \frac{\pi}{L}Nx\Big) $$
로 얻어진다(마지막 항은 $N$이 짝수일 때만 더해짐). 이 보간함수의 도함수는
$$y'(x)= \sum_{0 <k <N/2} \frac {2\pi i }{L} k \Big(Y_k e^{+\frac {2\pi i}{L} k x} - Y_{N-k} e^{ -\frac{2\pi i}{L} k x} \Big) - \frac{\pi}{L} N Y_{N/2} \sin \Big( \frac{\pi}{L} Nx\Big)$$
로 계산된다(마지막 항은 $N$이 짝수일 때만 더해짐). 샘플링 위치에서 미분계수는
$$ y'_n = y'(nL/N) = \sum_{0<k<N/2} \frac{2\pi i}{L} k \Big( Y_k e^{+\frac{2\pi i}{N} kn} - Y_{N-k} e^{ -\frac{2\pi i}{N} kn} \Big) = \sum_{k=0}^{N-1} Y'_k e^{+\frac{2\pi i}{N} kn }$$
처럼 쓸 수 있는데 이 식은 $Y'_k$의 IDFT 형태이다.
따라서 DFT를 이용해서 미분계수를 구하는 전체적인 알고리즘은
1. DFT로 $Y_k$을 구하고,
2. $Y'_k$는 $$ Y_k' = \left\{ \begin {array}{ll} \frac {2\pi i}{L} k \times Y_k & 0 \le k <N/2\\ \frac{2\pi i}{L} (k-N) \times Y_k & k > N/2 \\ 0 & k=N/2 ~\text {when } N=\text {even} \end {array} \right. $$
3. $Y'_k$에 IDFT를 해서 $y'_n$을 구한다.
위에서 구한 보간 함수를 이용하면 이차 미분 값도 DFT를 이용해서 구할 수 있다.
아래 그림은 6개의 주파수 $f_k = 2,5,9,11,21,29\text{Hz}$가 섞인 신호 함수의 미분을 FFT를 이용해서 구한 결과를 정확한 도함수(green)와 함께 보여준다.
$$f(x) = \sum_{k=0}^{5} (k+1) \cos (2\pi f_k x)$$
$0\le x < 1=L$ 구간에서 일정한 간격으로 샘플링한 $N=64$개의 데이터를 사용하였다.
void fftgrad1d() {
const double TWOPI = 8.0 * atan(1.0);
const int samples = 64;
double re[samples], im[samples];
const int nfreqs = 6;
const double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
for ( int i = 0; i < samples; i++ ) {
double ti = double(i) / samples;
double signal = 0;
for ( int k = 0; k < nfreqs; k++ )
signal += (k + 1) * cos ( TWOPI * freq[k] * ti);
re[i] = signal;
im[i] = 0.;
}
// int fft1d(int dir, int nn, double x[], double y[]);
fft1d(1, samples, re, im);
const double norm = TWOPI / 1.0;
for (int i = 0; i < samples; i++) {
double rv = re[i];
double iv = im[i];
double fac = norm * double (i > samples / 2 ? i - samples: i);
if (i == samples / 2) fac = 0;
re[i] = - fac * iv;
im[i] = fac * rv;
}
//invers FFT: re[] stores derivatives.
fft1d(-1, samples, re, im);
}이차원에서 DFT를 이용해서 gradient을 구하는 경우는 복소함수 $f = Re(x, y) + i Im(x, y)$의 gradient가
$$\Big(\frac{\partial}{ \partial x} +i \frac{\partial}{\partial y} \Big) (Re + i Im)= i\Big( \frac{\partial Re}{\partial y} + \frac{\partial Im}{\partial x} \Big) - \Big(- \frac{\partial Re}{\partial x} + \frac{\partial Im}{\partial y}\Big) $$
임을 이용하면 된다.
영상(실수 함수)의 IDFT 식에 적용하면 미분 연산자가 우변에 전체적으로 $i$를 곱하는 효과를 주게 되고, 우변식의 실수부와 허수가 각각 영상의 $x$-방향/ $y$-방향 미분에 해당한다.
void FFTGrad(CRaster& raster, CRaster& out) {
const double TWOPI = 8.0 * atan(1.0);
CSize sz = raster.GetSize();
int nx = sz.cx, ny = sz.cy;
int stride = nx << 1;
std::vector<double> data(stride * sz.cy);
for (int y = 0; y < ny; y++) {
BYTE *p = (BYTE *)raster.GetLinePtr(y);
double *re = &data[y * stride];
double *im = re + nx;
for (int x = 0; x < nx; x++) {
re[x] = *p++;
im[x] = 0;
}
}
// forward FFT(1): https://kipl.tistory.com/22
fft2d(&data[0], nx, ny, 1);
double normx = TWOPI / double(nx);
double normy = TWOPI / double(ny);
int hnx = nx >> 1;
int hny = ny >> 1;
for (int y = 0; y < ny; y++) {
double *re = &data[y * stride];
double *im = re + nx;
for (int x = 0; x < nx; x++) {
double rv = re[x];
double iv = im[x];
double xx = normx * double(x > hnx ? x - nx: x);
double yy = normy * double(y > hny ? y - ny: y);
if (x == hnx) xx = 0; // N = even;
if (y == hny) yy = 0;
re[x] = yy * rv + xx * iv;
im[x] = -xx * rv + yy * iv;
}
}
/* inverse FFT of {Y'_k} */
fft2d(&data[0], nx, ny, -1);
std::vector<double> mag(nx * ny);
double magmax = 0, magmin = 1.e+20;
for (int y = 0, k = 0; y < ny; y++) {
double *re = &data[y * stride];
double *im = re + nx;
for (int x = 0; x < nx; x++) {
double m = hypot(re[x], im[x]);
if (m > magmax) magmax = m;
if (m < magmin) magmin = m;
mag[k++] = m;
}
}
// stretch the magnitude;
for (int y = 0, k = 0; y < ny; y++)
for (int x = 0; x < nx; x++) {
double a = 255 * (mag[k++] - magmin) / (magmax - magmin);
out.SetPixel(x, y, 255-int(a));
}
}'Image Recognition > Fundamental' 카테고리의 다른 글
| Valley emphasis Otsu threshold (0) | 2022.02.23 |
|---|---|
| Image rotation by FFT (0) | 2022.02.18 |
| SVD Fitting (0) | 2022.02.07 |
| Color Histogram Equalization (4) | 2022.02.07 |
| Least Squares Fitting of Ellipses (0) | 2022.01.27 |
