1. 程式人生 > >【演算法】快速傅立葉變換(FFT)的遞迴實現

【演算法】快速傅立葉變換(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計算出的結果: