사용자 삽입 이미지

사용자 삽입 이미지

/* 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) {
    int  i,i1,j,k,i2,l,l1,l2;
    double c1,c2,tx,ty,t1,t2,u1,u2,z;
    /* Calculate the number of points */
    int m = 0;
    while(nn > 1) {
        nn >>= 1;  m++;
    }
    nn = 1 << m;

    /* Do the bit reversal */
    i2 = nn >> 1;
    j = 0;
    for (i=0;i<nn-1;i++) {
        if (i < j) {
            tx = x[i];    ty = y[i];
            x[i] = x[j]; y[i] = y[j];
            x[j] = tx;    y[j] = ty;
        }
        k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
   
    /* Compute the FFT */
    c1 = -1.0;
    c2 = 0.0;
    l2 = 1;
    for (l=0;l<m;l++) {
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0;
        u2 = 0.0;
        for (j=0;j<l1;j++) {
            for (i=j;i<nn;i+=l2) {
                i1 = i + l1;
                t1 = u1 * x[i1] - u2 * y[i1];
                t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            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 i, j, index1, index2;
   
    /* Transform by ROWS for forward transform */
    if (isign == 1)
    {
        index1 = 0;
        for (i = 0; i < mm; i++)
        {
            fft1d(isign, nn, &data[index1], &data[index1+nn]);
            index1 += (nn << 1);
        }
    }
   
    /* Transform by COLUMNS */
    for (j = 0; j < nn; j++)
    {
        /* Copy pixels into temp array */
        index1 = j ;
        index2 = 0;
        for (i = 0; i < mm; i++)
        {
            copy[index2] = data[index1];
            copy[index2+mm] = data[index1 + nn];
            index2++;
            index1 += (nn << 1);
        }
       
        /* 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 += (nn << 1);
        }
    }
   
    /* Transform by ROWS for inverse transform */
    if (isign == -1)
    {
        index1 = 0;
        for (i = 0; i < mm; i++)
        {
            fft1d(isign, nn, &data[index1], &data[index1+nn]);
            index1 += (nn << 1);
        }
    }
}
void butterworth( double * data, int Xdim, int Ydim,
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost)
{
    int x, y, x1, y1, halfx, halfy;
    double Filter;
    double CutOff2 = CutOff * CutOff;
   
    /* Prepare to filter */
    halfx = Xdim / 2;
    halfy = Ydim / 2;
    double *rdata = &data[0];
    double *idata = &data[Xdim];
    /* Loop over Y axis */
    for (y = 0; y < Ydim; y++) {
        if (y < halfy)  y1 = y;
        else             y1 = y - Ydim;
       
        /* Loop over X axis */
        for (x = 0; x < Xdim; x++) {
            if (x < halfx)  x1 = x;
            else             x1 = x - Xdim;
           
            /* Calculate value of Butterworth 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 += (Xdim << 1) ; //go to next-line;
        idata += (Xdim << 1) ;
    }
}
Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)

사용자 삽입 이미지

PowerSpectrum 변화
사용자 삽입 이미지

'Image Recognition > Fundamental' 카테고리의 다른 글

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk
TAG ,