简介一份MATLAB中手动实现快速傅里叶变换FFT的源代码主要面向希望深入理解FFT算法原理、学习数字信号处理或需要定制优化变换逻辑的MATLAB用户可用于课堂教学、课题实验和个人进阶。压缩包内共2个文件包括1个.m源文件和1个.doc说明文档整体大小仅14KB轻量且便于阅读与移植。源码采用分治策略实现基2FFT覆盖蝶形运算、复数运算、位反变换与递归求解等核心环节配套文档则补充了FFT基础知识、性能对比和调试测试建议方便读者对照MATLAB内置fft函数验证结果、排查实现细节。已有184人学习/下载对想摆脱黑盒调用、从底层掌握FFT计算流程的学习者而言是一份直观、完整且可直接运行的入门级参考。1. 手写FFT源代码不调fft()但要看透变换本身“matlab 实现fft变化的源代码这是源代码不是直接调用函数”——搜到这句话的人多半不是不知道fft(x)怎么用而是被“不能用内置函数”卡住了课程设计要求手写面试要你现场推导蝶形或者你需要把算法翻译成C语言、Vivado IP核之外的自有实现。fft(x)一行出结果但一旦要搬到别处位反转、旋转因子、蝶形运算这些概念一个都绕不开。这篇文章从朴素DFT的O(N²)复杂度讲起推导出基2时间抽取FFT的拆分逻辑给出完整的MATLAB手写源码不依赖fft()、bitrevorder等内置函数再讲透验证方法和典型错误最后补上IFFT复用、旋转因子预计算和基4思路。适合课程设计、面试复习以及准备把算法落地到其他平台的工程师新手能跟着走通老手可以跳过推导直接看代码和坑。2. 从DFT到FFT蝶形运算怎么把O(N²)降成O(N log N)2.1 先写一个朴素DFT拿到正确的基准输出N点离散傅里叶变换的教科书定义X(k) Σ x(n)·e^(-j2πnk/N)k和n都从0到N-1。按这个公式直接写两层循环每个频点要遍历全部N个时域采样总共N²次复数乘法。N1024时约100万次还能忍到N8192时约6700万次MATLAB里已经肉眼可见地卡顿。更关键的是大量旋转因子被重复计算FFT就是从这一点下手。先实现一个my_dft作为后续所有验证的对照基准只保证正确不追求效率function X my_dft(x) % MY_DFT 按定义直接计算DFTO(N^2)复杂度 % x: 行向量N点时域输入 % X: 行向量N点频域输出 N length(x); X zeros(1, N); for k 0:N-1 for n 0:N-1 X(k1) X(k1) x(n1) * exp(-1j*2*pi*k*n/N); end end end两个容易忽略的点exp(-1j2pikn/N)里的负号对应正变换方向换成正号就变成了逆变换MATLAB下标从1开始所以公式里的k和n都要做加1偏移。这个偏移在后面的FFT实现里是最大的错误源。2.2 旋转因子的三个性质是拆分的全部依据约定W_N e^(-j2π/N)DFT可以写成X(k) Σ x(n)·W_N^(nk)。W_N有三个基本性质性质表达式对算法的意义周期性W_N^(kN) W_N^k旋转因子只需存N个值对称性W_N^(kN/2) -W_N^k后半段频点直接复用前半段取负即可可约性W_N^(2k) W_(N/2)^k大点DFT可以拆成小点DFT周期性说明W_N只有N个独立取值值得预先算好存表对称性让DFT后半段频点可以直接复用前半段的中间结果可约性则是递归拆分的钥匙它把N点变换里的旋转因子映射到N/2点变换里。这三个性质合起来决定了蝶形运算只算前半段、后半段取负的结构。不理解这一点你只是在抄循环。2.3 奇偶分裂推导蝶形复杂度从N²降到N log N把x(n)按下标奇偶拆成两组偶序列g(r) x(2r)奇序列h(r) x(2r1)r从0到N/2-1。代入DFT定义利用可约性把W_N^(2rk)变成W_(N/2)^(rk)X(k) Σ g(r)·W_(N/2)^(rk) W_N^k·Σ h(r)·W_(N/2)^(rk) G(k) W_N^k·H(k)其中G(k)和H(k)是两个N/2点DFT都以N/2为周期。k取0到N/2-1时算完前半段再用对称性W_N^(kN/2) -W_N^k后半段频点直接得到X(kN/2) G(k) - W_N^k·H(k)这就是基2蝶形两个输入G(k)、H(k)一个旋转因子W_N^k产生两个输出。一分为二、二分为四递归到2点DFT级数深度log2(N)每级N/2个蝶形复数乘法数约(N/2)·log2(N)。N1024时约5120次和朴素DFT的100万次差两个数量级。2.4 位反转是DIT的必然结果不是额外步骤自顶向下看每一级按二进制最低位分奇偶递归到底后输入顺序自然就是位反转。以N8为例递归三次后的叶子序列是x(0), x(4), x(2), x(6), x(1), x(5), x(3), x(7)对应原始下标的二进制逐位反转原始下标二进制位反转重排后下标0000000010011004201001023011110641000011510110156110011371111117位反转不是额外的排序开销而是DIT结构下叶子节点的自然顺序。自底向上的迭代实现必须先按这个顺序排列输入再逐级合并蝶形。理解这一点你在写位反转代码时才不会觉得它是个莫名其妙的预处理。3. 手写基2 DIT-FFTMATLAB完整实现与参数说明3.1 手动位反转不调用bitrevorder的写法bitrevorder是MATLAB内置函数但既然目标是弄懂算法、方便移植最好自己实现。位反转的原理是把原始下标i的二进制逐位翻转最低位搬去变成结果的最高位次低位变次高位以此类推。function y my_bitrevorder(x) % MY_BITREVORDER 手动位反转重排N必须是2的幂 % 不依赖内置bitrevorder便于移植到C或其他语言 N length(x); y x; bits log2(N); for i 0:N-1 j 0; n i; for b 1:bits j j * 2 mod(n, 2); % 当前最低位追加到j n floor(n / 2); % n右移一位 end if j i % 每对位置只交换一次 tmp y(i1); y(i1) y(j1); y(j1) tmp; end end end内层循环每次取出n的最低位mod(n,2)追加到j的末尾j*2等价于左移一位给新位腾位置floor(n/2)则丢弃已处理的位。bitslog2(N)正好是二进制位数。if ji保证每对位置只交换一次同时也跳过了ji这种自反转情况这个条件省去了额外的visited数组。3.2 蝶形主循环级、块、蝶的三层结构蝶形循环需要三个层次的索引控制。外层len表示当前级的半跨度step2*len是完整蝶形跨度中间层k遍历0到len-1对应本级每个旋转因子内层i以step为步长定位所有应用同一旋转因子的蝶形起始点。len 1; while len N step len * 2; w exp(-1j * pi / len); % 本级基础旋转因子 for k 0:len-1 w_k w^k; for i k1:step:N % MATLAB从1开始所以加1 t w_k * x(ilen); x(ilen) x(i) - t; % 下半输出用旧x(i)计算 x(i) x(i) t; % 上半输出仍然基于旧x(i) end end len step; end注意更新顺序必须先算x(ilen)再更新x(i)。MATLAB的赋值是逐语句执行的第一行t用旧x(i)算好第二行x(ilen) x(i) - t里的x(i)还是旧值如果先改x(i)第三行就再也拿不到旧值了。这里的w exp(-1jpi/len)来自W_(2len)^1 e^(-j2π/(2*len)) e^(-jπ/len)与2.3节的推导一一对应。级数、跨度与旋转因子的对照关系级数lenstep蝶形总数k取值范围112N/20224N/20~1348N/20~3……………log2(N)N/2NN/20~N/2-1每级蝶形总数都是N/2区别只在旋转因子个数和每个旋转因子被复用的次数。第一级只有W^01理论上可以省掉乘法但为了结构统一先不优化。3.3 完整函数长度检查、方向兼容与调用方式把位反转和蝶形循环组装成完整函数function X my_fft_dit(x) % MY_FFT_DIT 基2时间抽取FFT完全手写实现 % 输入: x - 向量长度必须是2的幂 % 输出: X - 与输入同方向的频域向量 N length(x); if N 1 || mod(N, 2) ~ 0 || log2(N) ~ round(log2(N)) error(输入长度必须是2的幂当前N%d, N); end isCol iscolumn(x); x x(:).; % 统一为行向量 x my_bitrevorder(x); len 1; while len N step len * 2; w exp(-1j * pi / len); for k 0:len-1 w_k w^k; for i k1:step:N t w_k * x(ilen); x(ilen) x(i) - t; x(i) x(i) t; end end len step; end if isCol X x.; else X x; end endmod(N,2)~0排除奇数log2(N)~round(log2(N))排除非2的幂偶数iscolumn判断输入方向结束时再转回去避免调用方拿到不一致的形状。调用方式N 1024; x randn(1, N); % 任意实或复信号 X my_fft_dit(x); % 与内置fft(x)结果一致X是复数频谱幅度用abs(X)取模如果采样率是fs频率轴刻度按(0:N-1)*fs/N换算。N/2以上是镜像频率画单边谱时取前N/21个点。4. 验证与排错my_fft_dit和内置fft()对比、误差与性能4.1 三路对比N16的误差基准实现写完先别急着上大数据。N16时同时跑my_dft、my_fft_dit和内置fft用最大绝对误差判断正确性N 16; x randn(1, N) 1j*randn(1, N); X_dft my_dft(x); X_fft my_fft_dit(x); X_builtin fft(x); fprintf(dft vs builtin: %e\n, max(abs(X_dft - X_builtin))); fprintf(dit vs builtin: %e\n, max(abs(X_fft - X_builtin))); fprintf(dft vs dit: %e\n, max(abs(X_dft - X_fft)));三路两两对比比只跟内置比更稳两个独立实现同时错在同一个点的概率很低。N16时误差都应该在1e-13量级因为浮点累加顺序不同会有微小差异如果误差到1e-3或直接出现NaN一定是位反转、索引偏移或旋转因子方向的问题不是数值精度。误差量级结论排查方向~1e-14实现正确无需处理~1e-3 到 1e-6精度受损输入是否单精度、N是否过大明显错乱/NaN逻辑错误位反转、索引偏移、旋转因子符号提示误差在1e-13量级就说明算法正确不需要追求和内置fft()完全相同的位级结果。4.2 信号级验证多频率分量的合成信号更直观的验证用已知频率的合成信号。50Hz正弦加120Hz余弦采样率1000Hz取1024点保证两个频率都包含整数个完整周期fs 1000; t (0:1023) / fs; x 2 * sin(2*pi*50*t) 0.5 * cos(2*pi*120*t); X my_fft_dit(x); f (0:511) * fs / 1024; stem(f, abs(X(1:512)));幅度谱在50Hz和120Hz处各有一个峰幅度分别接近2和0.5。因为x长度恰好是整数周期峰值不会展宽。如果你把长度改成1000点同样两个频率会出现频谱泄漏峰变矮、底部拖宽——那是截断效应不是FFT实现的问题验证时别把这两件事搞混。4.3 性能对比手写比内置慢多少是正常的N 4096; x randn(1, N); t1 timeit(() fft(x)); t2 timeit(() my_fft_dit(x)); fprintf(内置fft: %.6f s, 手写: %.6f s, 慢 %.0f 倍\n, t1, t2, t2/t1);内置fft()底层是FFTW的高性能实现还可能多线程纯MATLAB脚本的for循环每轮都要经过解释器和边界检查慢一两个数量级是完全正常的。手写实现的价值在算法结构和参数语义不在速度。真要提速有三个方向把exp计算换成预生成的旋转因子表把最内层的i循环向量化或者改用MATLAB Coder生成C代码。这些都改完蝶形算法本身也没变变的只是运行环境。4.4 三个典型坑非2幂、索引偏移、旋转因子方向第一个坑是输入长度不是2的幂CSV导入的实测数据最常见。解决方法是补零到最近的2的幂data readmatrix(sensor_readings.csv); x data(:, 1); % 取第一列信号 N 2^nextpow2(length(x)); X my_fft_dit([x; zeros(N-length(x), 1)]);nextpow2返回能容纳length(x)的最小2的幂补零不改变谱峰位置只是让频点变密。如果CSV里有多列先确认哪一列是时间、哪一列是信号别把时间戳拿去变换。第二个坑是索引偏移。所有0-based的算法公式搬进MATLAB都要加1。常见错误是把内层循环写成for ik:step:N少了1的偏移N4时也许碰巧还能出结果N8开始错乱而且错误模式很隐蔽。定位办法把N8的输出和fft(x)逐点打印对比第一个不一致的频点通常就能告诉你错在哪一级。第三个坑是旋转因子方向。正变换必须用e^(-j2πnk/N)写成正号就是逆变换。最省事的自检是x[1 0 0 0]N4时正确FFT结果应该全为1如果出现共轭或虚部顺序不对先查exp里的符号。5. 往前一步IFFT复用、旋转因子预计算与基4思路5.1 利用共轭对称一行实现IFFT逆变换不必另写一套蝶形。DFT和IDFT的旋转因子互为共轭所以IDFT可以把正变换函数包一层function x my_ifft_dit(X) % MY_IFFT_DIT 基于my_fft_dit的逆变换 N length(X); x conj(my_fft_dit(conj(X))) / N; end验证方法先用my_fft_dit正向变换再用my_ifft_dit逆变换回来max(abs(x - x_orig))应该在1e-15量级。这一行代码说明正逆变换共用同一套蝶形核只是输入输出各取一次共轭移植到C语言时只需要写一个FFT函数IFFT通过包装实现。5.2 预计算旋转因子表消除重复exp蝶形循环里每级都要算exp(-1j*pi/len)级数多了之后这部分开销不小。常见做法是提前生成一张表每个旋转因子只算一次function w_table build_twiddle_table(N) % BUILD_TWIDDLE_TABLE 预计算各级旋转因子 levels log2(N); w_table cell(levels, 1); len 1; for level 1:levels w_table{level} exp(-1j * pi * (0:len-1) / len); len len * 2; end end使用时在蝶形循环里直接取w w_table{level}(k1)。注意MATLAB索引从1开始旋转因子指数k对应表内下标k1。exp的运算量从每级N/2次降到总共N-1次在纯循环实现里提速非常明显而且这个表后续做频域滤波也能复用。5.3 基4思路与实数优化的方向基4FFT的基本单元是4点DFT每级把N点拆成4个N/4点子序列一个蝶形用三个旋转因子、输出四个中间值复数乘法数从(N/2)·log2(N)降到(3N/8)·log2(N)64点以上大约省25%。代价是位反转规则变成基4反转索引复杂度明显上升只有当目标平台的乘法器特别贵或者需要严格对齐FFT IP核的吞吐时才值得从基2改基4。实数输入是另一个常见场景。因为实信号频谱满足共轭对称X(N-k)conj(X(k))可以把N点实数序列打包成N/2点复数序列做一次N/2点FFT再拆包恢复名义计算量减半。判断拆包公式是否写对的快速标准是拆包后的X(k)必须满足X(N-k)conj(X(k))不满足就一定是哪里取错了共轭或下标。学习阶段直接用my_fft_dit处理实数序列就够先确保主链路正确再考虑打包拆包这种优化。本文还有配套的精品资源点击获取
