[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