C#代码实现DFT与FFT

一、简介

该文章用于记录C#代码实现DFT和FFT的过程。

二、代码步骤

1. 横坐标换算

横坐标需要几个公式。
时域数据长度$T=(N-1)\times △T$，$△T$为采样间隔，频率域的频率间隔$△f=1/T$,频域中的第r个点就代表频率为$f=r\times △f$或者$f=r/T$。如果知道采样周期为一个周期则$T=W/2Pi（W=2\ast PI\ast f）$。

2. 求WN

首先$WN$的公式为：
$$W_N^K=e^{-j\frac{2Pi}{N}K}$$
并且：
$$e^{-j\frac{2Pi}{N}K}=\cos (\frac{2PiK}{N})-j\sin (\frac{2PiK}{N})$$
于是就可以预先求出全部$WN$的实数值与虚数值进行备用，函数代码如下。

public void MN_Init()
{
    Count = X.Count;
    List<double> MN_r = new List<double>();
    List<double> MN_i = new List<double>();
    Dictionary<string, List<double>> WN_placehold = new Dictionary<string, List<double>>();
    for (int i = 0; i < Count; i++)
    {
        MN_r.Add(Math.Cos(2 * Math.PI * i / Count));
        MN_i.Add(-Math.Sin(2 * Math.PI * i / Count));
    }
    WN_placehold.Add("R", MN_r);
    WN_placehold.Add("I", MN_i);
    WN = WN_placehold;
}

3. 进行码位倒序

代码实现，这里如果数组Count长度不为2的N次，便用0填充。

public void aaa(ref List<int> array)
{
    array = new List<int>();
    int iii = (int)Math.Log((double)Count,2);
    int N =0;
    if (iii % 1 > 0)
    {
        iii++;
        N = (int)Math.Pow(2, iii);
    }
    else N = (int)Math.Pow(2, iii);
    for(int ii=0;ii<Count;ii++)
    {
        array.Add(ii);
    }
    int i, j, k;
    int temp;
    j = 0;
    for (i = 0; i < N - 1; i++)
    {
        if (i < j)
        {
            temp = array[i];
            array[i] = array[j];
            array[j] = temp;
        }
        k = N >> 1;
        while (k <= j)
        {
            j = j - k;
            k >>= 1;
        }
        j = j + k;
    }
}

4. DFT

DFT公式如下：
$$X(k)=\sum _{n=0}^{N-1}X(n)\ast e^{-j\frac{2Pi}{N}kn}$$
所以直接计算每项的实数计算还有虚数计算然后最后做总和就可以了，这里也没什么好说的，直接套公式就好

public void output_2(int n, ref double output_r, ref double output_i)
{
    List<double> odd_r = new List<double>();
    List<double> odd_i = new List<double>();
    double odd_r_t = 0;
    double odd_i_t = 0;
    for (int i = 0; i < Count; i++)
    {
        double WN_R = Math.Cos(2 * Math.PI * (i * n) / Count);// WN["R"][(i * 2 * n) % Count];
        double WN_I = Math.Sin(2 * Math.PI * (i * n) / Count);

        odd_r_t += (X[i] * WN_R + Y[i] * WN_I);
        odd_i_t += (Y[i] * WN_R - X[i] * WN_I);
        odd_r .Add(X[i] * WN_R + Y[i] * WN_I);
        odd_i .Add( Y[i] * WN_R- X[i] * WN_I);
    }
    output_r = odd_r_t;
    output_i = odd_i_t;

    GC.Collect();
}

5. FFT

首先需要寻找蝶形计算的规律。

总体可以分三个循环步骤：
第一个循环为总体级数，如下图红框部分，级数为$k=lo{g}_{2}N$

第二个循环为每级中的蝶形数量，第一级可以看出有4个，计数每次加$2^k$ $(k=0,1,\mathrm{2...}N)$。

第三个循环为每个蝶形中的计数次数，这里需要区分上下区，上区为相加，下区为相减。计数为$2^{k+1}$ $(k=0,1,\mathrm{2...}N)$。

这样三个循环条件都知道，就可以开始写程序了。三个大体循环代码如下

 for (int k=0;k<Math.Log(Count,2);k++)
 {
      REAL_ = new List<double>(REAL.ToArray());
      IMAGIN_ = new List<double>(IMAGIN.ToArray());
      for (int j=0;j<Count;j+=(int)Math.Pow(2,k+1))
      {
          int pow_for = (int)Math.Pow(2, k+1);
          for (int i=0;i<pow_for;i++)
          {
          }
      }
}

这里还需要一条公式：
$$W_N^K=W_{N/2}^{K/2}$$
所以每级$W_N^K$中$K$的数量为$2^n$，其中$n$为目前的级数$0,1,2,...,N$，间隔为$N/2^{n+1}$，第一级就是$W_N^0$，第二级就是$W_N^0、W_N^2$，如此类推。
然后再注意一下上下区对应的数值就可以了。下面是整体循环的代码部分。

for (int k=0;k<Math.Log(Count,2);k++)
{
     REAL_ = new List<double>(REAL.ToArray());
     IMAGIN_ = new List<double>(IMAGIN.ToArray());
     for (int j=0;j<Count;j+=(int)Math.Pow(2,k+1))
     {
         int pow_for = (int)Math.Pow(2, k+1);
         int pow = (int)Math.Pow(2, k);
         for (int i=0;i<pow_for;i++)
         {
             double REAL_D, IMAGIN_D;
             if(i<pow)//上层使用加法
             {
                 REAL_D = REAL[j + i + pow] * WN["R"][Count / (pow_for) * i] - WN["I"][Count / (pow_for) * i] * IMAGIN[j + i + pow];
                 IMAGIN_D = IMAGIN[j + i + pow] * WN["R"][Count / (pow_for) * i] + WN["I"][Count / (pow_for) * i] * REAL[j + i + pow];

                 REAL_[j + i] += REAL_D;
                 IMAGIN_[j + i] += IMAGIN_D;
                 //X[j + i] = X[j + i] + X[j + i +1] * W[Count / (pow_for) * i];
             }
             else//下层使用减法
             {
                 REAL_D = -(REAL[j + i] * WN["R"][Count / (pow_for) * (i - pow)] - WN["I"][Count / (pow_for) * (i - pow)] * IMAGIN[j + i]);
                 IMAGIN_D =-( IMAGIN[j + i] * WN["R"][Count / (pow_for) * (i - pow)] + WN["I"][Count / (pow_for) * (i - pow)] * REAL[j + i ]);

                 REAL_[j + i] = REAL[j + i-pow]+ REAL_D;
                 IMAGIN_[j + i] = IMAGIN[j + i - pow]+IMAGIN_D;
                 // X[j + i] = X[j + i] - X[j + i] * W[Count / (pow_for) * i];
                 //X[j+i] = X[j + i - pow] - 
             }
         }
     }
     REAL = REAL_;
     IMAGIN = IMAGIN_;
 }