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
  6using namespace std;
  7
  8//************************************/
  9// Fushu
 10/**//**************************************/
 11class 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
 45FuShu operator+ (FuShu const & lhs, FuShu const & rhs)
 46{
 47    FuShu res = lhs;
 48    return res+=rhs;
 49}    
 50
 51FuShu operator- (FuShu const & lhs, FuShu const & rhs)
 52{
 53    FuShu res = lhs;
 54    return res-=rhs;
 55}        
 56
 57FuShu operator* (FuShu const & lhs, FuShu const & rhs)
 58{
 59    FuShu res = lhs;
 60    return res*=rhs;
 61}    
 62
 63double 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
 72void 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
 80void 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
 93template <typename T>
 94vector <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
105template <typename T>
106vector <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
118template <typename T>
119vector <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
129template <typename T>
130double 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
153template <typename T>
154vector <T> & operator<<(vector <T> & lhs, T const &rhs)
155{
156    lhs.push_back(rhs);
157    return lhs;
158}
159
160/**//***************************************************/
161// FFT
162/**//***************************************************/
163vector<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
180vector<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
211unsigned 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
223vector<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
257void 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}