주어진 함수 $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$-방향 미분에 해당한다.

spectrum
Canny edge

 

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));
        }
}
728x90

'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
Posted by helloktk
,