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";elsereturn a + "+" + b + "i";}elsereturn 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)=1if (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所示。

Java编程实现快速傅里叶变换FFT相关推荐

Java中实现快速傅里叶变换FFT
Java中实现快速傅里叶变换FFT 一.概述 1.傅里叶变换(FT) 2.离散傅里叶变换(DFT) 3.快速傅里叶变换(FFT) 1)单位根 2)快速傅里叶变换的思想 3)蝶形图 4)快速傅里叶变换的 ...
java fft_java实现快速傅里叶变换（FFT）
最近做音频信号处理的时候,需要对数据做fft变换.关于fft原理: 请参考:FFT算法讲解--麻麻我终于会FFT了! matlab实现fft函数很简单,直接调用fft即可.但java实现起来就有点难了 ...
基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...
基于python的快速傅里叶变换FFT（一）
基于python的快速傅里叶变换FFT(一) FFT可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了.这就是很多信号分析采用FFT变换的原因. ...
MIT 线性代数 Linear Algebra 26:复矩阵，傅里叶矩阵, 快速傅里叶变换 FFT
这一讲我们来讲一下复矩阵.线性代数中,复矩阵是避免不了的话题,因为一个简单实矩阵都有可能有复数特征值. 复矩阵我们着重看一下复矩阵和实矩阵在运算上的区别. 距离首先,一个复数向量的的距离求法发生了 ...
OpenCV快速傅里叶变换(FFT)用于图像和视讯流的模糊检测
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测翻译自[OpenCV Fast Fourier Transform (FFT) for blur detection in images ...
Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析
文章目录 1. 定义 2. 变换和处理 3. 函数 4. 实例演示例1:单频正弦信号(整数周期采样) 例2:单频正弦信号(非整数周期采样) 例3:含有直流分量的单频正弦信号例4:正弦复合信号例5 ...
快速傅里叶变换FFT进行频谱分析（matlab）
快速傅里叶变换FFT进行频谱分析(matlab) 本章摘要:FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了 ...
快速傅里叶变换FFT C语言实现可用于嵌入式系统进行模拟采样频谱分析
快速傅里叶变换C语言实现模拟采样进行频谱分析 FFT是DFT的快速算法用于分析确定信号(时间连续可积信号.不一定是周期信号)的频率(或相位.此处不研究相位)成分,且傅里叶变换对应的 ω \omega ...

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

快速傅里叶变换的时间复杂度分析