Fourier Interpolation

$0\le x <L$ 구간에서 일정한 간격($\Delta = L/N$)으로 샘플링된 데이터 $\{y_k|k=0,...,N-1\}$이 주어졌을 때 이들의  DFT 계수 $\{Y_k | k=0, N-1\}$을 이용하면 보간함수를 만들 수 있다.

$$ 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) $$

마지막 항(Nyquist term)은 $N$이 짝수일 때만 더해진다. 구체적인 유도과정은 https://kipl.tistory.com/406에 있다. 샘플링 데이터가 실수인 경우는 $Y_{N-k} = (Y_k)^*$가 성립하므로 표현이 더 간단해진다.

FFT를 이용한 영상의 미분

 

$x/L\to x$으로 하면 이 보간함수의 domain은 $[0,1]$로 제한된다.

// [x=0:1]
// Fourier coefficient (re[], im[]);
double FourierInterp(double x, std::vector<double>& re, std::vector<double>& im) {
    const int N = re.size();
    double fval = re[0];
    for (int k = N >> 1; k-->1;) {
        double ang = TWOPI * k * x;
        double cc = cos(ang);
        double ss = sin(ang);
        fval += re[k] * cc - im[k] * ss;     // 입력데이터가 실수인 경우:
        fval += re[N-k] * cc + im[N-k] * ss; // Y[N-k] = (Y[k])*
    }
    if (!(N&1))           //k = N/2 for even N;
        fval += re[N>>1] * cos(PI * N * x);
    return fval;
}

$f(x) = \sin^2(x)\cos(x)$을 $[0,3\pi]$ 구간에서 일정한 간격으로 얻은 30개의 데이터를 이용해서 보간함수를 만든다.

double f(double x) {
    return sin(x) * sin(x) * cos(x);
}
void Finterpolate() {
    const int nsamples = 30;
    std::vector<double> sample(nsamples);
    const double L = 3 * PI, dx = L / nsamples;
    for (int k = 0; k < nsamples; k++) sample[k] = f(k * dx);

    std::vector<double> re(sample.size(), 0), im(sample.size(), 0);
    // perform DFT;
    DFT(sample, re, im);

    FILE *fp = fopen("finp.txt", "w");
    const double N = 200;
    const double dy = dx * sample.size() / N;
    for (int k = 0; k < N; k++) { 
        double x = k * dy;
        fprintf(fp, "%f %f\n", x, FourierInterp(x / L, re, im));
    }
    fclose(fp);
}