Geometry & Recognition

주어진 함수 $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}nk},$$

로 정의된다. 역변환(IDFT)은 $Y_k$에서 $y_n$을 얻는 과정으로 다음 식으로 계산된다:

$$y_n=\sum _{k=0}^{N-1}{Y}_k{e}^{+\frac{2\pi i}{N}kn}.$$

이산적인 data $y_n$만을 가지고 $x_n=nL/N$에서 미분계수를 구하려면 어떻게 해야 할까? 샘플링 점 사이 구간에서 값을 추정할 수 있도록 보간 함수를 만들어 사용해야 한다. 간단하게는 IDFT 식에서 $n$번째 샘플링 위치인 $x_n=n(L/N)$을 실수 $x$로 대체하여 연속함수를 얻는 방법을 생각할 수 있다. 그리고 IDFT 식에서 $Y_k{e}^{+\frac{2\pi i}{N}kn}$ 대신에 추가적으로 위상 요소를 곱해서 $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^{\prime }(x){|}^{2}dx$$

그런데 $$\frac{1}{L}{\int }_0^L|y^{\prime }(x){|}^{2}dx=(\frac{2\pi }{L}{)}^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}\to \frac{1}{2}{Y}_{N/2}(e^{\frac{\pi i}{L}Nx}+e^{-\frac{\pi i}{L}Nx})=Y_{N/2}\cos (\frac{\pi }{L}Nx)$$ 따라서 $m_k$ 선택은 정리하면 다음과 같다:

$$m_k=\{\begin{array}{ll}0,& 0\le k<N/2\\ -1,& N/2<k<N\end{array}$$

따라서 주어진 샘플링 데이터 $\{y_n\}$의 DFT 계수 $\{Y_n\}$을 구하면 샘플링 데이터를 보간하는 함수(Fourier interpolation)는

$$y(x)=Y_0+\sum _{0<k<N/2}(Y_k{e}^{+\frac{2\pi i}{L}kx}+Y_{N-k}{e}^{-\frac{2\pi i}{L}kx})+Y_{N/2}\cos (\frac{\pi }{L}Nx)$$

로 얻어진다(마지막 항은 $N$이 짝수일 때만 더해짐). 이 보간함수의 도함수는 

$$y^{\prime }(x)=\sum _{0<k<N/2}\frac{2\pi i}{L}k(Y_k{e}^{+\frac{2\pi i}{L}kx}-Y_{N-k}{e}^{-\frac{2\pi i}{L}kx})-\frac{\pi }{L}N{Y}_{N/2}\sin (\frac{\pi }{L}Nx)$$

로 계산된다(마지막 항은 $N$이 짝수일 때만 더해짐). 샘플링 위치에서 미분계수는

$$y_n^{\prime }=y^{\prime }(nL/N)=\sum _{0<k<N/2}\frac{2\pi i}{L}k(Y_k{e}^{+\frac{2\pi i}{N}kn}-Y_{N-k}{e}^{-\frac{2\pi i}{N}kn})=\sum _{k=0}^{N-1}{Y}_k^{\prime }{e}^{+\frac{2\pi i}{N}kn}$$

처럼 쓸 수 있는데 이 식은 $Y_k^{\prime }$의 IDFT 형태이다. 

따라서 DFT를 이용해서 미분계수를 구하는 전체적인 알고리즘은

1. DFT로 $Y_k$을 구하고,

2. $Y_k^{\prime }$는 $$Y_k^{\prime }=\{\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}$$

3. $Y_k^{\prime }$에 IDFT를 해서 $y_n^{\prime }$을 구한다.

위에서 구한 보간 함수를 이용하면 이차 미분 값도 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)+iIm(x,y)$의 gradient가

$$(\frac{\partial }{\partial x}+i\frac{\partial }{\partial y})(Re+iIm)=i(\frac{\partial Re}{\partial y}+\frac{\partial Im}{\partial x})-(-\frac{\partial Re}{\partial x}+\frac{\partial Im}{\partial y})$$

임을 이용하면 된다.

영상(실수 함수)의 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));
        }
}