FFT算法

对于变换长度为N的序列x(n)其傅立叶变换可以表示如下：

	N	nk
X(k)=DFT[x(n)] =         Σx(n)W
	n="0"

　　　　　　　　　　　　　　　　　　　　　式（１）

其中，W="exp"(-2π/N)。

在下面的源码里,对应函数为 DFT

 

算法导论 里介绍了一种递归计算fft地方法, 函数为FFT_recursive

主要利用了 DFT(x) = DFT[0](x) + w*DFT[1](x)

继而导出迭代的fft计算方法,就是现在最常用的,函数为FFT_Iter

步骤为:

将输入数组元素反置变换

for 树的每一层:

    for 这一层上需要计算的每一对数据:

        根据这一对数据,及蝶形公式,计算上一层的数据

 另外,代码里写了一个简单的复数类和向量计算函数

 

  1 #include  < vector >
  2 #include  < math.h >
  3 #include  < assert.h >
  4 // #include <windows.h>
  5
  6 using   namespace  std;
  7
  8 // ************************************/
  9 //  Fushu
 10 /**/ /**************************************/
 11 class  FuShu
 12 {
 13    public:
 14        FuShu(double r=0, double i=0): x(r), y(i)
 15        {
 16        }
 17        
 18        FuShu & operator+= (FuShu const &rhs)
 19        {
 20            x += rhs.x;
 21            y += rhs.y;
 22            return *this;
 23        }
 24        
 25        FuShu & operator-= (FuShu const &rhs)
 26        {
 27            x -= rhs.x;
 28            y -= rhs.y;
 29            return *this;
 30        }
 31        
 32        FuShu & operator*= (FuShu const &rhs)
 33        {
 34            double x2 = x * rhs.x - y * rhs.y;
 35            double y2 =    x * rhs.y + y * rhs.x;
 36            x=x2;
 37            y=y2;
 38            return *this;
 39        }
 40        
 41    //private:
 42        double x,y;
 43} ;
 44
 45 FuShu  operator +  (FuShu  const   &  lhs, FuShu  const   &  rhs)
 46 {
 47    FuShu res = lhs;
 48    return res+=rhs;
 49}     
 50
 51 FuShu  operator -  (FuShu  const   &  lhs, FuShu  const   &  rhs)
 52 {
 53    FuShu res = lhs;
 54    return res-=rhs;
 55}         
 56
 57 FuShu  operator *  (FuShu  const   &  lhs, FuShu  const   &  rhs)
 58 {
 59    FuShu res = lhs;
 60    return res*=rhs;
 61}     
 62
 63 double  fabs(FuShu  const   &  lhs)
 64 {
 65    return fabs(lhs.x) + fabs(lhs.y);
 66}
 67 //  bool operator== (FuShu const & lhs, FuShu const & rhs)
 68 //  {
 69      //  return lhs.x == rhs.x && lhs.y == rhs.y;
 70 //  }    
 71
 72 void  Print(FuShu  const   &  lhs)
 73 {
 74    if(lhs.y>=0)
 75    printf("%g+%gj ", lhs.x, lhs.y);
 76    else
 77    printf("%g%gj ", lhs.x, lhs.y);
 78}
 79
 80 void  Print(vector  < FuShu >   const   &  lhs)
 81 {
 82    size_t n=lhs.size();
 83    for(size_t i=0; i<n; i++)
 84        Print(lhs[i]);
 85    printf("\n");
 86}
 87
 88
 89 // ************************************/
 90 //  vector
 91 /**/ /**************************************/
 92
 93 template  < typename T >
 94 vector  < T >   operator +  (vector  < T >   const   &  lhs, vector  < T >   const   &  rhs)
 95 {
 96    assert(lhs.size() == rhs.size());
 97    vector <T> res;
 98    size_t n = lhs.size();
 99    res.resize(n);
100    for(size_t i=0; i<n; i++)
101        res[i] = lhs[i] + rhs[i];
102    return res;
103}
104
105 template  < typename T >
106 vector  < T >   operator -  (vector  < T >   const   &  lhs, vector  < T >   const   &  rhs)
107 {
108    assert(lhs.size() == rhs.size());
109    vector <T> res;
110    size_t n = lhs.size();
111    res.resize(n);
112    for(size_t i=0; i<n; i++)
113        res[i] = lhs[i] - rhs[i];
114    return res;
115}
116
117
118 template  < typename T >
119 vector  < T >   operator *  (FuShu  const   &  lhs, vector  < T >   const   &  rhs)
120 {
121    vector <T> res;
122    size_t n = rhs.size();
123    res.resize(n);
124    for(size_t i=0; i<n; i++)
125        res[i] = lhs[i] * rhs[i];
126    return res;
127}
128
129 template  < typename T >
130 double  fabs(vector  < T >   const   &  lhs)
131 {
132    double res(0);
133    size_t n=lhs.size();
134    
135    for(size_t i=0; i<n; i++)
136        res += fabs(lhs[i]);
137    return res;
138}
139 //  template <typename T>
140 //  bool operator== (FuShu const & lhs, vector <T> const & rhs)
141 //  {
142      //  size_t n = lhs.size();
143      //  if(n != rhs.size())
144          //  return false;
145         
146      //  for(size_t i=0; i<n; i++)
147          //  if(lhs[i] != rhs[i])
148              //  return false;
149             
150      //  return true;
151 //  }
152
153 template  < typename T >
154 vector  < T >   &   operator << (vector  < T >   &  lhs, T  const   & rhs)
155 {
156    lhs.push_back(rhs);
157    return lhs;
158}
159
160 /**/ /***************************************************/
161 //  FFT
162 /**/ /***************************************************/
163 vector < FuShu >  DFT(vector < FuShu >   const   & a)
164 {
165    size_t n=a.size();
166    if(n==1) return a;
167    
168    vector<FuShu> res(n);
169    for(size_t k=0; k<n; k++)
170    {
171        FuShu wnk(cos(6.28*k/n), sin(6.28*k/n));
172        FuShu w = 1;
173        for(size_t j=0; j<n; ++j, w*=wnk)
174            res[k] += a[j] * w;
175    }
176    
177    return res;
178}
179
180 vector < FuShu >  FFT_recursive(vector < FuShu >   const   & a)
181 {
182    size_t n=a.size();
183    if(n==1) return a;
184    
185    FuShu wn(cos(6.28/n), sin(6.28/n));
186    FuShu w(1);
187    vector<FuShu> a0, a1;
188    a0.reserve(n/2);
189    a1.reserve(n/2);
190    for(size_t i=0; i<n/2; i++)
191    {
192        a0.push_back(a[i*2]);
193        a1.push_back(a[i*2+1]);
194    }
195    
196    vector<FuShu> y0 = FFT_recursive(a0);
197    vector<FuShu> y1 = FFT_recursive(a1);
198    
199    vector<FuShu> vy;
200    vy.resize(n);
201    for(size_t k=0; k<n/2; k++)
202    {
203        vy[k] = y0[k] + w * y1[k];
204        vy[k + n/2] = y0[k] - w * y1[k];
205        w = w*wn;
206    }
207    
208    return vy;
209}
210
211 unsigned  int  rev(unsigned  int  num, unsigned  int  nBits)
212 {
213    unsigned int r = 0;
214    for(unsigned int i=0; i<nBits; i++)
215    {
216        if(num&(1<<i))
217            r |= 1<<(nBits-i-1);
218    }
219    
220    return r;
221}
222
223 vector < FuShu >  FFT_Iter(vector < FuShu >   const   & a, unsigned  int  nBits)
224 {
225    size_t n=a.size();
226    if(n==1) return a;
227    
228    vector<FuShu>  A;
229    A.reserve(n);
230    for(size_t i=0; i<n; i++)
231        A<<a[rev(i,nBits)];
232    
233    size_t m=2;
234    for(size_t s=0; s<nBits; s++, m*=2)
235    {
236        FuShu wm(cos(6.28/m), sin(6.28/m));
237        for(size_t k=0; k<n; k+=m)
238        {
239            FuShu w=1;
240            for(size_t j=0; j<m/2; j++, w*=wm)
241            {
242                FuShu t = w * A[k+j+m/2];
243                FuShu u = A[k+j];
244                A[k+j] = u+t;
245                A[k+j+m/2] = u-t;
246            }
247        }
248    }
249    
250    return A;
251}
252
253 /**/ /***************************************************/
254 //  Main
255 /**/ /***************************************************/
256
257 void  main()
258 {
259
260    srand(clock());
261    for(int i=0; i<10; i++)
262    {
263        FuShu a(rand(),rand()), b(rand(),rand());
264
265        vector<FuShu> vA;
266        vA<<a <<b<<a <<b;
267        
268        printf("input:\n");
269        Print(vA);
270        printf("FFT_recursive result:\n");
271        vector<FuShu> vB = FFT_recursive(vA);
272        Print(vB);
273        printf("DFT result:\n");
274        vector<FuShu> vBDft = DFT(vA);
275        Print(vBDft);
276        
277        printf("FFT_Iter result:\n");
278        vector<FuShu> vBItr = FFT_Iter(vA,2);
279        Print(vBItr);
280        
281        printf("diff: %g\n", fabs(vB - vBItr));
282        assert( fabs(vB - vBItr)<1e-1);
283    }
284#if 0        
285    for(unsigned int i=0; i<8; i++)
286        printf("%d ", rev(i,3));
287#endif        
288}