自己寫了一個,和Matlab對比了一下,結果是一樣的,供各位參考吧

// ==============================================================================
// 快速離散傅立葉變換和功率譜
// 一維快速傅立葉變換FFT1和二維快速傅立葉變換FFT2
// 測試環境 C++ builder 2010
// MinGW (GCC) 4.5
// wuxuping 2010-08-10
// 一定要將檔案儲存為UTF-8的格式,這是MinGW的預設格式.
// ==============================================================================

#ifndef FFT_H
#define FFT_H
// --------------------
#include <vector>
#include<complex>
#include<bitset>
#include <iterator>
#include <algorithm>
#include <iostream>
#include <sstream>
using namespace std;

const double pai = 3.1415926535897932384626433832795028841971; // 圓周率

typedef vector<vector<complex<double> > >VVCD;
typedef vector<vector<double> >VVD;

// -------------------------------------------------
class FFT1 {
private:
	unsigned N; // 資料長度

	unsigned NFFT; // 使用快速傅立葉變化所計算的長度

	unsigned r; // 資料長度NFFT=2^r,且2^(r-1)<N<<NFFT

	vector<complex<double> >TS; // 時間域序列

	vector<complex<double> >FS; // 頻率域序列FrequencySeries

	vector<double>PSD; // PSD功率譜密度(離散譜)

	vector<double>PS; // 返回PS=PSD*Frequency
	int sign;

	// sign > 0為正變換否則為逆變換
	// ------------------------------------
	void Set_NFFT_r(unsigned n);
	unsigned RevBit(unsigned num);
	void GetIndex(unsigned L, unsigned ALindex, unsigned &AL_1index0,
		unsigned &AL_1index1, unsigned &Windex);

	// ------------------------------------

	void doFFT1();

	// 執行傅立葉變換
	// ------------------------------------
public:

	// ------實序列建構函式-------------
	FFT1(vector<double>TimeSeries, int sign_ = 1) {
		N = TimeSeries.size();
		if (sign_ > 0) {
			sign = 1.0; // 正變換
		}
		else {
			sign = -1.0; // 逆變換
		}
		Set_NFFT_r(N);

		complex<double>W0(0.0, 0.0); // 旋轉因子
		// --------------------------------
		for (unsigned i = 0; i < NFFT; i++) { // 賦初值
			if (i < N) {
				FS.push_back(W0);
				complex<double>Wi(TimeSeries[i], 0.0); // 旋轉因子
				TS.push_back(Wi);
			}
			else {
				FS.push_back(W0); // 補零
				TS.push_back(W0); // 補零
			}

		}
		// --------------------------------

		doFFT1(); // 執行傅立葉變換
		// --------------------------------
	};

	// -------復序列建構函式-------------
	FFT1(vector<complex<double> >TimeSeries, int sign_ = 1) { // 賦初值
		N = TimeSeries.size();
		if (sign_ > 0) {
			sign = 1.0; // 正變換
		}
		else {
			sign = -1.0; // 逆變換
		}
		Set_NFFT_r(N);

		complex<double>W0(0.0, 0.0); // 旋轉因子
		// --------------------------------
		for (unsigned i = 0; i < NFFT; i++) { // 賦初值
			if (i < N) {
				FS.push_back(W0);
				TS.push_back(TimeSeries[i]);
			}
			else {
				FS.push_back(W0); // 補零
				TS.push_back(W0); // 補零
			}

		}
		// --------------------------------

		doFFT1(); // 執行傅立葉變換
		// --------------------------------
	};

	// ----------成員函式------

	// ----------成員函式 -------------
	vector<complex<double> >GetFS() { // 返回頻率域序列FrequencySeries
		return FS;
	};

	// ----------成員函式 -------------
	vector<double>GetFrequency() { // 返回頻率
		vector<double>Frequency;
		// ----------------------
		for (unsigned int i = 0; i < FS.size(); i++) {
			Frequency.push_back(i);
		}
		return Frequency;
	};

	// ----------成員函式 -------------
	vector<double>GetPSD() { // 返回FS^2
		if (PSD.size() > 0) {
			PSD.clear();
		}
		// ----------------------
		for (unsigned int i = 0; i < FS.size(); i++) {
			PSD.push_back(real(FS[i] * conj(FS[i])));
		}
		return PSD;
	};

	vector<double>GetPS() { // 返回PS=PSD*Frequency
		if (PS.size() > 0) {
			PS.clear();
		}
		// ----------------------
		for (unsigned int i = 0; i < FS.size(); i++) {
			PS.push_back(i * 1.0 * real(FS[i] * conj(FS[i])));
		}
		return PS;
	};
	// ----------解構函式 -------------

	~FFT1() {
	};
	// -----------------------------------
};

// ------設定NFFT和r的值--------
void FFT1::Set_NFFT_r(unsigned n) {
	bitset<32>bsn(n);
	r = 0;
	while ((n >> r) > 0) {
		r++;
	}

	if ((bsn.count() == 1) && (r > 0)) {
		r--;
	}
	bitset<32>bsNFFT(0);
	bsNFFT.set(r);
	NFFT = (unsigned)bsNFFT.to_ulong();
};

// ---- 位元倒序函式
unsigned FFT1::RevBit(unsigned num) {
	bitset<32>bs(num);
	bitset<32>bsR;
	// ---------逆序
	for (unsigned i = 0; i < r; i++) {
		bsR.set(i, bs[r - i - 1]);
	}
	return(unsigned)bsR.to_ulong();
};

// ------------------------------------------------------------------------------
// 計算公式中的索引AL(ALindex)=AL-1(AL_1index0)+AL-1(AL_1index1)W(Windex)
// ------------------------------------------------------------------------------
void FFT1::GetIndex(unsigned L, unsigned ALindex, unsigned &AL_1index0,
	unsigned &AL_1index1, unsigned &Windex) {
	// ALindex為AL()的索引
	// AL_1index0 = 0; // AL-1()的第一項索引
	// AL_1index1 = 0; // AL-1()的第二項索引
	// Windex = 0; // 相位角中的索引
	bitset<32>bs(ALindex);
	bitset<32>bs1(ALindex);
	bs1.set(r - L, 0);
	AL_1index0 = bs1.to_ulong();
	// -----------------------------------
	bitset<32>bs2(ALindex);
	bs2.set(r - L, 1);
	AL_1index1 = bs2.to_ulong();
	// -----------------------------------
	bitset<32>bs3; // 左L位
	for (unsigned i = 0; i < L; i++) {
		bs3.set(r - L + i, bs[r - i - 1]);
	}
	Windex = bs3.to_ulong();
}

// ------------------------------------
void FFT1::doFFT1() { // 一維快速Fourier變換
	unsigned AL_1index0 = 0; // AL-1()的第一項索引
	unsigned AL_1index1 = 0; // AL-1()的第二項索引
	unsigned Windex = 0; // 相位角中的索引
	double alpha; // 相位角
	// ----------------------------
	vector<complex<double> >X = TS; // X中間臨時變數

	for (unsigned L = 1; L <= r; L++) { // 蝶形計算過程
		for (unsigned ALindex = 0; ALindex < NFFT; ALindex++) {
			// ALindex為AL()的索引
			GetIndex(L, ALindex, AL_1index0, AL_1index1, Windex);
			alpha = -2.0 * sign * pai * Windex / (1.0 * NFFT); // 旋轉因子的相位角
			complex<double>W(cos(alpha), sin(alpha)); // 旋轉因子
			FS[ALindex] = X[AL_1index0] + X[AL_1index1] * W;
		}
		X = FS; // 複製陣列
	}
	//
	if (sign > 0) { // 為正變換
		for (unsigned i = 0; i < NFFT; i++) {
			FS[i] = X[RevBit(i)]; // 索引按二進位制倒序
		}
	}
	else { // 為逆變換
		complex<double>WN(NFFT, 0.0); // 資料的個數
		for (unsigned i = 0; i < NFFT; i++) {
			FS[i] = X[RevBit(i)] / WN; // 索引按二進位制倒序
		}
	}
}
// ==============================================================================
// 計算二維快速傅立葉變換FFT2
// ==============================================================================

// ----------------------------------------------
class FFT2 {

private:

	// ----------------------------------------------
	unsigned N1; // 資料行長度 0<=k1<N1  0<=j1<N1

	unsigned N2; // 資料列長度 0<=k2<N2  0<=j2<N2

	int sign; // sign > 0為正變換否則為逆變換

	VVCD Tk1k2; // 二維時間域序列,行k1列k2

	VVCD Yj1k2; // 中間序列對 Ak1k2每一行(k1)做傅立葉變換的結果

	VVCD Fj1j2; // 頻率域序列FrequencySeries

	VVD PSD; // PSD功率譜密度(離散譜)

	VVD PS; // 返回PS=PSD*Frequency

	// ------------------------------------
	VVCD TranspositionVVCD(VVCD dv); // 矩陣轉置
	// ------------------------------------

public:

	// ------實序列建構函式-------------
	FFT2(VVD TK1K2, int sign_ = 1) { // 賦初值

		N1 = TK1K2.size(); // 獲取行數
		Tk1k2.resize(N1);
		for (unsigned k1 = 0; k1 < N1; k1++) {
			N2 = TK1K2[k1].size();
			Tk1k2[k1].resize(N2);
			for (unsigned k2 = 0; k2 < N2; k2++) {
				complex<double>W(TK1K2[k1][k2], 0.0);
				Tk1k2[k1][k2] = W;
			}

		}
		// -----------------------------------------
		for (unsigned k1 = 0; k1 < N1; k1++) {
			FFT1 fft1(Tk1k2[k1], sign_);
			Yj1k2.push_back(fft1.GetFS());
		}
		// -----------------------------------------
		Yj1k2 = TranspositionVVCD(Yj1k2); // 將矩陣轉置
		// -----------------------------------------
		N2 = Yj1k2.size(); // 獲取列數
		for (unsigned k2 = 0; k2 < N2; k2++) {
			FFT1 fft1(Yj1k2[k2], sign_);
			Fj1j2.push_back(fft1.GetFS());
		}
		// --------------
		Fj1j2 = TranspositionVVCD(Fj1j2); // 將矩陣轉置
		// -----------------------
	};

	// -------復序列建構函式-------------
	FFT2(VVCD TK1K2, int sign_ = 1) { // 賦初值
		Tk1k2 = TK1K2;
		N1 = Tk1k2.size(); // 獲取行數
		for (unsigned k1 = 0; k1 < N1; k1++) {
			FFT1 fft1(Tk1k2[k1], sign_);
			Yj1k2.push_back(fft1.GetFS());
		}
		// -----------------------------------------
		Yj1k2 = TranspositionVVCD(Yj1k2); // 將矩陣轉置
		// -----------------------------------------
		N2 = Yj1k2.size(); // 獲取列數
		for (unsigned k2 = 0; k2 < N2; k2++) {
			FFT1 fft1(Yj1k2[k2], sign_);
			Fj1j2.push_back(fft1.GetFS());
		}
		// --------------
		Fj1j2 = TranspositionVVCD(Fj1j2); // 將矩陣轉置
		// -----------------------
	};

	// ----------成員函式 -------------
	VVCD GetFS() { // 返回頻率域序列Fj1j2
		return Fj1j2;
	};

	// ----------成員函式 -------------
	VVD GetPSD() { // 返回FS^2
		PSD.resize(Fj1j2.size());
		for (unsigned i = 0; i < Fj1j2.size(); i++) {
			PSD[i].resize(Fj1j2[i].size());
		}
		// ----------------------
		for (unsigned i = 0; i < Fj1j2.size(); i++) {
			for (unsigned j = 0; j < Fj1j2[i].size(); j++) {
				PSD[i][j] = real(Fj1j2[i][j] * conj(Fj1j2[i][j]));
			}
		}
		return PSD;
	};

	// ----------成員函式 -------------
	VVD GetPS() { // 返回PS=PSD*Frequency
		PS.resize(Fj1j2.size());
		for (unsigned i = 0; i < Fj1j2.size(); i++) {
			PS[i].resize(Fj1j2[i].size());
		}
		// ----------------------
		for (unsigned i = 0; i < Fj1j2.size(); i++) {
			for (unsigned j = 0; j < Fj1j2[i].size(); j++) {
				PS[i][j] = (double)i * (double)j * real
					(Fj1j2[i][j] * conj(Fj1j2[i][j]));
			}
		}
		return PS;
	};
	// ----------解構函式 -------------

	~FFT2() {
	};
	// -----------------------------------
};

// ----------------------------------------
VVCD FFT2::TranspositionVVCD(VVCD dv) { // 將矩陣轉置
	unsigned dvRow = dv.size();
	unsigned maxcol = 0;
	unsigned mincol = 99999;
	VVCD temp;
	if (dv.size() > 0) {
		// ------------------------------
		for (unsigned i = 0; i < dvRow; i++) {
			if (maxcol < dv[i].size()) {
				maxcol = dv[i].size();
			}
			if (mincol > dv[i].size()) {
				mincol = dv[i].size();
			}
		}
		// ------------------------------
		if (mincol == maxcol && mincol > 0) {
			temp.resize(mincol);
			for (unsigned i = 0; i < mincol; i++) {
				temp[i].resize(dvRow);
			}
			for (unsigned i = 0; i < dvRow; i++) {
				for (unsigned j = 0; j < mincol; j++) {
					temp[j][i] = dv[i][j];
				}
			}
		}
	}
	return temp;
}

// ==============================================================================
// 此處定義的函式是為了控制檯輸出檢視方便
// ==============================================================================
// VVCD==vector<vector<complex<double> > >
// VVD== vector<vector<double> >
// ---------------------------------
void PrintVVCD(vector<vector<complex<double> > >Data) {
	cout << "Data : " << endl;
	for (unsigned i = 0; i < Data.size(); i++) {
		cout << "Row[" << i << "]=\t";
		vector<complex<double> >v = Data[i];
		copy(v.begin(), v.end(), ostream_iterator<complex<double> >(cout,
				"\t"));
		cout << endl;
	}
}

// ---------------------------------
void PrintVCD(vector<vector<double> >Data) {
	cout << "Data : " << endl;
	for (unsigned i = 0; i < Data.size(); i++) {
		cout << "Row[" << i << "]=\t";
		vector<double>v = Data[i];
		copy(v.begin(), v.end(), ostream_iterator<double>(cout, "\t"));
		cout << endl;
	}
}
// ---------------------------------

// --------------------------------------------------------------------------------
#endif