Geometry & Recognition

• Forward transformation:
  $$X_{k} =  \sum_{n = 0} ^{N-1} x_n e^{ -2i  \pi k n /N }, ~\quad k=0,1,..., N-1;$$
• Backward transformation:
  $$x_n = \frac{1}{N}\sum_{k = 0} ^{N-1} X_k  e^{ + 2i  \pi k n /N }, ~\quad n=0,1,..., N-1;$$
• Butterworth Filter:
  $$H(x, y) = \frac{1}{1+\left( \frac{x^2 +y^2}{H_c^2 }\right)^N} ,\quad N=\text{order},~H_c =\text{cut-off};$$
• 주어진 수가 $2^n$인가?(kipl.tistory.com/84). 입력 데이터 개수가 2의 거듭제곱이 아니면 부족한 부분은 0으로 채워서 $2^n \ge \text{데이터 개수}$ 크기가 되도록 만든다. 

Cooley–Tukey FFT algorithm:  

/* 1차원 FFT:: dir=1 for forward, -1 for backward;
** x = real-component, and y = imaginary componet.
** nn = length of data: 2의 거듭제곱이어야 한다.
*/
int fft1d(int dir, int nn, double x[], double y[]) {
    /* Calculate the number of points */
    int m = 0;
    while (nn > 1) {
        nn >>= 1;  m++;
    }
    nn = 1 << m;
    
    /* Do the bit reversal */
    int i2 = nn >> 1;
    int j = 0;
    for (int i = 0; i < nn - 1; i++) {
        if (i < j) { // swap(x[i], x[j]), swap(y[i], y[j]);
            double tx = x[i], ty = y[i];
            x[i] = x[j];  y[i] = y[j];
            x[j] = tx;    y[j] = ty;
        }
        int k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
    
    /* Compute the FFT */
    double c1 = -1.0;
    double c2 = 0.0;
    int l2 = 1;
    for (l = 0; l < m; l++) {
        int l1 = l2;
        l2 <<= 1;
        double u1 = 1.0;
        double u2 = 0.0;
        for (int j = 0; j < l1; j++) {
            for (int i = j; i < nn; i += l2) {
                int i1 = i + l1;
                double t1 = u1 * x[i1] - u2 * y[i1];
                double t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            double z = u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1) c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }
    
    /* Scaling for forward transform */
    if (dir == 1) {
        for (i = 0; i < nn; i++) {
            x[i] /= double(nn);
            y[i] /= double(nn);
        }
    }
    return 1;
};

/* 2차원 FFT. data = (2 * nn) x mm 크기의 버퍼;
** 각 행은 nn개의 실수 성분 다음에 nn개의 허수 성분이 온다.
** nn, mm은 각각 행과 열은 나타내며, 2의 거듭제곱이어야 한다.
** isign = 1 for forward, -1 for backward. 
** copy = buffer of length 2*mm ;
*/
void fft2d(double* data, int nn, int mm, int isign, double* copy) { 
    int stride = nn << 1; //row_stride;
    /* Transform by ROWS for forward transform */
    if (isign == 1) {
        int index1 = 0;
        for (int i = 0; i < mm; i++) {
            // real = &data[index1], imag = &data[index1 + nn];
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
    
    /* Transform by COLUMNS */
    for (int j = 0; j < nn; j++) {
        /* Copy pixels into temp array */
        int index1 = j ;
        int index2 = 0;
        for (i = 0; i < mm; i++) {
            copy[index2] = data[index1];
            copy[index2 + mm] = data[index1 + nn];
            index2++;
            index1 += stride;
        }
        
        /* Perform transform */
        fft1d(isign, mm, &copy[0], &copy[mm]);
        
        /* Copy pixels back into data array */
        index1 = j;
        index2 = 0;
        for (i = 0; i < mm; i++) {
            data[index1] = copy[index2];
            data[index1 + mm] = copy[index2 + mm];
            index2++;
            index1 += stride;
        }
    }
    
    /* Transform by ROWS for inverse transform */
    if (isign == -1) {
        int index1 = 0;
        for (i = 0; i < mm; i++) {
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
}

void butterworth( double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost);

// H * F;
void butterworth(double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost) 
{
    double CutOff2 = CutOff * CutOff;
    /* Prepare to filter */
    int stride = Xdim << 1;  // width of a sigle row(real + imag);    
    int halfx = Xdim >> 1;
    int halfy = Ydim >> 1;
    double *rdata = &data[0];   //real part;
    double *idata = &data[Xdim];//imag part;
    /* Loop over Y axis */
    for (int y = 0; y < Ydim; y++) {
        int y1 = (y < halfy) ? y: y - Ydim;
        /* Loop over X axis */
        for (int x = 0; x < Xdim; x++) {
            int x1 = (x < halfx) ? x: x - Xdim;
            /* Calculate value of Butterworth filter */
            double filter;
            if (LowPass)
                Filter = (1 / (1 + pow((x1 * x1 + y1 * y1) / CutOff2, Power)));
            else if ((x1 != 0) || (y1 != 0))
                Filter = (1 / (1 + pow( CutOff2 / (x1 * x1 + y1 * y1), Power)));
            else
                Filter = 0.0;
            if (Homomorph)
                Filter = Boost + (1 - Boost) * Filter;
            /* Do pointwise multiplication */
            rdata[x] *= Filter;
            idata[x] *= Filter;
        };
        rdata += stride; //go to next-line;
        idata += stride; 
    }
};

Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)

PowerSpectrum 변화