FFT---快速傅里叶变换
最近在学习dsp数字信号处理算法,FFT是避免不了要学习的,本文将用实例理解FFT的实现。
FFT是DFT的快速计算方法,所以先简单的复习一下DFT的基本原理。
DFT:
式中,,称为蝶形因子。上面的式子是N点的DFT,需要做的乘法是N^2,FFT算法是利用了蝶形因子WN的对称性和周期性,将乘法的运算量减少到N^2 /2。
蝶形因子的三个性质 :
周期性:
对称性: 和
可约性: 和
FFT原理推导
以N=8点的FFT计算为例。
注意:
1.第一层L=1时,N=2,;第二层L=2时,N=4;第三层L=3时,N=8.
2.红线表示“+”,蓝线表示“-”.
3. i表示从左边出发原始的序列标号,j表示第j层。比如我标记绿色三角形,
,表示x[1]的第二层。
4.最左边的x是小写的,右边的X是大写的.
实例-----以计算X[1]为例
DFT的算法:每次算N个点
FFT的算法:每次算N/2个点
第一层:
第二层:
第三层:
然后通过递归迭代,把第二层带入第三层,把第一层带入第二层,最终算的X[1]。相当于把问题分解成了众多子问题,用算法描述就是T[N]=2*T[N/2]+cN,它的时间复杂度下降了,从而实现FFT的快速计算。
C语言实现:
#include <stdio.h> #include <stdlib.h> #include <string.h> #include <math.h> #define PI 3.1415926535 void kfft(pr, pi, n, k, fr, fi) int n, k; double pr[], pi[], fr[], fi[]; { int it, m, is, i, j, nv, l0; double p, q, s, vr, vi, poddr, poddi; for (it = 0; it <= n - 1; it++) //将pr[0]和pi[0]循环赋值给fr[]和fi[] { m = it; is = 0; for (i = 0; i <= k - 1; i++) { j = m / 2; is = 2 * is + (m - 2 * j); m = j; } fr[it] = pr[is]; fi[it] = pi[is]; } pr[0] = 1.0; pi[0] = 0.0; p = 6.283185306 / (1.0 * n); pr[1] = cos(p); //将w=e^-j2pi/n用欧拉公式表示 pi[1] = -sin(p); for (i = 2; i <= n - 1; i++) //计算pr[] { p = pr[i - 1] * pr[1]; q = pi[i - 1] * pi[1]; s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]); pr[i] = p - q; pi[i] = s - p - q; } for (it = 0; it <= n - 2; it = it + 2) { vr = fr[it]; vi = fi[it]; fr[it] = vr + fr[it + 1]; fi[it] = vi + fi[it + 1]; fr[it + 1] = vr - fr[it + 1]; fi[it + 1] = vi - fi[it + 1]; } m = n / 2; nv = 2; for (l0 = k - 2; l0 >= 0; l0--) //蝴蝶操作 { m = m / 2; nv = 2 * nv; for (it = 0; it <= (m - 1) * nv; it = it + nv) for (j = 0; j <= (nv / 2) - 1; j++) { p = pr[m * j] * fr[it + j + nv / 2]; q = pi[m * j] * fi[it + j + nv / 2]; s = pr[m * j] + pi[m * j]; s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]); poddr = p - q; poddi = s - p - q; fr[it + j + nv / 2] = fr[it + j] - poddr; fi[it + j + nv / 2] = fi[it + j] - poddi; fr[it + j] = fr[it + j] + poddr; fi[it + j] = fi[it + j] + poddi; } } for (i = 0; i <= n - 1; i++) { pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]); //幅值计算 } return; } int main() { int i, j; int N = 16; double pr[16] = { 0.537667139546100, 1.83388501459509, -2.25884686100365, 0.862173320368121, 0.318765239858981, -1.30768829630527, -0.433592022305684, 0.342624466538650, 3.57839693972576, 2.76943702988488, -1.34988694015652, 3.03492346633185, 0.725404224946106, -0.0630548731896562, 0.714742903826096, -0.204966058299775}; double pi[16] = {0.0}, fr[N], fi[N], t[N]; // pr表示采样数据的实部,pi表示采样数据的虚部,fr表示离散傅里叶变换结果的实部,fi表示离散傅里叶变换结果的虚部。 /* for (i=0; i<N; i++) //生成输入信号 { pr[i]=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i); pi[i]=0.0; //依次将值赋值给pr,pi赋值为0 } */ kfft(pr, pi, N,4, fr, fi, 0, 1); //调用FFT函数 for (i = 0; i < N; i++) { // printf("%d\t%lf\n", i, pr[i]); //输出结果 printf("%lf+(%lf)j\n",fr[i],fi[i]); } }
用matlab验证C代码
1 2 | x = [0.537667139546100, 1.83388501459509, -2.25884686100365, 0.862173320368121, 0.318765239858981, -1.30768829630527, -0.433592022305684, 0.342624466538650, 3.57839693972576, 2.76943702988488, -1.34988694015652, 3.03492346633185, 0.725404224946106, -0.0630548731896562, 0.714742903826096, -0.204966058299775]; y = fft (x) |
发现结果一致。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· DeepSeek R1 简明指南:架构、训练、本地部署及硬件要求
· NetPad:一个.NET开源、跨平台的C#编辑器