Java中实现快速傅里叶变换FFT

  • 一.概述
    • 1.傅里叶变换(FT)
    • 2.离散傅里叶变换(DFT)
    • 3.快速傅里叶变换(FFT)
      • 1）单位根
      • 2）快速傅里叶变换的思想
      • 3）蝶形图
      • 4）快速傅里叶变换的逆变换(IFFT)
  • 二.代码实现
    • 1.复数的实现
      • 1）复数的表示
      • 2）复数的运算
      • 3）复数数组
    • 2.单位根的实现
      • 1）单位根的计算
      • 2）获取单位根时的处理
    • 3.快速傅里叶变换的实现
      • 1）索引插值
      • 2）数据奇偶分组
      • 3）确定循环次数
      • 4）核心计算
  • 三.获取全部代码

一.概述

1.傅里叶变换(FT)

傅里叶变换是将一个满足条件的函数表示成三角函数或其积分的线性组合。
      实际应用中，将一个随时间变化的信号，经过傅里叶变换，可以得到这个信号的频率组成，实现从时域到频域的变换，进而研究其频谱结构和变化规律。
      傅里叶变换表明：一个信号可以表示成多个不同频率的正弦波的叠加。

2.离散傅里叶变换(DFT)

对一个连续的时域信号以一定的频率进行采样，可以得到这个信号时域的一组采样点，通过对这些采样点进行离散傅里叶变换，可以得到这个信号对应频域的一组采样点。

由公式可以得出离散傅里叶变换的时间复杂度为O(n2)。

3.快速傅里叶变换(FFT)

快速傅里叶变换是对离散傅里叶变换进行快速计算的一种方法。可以将时间复杂度降低到O(nlog2n)，因此称为快速傅里叶变换。

1）单位根

在一个复平面内将单位圆平均分成n份，这n个点表示的复数称为n次单位根。

其中k = 0，1，2，…n-1，方向逆时针。

单位根的性质：

2）快速傅里叶变换的思想

将DFT公式的求和符号去掉，变成相加的多项式A。用单位根替换旋转因子，将展开式中的项分成奇数项展开式A1和偶数项展开式A0，通过利用单位根的性质，可以将奇数项展开式A1和偶数项展开式A0化简为相同的单位根和x(n)相乘的形式，通过将奇数项展开式A1和偶数项展开式A0中的x(n)分别和单位根相乘，最后再将两部分相加，可以求出X(K)的值。
      进一步，可以将奇数项展开式A1再分成奇数项展开式A11和偶数项展开式A10，偶数项展开式A0也再分成奇数项展开式A01和偶数项展开式A00，一直向下分，最后分成只包含两项的展开式，使问题的规模变小，最后向上计算，求出X(K)。
      因此快速傅里叶变换的时间复杂度为O(nlog2n)。
eg：设对于某个信号的一组16点采样值为 0，1，2，3，4，5，6，7，8，9，10，11，12，13，14，15。
第一次奇偶分组：0，2，4，6，8，10，12，14 | 1，3，5，7，9，11，13，15
第二次奇偶分组：0，4，8，12 | 2，6，10，14 | 1，5，9，13 | 3，7，11，15
第三次奇偶分组：0，8 | 4，12 | 2，10 | 6，14 | 1，9 | 5，13 | 3，11 | 7，15

设对于某个信号的一组8点采样值为 0，1，2，3，4，5，6，7，8。
第一次奇偶分组：0，2，4，6 | 1，3，5，7
第二次奇偶分组：0，4 | 2，6 | 1，5 | 3，7

通过观察可以发现，对采样点对应位置的索引的二进制形式进行镜像操作，再还原成十进制形式，可以获得采样点在奇偶分组后的位置。

原位置 二进制表示 镜像表示 分组后位置
0 000 000 0
1 001 100 4
2 010 010 2
3 011 110 6
4 100 001 1
5 101 101 5
6 110 011 3
7 111 111 7

3）蝶形图

将奇数项和偶数项的展开式从下向上计算的过程，可以通过画图来表示，由于该图形似蝴蝶，故称为蝶形图。

快速傅里叶变换：

4）快速傅里叶变换的逆变换(IFFT)

对信号频域的采样点进行快速傅里叶变换的逆变换，可以得到信号在时域对应的采样点。
      具体实现就是快速傅里叶变换蝶形计算的逆过程，但是其单位根为快速傅里叶变换单位根的共轭形式，同时，每次计算后结果要除以2。
快速傅里叶变换的逆变换：

二.代码实现

1.复数的实现

1）复数的表示

private float realPart;
private float imaginaryPart;public ComplexNumber(float realPart, float imaginaryPart) {this.realPart = realPart;this.imaginaryPart = imaginaryPart;
}

2）复数的运算

复数的运算包括加、减、乘、除和共轭运算。在实际应用中，两个复数进行运算后的结果可以存储在调用的对象中，但为了防止数据污染，也可能需要存储在一个新的对象中。其它的计算包括计算复数的幅值(模)和相位。

存储在新对象中：

//加法
public ComplexNumber addNew(ComplexNumber addend) {if(addend == null) return new ComplexNumber(this.realPart, this.imaginaryPart);return new ComplexNumber(this.realPart+addend.realPart, this.imaginaryPart+addend.imaginaryPart);
}

//减法
public ComplexNumber subtractNew(ComplexNumber subtrahend) {if(subtrahend == null) return new ComplexNumber(this.realPart, this.imaginaryPart);return new ComplexNumber(this.realPart-subtrahend.realPart, this.imaginaryPart-subtrahend.imaginaryPart);
}

//乘法
public ComplexNumber multiplyNew(ComplexNumber multiplier) {if(multiplier == null) return new ComplexNumber(this.realPart, this.imaginaryPart);float newRealPart = this.realPart*multiplier.realPart - this.imaginaryPart*multiplier.imaginaryPart;float newImaginaryPart = this.realPart*multiplier.imaginaryPart + this.imaginaryPart*multiplier.realPart;return new ComplexNumber(newRealPart, newImaginaryPart);
}

//除法
public ComplexNumber divideNew(ComplexNumber divisor) {if(divisor == null) return new ComplexNumber(this.realPart, this.imaginaryPart);float sumBase = (float) (Math.pow(divisor.realPart, 2) + Math.pow(divisor.imaginaryPart, 2));float newRealPart = (this.realPart*divisor.realPart + this.imaginaryPart*divisor.imaginaryPart)/sumBase;float newImaginaryPart = (this.realPart*divisor.imaginaryPart - this.imaginaryPart*divisor.realPart)/sumBase;return new ComplexNumber(newRealPart, newImaginaryPart);
}

//共轭
public ComplexNumber conjugateNew() {return new ComplexNumber(this.realPart, -this.imaginaryPart);
}

存储在调用的对象中：

//加法
public ComplexNumber add(ComplexNumber addend) {if(addend == null)return this;this.realPart += addend.realPart;this.imaginaryPart += addend.imaginaryPart;return this;
}

//减法
public ComplexNumber subtract(ComplexNumber subtrahend) {if(subtrahend == null)return this;this.realPart -= subtrahend.realPart;this.imaginaryPart -= subtrahend.imaginaryPart;return this;
}

//乘法
public ComplexNumber multiply(ComplexNumber multiplier) {if(multiplier == null) return this;float newRealPart = this.realPart*multiplier.realPart - this.imaginaryPart*multiplier.imaginaryPart;float newImaginaryPart = this.realPart*multiplier.imaginaryPart + this.imaginaryPart*multiplier.realPart;this.realPart = newRealPart;this.imaginaryPart = newImaginaryPart;return this;
}

//除法
public ComplexNumber divide(ComplexNumber divisor) {if(divisor == null) return this;float sumBase = (float) (Math.pow(divisor.realPart, 2) + Math.pow(divisor.imaginaryPart, 2));float newRealPart = (this.realPart*divisor.realPart + this.imaginaryPart*divisor.imaginaryPart)/sumBase;float newImaginaryPart = (this.imaginaryPart*divisor.realPart - this.realPart*divisor.imaginaryPart)/sumBase;this.realPart = newRealPart;this.imaginaryPart = newImaginaryPart;return this;
}

//共轭
public ComplexNumber conjugate() {this.imaginaryPart = -this.imaginaryPart;return this;
}

计算幅值和相位：

//幅度
public float getAmplitude() {return (float) Math.sqrt(Math.pow(realPart, 2)+Math.pow(imaginaryPart, 2));
}

//相位
public float getPhase() {if(realPart == 0) {if(imaginaryPart == 0)return 0;else if(imaginaryPart>0)return (float) (Math.PI/2);elsereturn (float) (-Math.PI/2);}else {return (float) Math.atan(imaginaryPart/realPart);}
}

3）复数数组

由于快速傅里叶变换涉及到大量数据的处理，若采用上述定义的复数的数组形式进行数据的存储，会消耗大量的内存，因此需要更高效的数据结构来存储多个复数。
      这里采用两个数组来存复数，一个存储复数的所有实部，另一个存储虚部。

private float[] realArray;
private float[] imaginaryArray;
private int size;public ComplexNumberArray(int size) {if(size<= 0) {this.size = 128;realArray = new float[this.size];imaginaryArray = new float[this.size];}else {this.size = size;realArray = new float[size];imaginaryArray = new float[size];}
}

同样的，我们也需要定义复数数组的相关的运算和计算。通过指定对应复数的索引来对数组中的复数进行操作。

2.单位根的实现

1）单位根的计算

private int n;
private List<ComplexNumber> rootList;private void calculate() {for(int k=0;k<=n;k++) {float unit = (float) (2 * Math.PI  * k / n);rootList.add(new ComplexNumber((float)Math.cos(unit), -(float)Math.sin(unit)));}
}

2）获取单位根时的处理

需要利用单位根的性质。

public ComplexNumber getUnitRoot(int n ,int k) {if(this.n == n) {return rootList.get(k);}else if(this.n>n) {int multiple = this.n/n;return rootList.get(k*multiple);}else if(this.n<n) {int multiple = n/this.n;return rootList.get(k/multiple);}return null;
}

3.快速傅里叶变换的实现

1）索引插值

对于两点采样点分组后位置为：0，1
      对于四点采样点分组后位置为：0，2，1，3
      对于八点采样点分组后位置为：0，4，2，6， 1，5， 3，7
      对于十六点采样点分组后位置为：0，8，4，12，2，10，6，14，1，9，5，13，3，11，7，15

通过观察可以发现：
      四点为两点在0后和1后进行插值，插的值为前一个值加上2的1次方。
      八点为四点在0、2、1、3后进行插值，插的值为前一个值加上2的2次方。
      十六点同理。
      即：n点采样点分组后位置为n/2点采样点分组后位置每个位置后插入前一个值和2的((√n)-1)次方的和。

//创建奇偶分类后的索引
private int[] createIndex(int n) {int[] index = new int[n];int left = 0;int middle = n/2;int right = n;index[left] = 0;index[middle] = 1;createIndexMerge(index,left,middle,1);createIndexMerge(index,middle,right,1);return index;
}
//递归调用
private void createIndexMerge(int[] index, int left ,int right, int multiple) {if(right - left<=1) return;int value = (int) Math.pow(2, multiple);int middle = (right+left)/2;index[middle]=(index[left]+value);createIndexMerge(index,left,middle,multiple+1);createIndexMerge(index,middle,right,multiple+1);
}

2）数据奇偶分组

根据上一步产生的位置信息，对数据进行排序。为了不污染原数据，用新的数组来进行排序。

//根据索引内容将数据进行重排序
//用于FFT
private  float[] sortDataByIndexFFT(float[] data) {clearArray(sortedFFTData);for(int i=0;i<index.length;i++) {int p = index[i];sortedFFTData[i] = data[p];}return sortedFFTData;
}

3）确定循环次数

以八点采样点的快速傅里叶变换为例子，观察蝶形图，需要进行3次大交叉计算，第一个交叉计算涉及2个点，计算2次，一共需要重复4次。第二个交叉计算涉及4个点，计算4次，一共需要重复八次。第三次交叉计算涉及8个点，计算8次，一共需要重复以此。
      即：n点采样点的快速傅里叶变换，需要进行(√n)次计算。k = 1,2,3…,√n。
      第k次交叉计算涉及2的k次方个点，计算2的k次方次，一共需要重复n除以(2的k次方)次。

4）核心计算

将每一次计算的数据平均分成上部分和下部分，每次计算分别从上部分和下部分相同次序处各取一个数据，上部分数据对下部分数据与对应单位根的乘积做和，结果存放在当前上部分数据的位置。上部分数据对下部分数据与对应单位根的乘积做差，结果存放在当前下部分数据的位置。

//快速傅里叶变换
public  ComplexNumberArray fft(float[] data) {sortedFFTData = sortDataByIndexFFT(data);ComplexNumberArray complexNumberArray = createComplexNumberArray(sortedFFTData);for(int i=1;i<=flyTime;i++) {//每次蝶形运算中数据的数量int teamSize = (int)(Math.pow(2, i));//大蝶形运算中需要计算的组数int teams = size/teamSize;for(int j=0;j<teams;j++) {int start1 = j*teamSize;int start2 = start1+teamSize/2;for(int k=0;k<teamSize/2;k++) {//核心计算complexNumberArray.multiply(start2, root.getUnitRoot(teamSize, k));ComplexNumber cNum = complexNumberArray.getComplexNumber(start2);complexNumberArray.setComplexNumber(start2, complexNumberArray.getRealPart(start1), complexNumberArray.getImaginaryPart(start1));complexNumberArray.add(start1, cNum);complexNumberArray.subtract(start2, cNum);start1++;start2++;}}}return complexNumberArray;
}

三.获取全部代码

https://github.com/HolyCrazy/FFT/

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

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

    快速傅里叶变换的时间复杂度分析 1 快速傅里叶变换FFT 1.1 理论分析 1.1.1 离散傅里叶变换 1.1.2 快速傅里叶变换 1.2 编程实现 1.2.1 算法思想 1.2.2 实验结果 1 快 ...

  2. MATLAB中的快速傅里叶变换FFT与IFFT

    背景 FFT (Fast Fourier Transform)是离散傅立叶变换的快速算法,可以将一个信号从时域变换到频域.同时与之对应的是IFFT(Inverse Fast Fourier Trans ...

  3. OpenCV快速傅里叶变换(FFT)用于图像和视讯流的模糊检测

    OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测 翻译自[OpenCV Fast Fourier Transform (FFT) for blur detection in images ...

  4. java fft_java实现快速傅里叶变换（FFT）

    最近做音频信号处理的时候,需要对数据做fft变换.关于fft原理: 请参考:FFT算法讲解--麻麻我终于会FFT了! matlab实现fft函数很简单,直接调用fft即可.但java实现起来就有点难了 ...

  5. matlab cftool光滑曲线导出为什么就不光滑了_快速傅里叶变换(FFT)中为什么要“补零”？...

    为了大家能够复现各个图中的结果,我附上了所有我编写的MATLAB代码. 创作不易,未经允许,禁止转载. 另外,说明一下,用MATLAB做FFT并不要求数据点个数必须为以2为基数的整数次方.之所以很多资 ...

  6. C语言二维数组范德蒙,浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理...

    浅谈范德蒙德(Vandermonde)方阵的逆矩阵与拉格朗日(Lagrange)插值的关系以及快速傅里叶变换(FFT)中IDFT的原理 标签: 行列式 矩阵 线性代数 FFT 拉格朗日插值 只要稍微看 ...

  7. 浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理...

    浅谈范德蒙德(Vandermonde)方阵的逆矩阵与拉格朗日(Lagrange)插值的关系以及快速傅里叶变换(FFT)中IDFT的原理 标签: 行列式 矩阵 线性代数 FFT 拉格朗日插值 只要稍微看 ...

  8. 基于python的快速傅里叶变换FFT（二）

    基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点   FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...

  9. 基于python的快速傅里叶变换FFT（一）

    基于python的快速傅里叶变换FFT(一) FFT可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了.这就是很多信号分析采用FFT变换的原因. ...

最新文章

  1. java tree degree_生成树计数-Matrix-Tree定理
  2. 学习LOWORD、 HIWORD、LOBYTE、HIBYTE
  3. mysql 执行计划详解_mysql explain执行计划详解
  4. 全球及中国散热市场应用格局与竞争态势研究报告2022-2027年
  5. 去掉日志服务器性能,日志服务器及性能监控
  6. 微软云架构服务器,微软云存储架构（Azure Cloud Storage）
  7. php 打印 域名ip_php如何获取域名IP地址代码函数
  8. spring boot配置ip_Zookeeper作为配置中心使用说明
  9. windows 快捷方式(.lnk)代码执行漏洞(CVE-2017-8464 )[附EXP生成工具]
  10. systemtap打点方法
  11. k3cloud6.0文件服务器,K3Cloud系统集成配置详解
  12. Android 闹钟设置最新版
  13. python合成gif动图
  14. web全栈工程师技能介绍
  15. css特效实例——纯css实现带边角卷边阴影的纸
  16. [评论送书 ]手撕源码，实现一个Koa。
  17. SQLserver2008R2详细安装教程
  18. 移动应用数据保护都需要哪些安全技术
  19. 三菱FX3U，单控气缸报警功能块，双控气缸报警功能块，真空报警功能块，伺服定位功能块，伺服手动操作功能块，产量节拍功能块
  20. 电脑安装java环境

热门文章

  1. contest18 CF788 div1 ooxxx oooox oooox
  2. 2021年全球圆锥破碎机收入大约1357.4百万美元，预计2028年达到1665.6百万美元
  3. 【食品化学与营养】第二章 水的化学与营养 笔记
  4. Python Tkinter——数字拼图游戏详解版
  5. C#编写的Word操作类,有换页,添加表格,文本功能
  6. 价值1K的微信二维码活码源码
  7. 查看linux设备Ran的大小,linux – 来自/ dev / zero和/ dev / urandom的不同文件大小
  8. 网络安全毕业设计选题题目大全
  9. manifestintert-filter详解
  10. @sequencegenerator oracle,SequenceGenerator注解的使用