fft.cpp

数字信号的分析处理应用中常用的程序：快速傅立叶变换的C＋＋实现

/*Parameter:
f:        point a source sequance f(x)
F:        point a target sequance F(x) by FFT procedure
Function:
the procedure through FFT transform f(x) to F(x)
*/
void CJPEG::FFT(const complex<double> *f, complex<double> *F, int N)
...{
    
    //将时域复制至频域
    memcpy(F,f,sizeof(complex<double>) * N);

    //需要多少级蝶形运算
    int r = log(N +1) / log(2);

    //进行倒位序运算
    BitReOrder(F,r);

    //计算Wr因子, 加权系数 
    complex<double> *W = new complex<double>[N/2];
    double angle;    
    for(int i = 0; i < N / 2; i++) 
    ...{ 
        angle = - i * PI * 2 / N; 
        W[i] = complex<double> (cos(angle), sin(angle)); 
        //W[i] = complex<double>(0, exp(-i*2*PI/N));
    } 

    //采用蝶形算法进行快速傅立叶变换
    int DFTn;//第DFTn级
    int k;//分级有多少个蝶形运算
    int d;//蝶形运算的偏移
    int p;//index
    
    //complex<double> *X,*X1,*X2;
    //X1 = new complex<double>[N]; 
    //X2 = new complex<double>[N]; 
    complex<double> X1,X2;
    for( DFTn = 0; DFTn < r; DFTn++)...{//第几级蝶形运算
        for( k= 0; k < (1 << (r - DFTn - 1)); ++k)...{//当前列所有的蝶形进行运算
            p = 2 * k *  (1 << DFTn);
            for ( d = 0; d < (1 << DFTn); ++d)...{//相邻蝶形运算
                //p = k * (1 << (k + 1)) + d;
                X1 = F[p+d];
                X2 = F[p+d+(1 << DFTn)];
                F[p+d] = X1 + X2 * W[d*(1<<(r - DFTn - 1))];
                F[p+d+(1 << DFTn)] = X1 - X2 * W[d*(1<<(r - DFTn - 1))];
            }
        }
    }
    //BitReOrder(F,r);
    for(i = 0; i < N; ++i)
        F[i] /= N;
    delete[] W;

    return;
}



void CJPEG::IFFT(const complex<double> *F, complex<double> *f, int N)
...{
    int i;
    complex<double>*T = new complex<double>[N];
    memcpy(T,F,sizeof(complex<double>) * N);
    for(i = 0; i < N; ++i)...{
        //取共轭
        T[i] = complex<double>(T[i].real(), -T[i].imag());        
    }
    CJPEG::FFT(T,f,N);
    for(i = 0; i < N; ++i)...{//取共轭，并乘以N
        f[i] = complex<double>(f[i].real(), -f[i].imag());    
        f[i] *= complex<double>(N,0) ;
    }
    delete[] T;
}

int CJPEG::ReadFraqData(const char *lpszFileName, int N, double *pData)
...{
    return 0;
}

int CJPEG::ReadTimeData(const char *lpszFileName, int N, double *pData)
...{
    if(pData == NULL) return -1;
    CTextFile ctf(lpszFileName);
    char sLine[64];
    double dData;
    int i = 0;
    while(!ctf.Eof() && i < N)...{
        ctf.GetLine(sLine, 64);
        dData = atof(sLine);
        pData[i++] = dData;
    }
    return i;
}

/**//////////////////////////////////////////// 
// Function name : BitReverse 
// Description : 二进制倒序操作 
// Return type : int 
// Argument : int src 待倒读的数 
// Argument : int size 二进制位数 
int CJPEG::BitReverse(int src, int size) 
...{ 
    int tmp = src; 
    int des = 0; 
    for (int i=size-1; i>=0; i--) 
    ...{ 
        des = ((tmp & 0x1) << i) | des; 
        tmp = tmp >> 1; 
    } 
    return des; 
} 

/**//*
Parameter:
sequ:    point a sequance that f(x)
N:        f(x), the x maxinum
Function:
the procedure change the sequance by this:
    // 101 -> 101
    // 110 -> 011
    // 1010 -> 0101
f(5) -> f(5)
f(6) -> f(3)
f(10) -> f(5)
and so on.
*/
void CJPEG::BitReOrder(complex<double> *sequ, int r)
...{
    complex<double> dTemp;
    int x;
    for(int i = 0; i <(1 << r); ++i)...{
        x = BitReverse(i,r);
        if (x > i)...{
            dTemp = sequ[i];
            sequ[i] = sequ[x];
            sequ[x] = dTemp;
        }
    }    
}