C＃的FFT函数

同学编的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 = complex.commulcc(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<n; i++)
                {
                        w[i] = complex.powcc(w1,i);
                }
                complex[] f = new complex[n];
                complex temp;  //w[i]的次方
                complex temp1; //f中单项的值
                for (i=0; i<n; i++)
                {
                        f[i].a = 0;
                        f[i].b = 0;
                        for (j=0; j<n; j++)
                        {
                                temp = complex.powcc(w[i],j);                           
                                temp1 = complex.commul(signal[j],temp);
                                f[i] = complex.comsum(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<n; 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);





--