C#實現FFT(遞迴法)的範例程式碼

2022-07-11 18:03:45

1. C#實現複數類

我們在進行訊號分析的時候，難免會使用到複數。但是遺憾的是，C#沒有自帶的複數類，以下提供了一種複數類的構建方法。

複數相比於實數，可以理解為一個二維數，構建複數類，我們需要實現以下這些內容：

複數實部與虛部的屬性
複數與複數的加減乘除運算
複數與實數的加減乘除運算
複數取模
複數取相位角
尤拉公式（即e^ix+y）

C#實現的程式碼如下：

 public class Complex
    {
        double real;
        double imag;
        public Complex(double x, double y)   //建構函式
        {
            this.real = x;
            this.imag = y;
        }
        //通過屬性實現對複數實部與虛部的單獨檢視和設定
        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));           //保留四位小數後輸出
            }
        }
      
        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));           //保留四位小數後輸出
            }
        }
      
        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));           //保留四位小數後輸出
            }
        }
        //建立一個取模的方法
        public static double Abs(Complex c)
        {
            return Math.Sqrt(c.imag * c.imag + c.real * c.real);
        }
        //建立一個取相位角的方法
        public static double Angle(Complex c)
        {
            return Math.Round(Math.Atan2(c.real, c.imag), 6);//保留6位小數輸出
        }
        //過載字串轉換方法,便於顯示覆數
        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));//保留四位小數輸出
        }
    }

2. 遞迴法實現FFT

以下的遞迴法是基於奇偶分解實現的。

奇偶分解的原理推導如下：

x(2r)和x(2r+1)都是長度為N/2−1的資料序列，不妨令

則原來的DFT就變成了：

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

但是，N/2點傅立葉變換隻能確定N/2個頻域資料，另外N/2個資料怎麼確定呢？

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

從而得到：

綜上，我們就可以得到遞迴法實現FFT的流程：

1.對於每組資料，按奇偶分解成兩組資料

2.兩組資料分別進行傅立葉變換，得到X₁(k)和X₂(k)

3.總體資料的X(k)由下式確定：

4.對上述過程進行遞迴

具體程式碼實現如下：

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. 補充：窗函數

順便提供幾個常用的窗函數：

Rectangle
Bartlett
Hamming
Hanning
Blackman

    public class WDSLib
    {
        //以下窗函數均為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#實現FFT(遞迴法)的範例程式碼的詳細內容，更多關於C# FFT遞迴法的資料請關注it145.com其它相關文章！

C#實現FFT(遞迴法)的範例程式碼

目錄

1. C#實現複數類

2. 遞迴法實現FFT

3. 補充：窗函數

熱門文章