Java编程实现快速傅里叶变换FFT

1 快速傅里叶变换FFT

1.1 理论分析

1.1.1 离散傅里叶变换

有限长序列x(n)的N点DFT为

考虑为x(n)复数序列的一般情况，对某一个k值，直接按照式(1)计算X(k)的1个值需要N次复数乘法和(N-1)次复数加法。因此，计算X(k)的所有N个值，共需N^2次复数乘法和N(N-1)次复数加法。当N>>1时， N(N-1)约等于N^2。由上述可见，N点DFT的乘法和加法运算次数均为 N*N，也即DFT的时间复杂度为O(N^2)。

1.1.2 快速傅里叶变换

由于DFT时间复杂度较高，在处理较大计算量时花费时间较长，所以JWCooley和JohnTukey在1965年的文章中提出了Cooley-TukeyFFT算法，FFT算法就是不断把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数。其周期性表现为：

其对称性表现为：

设序列 的长度为N，且满足N=2的M次方，M为自然数。经过第一级分解后N点DFT被分解为两个点DFT X1(k)和X2(k)以及式(4)和式(5)的运算。

经过上述一次分解后，使运算量减少将近一半，由于N=2的M次方，N/2仍然是偶数，故可以进行第二次分解，进一步将N/2 DFT分解为2个N/4点DFT。以此类推，经过M次分解，最后将N点DFT分解成N个1点DFT和M级蝶形运算。故FFT的运算流图应有M级蝶形，每一级都由N/2个蝶形运算构成。并且完成一个蝶形运算，需要一次复数乘法和两次复数加法运算。因此每一级运算都需要N/2次复数乘和N次复数加。所以，M级运算总的复数乘次数为 Cm=(N/2)M=(N/2)log2N，复数加次数为Ca=NM=Nlog2N ，也即FFT的时间复杂度为O(Nlog2N)。

1.2 编程实现

本次实验DFT算法和FFT算法实现是通过Java语言，然后通过Matlab进行验证。FFT本质上是采用分治算法来计算的DFT。一个初始长度为N的初始数据，最后被分治成很多长度为2的数据，分别对这些长度为2的数据做DFT，最终等价整个N长度的数据做傅里叶变换。分治是一种思想，具体实现其实就是使用递归。

1.2.1 算法思想

FFT运算过程中涉及到复数的运算，公式(4)和(5)包含了复数的加法和乘法。因此在Java程序中需要建立一个Complex复数类，除了real实部和imag虚部两个属性外，还需要给Complex类定义plus方法和multiple方法，为了复数类的完整性，本次实验将加减乘除四个方法全部定义。此时，Complex类不仅可以用作复数对象存储实部和虚部，还能实现复数的运算。
Complex类的定义如下：

public class Complex {
    private double a, b;
    public Complex(){
        this.a = 0.0;
        this.b = 0.0;
    }
    public Complex(double a, double b){
        this.a = a;
        this.b = b;
    }
    public Complex(Complex C){
        this.a = C.a;
        this.b = C.b;
    }
     // 获取实部   
    public double getRealPart(){ 
        return a;
    } 
    //获取虚部
    public double getImaginaryPart(){
        return b;
    } 
   // 字符串表示复数
    public String toString(){
        if(b != 0) {
            if (a == 0)
                return b + "i";
            else
                return a + "+" + b + "i";
        }
        else
            return a + "";
    }  
    // 加法
    public Complex plus(Complex C){
        if (C == null){
            System.out.println("对象为null");
            return new Complex();
        }
        return new Complex(this.a + C.a, this.b + C.b);
    } 
    // 减法
    public Complex minus(Complex C){
        if (C == null){
            System.out.println("对象为null");
            return new Complex();
        }
        return new Complex(this.a - C.a, this.b - C.b);
    }
    // 乘法
    public Complex multiple(Complex C){
        if (C == null){
            System.out.println("对象为null");
            return new Complex();
        }
        double newa = this.a * C.a - this.b * C.b;
        double newb = this.a * C.b + this.b * C.a;
        return new Complex(newa, newb);
    } 
    // 除法
    public Complex divide(Complex C){
        if (C == null){
            System.out.println("对象为null");
            return new Complex();
        }
        if(C.a == 0 && b == 0){
            System.out.println("除数不能为0！");
            return new Complex();
        }
        double newa = (this.a * C.a + this.b * C.b) / (C.a * C.a + C.b * C.b);
        double newb = (this.a * C.a - this.b * C.b) / (C.a * C.a + C.b * C.b);
        return new Complex(newa, newb);
    } 
}

接下来，根据傅里叶变换的定义实现DFT算法。首先判断N是否为1，若为1直接递归到原点，否则按照公式(1)进行双重N次循环。因此其时间复杂度为O(N^2)。
DFT函数定义如下：

public static Complex[] dft(Complex[] x) {
         int N = x.length;

            // exp(-2i*n*k*PI)=cos(-2*n*k*PI)+i*sin(-2*n*k*PI)=1
            if (N == 1)
               return new Complex[]{x[0]};

            Complex[] result = new Complex[N];
            for (int k = 0; k < N; k++) {
                result[k] = new Complex(0, 0);
                for (int n = 0; n < N; n++) {
              //使用欧拉公式e^(-i*2pi*k*n/N) = cos(-2pi*k*n/N) + i*sin(-2pi*k*n/N)
                    double p = -2 * k * n* Math.PI / N;
                    Complex m = new Complex(Math.cos(p), Math.sin(p));
                    result[k] = result[k].plus(x[n].multiple(m));
                }
            }
            return result;
           }

FFT函数与DFT类似，最开始也是先进行N值判断，若为1则返回原点x[0]的值。然后在判断信号数N是否为奇数，若为奇数，则无法使用FFT，转而使用DFT算法；若为偶信号数，则提取下标为偶数和下标为奇数的原始信号分别进行递归FFT计算，然后根据公式(4)和(5)进行N/2次蝶形计算。综合来看，在FFT函数中递归次数为log2N，循环次数为N/2，因此FFT算法时间复杂度为O(Nlog2N)。
FFT函数定义如下：

 public static Complex[] fft(Complex[] x) {
        int N = x.length;

        // 因为exp(-2i*n*k*PI)=1，N=1时递归原点
        if (N == 1){
            return new Complex[]{x[0]};
        }

        // 如果信号数为奇数，使用dft计算
        if (N % 2 != 0) {
            return dft(x);
        }

        // 提取下标为偶数的原始信号值进行递归fft计算
        Complex[] even = new Complex[N / 2];
        for (int k = 0; k < N / 2; k++) {
            even[k] = x[2 * k];
        }
        Complex[] evenValue = fft(even);

        // 提取下标为奇数的原始信号值进行fft计算
        // 节约内存
        Complex[] odd = even;
        for (int k = 0; k < N / 2; k++) {
            odd[k] = x[2 * k + 1];
        }
        Complex[] oddValue = fft(odd);

        // 偶数+奇数
        Complex[] result = new Complex[N];
        for (int k = 0; k < N / 2; k++) {
         //使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)
            double p = -2 * k * Math.PI / N;
            Complex m = new Complex(Math.cos(p), Math.sin(p));
            result[k] = evenValue[k].plus(m.multiple(oddValue[k]));
          //exp(-2*(k+n/2)*PI/n) 相当于 -exp(-2*k*PI/n)，其中exp(-n*PI)=-1(欧拉公式);
            result[k + N / 2] = evenValue[k].minus(m.multiple(oddValue[k]));
        }
        return result;
       }

1.2.2 实验结果

为了更好地对比DFT和FFT算法的时间复杂度，本次实验将产生幅值为0.7的100 Hz正弦量和幅值为1的220 Hz正弦量（一共1500个信号量），同时获取DFT和FFT运算时间。另外为了方便后续使用Matlab验证算法的正确性，本次实验将产生的FFT数据写入到文件中，然后在Matlab中调用FFT函数并画出频谱图进行对比。
产生信号函数如下：

public static void CreateSignal(int length, Complex[] x){
        double Fs = 1000; // 采样频率
        double T = 1 / Fs; // 采样周期
        System.out.println("产生幅值为 0.7 的 100 Hz正弦量和幅值为 1 的 220 Hz正弦量(一共"+length+"个)");
        for (int i = 0; i < length; i++) {
            // 产生幅值为0.7的100Hz正弦量和幅值为1的220Hz正弦量。
            double S = 0.7 * Math.sin(2 * Math.PI * 100 * i * T) + Math.sin(2 * Math.PI * 220 * i * T);
            x[i] = new Complex(S, 0);
        }
    }

生成txt文件的函数如下：

public static void writeText(String filename,Complex[]result) {
        try {
            File writeName = new File(filename); // 相对路径，如果没有则要建立一个新文件
            writeName.createNewFile(); // 创建新文件,有同名的文件的话直接覆盖
            try (FileWriter writer = new FileWriter(writeName);
                 BufferedWriter out = new BufferedWriter(writer)
            ) {
                for (int i = 0; i < result.length; i++) {
                    out.write(result[i].toString()+"\r\n"); // 即为换行
                }
            }
        } catch (IOException e) {
            e.printStackTrace();
        }
        System.out.println("数据写入成功！");
    }

主函数如下：

public static void main(String[] args) {
        int length = 1500;
        Complex[] x = new Complex[length];
        Complex[] result;
        CreateSignal(length, x);

        long start1 = System.currentTimeMillis();
        result = dft(x);
        long end1 = System.currentTimeMillis();

        long start2 = System.currentTimeMillis();
        result = fft(x);
        long end2 = System.currentTimeMillis();

        System.out.println("FFT结果：");
        for (int i = 0; i < length; i++) {
           System.out.println("第"+(i+1)+"个数据："+result[i].toString());
        }
        System.out.println("-----------------------------------------");

        System.out.println("DFT运行时间: " + (end1 - start1) + "ms");
        System.out.println("FFT运行时间: " + (end2 - start2) + "ms");

        String filename = "D:\\研一课程\\数据结构\\FFT.txt";
        writeText(filename, result);
    }

Matlab程序如下：

clc
clear
close all;
Fs = 1000;                   
T = 1/Fs;             
L = 1500;          
t = (0:L-1)*T; 
S = 0.7*sin(2*pi*100*t) + sin(2*pi*220*t);

%matlab fft函数
Y = fft(S);
%导入java产生的数据
file = fopen('FFT.txt');
p = textscan(file,'%s');
pp =[str2double(p{:})]';

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(1,2,1)
plot(f,P1) 
hold on;
title('FFT of Matlab')
xlabel('f (Hz)')
ylabel('|P1(f)|')

X2 = abs(pp/L);
X1 = X2(1:L/2+1);
X1(2:end-1) = 2*X1(2:end-1);
f1 = Fs*(0:(L/2))/L;
subplot(1,2,2);
plot(f1,X1) 
title('FFT of Java')
xlabel('f (Hz)')
ylabel('|X1(f)|')
hold off;

将1500个信号量为幅值为0.7的100 Hz正弦量和幅值为1的220 Hz正弦量FFT后的结果如图1所示。

FFT运行时间相较于DFT提升了5倍左右，如图2所示。

同样的，与Matlab中自带的FFT算法进行对比，得到的频谱图完全一致！算法验证成功！如图3所示。