使用java实现一维FFT

傅里叶变换基础

FFT是快速傅里叶变换，它是DFT的快速算法，但是计算结果与DFT等价。首先回顾一下傅里叶变换的计算公式：

首先，非周期性连续时间信号x(t)的傅里叶变换可以表示为：（公式1）

在离散的情况下，傅里叶变换的公式为：（公式2）

上面就是DFT的计算公式。DFT用代码实现是很简单的，因为，对于任意指定数值的K，只需要遍历N个采样点即可得到该K值下的频谱。不过DFT的复杂度为O(N^2)，在N较大的情况下，计算速度是很慢的。

FFT

FFT本身是DFT的变体，所以我们需要先了解DFT的理论。下面先简单介绍FFT的原理：
还是对x(n)信号求傅里叶变换。x(n)本身可以拆分成下标为奇数子集和偶数子集。

x1(n)和x2(n)的长度都是N／2，x1(n)是偶数序列，x2(n)是奇数序列，则：（公式3）

其中X1(k)和X2(k)分别为x1(n)和x2(n)的N／2点DFT。由于X1(k)和X2(k)均以N／2为周期，且WN k+N/2=-WN k，所以X(k)又可表示为：（公式4）

写到这里就不难发现，FFT本质上是采用分治算法来计算的DFT。一个初始长度为N的初始数据，最后会被分治为很多长度为2 的数据，然后分别对这些长度为2的数据做DFT，最终等价于对整个N长度的数据做傅里叶变换。
分治其实是一种思想，具体的实现其实就是使用递归。在写Java程序之前，我们需要做一些准备工作。FFT的运算过程中涉及到复数的运算。公式4中包含了复数的加法和乘法。因此我们在Java程序中需要建立一个Complex（复数）类，除了real（实部）和imag（虚部）两个属性以外，需要给Complex类定义add方法和multiply方法。这样的话，Complex类不仅仅可以用作复数对象存储实部和虚部，还能实现复数的运算。Java中不能实现运算符的重载，这是挺遗憾的事情，这样的话就不能直接用加减乘除号来进行Complex对象之间的运算了。
以下是Complex类的定义：

// Struct to hold two primary type
    public static class Complex2 {
        public double real;
        public double imag;

        Complex2() {
            this.real = 0;
            this.imag = 0;
        }

        Complex2(double real, double imag) {
            this.real = real;
            this.imag = imag;
        }

        Complex2 add(Complex2 another) {
            return new Complex2(this.real + another.real, this.imag + another.imag);
        }

        Complex2 scalarMultiply(Complex2 another) {

            double x1 = this.real;
            double x2 = another.real;
            double y1 = this.imag;
            double y2 = another.imag;

            return new Complex2(x1*x2-y1*y2, x1*y2+x2*y1);
        }
    }

复数的加法是实部加实部，虚部加虚部。复数的乘法是类似于多项式相乘。正如代码所示。
下面进入正题。实现FFT主函数。FFT主函数必须分成两个函数来写，因为其中一个是递归函数，需要分离出来才能做递归。代码如下：

public static Complex2[] fft(List<Double> dataList) {

        // Get the size of dataList
        final int N = dataList.size();

        Complex2[] fftResult = new Complex2[N];

        // Traverse dataList by index k from 0 to size of dataList
        for(int k=0 ; k<N ; k++) {
            fftResult[k] = fftRecursion(dataList, N, k);
        }

        return fftResult;
    }

这是fft的主函数，主函数中主要是做了一个循环，因为我们需要求K频域上K从0~N-1变化过程中的频谱。对于每一个不同的K值，我们分别调用fftRecursion函数来实现分治算法。下面列出fftRecursion函数的实现：

/**
     * Realize FFT by divide-and-conquer algorithm
     *
     * @param ALineUnit This is List of fragemnt of ALine Points. The termination of recursion
     *                  is sometime when size of ALineUnit is 2
     */
    private static Complex2 fftRecursion(List<Double> ALineUnit, int unitLength, int k) {

        if(unitLength == 2) {

            Complex2 unit2 = new Complex2();

            // Apply FFT to data list (size:2)
            double xn;

            // Get real value
            for(int k1=0 ; k1<2 ; k1++) {
                for(int n=0 ; n<2 ; n++) {

                    xn = ALineUnit.get(n);

                    unit2.real += xn * cos(2*PI*k1*n/2);
                    unit2.imag -= xn * sin(2*PI*k1*n/2);
                }
            }

            return unit2;

        }

        else {
            // Cut odd and even component of the current unit
            List<Double> oddList = new ArrayList<>();
            List<Double> evenList = new ArrayList<>();

            for(int i=0 ; i<unitLength/2 ; i++) {
                oddList.add(ALineUnit.get(2*i + 1));
                evenList.add(ALineUnit.get(2*i));
            }

            return fftRecursion(evenList, unitLength/2, k).add(new Complex2(cos(2*PI*k/unitLength), -sin(2*PI*k/unitLength)).scalarMultiply(fftRecursion(oddList, unitLength/2, k)));
        }

    }

需要给递归设置一个终止递归的条件。也就是说当递归进行到每个数据包的长度已经变成2的时候，就可以对每个数据包做DFT，最后整体实现一个FFT的效果。
在递归还没有进行到数据包长度等于2之前，我们需要依照公式4的第一个公式对当前的数据包进行拆分，拆成奇数项和偶数项，然后对照公式4中的第一个公式，把k代进去（k在递归的过程中是个常数），同时做复数的加法和乘法（这个是我们在Complex类中就已经实现了的）。
以上。