[C++]自编FFT(递归形式)

1、对C++的运算符重载不太熟悉，所以相关复数运算操作是通过定义静态函数来操作;

2、C++的类变量貌似和java定义类变量机制不一样，PI只能用宏来定义了。

复数类：

 1 #pragma once#
 2 #include <math.h>
 3 class CComplex/*复数运算类*/
 4 {
 5 private:
 6     double m_real;   /*实部*/
 7     double m_img ;   /*虚部*/
 8 
 9 public:
10     /*构造函数和析构函数*/
11     CComplex();
12     CComplex(CComplex &c);
13     CComplex(double real, double img);
14     ~CComplex();
15     /* setters & getters */
16     void    setReal ( double value   );
17     double   getReal (   );
18     void    setImg  ( double value   );
19     double   getImg  ( );
20 
21     /*加法*/
22     static void add(CComplex &c1 ,CComplex &c2 ,CComplex &c3)
23     {
24         double _real = c1.getReal() + c2.getReal();
25         double _img = c1.getImg() + c2.getImg();
26         c3.setReal(_real);
27         c3.setImg(_img);
28     }
29 
30     /*减法*/
31     static void minus(CComplex &c1, CComplex &c2, CComplex &c3)
32     {
33         double _real = c1.getReal() - c2.getReal();
34         double _img = c1.getImg() - c2.getImg();
35         c3.setReal(_real);
36         c3.setImg(_img);
37     }
38 
39     /*乘法*/
40     static void multi(CComplex &c1, CComplex &c2, CComplex &c3)
41     {
42 
43         double a = c1.getReal();
44         double b = c1.getImg();
45         double c = c2.getReal();
46         double d = c2.getImg();
47 
48         double _real = a * c - b * d;
49         double _img = a * d + b * c;
50         c3.setReal(_real);
51         c3.setImg(_img);
52     }
53 
54     double getModulus   ();    //模长计算
55     double getAngle     ();    //幅角计算
56 };
CComplex.h
/*
    复数类
*/
#include "stdafx.h"
#include "CComplex.h"


/*构造函数和析构函数*/
CComplex::CComplex()
{
    this->m_real = 0;
    this->m_img = 0;

}
CComplex::CComplex(CComplex &c)
{
    this->m_real = c.getReal();
    this->m_img = c.getImg();
}
CComplex::CComplex(double real, double img)
{
    this->m_real = real;
    this->m_img = img;
}
CComplex::~CComplex()
{
    this->m_real = 0;
    this->m_img = 0;
}
/* setters & getters */
void CComplex::setReal(double value)
{
    this->m_real = value;
}
double CComplex::getReal()
{
    return this->m_real;
}
void CComplex::setImg(double value)
{
    this->m_img = value;
}
double CComplex::getImg()
{
    return this->m_img;
}

double CComplex::getModulus() 
{
    return sqrt(m_real*m_real + m_img * m_img);
}
double CComplex::getAngle() 
{
    double modulus = getModulus();
    if (modulus == 0) {
        return 0;
    }
    else 
    {
        return acos(m_real / modulus);
    }
}
CComplex.cpp

fft计算类：

 1 #pragma once
 2 #include "CComplex.h"
 3 #include <math.h>
 4 #include <iostream>
 5 
 6 #define PI 3.141592654
 7 class CMyFFT
 8 {
 9 
10 public:
11     CMyFFT();
12     ~CMyFFT();
13     static void doFFT(int N, int m, CComplex *x, CComplex *X); /* FFT程序*/
14     static void W_N_nk(int N, int n, int k, CComplex& factor); /* 计算旋转因子*/
15 };
CMyFFT.h
#include "stdafx.h"
#include "CMyFFT.h"



CMyFFT::CMyFFT()
{
}


CMyFFT::~CMyFFT()
{
}

void CMyFFT::doFFT(int N, int m, CComplex *x, CComplex *X)
{
    if (N != int(exp2(m + 1))) //检查输入参数
    {   
        exit(0);
        return;
    }
    if ((N != 2) && (m != 0)) /* 没到基2*/
    {
        /*
        准备存储空间
        */
        CComplex *ou = new CComplex[ N / 2 ];
        CComplex *ji = new CComplex[ N / 2 ];
        CComplex *Ak = new CComplex[ N / 2 ];
        CComplex *Bk = new CComplex[ N / 2 ];
        /*
        *奇偶分离
        */
        for (int i = 0; i < N;i++)
        {
            if (i % 2 == 0)
            {
                ou[i / 2] = x[i];
            }
            else
            {
                ji[i / 2] = x[i];
            }
        }
        //求出Ak 和 Bk
        doFFT(N / 2, m - 1, ou, Ak);
        doFFT(N / 2, m - 1, ji, Bk);
        for (int k = 0; k < N / 2; k++)
        {
            CComplex factor;
            W_N_nk(N, 1, k, factor);
            CComplex tmp;
            CComplex::multi(factor, Bk[k], tmp);
            CComplex::add(Ak[k], tmp, X[k]);//X[k] = Ak[k] + (factor*Bk[k]);
            CComplex::minus(Ak[k], tmp, X[k + N / 2]);//X[ k + N / 2 ]  = Ak[k] - (factor*Bk[k]);

        }
        delete[] ou;
        delete[] ji;
        delete[] Ak;
        delete[] Bk;
    }
    else
    {
        int r = 0;
        int k = 0;
        CComplex factor1;
        CComplex factor2;
        W_N_nk(N / 2, r, k, factor1);
        W_N_nk(N, 1, k, factor2);
        CComplex A_0;
        CComplex::multi(x[0], factor1, A_0);
        CComplex B_0;
        CComplex::multi(x[1], factor1, B_0);
        CComplex tmp;
        CComplex::multi(B_0, factor2, tmp);
        CComplex::add(A_0, tmp, X[0]);      //X[ 0     ] = A_0 + B_0*factor2;
        CComplex::minus(A_0, tmp, X[1]);    //X[ 0 + 1 ] = A_0 - B_0*factor2;
    }
}


void CMyFFT::W_N_nk(int N, int n, int k, CComplex& factor)
{
    // W_N_nk = e^( -j*2*pi*n*k/N)
    // cos(-2*pi*n*k/N) + jsin(-2*pi*n*k/N)
    factor.setReal(cos(-2 * PI * n * k / N));
    factor.setImg(sin(-2 * PI * n * k / N));
}
CMyFFT.cpp

测试代码：

 1 #include "CComplex.h"
 2 #include "CMyFFT.h "
 3 using namespace std;
 4 int main()
 5 {
 6     int m = 10;
 7     int N = int(exp2(m));
 8     CComplex x[1024];
 9     CComplex X[1024];
10 
11     for (int i = 0; i < N;i++) {
12         x[i].setReal(i);
13     }
14     CMyFFT::doFFT(N, m-1, x, X);
15 
16 
17     cout << "输出幅度值:" << endl;
18     for (int i = 0; i < 20;i++) {
19       cout << X[i].getModulus() <<endl;
20     
21     }
22     cout << "输出相角值:" << endl;
23     for (int i = 0; i < 20;i++) {
24         cout << X[i].getAngle() << endl;
25     }
26     system("pause");
27     return 0;
28 }
main.cpp

 

测试结果与Matlab自带的fft函数一致.

测试工程：

VS2017 

https://files.cnblogs.com/files/alimy/MyFFT_Cpp.zip

posted @ 2018-05-29 11:39  Alimy  阅读(722)  评论(0)  编辑  收藏  举报