我是编程的初学者,目前正在尝试开发一个需要快速傅里叶变换实现的项目.
到目前为止,我设法实现了以下内容:
有没有人有任何替代方案和建议来提高程序的速度而不会失去准确性.
short FFTMethod::FFTcalc(short int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;
/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++)
n *= 2;
/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<n;i++) {
x[i] /= n;
y[i] /= n;
}
}
return(1);
}
Dr.*_*ABT 22
我最近在Eric Postpischil 的构建高性能FFT上发现了这篇优秀的PDF .自己开发了几个FFT后,我知道与商业图书馆竞争是多么困难.相信我,如果您的FFT比英特尔或FFTW慢4倍,而不是40倍,那么你做得很好!然而,你可以竞争,这是如何.
总结一篇文章,作者指出Radix2 FFT简单但效率低,最有效的结构是radix4 FFT.一个更有效的方法是Radix8,但是这通常不适合CPU上的寄存器,所以Radix4是首选.
FFT可以分阶段构建,因此要计算1024点FFT,您可以执行Radix2 FFT的10个阶段(2 ^ 10-1024),或Radix4 FFT的5个阶段(4 ^ 5 = 1024).如果您愿意,您甚至可以在8*4*4*4*2的阶段计算1024点FFT.较少的阶段意味着对内存的读取和写入更少(FFT性能的瓶颈是内存带宽)因此动态选择基数4,8或更高是必须的.Radix4阶段特别有效,因为所有权重都是1 + 0i,0 + 1i,-1 + 0i,0-1i,并且可以编写Radix4蝶形代码以完全适合缓存.
其次,FFT中的每个阶段都不相同.第一阶段的权重都等于1 + 0i.没有必要计算这个权重,甚至乘以它,因为它是一个复数乘以1,所以第一阶段可以在没有权重的情况下进行.最后阶段也可以区别对待,并可用于执行时间抽取(位反转).Eric Postpischil的文件涵盖了所有这些.
权重可以预先计算并存储在表格中.在x86硬件上进行Sin/cos计算大约需要100-150个周期,因此预计算这些可以节省10-20%的总计算时间,因为在这种情况下,内存访问比CPU计算更快.使用快速算法一次性计算sincos是特别有益的(注意cos等于sqrt(1.0 - 正弦*正弦),或者使用表查找,cos只是正弦的相移).
最后,一旦你拥有了超级简化的FFT实现,你就可以利用SIMD矢量化来计算蝶泳程序内每个循环的4x浮点或2x双浮点运算,从而提高速度100-300%.综合以上所有内容,您将拥有一个非常漂亮且快速的FFT!
为了更进一步,您可以通过提供针对特定处理器架构的FFT阶段的不同实现来动态执行优化.每个机器的高速缓存大小,寄存器计数,SSE/SSE2/3/4指令集等不同,因此选择一种适合所有方法的方法通常会被目标例程打败.例如,在FFTW中,许多较小尺寸的FFT是针对特定架构的高度优化的展开(无环路)实现.通过组合这些较小的构造(例如RadixN例程),您可以为手头的任务选择最快和最好的例程.