最近工作中有个需求,在C#环境中实现FFT算法,在网上找了些资料,最后实现了下面的两种方式,实际应用任选其一就好。
第一种方法:
不依赖C#中的Complex,需要实现计算过程的每一步详细步骤。
输入序列长度为2的N次幂,使用前需先定义序列长度:
FFT filter = new FFT(256);
filter.fft(x,y) 其中x为实部y为虚部,计算后x为FFT后的实部,y为FFT后的虚部。
class FFT
{
int n, m;
// Lookup tables. Only need to recompute when size of FFT changes.
double[] cos;
double[] sin;
public FFT(int n)
{
this.n = n;
this.m = (int)(Math.Log(n) / Math.Log(2));
// Make sure n is a power of 2
if (n != (1 << m))
Console.Out.WriteLine("FFT length must be power of 2");
// precompute tables
cos = new double[n / 2];
sin = new double[n / 2];
for (int i = 0; i < n / 2; i++)
{
cos[i] = Math.Cos(-2 * Math.PI * i / n);
sin[i] = Math.Sin(-2 * Math.PI * i / n);
}
}
///x为实部y为虚部///
private void fft(double[] x, double[] y)
{
int i, j, k, n1, n2, a;
double c, s, t1, t2;
// Bit-reverse
j = 0;
n2 = n / 2;
for (i = 1; i < n - 1; i++)
{
n1 = n2;
while (j >= n1)
{
j = j - n1;
n1 = n1 / 2;
}
j = j + n1;
if (i < j)
{
t1 = x[i];
x[i] = x[j];
x[j] = t1;
t1 = y[i];
y[i] = y[j];
y[j] = t1;
}
}
// FFT
n1 = 0;
n2 = 1;
for (i = 0; i < m; i++)
{
n1 = n2;
n2 = n2 + n2;
a = 0;
for (j = 0; j < n1; j++)
{
c = cos[a];
s = sin[a];
a += 1 << (m - i - 1);
for (k = j; k < n; k = k + n2)
{
t1 = c * x[k + n1] - s * y[k + n1];
t2 = s * x[k + n1] + c * y[k + n1];
x[k + n1] = x[k] - t1;
y[k + n1] = y[k] - t2;
x[k] = x[k] + t1;
y[k] = y[k] + t2;
}
}
}
}
}
第二种方式:
傅里叶变换的计算中需要用到复数,无意中发现C#中有复数的定义,考虑使用复数直接计算的方式实现,这样可以利用C#中的复数计算。C#中使用复数需要引用using System.Numerics;
FFT filter = new FFT(256);
filter.DoFFT(x) ; x为输入参数(Complex[]类型),计算后x为FFT后的结果。
突然发现C#好好用啊。。。。。。
class FFT
{
private Complex[] W;//定义旋转因子
/// <summary>
/// 构造函数用于初始化旋转因子。
/// </summary>
public FFT(int N)
{
W = initW(N);
}
/// <summary>
/// 快速傅立叶正变换函数
/// W:输入旋转因子序列数组名
/// x: 输入信号序列数组名
/// </summary>
public void DoFFT(Complex[] x)
{
int i = 0, j = 0, k = 0, l = 0;
Complex up, down, product;
ReverseOder(x); ///对输入序列进行倒序
for (i = 0; i < Math.Log(N, 2); i++) ///外循环,变换级数循环
{
l = 1 << i;
for (j = 0; j < N; j += 2 * l) ///中间循环,旋转因子循环
{
for (k = 0; k < l; k++) ///内循环,序列循环
{
product = x[j + k + l] * W[N * k / 2 / l];
up = x[j + k] + product;
down = x[j + k] - product;
x[j + k] = up;
x[j + k + l] = down;
}
}
}
}
/// <summary>
/// 倒序函数
/// W:输入旋转因子序列长度
/// </summary>
private Complex[] initW(int N)
{
Complex[] W = new Complex[N];
for (int i = 0; i < N; i++)
{
double real = Math.Cos(2 * Math.PI * i / N); //旋转因子实部展开
double imag = -1 * Math.Sin(2 * Math.PI * i / N); //旋转因子虚部展开
W[i] = new Complex(real, imag);
}
return W;
}
/// <summary>
/// ReverseOder()倒序变换
/// </summary>
/// <param name="x"></param>
private void ReverseOder(Complex[] x)
{
Complex temp;
int i = 0, j = 0, k = 0;
double t;
for (i = 0; i < N; i++)
{
k = i; j = 0;
t = Math.Log(N, 2);
while ((t--) > 0)
{
j = j << 1;
j |= (k & 1);
k = k >> 1;
}
if (j > i)
{
temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
}