FFT算法理解与实现

快速傅里叶变换（一）
参考这本书理解FFT基本概念
FFT是DFT（discrete Fourier Transform）一种快速算法。
傅里叶变换： 任何连续的周期信号可以分解成许多正弦波与余弦波。

为什么分解成正余弦波？
因为正余弦波输入一个系统，输出仍然是正余弦波，只有幅度与相位发生变化，频率与波形保持不变，更容易分析信号特征。

傅里叶变换可以分为四种：
1.非周期连续信号———->傅里叶变换（FourierTransform）
例如：高斯分布曲线
2. 周期连续信号———–>傅里叶级数 （Fourier Series）
例如：正弦波，方波
3. 非周期离散信号———>离散时间域的傅里叶变换（Discrete TimeFourierTransform）DTFT
4. 周期离散信号———->离散傅里叶级数（DiscreteFourierSeries）但一般都叫（Discrete Fourier Transform）DFT

傅里叶变换分类比较容易弄混，记住就好了。
实际用于DSP的只有DFT（离散傅里叶变换）
DFT 基本方程：
$c_K[i]=cos(2πKi/N),(i\in[0,N-1]);$
$s_K[i]=sin(2πKi/N),(i\in[0,N-1]);$
$k\in[0,N/2];$
$k$代表频率，$c_K[i],s_K[i]$分别表示余弦波、正弦波。
由下面公式

$$ReX[k]=\sum_{i=0}^{N-1}x[i]cos(2πki/N)$$
$$ImX[k]=-\sum_{i=0}^{N-1}x[i]sin(2πki/N)$$
$k\in[0,N/2];$
可以实现DFT算法

        //1 real inverse_DFT ,[in]frequency domain realARR;
        //2. [in] frequency domain imagine maguide
        //3. [out] y synthesis time domain signal
        public static void I_DFT(double[] X_R, double[] X_I, out double[] y)
        {
            int N = (X_R.Length - 1) * 2;/////N 点时域信号。
            double[] yy = new double[N];//X_R[] Frequency  real Part, X_I[] Frequency Im part.
            for (int i = 0; i <= N / 2; i++)
            {
                X_R[i] = X_R[i] / (N / 2);
                X_I[i] = -X_I[i] / (N / 2);

            }
            X_R[0] = X_R[0] / 2;
            X_R[N / 2] = X_R[N / 2] / 2;

            for (int k = 0; k <= N / 2; k++)
            {
                for (int i = 0; i < N; i++)
                {
                    yy[i] = yy[i] + X_R[k] * Math.Cos(2 * Math.PI * k * i / N);
                    yy[i] = yy[i] + X_I[k] * Math.Sin(2 * Math.PI * k * i / N);
                }
            }
            /////////////////////////////////////////////////
            y = yy;
        }
        public static void DFT(double[] xx, out double[] REX, out double[] IMX)
        {
            int N = xx.Length;
            double[] temp_re = new double[N / 2 + 1];
            double[] temp_im = new double[N / 2 + 1];

            for (int k = 0; k <= N / 2; k++)
            {
                for (int i = 0; i < N; i++)
                {
                    temp_re[k] = temp_re[k] + xx[i] * Math.Cos(2 * Math.PI * k * i / N);
                    temp_im[k] = temp_im[k] - xx[i] * Math.Sin(2 * Math.PI * k * i / N);

                }
            }
            REX = temp_re;
            IMX = temp_im;
        }