Savitzky-Golay smoothing

   스무딩은 자료를 매끄럽게 하는 것이라 할 수 있다. 이러한 작업은 데이터를 다루는 것에 있어 거의 필수적이다. 왜냐 하면, 관측되는 거의 모든 데이터는 여러 요인으로 인해 '오차'를 포함하고 있기 때문이다. 그래서 그와 같은 오차를 없애기 위해 여러 통계적 기법을 사용하는데, 스무딩은 통계적이진 않지만 오차를 없애고 원래의 데이터를 추정해 내기 위한 단순한 기법이면서도 효과적인 방법이다.

   가장 쉽게 생각할 수 있는 방법은 이동평균(moving average)로, 주어진 데이터 전/후의 일정 개수의 데이터의 평균을 그 데이터의 값으로 추정하는 방법이다. 만약 주어진 데이터에서 멀어지는 점일수록 중요도가 떨어진다면 중요도를 낮추어 주면 된다. 즉, 이동평균은 주어진 데이터와의 거리에 상관없이 모두 동일한 가중치(1)를 두는 것으로 생각할 수 있는 것이다. 이와 같은 개념을 좀 더 일반화시킨 다음 다항식으로까지 전개시켜서 smoothing 하는 알고리즘이 Saivtzky-Golay 알고리즘이다.

   이 글은 Smoothing and Differentiation of Data by Simplified Least Squares Procedures, Abraham Savitzky and Marcel J. E. Golay, 1994, Analytical Chemistry 논문에 근거하여 작성한다. 

이동평균은 y=1 인 함수와 원 함수와의 convolution 이라고 할 수 있다. 그 이외에도 여러 convolution 함수를 사용할 수 있다. 위 그림 참조. Savitzky-Golay 알고리즘은 데이터와 양옆의 데이터들을 다항식으로 회귀시켰을 때 만들어지는 다항식으로 현재의 데이터 값을 설정하자는 것이다. 핵심은 각각의 데이터에 대해 매번 회귀시킨 다항식을 계산할 필요가 없이 특정 convolution 을 사용하면 다항식을 회귀시킨 것을 '정확히' 찾아낼 수 있다는 것이다^[각주:1]. 이와 같은 것을 논문에서는 다음과 같이 표현하고 있다.

즉, xi의 새로운 값은 양쪽 근방 2n+1 개의 점으로 다항식 회귀한 식으로부터 다시 추정해 내는 것이다. 즉, xi 를 포함하고 있는 양쪽 근방 2n+1 개의 점으로 다항식 회귀한 식 Si을 찾고, Si(xi)로 xi 를 대체하는 것이다. 문제는 Si를 찾아 내는 것인데, 각 xi에 대해서 Si를 매번 찾아내야 한다면 계산횟수가 많아질텐데, Si를 3차 다항식으로 식을 만든 후, 각 점에 대해서 다항연립방정식을 만들면, 놀랍게도, 각 계수가 정확히 계산이 된다! 이렇게 계산된 계수로부터 추정된 값, 즉 Si(xi)을 계산하는 것은 결국 convolution 형태로 나타나게 된다. 쉽게 말해, 이동평균을 구할 때는 앞/뒤 구간에 있는 값들에 1씩 곱해서 더한 후 평균을 구하는 것과 비슷하게, Savitzky-Golay smoothing 은 정해진 weight로 이동평균을 구하는 것과 비슷하게 하는 것이다. 더욱 놀라운 것은, 원래의 데이터를 smoothing 하는 것 뿐만이 아니라, 그 데이터를 미분한 그래프가 어떻게 생겼는지까지 이렇한 convolution 형태로 쉽게 계산할 수 있다는 것이다. 데이터가 급격하게 변하는 곳은 고차 다항식으로, 그렇지 않은 곳은 저차 다항식으로 smoothing 할 수도 있는 것이다. 이와 같은 테이블의 예는 다음과 같다. 

이와 같은 알고리즘을 C++ 코드로 나타내면 다음과 같다. 

bool Savitzky_Golay_smoothing(std::vector<double>* x_series, std::vector<double>* y_series, std::vector<double> *destX, std::vector<double>* destY)

{
        int A[] = { -2, 3, 6, 7, 6, 3, -2 };

        int n = 3;
        int k = 0;
        int point_number = (int)x_series->size();
        for(k = 0; k<n; k++){
                if(destX != NULL)     destX->push_back(x_series->operator [](k));
                if(destY != NULL)     destY->push_back(y_series->operator [](k));
        }
        for(k = n; k<point_number-n; k++){
                double x = x_series->operator[](k);
                int i = 0;
                double nominator = 0;
                double denominator = 0;
                for(i = -n; i <= n; i++){
                        nominator += (A[n+i]*y_series->operator[](k+i));
                        denominator += A[n+i];
                }
                double y = nominator / denominator;
                if(destX != NULL)             destX->push_back(x);
                if(destY != NULL)             destY->push_back(y);
        }
        for(k = point_number-n; k<point_number; k++){
                if(destX != NULL)     destX->push_back(x_series->operator [](k));
                if(destY != NULL)     destY->push_back(y_series->operator [](k));
        }

        return true;

  }

      

위 코드는 3차 다항식으로 regression 하는 것을 사용하고 있으며, 위 테이블에 NORM으로 되어 있는 normalizing factor 로 나누어 주고 있지 않다. 내가 사용한 코드에서는 특성상 NORM으로 나눌 필요가 없었기 때문이며, 원래는 나누어 주어야지 xi와 Si(xi)의 scale이 같아진다. 또한 이런 식으로 window를 사용하는 알고리즘들은 가장 앞/뒤의 window size /  2 크기만큼의 데이터가 문제가 되는데, 위 코드에서는 이 값을 그냥 원래의 값을 사용하고 있다. 이것도 나는 문제의 특성상 가장 앞/뒤 몇 개의 데이터는 필요 없어서 그냥 저렇게 한 것이므로 그렇지 않은 경우엔 적당한 수정이 필요하다. 

1. 원 논문을 읽어 보면, 'It is not approximate.' 란 구절이 나온다. 감동!! [본문으로]