zoukankan html css js c++ java

快速傅里叶变换（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());
    }



}

####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());
    }

}

查看全文

相关阅读:
FTP 协议和 HTTP 协议的比较
 HttpURLConnection的post请求，什么时候发出，writeData存在什么地方
 装饰器
 函数参数以及名称空间作用域
 函数的调用
 函数的返回值
 定义函数的三种方式
 函数
 day05
day04

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