最近在看快速傅立叶变换的原理,并尝试着用c++来实现 。
要实现fft,肯定要处理复数的运算,本来想着自己定义一个复数的结构体,我同学说c++自带复数类,我就直接用了。
可是当我写完程序后进行测试的时候,发现程序运行速度特别慢,很短的测试数据进行fft要花上好长时间。
我在网上搜了下,说是用c++自带的complex类运行速度在某些方面会很慢,看来还是自己定义个复数结构体或者写个复数类比较好。
以上的结论我并没有进行具体的速度上的分析,速度慢也有可能是我写的代码的效率不高,具体分析请见:
具体的FFT实现代码如下:
头文件:
#ifndef FFT_H_
#define FFT_H_
#include "stdafx.h"
#include <complex>
#define PI 3.1415926
using namespace std;
class SignalProcess
{
public:
SignalProcess();
~SignalProcess();
complex<double> multiply(const complex<double> &first, const complex<double> &second);
bool doFft(complex<double> *source, int length, bool inverse);
private:
bool mIsInverse;
int mLength;
complex<double> *mSource;
bool checkLength();
void indexReverse();
};
#endif
具体实现:
#include "stdafx.h"
#include "SignalProcess.h"
#include <xcomplex>
using namespace std;
SignalProcess::SignalProcess()
:mLength(0),
mIsInverse(false),
mSource(NULL),
{
}
SignalProcess::~SignalProcess()
{
mSource = NULL;
}
//检查输入的长度是否为2的幂次方
bool SignalProcess::checkLength()
{
int tmp = mLength;
if (tmp <= 1)
{
return false;
}
else
{
int N;
for (N = 1; tmp != 1; tmp /= 2)
{
N *= 2;
}
if (N != mLength)
{
return false;
}
else
return true;
}
}
/*
*将输入的数据进行倒序处理,i为正常序列的索引,j为倒序后的索引,
*为了避免重复交换,规定当i<j时才进行交换操作。
*如果倒序序列的索引值大于等于mLength/2,说明索引值的二进制形式
*的最高位为1,这时候要将索引值减去mLength/2,并判断减后的值与
*mLength/4的大小,以此进行。
*调用该函数之前必须调用checkLength函数,否则程序会崩溃。
*/
void SignalProcess::indexReverse()
{
int j = mLength / 2;
complex<double> tmp;
int k;
for (int i = 1; i < mLength - 2; ++i)
{
if (i < j)
{
tmp = mSource[j];
mSource[j] = mSource[i];
mSource[i] = tmp;
}
k = mLength / 2;
while(j >= k)
{
j -= k;
k /= 2;
}
j += k;
}
}
//复数相乘函数
complex<double> SignalProcess::multiply(const complex<double> &first, const complex<double> &second)
{
complex<double> tmp;
tmp.real(first.real() * second.real() - first.imag() * second.imag());
tmp.imag(first.real() * second.imag() + first.imag() * second.real());
return tmp;
}
/*
*快速傅立叶变换具体实现函数
*输入参数:
*source 复数类类型指针,指向要进行fft变换的数据
*length fft变换的长度,必须为2的幂次方
*inverse inverse为false进行正变换,为true进行反变化
*
*输出参数:
*source 将变换后的数据放在原数据的存储空间
*
*返回值类型:
*bool fft变换正确返回true,有错误返回false
*/
bool SignalProcess::doFft(complex<double> *source, int length, bool inverse)
{
mSource = source;
mLength = length;
mIsInverse = inverse;
//检查输入的fft变换长度是否为2的幂次方
if (checkLength() == false)
{
return false;
}
//进行倒序
indexReverse();
int number = mLength;
int m;
//m是蝶形运算的级数
for (m = 1; (number /= 2) != 1; ++m)
{
}
int level, num;
int i, ip;
complex<double> v, w, tmp;
//最外层循环控制蝶形的级数
for (int k = 1; k <= m; ++k)
{
level = (int)pow((float)2.0, k);
num = level / 2;
v.real(1.0);
v.imag(0.0);
if (mIsInverse)
{
w.real(cos(PI / num));
w.imag(sin(PI / num));
}
else
{
w.real(cos(PI / num));
w.imag(-sin(PI / num));
}
//中层循环控制旋转因子的个数,也就是不同的Wn的个数
for (int j = 0; j < num; ++j)
{
//内层循环控制特定的Wn对应的运算
for (i = j; i < mLength; i += level)
{
ip = i + num;
if (mIsInverse)
{
mSource[i].real(mSource[i].real() * 0.5) ;
mSource[i].imag(mSource[i].imag() * 0.5) ;
mSource[ip].real(mSource[ip].real() * 0.5) ;
mSource[ip].imag(mSource[ip].imag() * 0.5) ;
}
tmp = multiply(v, mSource[ip]);
//蝶形运算
mSource[ip] = mSource[i] - tmp;
mSource[i] += tmp;
}
v = multiply(v, w);
}
}
}