声明:
本人在校期间主修过《数字信号处理》这门课程
对离散傅里叶变换(DFT)和快速傅里叶变换(FFT)深有了解
现编写了基于C语言的FFT算法,已完成对抽样序列的FFT变换并通过窗口输出。
编写思路:
由于FFT变换里面含有对虚数的运算,现将输出序列定义成一个结构体函数
结构体函数里面包含输出序列的实部和虚部,并定义处理复数运算的几种函数
快速傅里叶变换的核心在于倒序运算和蝶形运算
现将倒序运算和蝶形运算的核心思想用代码的形式表现出来
代码部分
## 倒序代码:
for (int I = 1; I < N1; I++) //定义倒序序列函数
{
if (I < J)
{
T = x[I];
x[I] = x[J];
x[J] = T;
}
K = LH;
if (J >= K)
{
do
{
J = J - K;
K = K / 2;
} while (J >= K);
}
J = J + K;
}
## 倒序代码解释:
用J表示当前倒序数的十进制数值。
对于N=2^M,M位二进制数最高位的十进制权值为N/2,且从左向右二进制位的权值依次为N/4,N/8,……,2,1。
因此,最高位加一相当于十进制运算J+N/2。如果最高位是0(J<N/2),则直接由J+N/2得下一个倒序值;
如果最高位是1(J>=N/2),则先将最高位变成0(J=J-N/2),然后次高位加1(J+N/4)。
但次高位加1时,同样需要判断0、1值,如果为0(J<N/4),则直接加1(J=J+N/4),否则将次高位变成0(J=J-N/4),再判断下一位;
依次类推,直到完成最高位加1,逢2向右进位的运算
## FFT代码:
for (L = 1; L <= M; L++) //FFT运算
{
B = (int)(pow(2, L - 1));
for (J = 0; J < B; J++)
{
P = (int)(J*pow(2, M - L));
for (k = J; k < N; k = (int)(k + pow(2, L)))
{
K1 = k + B;
complex wn, t;
Wn_i(N, P, &wn);
c_mul(f[K1], wn, &t); //。。。。。。。。。。。。
c_sub(f[k], t, &(f[K1])); //蝶形运算
c_plus(f[k], t, &(f[k])); //。。。。。。。。。。。。
}
}
}
## FFT代码解释:
第L级中,每个蝶形的两个输入数据相距B=2^(L-1)个点
每级有B个不同的旋转因子;同一旋转因子对应着间隔为2^L点的2^(M-1)个蝶形
先从输入端(第1级)开始,逐级进行,共进行M级运算
在进行第L级运算时,依次求出B个不同的旋转因子,每求出一个旋转因子,就计算完它对应的所有2^(M-L)个蝶形
## 复数函数定义与运算代码部分:
void c_plus(complex a, complex b, complex *c) //复数加法
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
void c_sub(complex a, complex b, complex *c) //复数减法
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
void c_mul(complex a, complex b, complex *c) //复数乘法
{
c->real = a.real*b.real - a.imag*b.imag;
c->imag = a.real*b.imag + a.imag*b.real;
}
void Wn_i(int n1, int i, complex *Wn) //定义FFT旋转因子
{
Wn->real = cos(2 * PI*i / n1);
Wn->imag = -sin(2 * PI*i / n1);
}
## 完整代码实现的功能:
1.基础功能:
- 可自定义序列的抽样点数
- 对序列进行FFT和DFT运算输出
2.拓展功能:
- 在定义抽样点的基础上比较DFT运算和FFT运算的时间