Java中实现快速傅里叶变换FFT
傅里叶变换是将一个满足条件的函数表示成三角函数或其积分的线性组合。
实际应用中，将一个随时间变化的信号，经过傅里叶变换，可以得到这个信号的频率组成，实现从时域到频域的变换，进而研究其频谱结构和变化规律。
傅里叶变换表明：一个信号可以表示成多个不同频率的正弦波的叠加。
对一个连续的时域信号以一定的频率进行采样，可以得到这个信号时域的一组采样点，通过对这些采样点进行离散傅里叶变换，可以得到这个信号对应频域的一组采样点。
由公式可以得出离散傅里叶变换的时间复杂度为O(n2)。
快速傅里叶变换是对离散傅里叶变换进行快速计算的一种方法。可以将时间复杂度降低到O(nlog2n)，因此称为快速傅里叶变换。
1）单位根 在一个复平面内将单位圆平均分成n份，这n个点表示的复数称为n次单位根。
其中k = 0，1，2，…n-1，方向逆时针。
单位根的性质：
将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 |
将奇数项和偶数项的展开式从下向上计算的过程，可以通过画图来表示，由于该图形似蝴蝶，故称为蝶形图。
快速傅里叶变换：
对信号频域的采样点进行快速傅里叶变换的逆变换，可以得到信号在时域对应的采样点。
具体实现就是快速傅里叶变换蝶形计算的逆过程，但是其单位根为快速傅里叶变换单位根的共轭形式，同时，每次计算后结果要除以2。
快速傅里叶变换的逆变换：
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);
else
return (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 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次方)次。
将每一次计算的数据平均分成上部分和下部分，每次计算分别从上部分和下部分相同次序处各取一个数据，上部分数据对下部分数据与对应单位根的乘积做和，结果存放在当前上部分数据的位置。上部分数据对下部分数据与对应单位根的乘积做差，结果存放在当前下部分数据的位置。
//快速傅里叶变换
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/
作者：LeeDuo.
