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工作目录。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）