快速傅里叶变换及java实现

傅里叶变换：

    傅里叶变换是一种线性的积分变换。它的理论依据是：任何连续周期信号都可以由一组适当的正弦曲线组合而成，即使用简单的正弦、余弦函数（如sinx,Acos(ωx+θ)）,可以拟合复杂函数。

    使用正弦曲线的原因：在信号处理中，正弦曲线可以更简单地处理信号，且一个正弦曲线信号经过处理仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波形不变。

   

    在信号处理中，傅里叶变换（连续）是将时域信号积分，得到频域上的信息（将一条曲线拆分成正弦曲线后，各正弦曲线的振幅，图中红色部分）：

    

    傅里叶逆变换（连续）是将频域信号积分，得到时域上的信息（将各个正弦曲线合成后的曲线，图中蓝色部分）：

其中，e^iwt = cos(wx)+ i * sin(wx)，（欧拉公式），表示复平面上的一个点。

傅里叶变换类型：

    1.   连续傅立叶变换：非周期性连续信号；

    2.   傅立叶级数：周期性连续信号；

    3.   离散时域傅立叶变换：非周期性离散信号；

    4.   离散傅立叶变换：周期性离散信号。

 

离散傅立叶变换（DFT）：

    对于计算机来说只有离散的和有限长度的数据才能被处理，所以计算机只能处理离散傅立叶变换，而其它的变换类型只有在数学演算中才能用到。关于离散傅里叶变换的公式：

    

离散傅里叶变换的算法复杂度为O( N²)，其中x_n表示合成的离散信号，X_k表示x_n在时域上的变换（这里k没有固定值，但k越大，对x_n的拟合曲线就越多，拟合更精确）。以下是一个例子：

 

图中原始信号有4个采样点，x轴表示时间（可以看出是等距采样），y轴表示振幅。而后面两个坐标，x轴表示频率个数（时间被积分了），y轴表示振幅。

离散傅里叶变换java实现如下，输入：离散信号值（实数），输出：经过傅里叶变换得到的一组频率（实数）。（Complex是一个复数类，为了方便阅读，省略Complex的实现）

[java]  view plain  copy

1. <span style="font-size:12px;">public class DFT {  
2.   
3.     public Complex[] dft(Complex[] x) {  
4.         int n = x.length;  
5.   
6.         // exp(-2i*n*PI)=cos(-2*n*PI)+i*sin(-2*n*PI)=1  
7.         if (n == 1)  
8.             return x;  
9.   
10.         Complex[] result = new Complex[n];  
11.         for (int i = 0; i < n; i++) {  
12.             result[i] = new Complex(0, 0);  
13.             for (int k = 0; k < n; k++) {  
14.                 //使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)  
15.                 double p = -2 * k * Math.PI / n;  
16.                 Complex m = new Complex(Math.cos(p), Math.sin(p));  
17.                 result[i].plus(x[k].multiple(m));  
18.             }  
19.         }  
20.         return result;  
21.     }  
22. }  
23. </span>  

计算1000个采样信号用时：

 

快速傅里叶变换（FFT）：

    由于离散傅里叶变换的算法复杂度为O( N² )，在处理较大计算量时花费时间较长，所以JWCooley和JohnTukey在1965年的文章中提出了Cooley-TukeyFFT算法，利用X_N+k = X_k这一对称性（证明略），将DFT算法的时间复杂度降低到了O(NlogN )：

      

可以看到，式子被分解为两个更小的DFT（每个含有N/2个元素）。在这两个更小的DFT中，可以使用X_N/2+m= X_m这一性质，因为k < N, m < N/2，所以m + N/2 < N，所以X_N/2+m = X_m = X_mn + e^{-i2π m/N} * X_mm，即可以递归。又因为X₁ = x₀（exp(-2*n*PI)=1），所以n=1为原点。

需要注意的是，N必须是2的倍数。在递归中，如果出现N是奇数，则等式不成立，需要转到DFT进行计算。

快速傅里叶变换java实现如下，输入：离散信号值（实数），输出：经过傅里叶变换得到的一组频率（实数）。

[java]  view plain  copy

1. <span style="font-size:12px;">public class FFT {  
2.   
3.     public Complex[] fft(Complex[] x) {  
4.         int n = x.length;  
5.   
6.         // 因为exp(-2i*n*PI)=1，n=1时递归原点  
7.         if (n == 1){  
8.             return x;  
9.         }  
10.   
11.         // 如果信号数为奇数，使用dft计算  
12.         if (n % 2 != 0) {  
13.             return dft(x);  
14.         }  
15.   
16.         // 提取下标为偶数的原始信号值进行递归fft计算  
17.         Complex[] even = new Complex[n / 2];  
18.         for (int k = 0; k < n / 2; k++) {  
19.             even[k] = x[2 * k];  
20.         }  
21.         Complex[] evenValue = fft(even);  
22.   
23.         // 提取下标为奇数的原始信号值进行fft计算  
24.         // 节约内存  
25.         Complex[] odd = even;  
26.         for (int k = 0; k < n / 2; k++) {  
27.             odd[k] = x[2 * k + 1];  
28.         }  
29.         Complex[] oddValue = fft(odd);  
30.   
31.         // 偶数+奇数  
32.         Complex[] result = new Complex[n];  
33.         for (int k = 0; k < n / 2; k++) {  
34.             // 使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)  
35.             double p = -2 * k * Math.PI / n;  
36.             Complex m = new Complex(Math.cos(p), Math.sin(p));  
37.             result[k] = evenValue[k].plus(m.multiple(oddValue[k]));  
38.             // exp(-2*(k+n/2)*PI/n) 相当于 -exp(-2*k*PI/n)，其中exp(-n*PI)=-1(欧拉公式);  
39.             result[k + n / 2] = evenValue[k].minus(m.multiple(oddValue[k]));  
40.         }  
41.         return result;  
42.     }  
43.       
44.     public Complex[] dft(Complex[] x) {  
45.         int n = x.length;  
46.   
47.         // 1个信号exp(-2i*n*PI)=1  
48.         if (n == 1)  
49.             return x;  
50.   
51.         Complex[] result = new Complex[n];  
52.         for (int i = 0; i < n; i++) {  
53.             result[i] = new Complex(0, 0);  
54.             for (int k = 0; k < n; k++) {  
55.                 //使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)  
56.                 double p = -2 * k * Math.PI / n;  
57.                 Complex m = new Complex(Math.cos(p), Math.sin(p));  
58.                 result[i].plus(x[k].multiple(m));  
59.             }  
60.         }  
61.         return result;  
62.     }  
63. }</span>  

计算1000个采样信号用时：

 

与DFT算法相比，速度有明显提高。