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 **fftok = false; // 判断计算点数是否为2的整数次幂 if (n != n1) { // 不是2的整数次幂，直接计算dft **fftok = [self dft:yvector andoriginallength:length andfftcount:n andfftreal:fftyreal andfftyimg:fftyimg]; // 返回成功标志 return **fftok; } // 如果计算点数位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 d**tanceline = 1; // 行间距，第level级运算每个蝶形的两个节点距离为d**tanceline = 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 = d**tanceline - 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) { d**tanceline = powf(2, level - 1); // 初始条件 d**tanceline = 2^(level-1) j = 0; nsinteger pow2_level = d**tanceline * 2; // 2^level nsinteger pow2_nsubl = n1/pow2_level; // 2^(m-l) // 第二层循环 while (j < d**tanceline) { 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 + d**tanceline; /*--------------------------------------- * 用复数运算计算每级中各行的蝶形运算结果 * 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++; } **fftok = true; return **fftok; } 20210311