[C++] 頻譜圖中 FFT快速傅里葉變換C++實現

在項目中，需要畫波形頻譜圖，因此進行查找，不是很懂相關知識，下列代碼主要是針對這篇文章。

http://blog.csdn.net/xcgspring/article/details/4749075

//快速傅里葉變換
/*
入口參數：
inv： =1,傅里葉變換; =-1,逆傅里葉變換
N：輸入的點數，為偶數，一般為2的冪次級，2,4,8,16...
k: 滿足N=2^k（k>0）,實質上k是N個采樣數據可以分解為偶次冪和奇次冪的次數
real[]: inv=1時，存放N個采樣數據的實部，inv=-1時，存放傅里葉變換的N個實部
imag[]: inv=1時，存放N個采樣數據的虛部，inv=-1時，存放傅里葉變換的N個虛部

出口參數：
real[]: inv=1時，返回傅里葉變換的實部，inv=-1時，返回逆傅里葉變換的實部
imag[]: inv=1時，返回傅里葉變換的虛部，inv=-1時，返回逆傅里葉變換的虛部
*/
void  FFT::dealFFT(double real[], double imag[], double dSp[], int N, int k, int inv)
{
    int i, j, k1, k2, m, step, factor_step;
    double temp_real, temp_imag, factor_real, factor_imag;
    if (inv != 1 && inv != -1)
        return;

    //double *real = new double[N];
    //double *imag = new double[N];
    //倒序
    j = 0;
    for (i = 0; i < N; i++)
    {
        if (j>i)
        {
            temp_real = real[j];
            real[j] = real[i];
            real[i] = temp_real;

            temp_imag = imag[j];
            imag[j] = imag[i];
            imag[i] = temp_imag;
        }
        m = N / 2;
        while (j >= m&&m != 0)
        {
            j -= m;
            m >>= 1;
        }
        j += m;
    }

    //蝶形運算
    for (i = 0; i < k; i++)
    {
        step = 1 << (i + 1);
        factor_step = N >> (i + 1); //旋轉因數變化速度

        //初始化旋轉因子
        factor_real = 1.0;
        factor_imag = 0.0;

        for (j = 0; j < step / 2; j++)
        {
            for (k1 = j; k1 < N; k1 += step)
            {
                k2 = k1 + step / 2; //蝶形運算的兩個輸入

        /*        temp_real = real[k1] + real[k2] * factor_real - imag[k2] * factor_imag;
                temp_imag = imag[k1] + real[k2] * factor_imag + imag[k2] * factor_real;
                real[k2] = real[k1] - (real[k2] * factor_real - imag[k2] * factor_imag);
                imag[k2] = imag[k1] - (real[k2] * factor_imag + imag[k2] * factor_real);
                real[k1] = temp_real;
                imag[k1] = temp_imag;*/
                //上面方法雖然直白，但效率太低，稍微改變結構如下：
                temp_real = real[k2] * factor_real - imag[k2] * factor_imag;
                temp_imag = real[k2] * factor_imag + imag[k2] * factor_real;
                real[k2] = real[k1] - temp_real;
                imag[k2] = imag[k1] - temp_imag;
                real[k1] = real[k1] + temp_real;
                imag[k1] = imag[k1] + temp_imag;
            }

            factor_real = inv*cos(-2 * PI*(j + 1)*factor_step / N);
            factor_imag = inv*sin(-2 * PI*(j + 1)*factor_step / N);
        }
    }
     
    if (inv == -1)
    {
        for (i = 0; i <= N - 1; i++)
        {
            real[i] = real[i] / N;
            imag[i] = imag[i] / N;
        }
    }
    for (i = 0; i<N;i++)
    {
        dSp[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);
    }
}

一般好像需要進行下轉換，即后半部分和前半部分置換，即1234變成3412.

void  FFT::FFTShift(double dp[], int len)
{
    for (int i = 0; i < len / 2; i++)
    {
        double tmp = dp[i];
        dp[i] = dp[i + len / 2];
        dp[i + len / 2] = tmp;
    }
}

此時得到的應該是實部和虛部解出來的頻譜圖的Y軸電壓值，一般頻譜圖Y軸為dB，因此需要進行轉換

void FFT::getFFT(double dRe[], double dIm[], double dSp[], int len, int nBits, double dWorkingImpedance)
{
    dealFFT(dRe, dIm, dSp, len, nBits, 1); 
    FFTShift(dSp,len); //此時得到的應該是實部和虛部解出來的頻譜圖的Y軸電壓值，還需要轉換 ////dBW = 10lg(電壓^2/阻抗);dBm =dBW+30，注意電壓單位是V
    for (int i = 0; i<len; i++)
    {
        dSp[i] = 10 * log10(dSp[i] * dSp[i] / dWorkingImpedance)+30;
    }
}

getFFT()輸出之后的dp才是要的頻譜圖Y軸值，頻譜圖X軸的坐標得到通過以下方式：

//X軸精確度，采樣頻率/數據個數 = 步長
m_DeltaX_S = m_dataPara.nSampleFrequency / nDataNumOfPage_S;

data_SX[i / 2] = m_dataPara.nCenterFrequency + count*m_DeltaX_S - m_dataPara.nWorkingBandWidth/2;//中心頻率+當前點*步長-帶寬/2

在項目中，實際代碼如下：

int count = 0;
    for (int i = 0; i < nDataNumOfPage_S * 2; i++)
    {
        if (i % 2 == 0)
            data_SQ[i / 2] = data_S[i] * m_DeltaY_S;
        else
            data_SI[i / 2] = data_S[i] * m_DeltaY_S;
         
        if (i % 2 == 0)
        {
            count++;
            data_SX[i / 2] = m_dataPara.nCenterFrequency + count*m_DeltaX_S - m_dataPara.nWorkingBandWidth/2;
        }
    }
    m_dataPara.nWorkingImpedance = 50;
    FFT fft;
    int nBits = log10(nDataNumOfPage_S) / log10(2);//因為參數需要是2的N次方

fft.getFFT(data_SQ, data_SI, data_SS, nDataNumOfPage_S, nBits, m_dataPara.nWorkingImpedance); LoadData_S(data_SX, data_SS, nDataNumOfPage_S);

。。。

 

 

其他參考文章：

http://blog.sina.com.cn/s/blog_65d639d50101buo1.html

http://blog.csdn.net/hippig/article/details/8778753

http://www.makaidong.com/%E5%8D%9A%E5%AE%A2%E5%9B%AD%E6%8E%92%E8%A1%8C%E6%A6%9C/20151025/365773.html