public unsafe static class FourierDNR
{
public static void ForwardTransform(Complex* samples, int length)
{
int nm1 = length - 1; // 变换的点数，也就是算法中的 N ，nm1=N-1,
int nd2 = length / 2; // 算法中N/2
int i, j, jm1, k, l, m, le, le2, ip;
float ur, ui, sr, si, tr, ti;
m = 0;
i = length;
// 功能：通过移位次数计算蝶形算法的阶数m，
// 条件：i是2的m次方，
// 例 ：i=2^3=8,循环3次,所以m=3，三阶蝶形算法
while (i > 1)
{
++m;
i = (i >> 1);
}
j = nd2;
// 该for 循环实现倒位序；
// 时域抽取算法中，输入序列需要倒位序，输出序列是正序
// 频域抽取算法中，输入序列是正序，输出序列需倒位序
// tr、ti : 将输入序列看做复数，tr、为实部，ti是虚部；
// samples : 是输入序列缓冲区，在数组内部实现倒位序
// 注意：从循环可以看出，i=0,i=nm1没有参与操作，也就是第一点和最后一点不参与倒序
for (i = 1; i < nm1; ++i)
{
if (i < j) // j=N/2
{
tr = samples[j].Real;
ti = samples[j].Imag;
samples[j].Real = samples[i].Real;
samples[j].Imag = samples[i].Imag;
samples[i].Real = tr;
samples[i].Imag = ti;
}
k = nd2; // k=N/2
while (k <= j)
{
j = j - k;
k = k / 2;
}
j += k;
}
// 以下为蝶形算法过程
// l : 表示第几阶蝶形运算
// m : 上面求出的总共的级数或阶数，若N=8,则m=3
for (l = 1; l <= m; ++l) // 第一级循环，代表当前是第几级蝶运算；若N=8,则m=3,=>总共循环3次
{
le = 1 << l;
le2 = le / 2; // le2：节点间距离，也就是本级旋转因子的个数
ur = 1; // Wn0实部，也是中间变量
ui = 0; // Wn0虚部，也是中间变量
sr = (float) Math.Cos(Math.PI / le2); // 旋转因子实部，每一级的旋转解度不同，le2决定了旋转角度
si = (float) -Math.Sin(Math.PI / le2); // 旋转因子虚部，每一级的旋转解度不同，le2决定了旋转角度
for (j = 1; j <= le2; ++j) // 第二级循环，代表当前是第几个旋转因子，
{
jm1 = j - 1;
for (i = jm1; i <= nm1; i += le) // 第三级循环，代表同一级同一个旋转因子作的运算
{
ip = i + le2; // i 是参与蝶形运算上方节点的序号，+le2（距离）后，=ip,即下方节点的序号
tr = samples[ip].Real * ur - samples[ip].Imag * ui; // 蝶形运算下方的值*旋转因子
ti = samples[ip].Real * ui + samples[ip].Imag * ur;
samples[ip].Real = samples[i].Real - tr; // 下方为减
samples[ip].Imag = samples[i].Imag - ti;
samples[i].Real = samples[i].Real + tr; // 上方为加
samples[i].Imag = samples[i].Imag + ti;
}
// 以下三名是从当前旋转因子，以旋转生成下一个旋转因子
tr = ur;
ur = tr * sr - ui * si;
ui = tr * si + ui * sr;
}
}
}
public static void InverseTransform(Complex* samples, int length)
{
for (int i = 0; i < length; i++)
{
samples[i].Imag = -samples[i].Imag; // 输入序列虚部取反，相当于对这个复数序列取共轭
}
ForwardTransform(samples, length); // 做一次正变换
var factor = 1.0f / length; // 计算缩放因子
for (int i = 0; i < length; i++)
{
samples[i].Real = samples[i].Real * factor; // 实部缩放
samples[i].Imag = -samples[i].Imag * factor; // 虚部缩放，并取共轭
}
}
}
评论0