FFT(Fast Fourier Transform)

Fast Fourier Transform 코드입니다.

void FourierTransform::revertBit(float* real, float* img, int size, float* realDst, float* imgDst, bool bHori)
{
    int reverseIdx = 0;
    int n = (int)log2((double)size);
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (i & (1 << j))
                reverseIdx |= (size >> (j + 1));
        }

        if (!bHori) {
            realDst[reverseIdx * size] = real[i * size];
             imgDst[reverseIdx * size] =  img[i * size];
        }
        else {
            realDst[reverseIdx] = real[i];
             imgDst[reverseIdx] = img[i];
        }

        reverseIdx = 0;
    }
}

1D FFT 코드입니다.

void FourierTransform::FFT1D(float* pSrcReal, float* pSrcImg, int size, float* pDstReal, float* pDstImg, bool bHori, bool bInverse)
{
    /* only support 2^n */
    revertBit(pSrcReal, pSrcImg, size, pDstReal, pDstImg, bHori);

    /* 1D FFT */
    int n = (int)log2((double)size);
    int bias = 1;
    int step = size >> 1;
    int idx = 0;
    float realOdd = 0.f, imgOdd = 0.f, realEven = 0.f, imgEven = 0.f;
    float term0 = 0.f, term1 = 0.f, term2 = 0.f;
    uint evenIdx = 0, oddIdx = 0;
    int realSign = +1, imgSign = -1;
    
    if (bInverse)
        realSign = -1, imgSign = +1;

    for (int s = 0; s < n; s++)
    {
        for (int i = 0; i < step; i++)
        {
            for (int k = 0; k < bias; k++)
            {
                evenIdx = idx + k;
                oddIdx = idx + k + bias;

                if (!bHori) {
                    evenIdx = (idx + k) * size;
                    oddIdx = (idx + (k + bias)) * size;
                }

                /* only set by index calc */
                realEven = pDstReal[evenIdx];
                imgEven = pDstImg[evenIdx];

                realOdd = pDstReal[oddIdx];
                imgOdd = pDstImg[oddIdx];

                // j = n, k = idx + j 
                term0 = (2.f * PI * (float)k) / (float)(bias << 1);
                term1 = (realOdd * cosf(term0)) + realSign * (imgOdd  * sinf(term0)); // real
                term2 =  (imgOdd * cosf(term0)) +  imgSign * (realOdd * sinf(term0)); // img

                pDstReal[evenIdx] = realEven + term1;
                pDstImg[evenIdx] = imgEven + term2;

                pDstReal[oddIdx] = realEven - term1;
                pDstImg[oddIdx] = imgEven - term2;
            }

            idx += bias << 1;
        }

        idx = 0;
        bias <<= 1;
        step >>= 1;
    }
}

2D FFT 코드입니다.

void FourierTransform::FFT(Mat* pSrcReal, Mat* pSrcImg, Mat* pDstReal, Mat* pDstImg, bool bInverse)
{
    float* srcReal = nullptr;
    float* srcImg  = nullptr;
    float* dstReal = nullptr;
    float* dstImg  = nullptr;

    for (int u = 0; u < pSrcReal->rows; u++) 
    {
        srcReal = pSrcReal->ptr<float>(u);
        srcImg  = pSrcImg ->ptr<float>(u);
        dstReal = _pTemp1ChnImage[2]->ptr<float>(u);
        dstImg  = _pTemp1ChnImage[3]->ptr<float>(u);

        FFT1D(srcReal, srcImg, pSrcReal->cols, dstReal, dstImg, true, bInverse);
    }

    srcReal = _pTemp1ChnImage[2]->ptr<float>(0);
    srcImg  = _pTemp1ChnImage[3]->ptr<float>(0);
    dstReal = pDstReal->ptr<float>(0);
    dstImg  = pDstImg->ptr<float>(0);

    for (int v = 0; v < pSrcReal->cols; v++) 
    {
        FFT1D(srcReal, srcImg, pSrcReal->rows, dstReal, dstImg, false, bInverse);
        srcReal++, srcImg++, dstReal++, dstImg++;
    }
}

DFT와 결과는 동일합니다.

왼쪽: box aperture / 오른쪽 DFT 결과 가시화(shiftImage & spectrum 함수 적용)

왼쪽: 육각 aperture / 오른쪽 DFT 결과 가시화(shiftImage & spectrum 함수 적용)

맨왼쪽: 입력이미지 / 중앙: DFT 결과 / 맨오른쪽: DFT 결과를 IDFT 한 결과