16点DFT的FFT算法

FFT（快速傅里叶变换）是DFT的一种特殊情况，就是当运算点的个数是2的整数次幂的时候进行的运算（不够用0补齐）。

FFT计算原理及流程图：

原理：FFT的计算要求点数必须为2的整数次幂，如果点数不够用0补齐。例如计算{2，3，5，8，4}的16点FFT，需要补11个0后进行计算。FFT计算运用蝶形运算，在蝶形运算中变化规律由W(N, p)推导，其中N为FFT计算点数，J为下角标的值。

L = 1时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0

L = 2时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1

L = 3时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3

所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1

又因为2^L = 2^M*2^(L-M) = N*2^(L-M)，这里N为2的整数次幂，即N=2^M，

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以，p = J*2^(M-L)，此处J = 0, 1, ..., 2^(L-1)-1，当J遍历结束但计算点数不够N时，J=J+2^L，后继续遍历，直到计算点数为N时不再循环。

流程图：

实现代码：

/*====================================================================== 

 * 方法名：  fft 

 * 方法功能：计算数组的FFT，运用蝶形运算 

 * 

 * 变量名称： 

 *          yVector   - 原始数据 

 *          length    - 原始数据长度 

 *          N         - FFT计算点数 

 *          fftYreal  - FFT后的实部 

 *          fftYImg   - FFT后的虚部 

 * 

 * 返回值：  是否成功的标志，若成功返回true，否则返回false 

 *=====================================================================*/  

  

+ (BOOL)fft:(floatfloat *)yVector andOriginalLength:(NSInteger)length andFFTCount:(NSInteger)N andFFTReal:(floatfloat *)fftYReal andFFTYImg:(floatfloat *)fftYImg  

{  

    // 确保计算时时2的整数幂点数计算  

    NSInteger N1 = [self nextNumOfPow2:N]  

      

    // 定义FFT运算是否成功的标志  

    BOOL isFFTOK = false  

      

    // 判断计算点数是否为2的整数次幂  

    if (N != N1)  

    {  

        // 不是2的整数次幂，直接计算DFT  

        isFFTOK = [self dft:yVector andOriginalLength:length andFFTCount:N andFFTReal:fftYReal andFFTYImg:fftYImg]  

          

        // 返回成功标志  

        return isFFTOK  

    }  

      

      

    // 如果计算点数位2的整数次幂，用FFT计算，如下  

    // 定义变量  

    float yVectorN[N1]             // N点运算的原始数据  

    NSInteger powOfN = log2(N1)    // N = 2^powOfN，用于标记最大运算级数(公式中表示为：M)  

    NSInteger level = 1            // 运算级数（第几次运算），最大为powOfN，初值为第一级运算(公式中表示为：L)  

    NSInteger lineNum              // 行号，倒序排列后的蝶形运算行号(公式中表示为：k)  

    float inverseOrderY[N1]        // yVector倒序x  

    NSInteger distanceLine = 1     // 行间距，第level级运算每个蝶形的两个节点距离为distanceLine = 2^(L-1)(公式中表示为：B)  

    NSInteger p                    // 旋转因子的阶数，旋转因子表示为 W(N, p)，p = J*2^(M-L)  

    NSInteger J                    // 旋转因子的阶数，旋转因子表示为 W(2^L, J)，J = 0, 1, 2,..., 2^(L-1) - 1 = distanceLine - 1  

    float realTemp, imgTemp, twiddleReal, twiddleImg, twiddleTheta, twiddleTemp = PI_x_2/N1  

    NSInteger N_4 = N1/4  

      

    // 判断点数是否够FFT运算点数  

    if (length <= N1)  

    {  

        // 如果N至少为length，先把yVector全部赋值  

        for (NSInteger i = 0 i < length i++)  

        {  

            yVectorN[i] = yVector[i]  

        }  

          

        if (length < N1)  

        {  

            // 如果 N > length 后面补零  

            for (NSInteger i = length i < N1 i++)  

            {  

                yVectorN[i] = 0.0  

            }  

        }  

    }  

    else  

    {  

        // 如果 N < length 截取相应长度的数据进行运算  

        for (NSInteger i = 0 i < N1 i++)  

        {  

            yVectorN[i] = yVector[i]  

        }  

    }  

      

    // 调用倒序方法  

    [self inverseOrder:yVectorN andN:N1 andInverseOrderVector:inverseOrderY]  

  

    // 初始值  

    for (NSInteger i = 0 i < N1 i++)  

    {  

        fftYReal[i] = inverseOrderY[i]  

        fftYImg[i] = 0.0  

    }  

      

    // 三层循环  

    //     第三层（最里）：完成相同旋转因子的蝶形运算  

    //     第二层（中间）：完成旋转因子的变化，步进为2^level  

    //     第一层（最外）：完成M次迭代过程，即计算出x(k) = A0(k), A1(k),...,Am(k) = X(k)  

      

    // 第一层循环  

    while (level <= powOfN)  

    {  

        distanceLine = powf(2, level - 1) // 初始条件 distanceLine = 2^(level-1)  

        J = 0  

        NSInteger pow2_Level = distanceLine * 2   // 2^level  

        NSInteger pow2_NSubL = N1/pow2_Level      // 2^(M-L)  

          

        // 第二层循环  

        while (J < distanceLine)  

        {  

            p = J * pow2_NSubL  

            lineNum = J  

            NSInteger stepCount = 0 // J运算的步进计数  

              

            // 求旋转因子  

            if (p==0)  

            {  

                twiddleReal = 1.0  

                twiddleImg = 0.0  

            }  

            else if (p == N_4)  

            {  

                twiddleReal = 0.0  

                twiddleImg = -1.0  

            }  

            else  

            {  

                // 计算尤拉公式中的θ  

                twiddleTheta = twiddleTemp * p  

                  

                // 计算复数的实部与虚部  

                twiddleReal = cos(twiddleTheta)  

                twiddleImg = -11 * sin(twiddleTheta)  

            }  

              

            // 第三层循环  

            while (lineNum < N1)  

            {  

                // 计算下角标  

                NSInteger footNum = lineNum + distanceLine  

                  

                /*--------------------------------------- 

                 * 用复数运算计算每级中各行的蝶形运算结果 

                 * X(k) = X(k) + X(k+B)*W(N,p) 

                 * X(k+B) = X(k) - X(k+B)*W(N,p) 

                 *---------------------------------------*/  

                realTemp = fftYReal[footNum] * twiddleReal - fftYImg[footNum] * twiddleImg  

                imgTemp  = fftYReal[footNum] * twiddleImg + fftYImg[footNum] * twiddleReal  

                  

                // 将计算后的实部和虚部分别存放在返回数组中  

                fftYReal[footNum] = fftYReal[lineNum] - realTemp  

                fftYImg[footNum]  = fftYImg[lineNum] - imgTemp  

                fftYReal[lineNum] = fftYReal[lineNum] + realTemp  

                fftYImg[lineNum]  = fftYImg[lineNum] + imgTemp  

                  

                stepCount += pow2_Level  

                  

                // 行号改变  

                lineNum = J + stepCount  

            }  

              

            // 旋转因子的阶数变换，达到旋转因子改变的效果  

            J++  

        }  

          

        // 运算级数加一  

        level++  

    }  

      

    isFFTOK = true  

    return isFFTOK  

}

Attempt to execute SCRIPT lab_3_fft_4_16 as a function.

你把脚本用成了函数，

你把文件名改一下，并且文件要放在MATLAB工作目录。

---