C#實現(xiàn)FFT(遞歸法)的示例代碼

更新時間：2022年07月11日 14:05:34 作者：Mokera

FFT是數(shù)字信號處理中的重要算法。這篇文章將為大家詳細介紹一下如何利用C#語言實現(xiàn)FFT(遞歸法)，文中的示例代碼講解詳細，感興趣的可以了解一下

1. C#實現(xiàn)復(fù)數(shù)類

我們在進行信號分析的時候，難免會使用到復(fù)數(shù)。但是遺憾的是，C#沒有自帶的復(fù)數(shù)類，以下提供了一種復(fù)數(shù)類的構(gòu)建方法。

復(fù)數(shù)相比于實數(shù)，可以理解為一個二維數(shù)，構(gòu)建復(fù)數(shù)類，我們需要實現(xiàn)以下這些內(nèi)容：

復(fù)數(shù)實部與虛部的屬性
復(fù)數(shù)與復(fù)數(shù)的加減乘除運算
復(fù)數(shù)與實數(shù)的加減乘除運算
復(fù)數(shù)取模
復(fù)數(shù)取相位角
歐拉公式（即e^ix+y）

C#實現(xiàn)的代碼如下：

 public class Complex
    {
        double real;
        double imag;
        public Complex(double x, double y)   //構(gòu)造函數(shù)
        {
            this.real = x;
            this.imag = y;
        }
        //通過屬性實現(xiàn)對復(fù)數(shù)實部與虛部的單獨查看和設(shè)置
        public double Real
        {
            set { this.real = value; }
            get { return this.real; }
        }
        public double Imag
        {
            set { this.imag = value; }
            get { return this.imag; }
        }
        //重載加法
        public static Complex operator +(Complex c1, Complex c2)
        {
            return new Complex(c1.real + c2.real, c1.imag + c2.imag);
        }
        public static Complex operator +(double c1, Complex c2)
        {
            return new Complex(c1 + c2.real, c2.imag);
        }
        public static Complex operator +(Complex c1, double c2)
        {
            return new Complex(c1.Real + c2, c1.imag);
        }
        //重載減法
        public static Complex operator -(Complex c1, Complex c2)
        {
            return new Complex(c1.real - c2.real, c1.imag - c2.imag);
        }
        public static Complex operator -(double c1, Complex c2)
        {
            return new Complex(c1 - c2.real, -c2.imag);
        }
        public static Complex operator -(Complex c1, double c2)
        {
            return new Complex(c1.real - c2, c1.imag);
        }
        //重載乘法
        public static Complex operator *(Complex c1, Complex c2)
        {
            double cr = c1.real * c2.real - c1.imag * c2.imag;
            double ci = c1.imag * c2.real + c2.imag * c1.real;
            return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
        }
        public static Complex operator *(double c1, Complex c2)
        {
            double cr = c1 * c2.real;
            double ci = c1 * c2.imag;
            return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
        }
        public static Complex operator *(Complex c1, double c2)
        {
            double cr = c1.Real * c2;
            double ci = c1.Imag * c2;
            return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
        }

        //重載除法
        public static Complex operator /(Complex c1, Complex c2)
        {
            if (c2.real == 0 && c2.imag == 0)
            {
                return new Complex(double.NaN, double.NaN);
            }
            else
            {
                double cr = (c1.imag * c2.imag + c2.real * c1.real) / (c2.imag * c2.imag + c2.real * c2.real);
                double ci = (c1.imag * c2.real - c2.imag * c1.real) / (c2.imag * c2.imag + c2.real * c2.real);
                return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));           //保留四位小數(shù)后輸出
            }
        }
      
        public static Complex operator /(double c1, Complex c2)
        {
            if (c2.real == 0 && c2.imag == 0)
            {
                return new Complex(double.NaN, double.NaN);
            }
            else
            {
                double cr = c1 * c2.Real / (c2.imag * c2.imag + c2.real * c2.real);
                double ci = -c1 * c2.imag / (c2.imag * c2.imag + c2.real * c2.real);
                return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));           //保留四位小數(shù)后輸出
            }
        }
      
        public static Complex operator /(Complex c1, double c2)
        {
            if (c2 == 0)
            {
                return new Complex(double.NaN, double.NaN);
            }
            else
            {
                double cr = c1.Real / c2;
                double ci = c1.imag / c2;
                return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));           //保留四位小數(shù)后輸出
            }
        }
        //創(chuàng)建一個取模的方法
        public static double Abs(Complex c)
        {
            return Math.Sqrt(c.imag * c.imag + c.real * c.real);
        }
        //創(chuàng)建一個取相位角的方法
        public static double Angle(Complex c)
        {
            return Math.Round(Math.Atan2(c.real, c.imag), 6);//保留6位小數(shù)輸出
        }
        //重載字符串轉(zhuǎn)換方法,便于顯示復(fù)數(shù)
        public override string ToString()
        {
            if (imag >= 0)
                return string.Format("{0}+i{1}", real, imag);
            else
                return string.Format("{0}-i{1}", real, -imag);
        }
        //歐拉公式
        public static Complex Exp(Complex c)
        {
            double amplitude = Math.Exp(c.real);
            double cr = amplitude * Math.Cos(c.imag);
            double ci = amplitude * Math.Sin(c.imag);
            return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));//保留四位小數(shù)輸出
        }
    }

2. 遞歸法實現(xiàn)FFT

以下的遞歸法是基于奇偶分解實現(xiàn)的。

奇偶分解的原理推導(dǎo)如下：

x(2r)和x(2r+1)都是長度為N/2−1的數(shù)據(jù)序列，不妨令

則原來的DFT就變成了：

于是，將原來的N點傅里葉變換變成了兩個N/2點傅里葉變換的線性組合。

但是，N/2點傅里葉變換只能確定N/2個頻域數(shù)據(jù)，另外N/2個數(shù)據(jù)怎么確定呢？

因為X₁(k)和X₂(k)周期都是N/2，所以有

從而得到：

綜上，我們就可以得到遞歸法實現(xiàn)FFT的流程：

1.對于每組數(shù)據(jù)，按奇偶分解成兩組數(shù)據(jù)

2.兩組數(shù)據(jù)分別進行傅里葉變換，得到X₁(k)和X₂(k)

3.總體數(shù)據(jù)的X(k)由下式確定：

4.對上述過程進行遞歸

具體代碼實現(xiàn)如下：

public Complex[] FFTre(Complex[] c)
{
    int n = c.Length;
    Complex[] cout = new Complex[n];
    if (n == 1)
    {
        cout[0] = c[0];
        return cout;
    }
    else
    {
        double n_2_f = n / 2;
        int n_2 = (int)Math.Floor(n_2_f);
        Complex[] c1 = new Complex[n / 2];
        Complex[] c2 = new Complex[n / 2];
        for (int i = 0; i < n_2; i++)
        {
            c1[i] = c[2 * i];
            c2[i] = c[2 * i + 1];
        }
        Complex[] c1out = FFTre(c1);
        Complex[] c2out = FFTre(c2);
        Complex[] c3 = new Complex[n / 2];
        for (int i = 0; i < n / 2; i++)
        {
            c3[i] = new Complex(0, -2 * Math.PI * i / n);
        }
        for (int i = 0; i < n / 2; i++)
        {
            c2out[i] = c2out[i] * Complex.Exp(c3[i]);
        }

        for (int i = 0; i < n / 2; i++)
        {
            cout[i] = c1out[i] + c2out[i];
            cout[i + n / 2] = c1out[i] - c2out[i];
        }
        return cout;
    }
}

3. 補充：窗函數(shù)

順便提供幾個常用的窗函數(shù)：

Rectangle
Bartlett
Hamming
Hanning
Blackman

    public class WDSLib
    {
        //以下窗函數(shù)均為periodic
        public double[] Rectangle(int len)
        {
            double[] win = new double[len];
            for (int i = 0; i < len; i++)
            {
                win[i] = 1;
            }
            return win;
        }

        public double[] Bartlett(int len)
        {
            double length = (double)len - 1;
            double[] win = new double[len];
            for (int i = 0; i < len; i++)
            {
                if (i < len / 2) { win[i] = 2 * i / length; }
                else { win[i] = 2 - 2 * i / length; }
            }
            return win;
        }

        public double[] Hamming(int len)
        {
            double[] win = new double[len];
            for (int i = 0; i < len; i++)
            {
                win[i] = 0.54 - 0.46 * Math.Cos(Math.PI * 2 * i / len);
            }
            return win;
        }

        public double[] Hanning(int len)
        {
            double[] win = new double[len];
            for (int i = 0; i < len; i++)
            {
                win[i] = 0.5 * (1 - Math.Cos(2 * Math.PI * i / len));
            }
            return win;
        }

        public double[] Blackman(int len)
        {
            double[] win = new double[len];
            for (int i = 0; i < len; i++)
            {
                win[i] = 0.42 - 0.5 * Math.Cos(Math.PI * 2 * (double)i / len) + 0.08 * Math.Cos(Math.PI * 4 * (double)i / len);
            }
            return win;
        }
    }

以上就是C#實現(xiàn)FFT(遞歸法)的示例代碼的詳細內(nèi)容，更多關(guān)于C# FFT遞歸法的資料請關(guān)注腳本之家其它相關(guān)文章！

您可能感興趣的文章: