自己写了一个,和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