N点复数FFT计算2N点实数FFT

最新推荐文章于 2024-07-14 21:59:41 发布

witcherlucas

最新推荐文章于 2024-07-14 21:59:41 发布

阅读量3.3k

点赞数 4

分类专栏：数字信号处理文章标签：算法 c++

本文链接： https://blog.csdn.net/witcherlucas/article/details/128419459

版权

数字信号处理专栏收录该内容

2 篇文章 1 订阅

订阅专栏

概述

一般来讲，不考虑效率的时候，实数作2N点FFT时可以将虚部设为0，直接调用2N点FFT的程序。但最近在公司读之前的代码，读了快一个月没读懂，最后还是前辈讲了之后才发现利用了N点FFT计算2N点实数FFT这个技巧，在这里做个简单的分享。

算法原理

首先贴上教材里对这个方法的说明：

注意，教材这里讲对称性时认为Ｘ（Ｎ）＝Ｘ（０）

代码实现

首先将输入的实数序列分奇偶项，奇数项作实部，偶数项作虚部构造一个复数序列。（Ｃｏｍｐｌｅｘ是自己写的复数类）。

    //将2N点的实数变换为N点的复数
    vector<Complex> vin(input_len);
    for (i = 0; i < input_len; i++)
    {
        
        vin[i].real = realsig[2 * i];
        vin[i].image = realsig[2 * i + 1];
    }

这里realsig就是输入的实数信号，S１６表示１６位有符号数，换成ｉｎｔ＼ｄｏｕｂｌｅ之类的都可以。

做完FFT按教材上的公式进行转换，注意下面的三目运算符的目的是使得ｘ（Ｎ）＝ｘ（０）

    FFT(vin);
    vector<Complex> vout(input_len);
    for (int i = 0; i < input_len; i++)
    {
        Complex XN_K = i ? vin[input_len - i] : vin[0];
        XN_K.conj();
        double er = 2;
        Complex fuj(0, -1 / 2);
        Complex X1 = (vin[i]+XN_K) ;
        X1 = X1 / er;
        Complex X2 =  (vin[i] - XN_K)* fuj;
        Complex W2N_K(cos(-pi / input_len * i), sin(-pi / input_len * i));
        vout[i] =  W2N_K * X2+X1;

    }

这样，ｖｏｕｔ中就是实数２N点FFT的前N个结果。至于后N个值可用上面书上的公式得到。

将结果在matlab中对比，当N＝１０２４时结果如下：

可见C中结果正好是实数２０４８点FFT的前半部分。

完整代码

为方便大家使用，我把整个完整代码贴在下面。

#include<iostream>
using namespace std;

#include<vector>

#include<fstream>
#include<string>
class Complex
{
public:
    double real;
    double image;
    Complex operator+(Complex& x);
    Complex operator-(Complex& x);
    Complex operator*(Complex& x);
    
    friend Complex operator/(Complex&c,double& x);
    Complex(double r = 0, double i = 0);
    void conj();
    double amp();
};

Complex Complex::operator+(Complex& x)
{
    Complex ans = *this;
    ans.real += x.real;
    ans.image += x.image;
    return ans;


}

Complex Complex::operator-(Complex& x)
{
    Complex ans = *this;
    ans.real -= x.real;
    ans.image -= x.image;
    return ans;
}

Complex Complex::operator*(Complex& x)
{
    Complex ans;
    ans.real = real * x.real - image * x.image;
    ans.image = real * x.image + image * x.real;
    return ans;
}

Complex::Complex(double r, double i)
{
    real = r;
    image = i;
}

double Complex::amp()
{
    return sqrt(real * real + image * image);
};



void Complex::conj()
{
    this->image = -this->image;
}

Complex operator/(Complex& c, double& x)
{
    Complex ans;
    ans.real = c.real / x;
    ans.image = c.image / x;
    return ans;
}

void FFT(vector<Complex>& A)
{
    int m = A.size();
    if (m <= 1)
        return;
    int i = 0;
    vector<Complex> even(m / 2);
    vector<Complex> odd(m / 2);
    for (i = 0; i < m - 1; i += 2)
    {
        odd[i >> 1] = A[i + 1];
        even[i >> 1] = A[i];
    }
    FFT(odd);
    FFT(even);
    for (int k = 0; k < m / 2; k++)
    {
        Complex t = Polar(1, -2 * pi * k / m) * odd[k];
        A[k] = even[k] + t;
        A[k + m / 2] = even[k] - t;
    }
}

void fin(string path, vector<S16>& v)
{
    //从文件中读取数据到v中
    //v为空
    ifstream infile;   //输入流
    string path1 = path;
    infile.open(path1, ios::in);
    cout << "isopen:" << infile.is_open() << endl;

    int m = 0;
    if (infile.is_open())
    {
        string s;
        while (getline(infile, s))
        {
            int num = stod(s);  //放缩
            v.push_back(num);
            m++;
            if (m > 2048)
                break;
        }

    }
    infile.close();
}


int main()
{

    write_twiddle(9);

    int i;

    S32 Stage_N = 5;
    int realN = Stage_N * 2;
    ifstream infile;   //输入流
    string path1 = "../FFTx.txt";
    string pathin2 = "../matlab_fftimagin.txt";
    vector<ｉｎｔ> realsig;  //实数信号

    fin(path1, realsig);
    
    int input_len = 1<<(2* Stage_N);
    //将2N点的实数变换为N点的复数(即1024点的实数切换成512点复数）
    vector<Complex> vin(input_len);
    for (i = 0; i < input_len; i++)
    {
        
        vin[i].real = realsig[2 * i];
        vin[i].image = realsig[2 * i + 1];
    }
    FFT(vin);
    vector<Complex> vout(input_len);
    for (int i = 0; i < input_len; i++)
    {
        Complex XN_K = i ? vin[input_len - i] : vin[0];
        XN_K.conj();
        double er = 2;
        Complex fuj(0, -1 / 2);
        Complex X1 = (vin[i]+XN_K) ;
        X1 = X1 / er;
        Complex X2 =  (vin[i] - XN_K)* fuj;
        Complex W2N_K(cos(-pi / input_len * i), sin(-pi / input_len * i));
        vout[i] =  W2N_K * X2+X1;

    }
｝

其实公司代码不长这样，这里只是为了方便理解写了复数类，FFT也是自己写的一个极低效率版本。如有错误还请指出。