목적

화상 처리 100개 노크에 도전중입니다. Q.27은 바이큐빅 보간입니다만, 아무래도 샘플링 정리의 응용인 것 같습니다. 여기 편을 제대로 이해하기 위해 샘플링 정리를 복습합니다.

샘플링 정리

관측 대상은 연속 실함수 $f(x)$입니다. 이산 샘플링은 com 함수로 표현됩니다.

s(x) = \sum_n \delta(x-nL)

$L$는 샘플링 간격입니다. 샘플링 후 신호는

\tilde f(x) = f(x)\sum_n \delta(x-nL)

그리고 나타납니다.

이것을 푸리에 변환합니다. 컨벌루션 정리보다

\mathcal F[\tilde f](\nu) = \mathcal F[f] \otimes \mathcal F[\sum_n\delta(x-nL)]

됩니다. 콤 함수의 푸리에 변환은

\mathcal F[\sum_n\delta(x-nL)] = \frac{1}{L}\sum_n\delta(\nu- n/L)

입니다. 이제부터

\mathcal F[\tilde f](\nu) = \frac{1}{L}\sum_n F(\nu- n/L)

됩니다. 원래의 연속 신호 $f$의 푸리에 변환 $F(\nu)$를 $1/L$ 주기로 늘어놓은 것이 되는 것입니다. 샘플링 간격을 줄이면 이 주기가 길어집니다. 즉, 충분히 작은 샘플링 간격을 취하면, 인접한 $F(\nu)$끼리가 겹치지 않게 됩니다. 실제 함수의 푸리에 변환이 짝수 함수가된다는 점에 유의하면,

\nu_{\rm max} < \frac{1}{2L}

를 채우도록 샘플링하면 됩니다. 이것이 샘플링 정리입니다. 신호가 포함하는 최대 주파수의 2배 이상의 샘플링 주파수가 필요하다는 것입니다. 2배인 것은 실함수의 푸리에 변환이 짝수함인 것에 유래하는 것입니다.

그런데, 이 조건을 만족하고 있는 경우, 인접하는 $F$끼리는 겹치지 않기 때문에, 주파수 공간상에서 폭 $1/L$의 구형 함수($\times L$)를 걸어 주면 원래의 $ F(\nu)$를 얻을 수 있습니다. 실제 공간에서 말하면,

\mathcal F[L\times {\rm rect}_{1/L}] = L\times \frac{\sin(\pi x/L)}{\pi x}

보다, sinc 함수를 이산 샘플링 된 신호 $\tilde f$에 컨벌루어 하면 원래의 연속 신호 $f$를 복원할 수 있게 됩니다. 이 sinc 함수의 주기는 $2L$이므로, 컨볼루션해도, 샘플링한 점에 있어서 기여하는 것은 그 점만이며 값이 변하지 않는 것을 알 수 있습니다. 잘 되었습니다.

전기 신호 등이라면, 샘플링된 이산 펄스열 $\tilde f$에 이상적인 로우 패스 필터를 걸어 주면 연속 신호 $f$를 생성할 수 있는 것입니다.
화상 처리와 같이 계산기상에서 디지털적인 계산을 하는 경우는, 실공간에서의 컨벌루션에 의해 임의의 위치 $x$에서의 아날로그 신호치를 얻을 수 있습니다.

이산 신호에서 연속 신호를 복원해 봅니다.

원래 신호를 $f(x) =\exp(-\frac{x^2}{2\sigma^2})$로 설정합니다. 샘플링 간격은 $L=1$로 합니다.

#include <iostream>
#include <vector>
const double pi = 3.14159265358979323846;

int main()
{
    int N = 1024;
    auto gauss = [&](double x)
    {
        double m = N / 2;
        double s = 100;
        return exp(-(x - m) * (x - m)/2./s/s);
    };
    std::vector<double> data;
    for (int i = 0; i < N; i++)
    {
        data.push_back(gauss(i));
    }

    auto f = [&](double x)
    {
        double ret = 0;
        for (int i = 0; i < N; i++)
        {
            double d = x - i;
            if (std::abs(d) > 1e-5)ret += data[i] * sin(pi * d) / pi / d;
            else ret += data[i];
        }
        return ret;
    };

    for (double x = N/2-10; x < N/2+10; x += 0.2)
    {
        std::cout << x << " " << f(x) << " " << gauss(x) << std::endl;
    }


    return 0;
}

실선이 참 값, 점이 복원한 값입니다. 정수점만 측정하고 있는데 제대로 복원할 수 있습니다.

바이큐빅 보간과의 관계

바이큐빅 보간은 sinc 함수를 3차 근사한 식을 이용하는 것 같습니다. 보간이라기보다는 "회복"이라고 말하는 것이 확고합니다. 물론, 단일 화소내에서의 급격한 변화가 없는 전제가 아니면 올바르게 회복할 수 있는 것은 아닙니다만.
자세한 계산은 향후 쓸 예정인 【이미지 처리 100개 노크에 도전】기사중에서 해 볼 생각입니다.