Fourier Interpolation

$0\le x<L$ 구간에서 일정한 간격($\mathrm{\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}(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)$$

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

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