【演算法】快速傅立葉變換(FFT)的遞迴實現
阿新 • • 發佈:2018-12-17
FFT是數字訊號處理中的重要演算法,在matlab中可以直接呼叫fft()函式,本文是C++版的FFT演算法,用遞迴方式進行實現。
//================================== //Signal_Process_FFT_v1.cpp //================================== #include<iostream> #include<vector> #include<complex> #include<cmath> using namespace std; //---------------------------------- class Signal { vector<complex<double> > xt; void init() {}; const double PI=3.141592653; complex<double> W(int N, int k); void fft(vector<complex<double> >& x,int N); public: Signal(complex<double> x[],int xLength); void display(); void FFT(); };//-------------------------------- Signal::Signal(complex<double> x[], int xLength) { if (nullptr==x || xLength <= 0) { cerr << "Input signal illegal."; exit(0); } xt.assign(x, x + xLength); init();//若訊號長度不是2的冪次,需要進行補0,此處略 }//--------------------------------- void Signal::display() { cout << "["; for (int i = 0; i < xt.size() - 1; i++) cout << xt[i] << ","; cout << xt[xt.size() - 1] << "]\n"; }//---------------------------------- complex<double> Signal::W(int N, int k) { return exp(complex<double>(0, -2 * PI*k / N)); }//---------------------------------- void Signal::fft(vector<complex<double> >& x,int N) { if (2 == x.size()) { complex<double> x0 = x[0]; complex<double> x1 = x[1]; x[0] = x0 + x1*W(N, 0); x[1] = x0 - x1*W(N, 0); } else { vector<complex<double> > y; for (int i = 0; i < x.size(); i++) y.push_back(x[i]); vector<complex<double> > g,h; for (int r = 0; r < y.size() / 2; r++) { g.push_back(y[2*r]); h.push_back(y[2*r+1]); } fft(g, N); fft(h, N); for (int k = 0; k < y.size() / 2; k++) { x[k] = g[k] + h[k] * W(N,k*N/y.size()); x[k+y.size()/2] = g[k] - h[k] * W(N,k*N/y.size()); } } }//---------------------------------- void Signal::FFT() { int N = xt.size(); fft(xt,N); }//---------------------------------- int main() { complex<double> x[] = {1,2,3,4,5,6,7,complex<double>(10.0,1.0) }; Signal s(x,sizeof(x)/sizeof(x[0])); s.display(); s.FFT(); s.display(); return 0; }//=================================
程式執行結果:
matlab計算出的結果: