快速傅裡葉變換C++遞迴演算法實現

快速傅裡葉變換C++遞迴演算法實現

 

    網上有些演算法資料經測試回合結果是錯誤的，雖然代碼的使用的是非遞迴形式。為了方便驗證快速傅裡葉變換的準確性，我提供了自己設計的遞迴演算法。

  基於時域抽取的“基2”快速傅裡葉變換演算法代碼：

    Fouier.h檔案:

#pragma once
#include"Complex.h"
class Fouier
{
    Complex * data;
    void fft(int start,int step,int len);
    Complex W(int k,int n);//e^(-i*2*pi*k/n)
public:
    Fouier(void);
    ~Fouier(void);
    
    void fft();
};

    Fouier.c檔案:

#include "Fouier.h"
#include<iostream>
using namespace std;
#include<cmath>
#include<ctime>
#define DATALEN 32
#define KEYVALUE 10000 //產生隨機浮點數的值，保證分子和分母在這個值之內
#define PI 3.14159265358979323846
Fouier::Fouier(void)
{
    data=new Complex[DATALEN];
    srand(unsigned int(time(0)));
    cout<<"來源資料："<<endl;
    for(int i=0;i<DATALEN;i++)
    {
        data[i]=(rand()%(KEYVALUE))/(double)(rand()%(KEYVALUE)+1);
        if(i%5==0&&i!=0)
            cout<<endl;
        cout<<data[i]<<" ";
    }
    cout<<endl;
}

Fouier::~Fouier(void)
{
    delete [] data;
}
Complex Fouier:: W(int k,int n)//歐拉公式
{
    double alpha=-2*PI*k/n;
    return Complex(cos(alpha),sin(alpha));
}
void Fouier::fft(int start,int step,int len)
{
    if(len==1)//一個元素
    {
        //一個元素不需要變換
        return ;
    }
    fft(start,step*2,len/2);//X1(k)
    fft(start+step,step*2,len/2);//X2(k)
    Complex X1,X2;
    for(int i=0;i<len/2;i++)
    {
        X1=data[start+step*i*2];
        X2=data[start+step*(i*2+1)];
        //計算X(k):k=0~N/2-1
        data[start+step*i]=X1+W(i,len)*X2;
        //計算X(k):k=N/2~N-1
        data[start+step*(i+len/2)]=X1-W(i,len)*X2;
    }
}
void Fouier::fft()
{
    fft(0,1,DATALEN);
    cout<<"變換後資料："<<endl;
    for(int i=0;i<DATALEN;i++)
    {
        if(i%5==0&&i!=0)
            cout<<endl;
        cout<<data[i]<<" ";
    }
}

  Complex.h檔案：

#pragma once
#include<iostream>
using namespace std;
class Complex//a+b*i
{
    double a;//實數部分
    double b;//虛數部分
public:
    Complex(double a=0,double b=0);
    //+操作
    friend Complex operator +(Complex &x,Complex &y);
    friend Complex operator +(double x,Complex &y);
    friend Complex operator +(Complex &x,double y);
    //-操作
    friend Complex operator -(Complex &x,Complex &y);
    friend Complex operator -(double x,Complex &y);
    friend Complex operator -(Complex &x,double y);
    //*操作
    friend Complex operator *(Complex &x,Complex &y);
    friend Complex operator *(double x,Complex &y);
    friend Complex operator *(Complex &x,double y);
    //=操作
    Complex operator =(Complex &x);
    Complex operator =(double x);
    //<<操作
    friend ostream & operator<<(ostream&out,Complex&c);

    ~Complex(void);
};

    Complex.c檔案:

#include "Complex.h"

Complex::Complex(double a,double b)//虛部預設是0
{
    this->a=a;
    this->b=b;
}

Complex::~Complex(void)
{
}

Complex operator +(Complex &x,Complex &y)
{
    return Complex(x.a+y.a,x.b+y.b);
}
Complex operator +(double x,Complex &y)
{
    return Complex(x+y.a,y.b);
}
Complex operator +(Complex &x,double y)
{
    return Complex(x.a+y,x.b);
}

Complex operator -(Complex &x,Complex &y)
{
    return Complex(x.a-y.a,x.b-y.b);
}
Complex operator -(double x,Complex &y)
{
    return Complex(x-y.a,-y.b);
}
Complex operator -(Complex &x,double y)
{
    return Complex(x.a-y,x.b);
}

Complex operator *(Complex &x,Complex &y)
{
    return Complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
}
Complex operator *(double x,Complex &y)
{
    return Complex(x*y.a,x*y.b);
}
Complex operator *(Complex &x,double y)
{
    return Complex(x.a*y,x.b*y);
}

Complex Complex::operator =(Complex &x)
{
    a=x.a;
    b=x.b;
    return *this;
}
Complex Complex::operator =(double x)
{
    a=x;
    b=0;
    return *this;
}
ostream & operator<<(ostream&out,Complex&c)
{
    if(c.a!=0||c.a==0&&c.b==0)
        out<<c.a;
    if(c.b!=0)
    {
        if(c.b>0)
            out<<"+";
        if(c.b!=1)
            out<<c.b;
        out<<"i";
    }
    return out;
}

    main.c檔案：

#include<iostream>
using namespace std;
#include"Fouier.h"

int main()
{
    Fouier f;
    f.fft();
    return 0;
}

    如有錯誤，歡迎批評指正！

    參考資料：http://zhoufazhe2008.blog.163.com/blog/static/63326869200971010421361/

       維基百科：http://zh.wikipedia.org/wiki/%E5%BF%AB%E9%80%9F%E5%82%85%E7%AB%8B%E5%8F%B6%E5%8F%98%E6%8D%A2