java FFT的实现及部分个人理解

0.FFT拿来干啥

简单来说：将时域波形转换到频域，更容易观察信号的规律

模拟信号f(t)->Fs时域采样频率采样->离散信号f(nT)->离散序列->f(k)(时间轴变为T）->进行FFt->变换到频域得到F(n)

1.复数

复数分为实部和虚部
复数的模=实部和虚部的平方之和
（实部，虚部类比x，y轴，模为到原点的距离）
复数运算：
加减法：实部虚部分开进行 加减
乘法：
(a1+b1j) * (a2+b2j)=(a1a2 -b1b2) + (a1b2j+b1a2j)
(j^2=-1)
类比乘法
(a1+b1)*( a2+b2)

2.DFT因子（也有叫旋转因子的）

如何理解：
将一个单位圆N等分，从x轴正半轴开始逆时针转一圈的第k个分隔点

有如下性质：
1.共轭复数对称性

2.周期性

3.不知道叫啥的特性（摘自书本信号与线性系统下册9.6节）

前两个公式懒得自己打了，从网上直接拿的别人的图
虽然这几个特性都蛮重要的，但是对于写代码的我们来说，貌似暂时用不上
（可能还是用的上的，在数据量非常大的情况下，需要对每个DFT因子的结果进行缓存）

3.欧拉公式（Euler）

如何理解：
结合单位圆一看就很清楚了。
虚数轴上的值等于单位圆半径1xsinx
实数轴上的值等于单位圆半径1xcosx
这个公式的作用：拿来计算上述的DFT因子

4.蝶形运算

蝶形运算分为：按时间抽取和按频率抽取


摘自万能的信号与线性系统分析下册

这两种算法我个人理解的区别：
时间抽取：
从输出出发，对输入进行排序归类，输入排序是位倒置的
频率抽取：
从输入出发，对输出进行排序归类，输入排序是位倒置的
计算复杂程度感觉是相同的，但是实际写的过程中，感觉第一种更容易被人理解

位倒置的意思：
首先做N点FFt，N是二的倍数，N=2^k次方
则k就是这个位数
比如：
8点fft
位数为3

输入	二进制	倒置	输入
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

人为进行计算时需要先将输入按照这种方式进行排序，然后两两做蝶形运算一层一层直到输出（按频率抽取不需要，直接输入层做蝶形运算直到输出层，再对输出层进行位倒置的排序）

实际我们写代码时，从输出倒推递归（针对按时间抽取这种方式），然后再对其通过奇偶的抽取分组会更方便，并且这样省去了对其位倒叙的操作

以下demo是按时间抽取的：

复数类的实现：

	public static class Complex {
		public double i;
		public double j;// 虚数部分

		public Complex(double i, double j) {
			this.i = i;
			this.j = j;
		}

		public double getMod() {// 求复数的模
			return Math.sqrt(i * i + j * j);
		}

		public static Complex Add(Complex a, Complex b) {
			return new Complex(a.i + b.i, a.j + b.j);
		}

		public static Complex Subtract(Complex a, Complex b) {
			return new Complex(a.i - b.i, a.j - b.j);
		}

		public static Complex Mul(Complex a, Complex b) {// 乘法
			return new Complex(a.i * b.i - a.j * b.j, a.i * b.j + a.j * b.i);
		}

		public static Complex GetW(int k, int N) {
			// 欧拉公式求DFT因子
			return new Complex(Math.cos(-2 * Math.PI * k / N), Math.sin(-2 * Math.PI * k / N));
		}

		public static Complex[] butterfly(Complex a, Complex b, Complex w) {
			//蝶形运算
			return new Complex[] { Add(a, Mul(w, b)), Subtract(a, Mul(w, b)) };
		}

		public static Double[] toModArray(Complex[] complex) {
			//求复数数组的各个模
			Double[] res = new Double[complex.length];
			for (int i = 0; i < complex.length; i++) {
				res[i] = complex[i].getMod();
			}
			return res;
		}
	}

FFt：

	private static Complex[] FFT(Complex[] input, int N) {
		if ((N / 2) % 2 == 0) {
			Complex[] even = new Complex[N / 2];// 偶数
			Complex[] odd = new Complex[N / 2];// 奇数
			for (int i = 0; i < N / 2; i++) {
				even[i] = input[2 * i];
				odd[i] = input[2 * i + 1];
			}
			even = FFT(even, N / 2);
			odd = FFT(odd, N / 2);
			for (int i = 0; i < N / 2; i++) {
				var res = Complex.butterfly(even[i], odd[i], Complex.GetW(i, N));
				input[i] = res[0];
				input[i + N / 2] = res[1];
			}
			return input;
		} else {// 两点DFT,直接进行碟形运算
			var res = Complex.butterfly(input[0], input[1], Complex.GetW(0, N));
			return res;
		}
	}

其中

			even = FFT(even, N / 2);
			odd = FFT(odd, N / 2);

就是递归，然后得到递归结果之后再对奇偶两个数组进行两两蝶形运算

摘自信号与线性系统分析下册：
如果N/2也是一个偶数，在计算N/2点DFT时，同样可以将其分为两个长度为N/4点的DFT进行计算， 这样的拆解过程可以不断地进行下去，直到DFT的点数不是一个偶数为止（对于2次幂DFT就是1）
这里就直接拍书本里的这个图了，很好理解，每次dft时按奇偶对输入分组（注以分组后的序列id为准）

验证结果：
利用FFT对f(n)=｛1，2，3，4｝做4点DFT

		Complex[] input = new Complex[4];
		for (int i = 1; i <= 4; i++)
			input[i - 1] = new Complex(i, 0);
		input = FFT(input, 4);
		for (Complex c : input) {
			System.out.println(c.i + "+" + c.j + "j");
		}

结果：

FFT后的物理意义：

1.频率分辨率：
Δf=Fs/N
即采样频率除以点数
2.得到复数数组后每个点的频率如何计算？
F=Δf*n
如：
第0个点是直流分量，第1个点的频率是Δf
第二个点是2Δf。。。。第n个点是(n-1)Δf
3.幅值
A=
Mod(complex)/N --直流分量
Mod(complex)/(N/2)–交流分量

还有相位频谱之类的分析。暂时还没看懂，以后用上了再去看，然后补充上来

最后放个动图，方便理解fft：