【算法】快速傅里叶变换(FFT)的递归实现
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计算出的结果: