「快速傅里叶变换」快速傅里叶变换及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

<span style="font-size:12px;">public class DFT {
public Complex[] dft(Complex[] x) {
int n = x.length;
// exp(-2i*n*PI)=cos(-2*n*PI)+i*sin(-2*n*PI)=1
if (n == 1)
return x;
Complex[] result = new Complex[n];
for (int i = 0; i < n; i++) {
result[i] = new Complex(0, 0);
for (int k = 0; k < n; 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[i].plus(x[k].multiple(m));
}
}
return result;
}
}
</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

<span style="font-size:12px;">public class FFT {
public Complex[] fft(Complex[] x) {
int n = x.length;
// 因为exp(-2i*n*PI)=1，n=1时递归原点
if (n == 1){
return x;
}
// 如果信号数为奇数，使用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;
}
public Complex[] dft(Complex[] x) {
int n = x.length;
// 1个信号exp(-2i*n*PI)=1
if (n == 1)
return x;
Complex[] result = new Complex[n];
for (int i = 0; i < n; i++) {
result[i] = new Complex(0, 0);
for (int k = 0; k < n; 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[i].plus(x[k].multiple(m));
}
}
return result;
}
}</span>

计算1000个采样信号用时：

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

快速傅里叶变换及java实现

快速傅里叶变换

相关阅读

栏目导航

推荐阅读

热门阅读