快速傅里叶算法 C语言实现

考研阶段学习过傅里叶级数，而最近的项目正好是用C语言编写傅里叶变换，于是很认真的复习了傅里叶级数。可是无奈，看来看去反而晕晕乎乎的。后经师兄师姐的指教，才得知对于工程中的信号处理，研究周期性的傅里叶变换都没有现实意义，而傅里叶级数更没有什么关系。
工程中待处理的信号，通常具有非周期性，故我们需要对离散傅里叶变换进行研究。离散公式：

【x(n)是采样的时域信号，X(k)是对于不同频率k的频域信号。】

而快速傅里叶变换又是对离散傅里叶变换的改进，通过蝶形运算（网上的图片如下），计算速度可大大提升，使计算量呈指数型下降。

【最左的x(n)是采样的时域信号，最右的X(k)是算出的频域信号。可以看到左边的x(n)中，n 的序列并非正常的递增，这是为了使得输出的X(k)中的k频率呈递增序列。[n]的序列是通过将十进制n转化成的二进制数字后，倒序排列生成的，如4的二进制码为100，倒序为001，故排在第二位。在上图中，共8个信号，可分成3级。第一级，每个分组两个节点，一个蝶形单元。第二级，每个分组四个节点，两个蝶形单元。第三级，每个分组八个节点，四个蝶形单元。】

核心代码：

[html] view plain copy

print?

<span style="font-size:12px;">void fft(complex f[], int N)
{
complex t, wn;//中间变量
int i, j, k, la, lb, lc, l, r, m, n;//中间变量
//M:分解运算级数
int M;
/*----计算分解的级数M=nextpow2(N)----*/
for (i = N, M = 1; (i = i / 2) != 1; M++);
/*----输入信号二进制码位倒序*/
for (i = 1, j = N / 2; i <= N - 2; i++)
{
if (i<j)
{
t = f[j];
f[j] = f[i];
f[i] = t;
}
k = N / 2;
while (k <= j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
/*----FFT算法----*/
for (m = 1; m <= M; m++)
{
la = pow(2.0, m); //la=2^m代表第m级每个分组所含节点数
lb = la / 2; //lb代表第m级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离
/*----碟形运算----*/
for (l = 1; l <= lb; l++)
{
r = (l - 1)*pow(double(2), M - m);
for (n = l - 1; n<N - 1; n = n + la) //遍历每个分组，分组总数为N/la
{
lc = n + lb; //n,lc分别代表一个碟形单元的上、下节点编号
Wn_i(N, r, &wn, 1);//wn=Wnr
c_mul(f[lc], wn, &t);//t = f[lc] * wn复数运算
c_sub(f[n], t, &(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
c_plus(f[n], t, &(f[n]));//f[n] = f[n] + f[lc] * Wnr
}
}
}
}</span>

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。