快速傅里叶变换（FFT）

2017年12月09日 20:09:51 爱玲姐姐阅读数 2747更多

分类专栏： ACM 算法

本文链接：https://blog.csdn.net/jal517486222/article/details/78761188

快速傅里叶变换不是一种新的变换，而是傅里叶变换的一种快速算法，这个算法可以将普通的离散傅里叶变换（DFT）的时间复杂度O（n*n）降到O（n
log n），大大提高了傅里叶变换的速度。快速傅里叶变换算法的提出，使傅里叶变换在通信领域得到了极大地运用和发展。
在ACM中，快速傅里叶变换通常用于大数的乘法。当两个数大到连JAVA的BI都无法承受时，就该使用快速傅里叶算法了。该算法是将两个数，都写成多项式矩阵，然后对其进行离散傅里叶变换。通信领域有句很经典的话，叫做时域卷积，频域相乘。所以将两个数的多项式系数矩阵进行离散傅里叶变换后，就可以对其进行直接相乘了。将相乘后得结果再进行傅里叶逆变换（DTFT）。

######我之前怎么也无法理解快速傅里叶算法，因为当时觉得太复杂了太难了，看都看不进去。但是，由于我的专业是电子信息工程，所以无可避免地会在各个专业课上与傅里叶变换碰面。就在上一周，我们的数字信号处理老师，花了好几节课的时间，给我们讲解了快速傅里叶算法，我当时听得很认真，并自认为听懂了。但是，当我想用自己脑子里的快速傅里叶算法求解大数相乘的题目时，却发现无从下手。我根本就不知道怎么去运用它。当时我就在想，这个算法我必须掌握它，不仅仅是因为它会出现在ACM中，更是因为它在我的专业上也占据了极其重要的地位。如果我连傅里叶变换都不会，那以后还有什么脸说自己是电子信息专业的。于是乎，我花了两天的时间，查资料，终于搞懂了这个强大的算法。我发现有些大神写的关于FFT的文章确实非常不错，清晰易懂，我自己也收藏了，推荐给大家，大家一起学习。

超级易懂的傅里叶变换

快速傅里叶变换（FFT）

####JAVA代码

import java.util.Arrays;
import java.util.Scanner;
import java.util.Vector;

/**
 * Created by jal on 2017/12/6 0006.
 */
public class DIT_FFT {
    private static int maxn;
    public static void main(String[] args) {
        Scanner scanner = new Scanner(System.in);
        String numstr1 = scanner.next();
        String numstr2 = scanner.next();
        int len =numstr1.length()+numstr2.length();
        maxn = 0;
        double temp = log2(len);
        double floortemp = Math.floor(temp);
        if(floortemp == temp){
            maxn = len;
        }else {
            maxn = (int)Math.pow(2,floortemp+1);
        }
        Complex []arra = createArray(maxn);
        Complex []arrb =createArray(maxn);
        Complex []arrc =createArray(maxn);
        Complex []arrA = createArray(maxn);
        Complex []arrB = createArray(maxn);
        Complex []arrC = createArray(maxn);

        char []a = numstr1.toCharArray();
        int j = 0;
        for(int i = a.length - 1; i >= 0; i--){
            arra[j++] = new Complex(a[i] - '0');
        }
        j = 0;
        char []b = numstr2.toCharArray();
        for(int i = b.length - 1; i >= 0; i--){
            arrb[j++] = new Complex(b[i] - '0');
        }
        //System.out.println(Arrays.toString(arra));

        //invert(arra);



        arrA = fft(arra);
        System.out.println(Arrays.toString(arrA));
        //invert(arrb);

        arrB = fft(arrb);
        //System.out.println(Arrays.toString(arrB));
        for(int i = 0; i < arrC.length; i++){
            arrC[i] = arrA[i].times(arrB[i]);
        }
        //System.out.println(Arrays.toString(arrC));
        arrc = ifft(arrC);
        Vector<Integer> vector = new Vector<>();
        vector = toIntOfString(arrc);
        String str = "";
        char vectorCHar[] = new char[vector.size()];
        for(int i = 0; i < vectorCHar.length; i++){
            vectorCHar[i] = (char)(vector.get(i) + '0');
        }
        str =String.valueOf(vectorCHar);
        str = new StringBuffer(str).reverse().toString();
        //System.out.println("str:"+str);
        str = trim(str);
        System.out.println(str);

    }
    private static String trim(String value) {

        int len = value.length();
        int st = 0;
        char[] val = value.toCharArray();    /* avoid getfield opcode */

        while ((st < len) && (val[st] <= '0')) {
            st++;
        }
//            while ((st < len) && (val[len - 1] <= ' ')) {
//                len--;
//            }
        return value.substring(st, len);

    }
    private static Vector<Integer> toIntOfString(Complex[] arrc) {
        Vector<Integer> result = new Vector<>();
        int n = arrc.length;
        int arrayTemp [] = new int[n+1];
        for(int i = 0; i < arrc.length; i++){
            arrayTemp[i] = (int)arrc[i].re();
        }
        int i;
        arrayTemp[n] = 0;
        for(i = 0; i < n; i++){
            result.add(arrayTemp[i]%10);
            arrayTemp[i+1] += arrayTemp[i] / 10;
        }
        int t = arrayTemp[n];
        while(t > 0){
            result.add(t % 10);
            t/= 10;
        }
        return result;
    }

    private static Complex[] ifft(Complex[] arrC) {
        int n = arrC.length;
        Complex arrayResult[] = new Complex[n];
        for(int i = 0; i < arrC.length; i++){
            arrC[i] = arrC[i].conjugate();
        }
        arrayResult= fft(arrC);
        for(int i = 0; i < arrayResult.length; i++){
            arrayResult[i] = arrayResult[i].conjugate().divides(new Complex(n));
        }
        return arrayResult;
    }

    private static Complex[] fft(Complex[] arrA) {
        int N = arrA.length;
        if(N == 1){
            return arrA;
        }
        Complex [] arrayEven = new Complex[N/2];
        Complex [] arrayOdd = new Complex[N/2];
        for(int i = 0; i < arrayEven.length; i++){
            arrayEven[i] = arrA[2*i];
            arrayOdd[i] = arrA[2*i+1];
        }
        arrayEven = fft(arrayEven);
        arrayOdd = fft(arrayOdd);
        Complex[]arrayResult = new Complex[N];
        for(int i = 0; i < N/2; i++){
            Complex W = new Complex(0,-2*Math.PI*i/N).exp();
            arrayResult[i] = arrayEven[i].plus(arrayOdd[i].times(W));
            arrayResult[i+N/2] = arrayEven[i].minus(arrayOdd[i].times(W));
        }
        return arrayResult;
    }

    private static  void bit_reverse_swap(Complex [] a) {
        int n = a.length;
        for (int i = 1, j = n >> 1, k; i < n - 1; ++i) {
            if (i < j) swap(a,i,j);
            // tricky
            for (k = n >> 1; j >= k; j -= k, k >>= 1)  // inspect the highest "1"
                ;
            j += k;
        }
    }

    private static void invert(Complex[] arra) {
        int n = (int)log2(maxn);
        for(int i = 0; i < arra.length; i++){
            String temp = Integer.toBinaryString(i);
            temp = new StringBuffer(temp).reverse().toString();
            int len = n - temp.length();
            char [] zeros = new char[len];
            Arrays.fill(zeros,'0');
            temp+=String.valueOf(zeros);

            int j = Integer.valueOf(temp,2);
            System.out.println("i:" + i +"j:" + j);
            if(j>i){
                swap(arra,i,j);
            }
        }
    }

    private static void swap(Complex[] arra, int i, int j) {
        Complex tmp = arra[i];
        arra[i] = arra[j];
        arra[j] = tmp;
    }
    private static Complex[] createArray(int maxn) {
        Complex[] result =new Complex[maxn];
        for(int i=0;i<maxn;i++)
            result[i]=new Complex(0,0);
        return result;
    }

    public static double log2(int n){
        return Math.log(n)/Math.log(2);
    }
}


/******************************************************************************
 *  Compilation:  javac Complex.java
 *  Execution:    java Complex
 *  Dependencies: StdOut.java
 *
 *  Data type for complex numbers.
 *
 *  The data type is "immutable" so once you create and initialize
 *  a Complex object, you cannot change it. The "final" keyword
 *  when declaring re and im enforces this rule, making it a
 *  compile-time error to change the .re or .im fields after
 *  they've been initialized.
 *
 *  % java Complex
 *  a            = 5.0 + 6.0i
 *  b            = -3.0 + 4.0i
 *  Re(a)        = 5.0
 *  Im(a)        = 6.0
 *  b + a        = 2.0 + 10.0i
 *  a - b        = 8.0 + 2.0i
 *  a * b        = -39.0 + 2.0i
 *  b * a        = -39.0 + 2.0i
 *  a / b        = 0.36 - 1.52i
 *  (a / b) * b  = 5.0 + 6.0i
 *  conj(a)      = 5.0 - 6.0i
 *  |a|          = 7.810249675906654
 *  tan(a)       = -6.685231390246571E-6 + 1.0000103108981198i
 *
 ******************************************************************************/

/**
 *  The {@code Complex} class represents a complex number.
 *  Complex numbers are immutable: their values cannot be changed after they
 *  are created.
 *  It includes methods for addition, subtraction, multiplication, division,
 *  conjugation, and other common functions on complex numbers.
 *  <p>
 *  For additional documentation, see <a href="https://algs4.cs.princeton.edu/99scientific">Section 9.9</a> of
 *  <i>Algorithms, 4th Edition</i> by Robert Sedgewick and Kevin Wayne.
 *
 *  @author Robert Sedgewick
 *  @author Kevin Wayne
 */
class Complex {
    private double re;   // the real part
    private  double im;   // the imaginary part

    /**
     * Initializes a complex number from the specified real and imaginary parts.
     *
     * @param real the real part
     * @param imag the imaginary part
     */
    public Complex(double real, double imag) {
        re = real;
        im = imag;
    }

    public Complex(double re) {
        this.re = re;
        this.im = 0;
    }

    /**
     * Returns a string representation of this complex number.
     *
     * @return a string representation of this complex number,
     *         of the form 34 - 56i.
     */
    public String toString() {
        if (im == 0) return re + "";
        if (re == 0) return im + "i";
        if (im <  0) return re + " - " + (-im) + "i";
        return re + " + " + im + "i";
    }

    /**
     * Returns the absolute value of this complex number.
     * This quantity is also known as the <em>modulus</em> or <em>magnitude</em>.
     *
     * @return the absolute value of this complex number
     */
    public double abs() {
        return Math.hypot(re, im);
    }

    /**
     * Returns the phase of this complex number.
     * This quantity is also known as the <em>angle</em> or <em>argument</em>.
     *
     * @return the phase of this complex number, a real number between -pi and pi
     */
    public double phase() {
        return Math.atan2(im, re);
    }

    /**
     * Returns the sum of this complex number and the specified complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this + that)}
     */
    public Complex plus(Complex that) {
        double real = this.re + that.re;
        double imag = this.im + that.im;
        return new Complex(real, imag);
    }

    /**
     * Returns the result of subtracting the specified complex number from
     * this complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this - that)}
     */
    public Complex minus(Complex that) {
        double real = this.re - that.re;
        double imag = this.im - that.im;
        return new Complex(real, imag);
    }

    /**
     * Returns the product of this complex number and the specified complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this * that)}
     */
    public Complex times(Complex that) {
        double real = this.re * that.re - this.im * that.im;
        double imag = this.re * that.im + this.im * that.re;
        return new Complex(real, imag);
    }

    /**
     * Returns the product of this complex number and the specified scalar.
     *
     * @param  alpha the scalar
     * @return the complex number whose value is {@code (alpha * this)}
     */
    public Complex scale(double alpha) {
        return new Complex(alpha * re, alpha * im);
    }

    /**
     * Returns the product of this complex number and the specified scalar.
     *
     * @param  alpha the scalar
     * @return the complex number whose value is {@code (alpha * this)}
     * @deprecated Replaced by {@link #scale(double)}.
     */
    @Deprecated
    public Complex times(double alpha) {
        return new Complex(alpha * re, alpha * im);
    }

    /**
     * Returns the complex conjugate of this complex number.
     *
     * @return the complex conjugate of this complex number
     */
    public Complex conjugate() {
        return new Complex(re, -im);
    }

    /**
     * Returns the reciprocal of this complex number.
     *
     * @return the complex number whose value is {@code (1 / this)}
     */
    public Complex reciprocal() {
        double scale = re*re + im*im;
        return new Complex(re / scale, -im / scale);
    }

    /**
     * Returns the real part of this complex number.
     *
     * @return the real part of this complex number
     */
    public double re() {
        return re;
    }

    /**
     * Returns the imaginary part of this complex number.
     *
     * @return the imaginary part of this complex number
     */
    public double im() {
        return im;
    }

    /**
     * Returns the result of dividing the specified complex number into
     * this complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this / that)}
     */
    public Complex divides(Complex that) {
        return this.times(that.reciprocal());
    }

    /**
     * Returns the complex exponential of this complex number.
     *
     * @return the complex exponential of this complex number
     */
    public Complex exp() {
        return new Complex(Math.exp(re) * Math.cos(im), Math.exp(re) * Math.sin(im));
    }

    /**
     * Returns the complex sine of this complex number.
     *
     * @return the complex sine of this complex number
     */
    public Complex sin() {
        return new Complex(Math.sin(re) * Math.cosh(im), Math.cos(re) * Math.sinh(im));
    }

    /**
     * Returns the complex cosine of this complex number.
     *
     * @return the complex cosine of this complex number
     */
    public Complex cos() {
        return new Complex(Math.cos(re) * Math.cosh(im), -Math.sin(re) * Math.sinh(im));
    }

    /**
     * Returns the complex tangent of this complex number.
     *
     * @return the complex tangent of this complex number
     */
    public Complex tan() {
        return sin().divides(cos());
    }



}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430

####C++代码

#include <complex>
#include <vector>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <string>
#include <cstdio>
using namespace std;
double log2(int n){
    return log(n)/log(2);
}
const int MAXN = 1 << 20;
complex<double> arra[MAXN],arrb[MAXN],arrC[MAXN];
complex<double>* arrA;
complex<double>* arrB;//,arrC[MAXN];
complex<double>* arrc;
const double PI = 3.141592653;
complex<double>* DIT_FFT(complex<double>* arr,int len){
    if(len == 1)return arr;
    complex<double> *arrayEven = new complex<double>[len/2];
    complex<double> *arrayOdd = new complex<double>[len/2];
    for(int i = 0; i < len/2; i++){
        arrayEven[i] = arr[2*i];
        arrayOdd[i] = arr[2*i+1];
    }
    arrayEven = DIT_FFT(arrayEven,len/2);
    arrayOdd = DIT_FFT(arrayOdd,len/2);
    complex<double> *arrayResult = new complex<double>[len];
    for(int i = 0; i < len/2; i++){
        //TODO
        //help me
        //help me
        // what is 'W' ?
        complex<double> W = exp(complex<double>(0,-2*PI *i / len));//ecomplex<double>(cos(-2*PI *i / len),sin(-2*PI *i / len));
        cout<<"W:"<<W<<endl;
        arrayResult[i] = arrayEven[i] + arrayOdd[i] * W;
        arrayOdd[i + len/2] = arrayEven[i] - arrayOdd[i] * W;
    }
    return arrayResult;
}

complex<double>* IFFT(complex<double>*arrC,int len){
    int n = len;
    complex<double>* arrayResult = new complex<double>[len];
    for(int i = 0; i < len; i++){
        arrC[i] = conj(arrC[i]);
    }
    arrayResult= DIT_FFT(arrC,len);
    for(int i = 0; i < len; i++){
        arrayResult[i] = conj(arrayResult[i]);
        arrayResult[i]/=n;
    }
    return arrayResult;

}
int main(){

    string str1,str2;
    cin>>str1>>str2;
    int len = str1.size() + str2.size();
    int maxn = 0;
    double floortemp =log2(len);
    if(floortemp == floor(floortemp)){
        maxn = len;
    }
    else{
        maxn = pow(2,(int)(floor(floortemp) + 1) );
    }

    for(int i = str1.size() - 1; i >= 0; i--){
        arra[i] = complex<double>(str1[i]-'0',0);
    }
    for(int i = str1.size() - 1; i >= 0; i--){
        arrb[i] = complex<double>(str2[i]-'0',0);
    }
    arrA = DIT_FFT(arra,maxn);
    arrB = DIT_FFT(arrb,maxn);
    for(int i= 0; i < maxn; i++){
        arrC[i] = arrA[i]*arrB[i];
    }
    arrc = IFFT(arrC,len);
    vector<int>v;
    for(int i = 0; i < maxn;i++){
        v.push_back((int)arrc[i].real());
    }

}

相关阅读:
背包九讲
最小生成树 prime + 队列优化
最小生成树 prime poj1287
树状数组--转载
O(n)求1-n的逆元
数据结构--线段树
博弈
RMQ 数据结构
BZOJ3687 计算子集和的异或和
Intel Code Challenge Final Round (Div. 1 + Div. 2, Combined)D Dense Subsequence

原文地址：https://www.cnblogs.com/grj001/p/12223990.html