谁人会FFT算法?高分!!!!

Question

1、MD5算法是对输入的数据进行补位，使得如果数据位长度LEN对512求余的结果是448。即数据扩展至K*512+448位。即K*64+56个字节，K为整数。具体补位操作：补一个1，然后补0至满足上述要求 
　　2、补数据长度：
　　用一个64位的数字表示数据的原始长度B，把B用两个32位数表示。这时，数据就被填补成长度为512位的倍数。 
　　3. 初始化MD5参数
　　四个32位整数 (A,B,C,D) 用来计算信息摘要，初始化使用的是十六进制表示的数字 A=0X01234567
B=0X89abcdef
C=0Xfedcba98
D=0X76543210 
　　4、处理位操作函数 
X，Y，Z为32位整数。
F(X,Y,Z) = X&Y|NOT(X)&Z
G(X,Y,Z) = X&Z|Y?(Z)
H(X,Y,Z) = X xor Y xor Z
I(X,Y,Z) = Y xor (X|not(Z)) 
　　5、主要变换过程：
　　使用常数组T[1 ... 64]， T[i]为32位整数用16进制表示，数据用16个32位的整数数组M[]表示。 
　　具体过程如下： 
/* 处理数据原文 */
For i = 0 to N/16-1 do 
/*每一次，把数据原文存放在16个元素的数组X中. */
For j = 0 to 15 do
Set X[j] to M[i*16+j].
end /结束对J的循环 
/* Save A as AA, B as BB, C as CC, and D as DD. */
AA = A
BB = B
CC = C
DD = D 
/* 第1轮*/
/* 以 [abcd k s i]表示如下操作
a = b + ((a + F(b,c,d) + X[k] + T[i]) <<< s). */ 
/* Do the following 16 operations. */
[ABCD 0 7 1] [DABC 1 12 2] [CDAB 2 17 3] [BCDA 3 22 4]
[ABCD 4 7 5] [DABC 5 12 6] [CDAB 6 17 7] [BCDA 7 22 8]
[ABCD 8 7 9] [DABC 9 12 10] [CDAB 10 17 11] [BCDA 11 22 12]
[ABCD 12 7 13] [DABC 13 12 14] [CDAB 14 17 15] [BCDA 15 22 16] 
/* 第2轮* */
/* 以 [abcd k s i]表示如下操作
a = b + ((a + G(b,c,d) + X[k] + T[i]) <<< s). */
/* Do the following 16 operations. */
[ABCD 1 5 17] [DABC 6 9 18] [CDAB 11 14 19] [BCDA 0 20 20]
[ABCD 5 5 21] [DABC 10 9 22] [CDAB 15 14 23] [BCDA 4 20 24]
[ABCD 9 5 25] [DABC 14 9 26] [CDAB 3 14 27] [BCDA 8 20 28]
[ABCD 13 5 29] [DABC 2 9 30] [CDAB 7 14 31] [BCDA 12 20 32] 
/* 第3轮*/
/* 以 [abcd k s i]表示如下操作
a = b + ((a + H(b,c,d) + X[k] + T[i]) <<< s). */
/* Do the following 16 operations. */
[ABCD 5 4 33] [DABC 8 11 34] [CDAB 11 16 35] [BCDA 14 23 36]
[ABCD 1 4 37] [DABC 4 11 38] [CDAB 7 16 39] [BCDA 10 23 40]
[ABCD 13 4 41] [DABC 0 11 42] [CDAB 3 16 43] [BCDA 6 23 44]
[ABCD 9 4 45] [DABC 12 11 46] [CDAB 15 16 47] [BCDA 2 23 48] 
/* 第4轮*/
/* 以 [abcd k s i]表示如下操作
a = b + ((a + I(b,c,d) + X[k] + T[i]) <<< s). */
/* Do the following 16 operations. */
[ABCD 0 6 49] [DABC 7 10 50] [CDAB 14 15 51] [BCDA 5 21 52]
[ABCD 12 6 53] [DABC 3 10 54] [CDAB 10 15 55] [BCDA 1 21 56]
[ABCD 8 6 57] [DABC 15 10 58] [CDAB 6 15 59] [BCDA 13 21 60]
[ABCD 4 6 61] [DABC 11 10 62] [CDAB 2 15 63] [BCDA 9 21 64] 
/* 然后进行如下操作 */
A = A + AA
B = B + BB
C = C + CC
D = D + DD 
end /* 结束对I的循环*/ 
6、输出结果。

Answer 1

VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
LONG count;

// 循环变量
int i,j,k;

// 中间变量
int bfsize,p;

// 角度
double angle;

complex<double> *W,*X1,*X2,*X;

// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W  = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X  = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}
这肯定行，我刚用完

Answer 2

///////////////////////////////////////////////////////////////////////
//  function name : fft()
//  purpose       : fourier transfer
//  input parameter:
//   xr ---pointer to input and output data real part;
//  xi ---pointer to input and output data imaginary part
//  l  ---data length;
//  inv---inv=1 means IFFT,otherwise FFT.
//  This is a standard programm.
///////////////////////////////////////////////////////////////////////
void fft(float *xr,float *xi,int l,int inv)
{
 int i,j,k,a,b,ierr;
 float ha,dr,di,pi,x;
 int cn,ar,p0,p1,po;
 float *co,*si;
 co=(float *)calloc(l,sizeof(float));
 if(co==NULL)
  {
   printf("\n Memory Error ! ");
   exit(1);
  }
 si=(float *)calloc(l,sizeof(float));
 if(si==NULL)
  {
   printf("\n Memory Error ! ");
         free(co);
   exit(1);
  }
 if(inv!=0)
  {
   for(i=0;i<l;i++)
    {
     xr[i]=xr[i]/l;
     xi[i]=xi[i]/l;
    }
  }
 po=0.0;
 i=l;
 while(i!=1){po=po+1;i=i/2;}
 pi=3.141592635897932;
 j=0;
 for(i=0;i<l-1;i++)
  {
   if(i<j)
    {
     x=xr[i];
     xr[i]=xr[j];
     xr[j]=x;
     x=xi[i];
     xi[i]=xi[j];
     xi[j]=x;
    }
   k=l/2;
   while(k<=j)
    {
     j=j-k;
     k=k/2;
    }
   j=j+k;
  }
 ha=2.0*pi/l;
 for(i=0;i<l;i++)
  {
   si[i]=sin(ha*i);
   co[i]=cos(ha*i);
  }
 if(inv==0)
  {
   for(i=0;i<l;i++)si[i]=-si[i];
  }
 a=2;
 b=1;
 for(cn=1;cn<=po;cn++)
  {
   p0=l/a;
   p1=0;
   for(k=0;k<b;k++)
    {
     i=k;
     while(i<l)
      {
       ar=i+b;
       if(k==0)
        {
         dr=xr[ar];
         di=xi[ar];
        }
       else
        {
         dr=xr[ar]*co[p1]-xi[ar]*si[p1];
         di=xi[ar]*co[p1]+xr[ar]*si[p1];
        }
       xr[ar]=xr[i]-dr;
       xi[ar]=xi[i]-di;
       xr[i]=xr[i]+dr;
       xi[i]=xi[i]+di;
       i=i+a;
      }
      p1=p1+p0;
    }
   a=2*a;
   b=b*2;
  }
 free(co);
 free(si);
 return;
}