C#的FFT函数(傅立叶函数)算法

发信人: jinnxjtu (Jinn), 信区: DotNET
标 题: 同学编的C＃的FFT函数拿出来共享
发信站: BBS 水木清华站 (Mon Jun 16 23:41:25 2003), 转信

struct complex//定义复数
{

//复数a+bi中
//a为实部，b为虚部
public double a;
public double b;
public static complex commul(double x,complex y)
{
complex c = new complex();
c.a = x*y.a;
c.b = x*y.b;

return c;
}
public static complex commulcc(complex x,complex y)
{
complex c = new complex();
c.a = x.a*y.a - x.b*y.b;
c.b = x.a*y.b + x.b*y.a;

return c;
}
public static complex comsum(complex x,complex y)
{
complex c = new complex();
c.a = x.a + y.a;
c.b = x.b + y.b;

return c;
}
public static complex comsum1(double x,complex y)
{
complex c = new complex();
c.a = x + y.a;
c.b = y.b;

return c;
}
public static complex decrease(complex x,complex y)
{
complex c = new complex();
c.a = x.a*y.a - x.b*y.b;
c.b = x.a*y.b + x.b*y.a;

return c;
}
public static complex powcc(complex x,double n)
{
int k;
complex xout;
xout.a = 1;
xout.b = 0;
for(k=1; k<=n; k++)
{
xout = mulcc(xout,x);
}
return xout;
}
public void show()
{
Console.Write(a+" + "+b+"i ");
}
}
class airthm
{
//计算ω=exp（j*2*pi/n）
public static complex omega( int n)
{
complex x;
x.a = Math.Cos(0 - 2*Math.PI/n);
x.b = Math.Sin(0 - 2*Math.PI/n);
return x;
}
public static complex[] dft(double[] signal, int n) //(信号，信号长度)
{
int i,j;
complex w1;
w1 = omega(n);
complex[] w = new complex[n];
for(i=0; i{
w[i] = complex.powcc(w1,i);
}
complex[] f = new complex[n];
complex temp; //w[i]的次方
complex temp1; //f中单项的值
for (i=0; i{
f[i].a = 0;
f[i].b = 0;
for (j=0; j{

temp = complex.powcc(w[i],j);
temp1 = mul(signal[j],temp);
f[i] = sum(f[i],temp1) ;
}
}
return f;
}
//求幅值 x 信号 n 信号长度 返回 幅值数组
public static double[] amplitude (complex[] x,int n)
{
int i;
double temp ;
double[] amp = new double[n];
for (i=0; i{
temp = x[i].a*x[i].a + x[i].b*x[i].b;
amp[i] = Math.Sqrt(temp);
}
return amp;
}
}
/////////////////////////////////////////////////////////////////////
//调用
int n=6;
double[] x = new double[]{1,2,3,4,5,6};//被计算的数组
complex[] y = new complex[n];//接收复数结果的数组
double[] z = new Double[n];//接收幅值结果的数组
y = airthm.dft(x,n);
z = airthm.amplitude(y,n);