简介这是一份基于C编写的二维快速傅里叶变换程序源码面向数字信号处理与图像处理初学者、相关课程设计与算法验证开发者。程序采用对图像逐行、逐列进行一维傅里叶变换的经典实现方式将空间域数据转换到频域低频成分对应图像整体变化高频成分对应边缘与细节可在此基础上开展低通滤波、高通滤波、去噪、频谱分析等实验。压缩包内仅含一个C源文件压缩后体积约1KB大小代码轻量易读既可直接编译运行也可提取核心函数用于二次开发便于深入理解傅里叶变换从理论到编码的落地过程。读者可自行准备图像或矩阵数据运行后观察二维频谱的幅度与相位信息比较不同频段处理前后的差异从而掌握频域分析的基本方法。目前已有299人学习/下载适合作为图像频域处理入门或课堂演示的参考资料。1. 二维傅里叶变换到底在算什么解压 FFT 包前先想清楚的四件事不管你是从图像处理转过来还是做光学测量、雷达回波处理碰到“二维傅里叶变换”这个词的频率不会低。和一位按列推进的一维 FFT 不同二维 FFT 处理的是二维采样网格输入是一个 M×N 矩阵输出是一个同样大小的复数矩阵每个输出点代表一个二维空间频率分量。标题里的 FFT.rar 通常就属于这类集合里面可能是一段 MATLAB 仿真也可能是给 FPGA 用的 C 代码关键点都落在 fft2 怎么拆、频谱怎么搬移、频域滤波怎么做这三件事上。这篇文章按我实际复现这类工程的顺序来写先立住二维 FFT 的数理骨架再用 MATLAB 把 CSV 数据导入后跑通最小样例接着讲窗函数、补零、谱域滤波器的参数边界最后给一套能放进嵌入式或 Vivado FFT 核流程的自检方法。适合两种人第一次接触二维频谱图、需要快速出结果的新手以及想在现有工程里减少无效调试的老手。2. 二维 FFT 的数学原理与频谱布局从 fft 到 fft2二维离散傅里叶变换的定义式并不复杂F(u,v) Σ_{x0}^{M-1} Σ_{y0}^{N-1} f(x,y) · exp[-2πj(ux/M vy/N)]难点在于怎么把它转换成高效计算以及算完之后怎么读懂输出矩阵。下面两节分别解决这两个问题。2.1 可分离性是把二维 FFT 落地的第一把钥匙观察定义式里指数项它能拆成 exp(-2πjux/M) 乘以 exp(-2πjvy/N)。这意味着先对每一列做一维 FFT再对得到的结果按每一行做一维 FFT和直接算二维 FFT 完全一致。这就是可分离性也是 FFT.rar 里那些 C 实现最常见的套路一个一维蝴蝶运算函数两个方向的循环调用就能完成二维变换。% 先对每列做一维 FFT再对每行做一维 FFT与内置 fft2 等价 colTransform fft(data, [], 1); % 沿第 1 维列变换 rowTransform fft(colTransform, [], 2); % 沿第 2 维行变换 isequal(rowTransform, fft2(data)) % 理论为 true浮点误差可忽略这段代码里第一个参数[]表示变换长度沿用该维原始长度第二个参数指定方向。fft2(data)实际内部顺序也是先沿一维再沿另一维只是对用户隐藏了细节。理解这层拆解的意义在于排错当二维结果出现“只有横条纹正常、纵条纹不对”的现象时问题多半出在某一维的数据读取 stride 上而不是算法本身。算法复杂度是 O(M·N·log(M·N))。朴素二维 DFT 是 O(M²·N²)像素为 1024×1024 时两者差距在万倍以上所以二维 FFT 在图像和相控阵处理里才具备实时性。2.2 频谱布局与 fftshift为什么零频不在矩阵中心一维 FFT 的输出前两点是直流和最低正频率分量二维 fft2 沿用同样约定F(0,0) 位于输出矩阵的左上角。但显示频谱图时大家习惯把低频放中间、高频放外侧于是必须有一步“频谱搬移”。MATLAB 函数fftshift把四象限对调让零频出现在矩阵中心。F fft2(data); Fshifted fftshift(F); % 把低频移到中心这里有一个常年出现的误区fftshift不是 FFT 的一部分它只改变数据的陈列顺序。如果你从 FFT.rar 里拿到一段自行实现的代码大概率会在输出后做同样操作对应到 C 语言里就是下标按 N/2 做环形移位。另一个相关函数是ifftshift它的作用恰好相反用于把中心化的频域数据恢复成 FFT 的输出顺序。注意fftshift和ifftshift在空间点数 N 为偶数时行为一致N 为奇数时函数实现不同所以不要混用。构造物理频率轴时我通常用下面这段代码它同时处理奇偶尺寸dx 0.1; % 空间采样间隔单位自定 dfx 1 / (N * dx); % 频率分辨率 dfy 1 / (M * dx); % 假设两个方向采样间隔相同 fx (-floor(N/2) : ceil(N/2)-1) * dfx; fy (-floor(M/2) : ceil(M/2)-1) * dfy;参数dfx是频谱图上相邻两个频率点的间隔它只由采样点数 N 和采样间隔 dx 决定与窗函数、补零长度都无关。这个关系在第四章还会用到。下面表格列出不同尺寸下频率轴的行列对应关系帮助你快速核对。矩阵尺寸未搬移的输出顺序fftshift 后中心位置频率轴范围N 为偶数DC 在 (0,0)Nyquist 在 (N/2)DC 在 (N/2)Nyquist 在 (N-1)-Fs/2 到 Fs/2-Fs/NN 为奇数DC 在 (0,0)无准确 NyquistDC 在 (N-1)/21-(N-1)/2·df 到 (N-1)/2·dfFFT 的 Xilinx IP 核、Stm32 DSP 库返回的都是一维输出搬到二维时同样要按这个规则计算每个 bin 对应的物理频率。3. 从 CSV 到频谱图在 MATLAB 里把二维 FFT 仿真跑通这一章给出可以直接复制运行的代码路径覆盖最小样例、真实 CSV 数据导入和复数数据三种常见输入形式。我默认你用的是 MATLAB R2019 之后版本因为readmatrix在这个版本开始稳定支持自动列类型识别。3.1 MATLAB 做二维 FFT 的最小样例最小验证不需要真实数据用单点冲激和余弦光栅就够了。先看冲激响应它能帮你确认频谱搬移和坐标映射是否正确。data zeros(128, 128); data(64, 64) 1; % 空间域的冲激位于矩阵中心 F fftshift(fft2(data)); imagesc(abs(F).^2); % 功率谱 axis image; colorbar; title(Point Source Spectrum);理论上冲激函数的傅里叶变换是一个常数幅度谱因此上面代码画出的abs(F).^2应该是全平面均匀的。实际画面里如果出现斜向条纹通常是因为fftshift用了一次但冲激没有放对位置或者数据被circshift干扰过。冲激验证通过后再构造一个已知频率的二维正弦用来验证频率轴刻度t (0:127) / 128; % 归一化坐标共 128 个点 [X, Y] meshgrid(t); % 生成二维网格 image2D cos(2*pi*8*X 2*pi*3*Y); % 水平方向 8 个周期垂直方向 3 个周期 F fftshift(fft2(image2D)); imagesc(abs(F)); axis image; colorbar;这个信号会在频谱图上产生两个亮点。亮点的坐标应该分别对应 (8,3) 和 (-8,-3) 的位置对应频率为 8/(128·dx) 和 3/(128·dx)。如果亮点实际坐标距离这个值差一位说明你的频率轴公式里用了 N 还是 N-1 要重新核对。3.2 将 CSV 数据导入到 MATLAB 中进行二维 FFT 仿真CSV 导入看起来简单实际坑经常在“非矩形网格”上。二维 FFT 要求输入矩阵在空间上等间距排列而 CSV 文件经常是测量仪器输出的散点列表第一列 x 坐标、第二列 y 坐标、第三列数值。这种数据不能直接fft2必须先重排到规则网格。我推荐的处理流程是先用readmatrix读入检查坐标是否等间隔再决定是否插值。raw readmatrix(surface_scan.csv); % 自动识别数值列 coordsX raw(:, 1); coordsY raw(:, 2); values raw(:, 3); % 判断 X 坐标是否为等间隔扫描 dx unique(diff(coordsX)); if length(dx) 1 warning(X 坐标非等间距建议先做 griddata 插值); end % 重排成矩阵按 y 行、x 列 [M, N] [numel(unique(coordsY)), numel(unique(coordsX))]; dataGrid reshape(values, [N, M]).; % 先按列填充再转置 % 缺失值处理此处直接置零不等同于真实值慎重 dataGrid(isnan(dataGrid)) 0; F fftshift(fft2(dataGrid)); imagesc(abs(F).^2); ...reshape(values, [N, M]).这行的逻辑是CSV 扫描顺序通常是 x 变化最快、y 变化最慢reshape默认按列填充所以先把数据排成 N 行再转置成 M 行 N 列的矩阵。isnan置零会把缺失点当成零值样本来处理在频谱上表现为高频噪声抬升。如果你的应用对频谱精度要求高应该用griddata插值补点或者用掩膜在频域剔除掉这些位置的贡献。很多 FFT 仿真结果不准根源不是 FFT 本身而是数据进 FFT 之前的网格重排出了问题。3.3 实数数据与复数数据二维 FFT 结果的共轭对称性真实图像、温度场这类数据全是实数经 fft2 后的频谱具有共轭对称性F(u,v) 和 F(-u,-v) 互为共轭所以频谱图只需要显示一半即可。但光学测量里经常遇到复数输入例如干涉图乘以一个复相位因子此时频谱不再对称幅度和相位都包含独立信息。[X, Y] meshgrid((0:255)/256, (0:255)/256); phaseMap exp(1j * 2*pi * (0.05*X 0.02*Y)); % 倾斜波前 Fc fftshift(fft2(phaseMap)); figure; subplot(2,1,1); imagesc(abs(Fc)); axis image; subplot(2,1,2); imagesc(angle(Fc)); axis image;这段代码里phaseMap是单位幅度的纯相位信号它的频谱会有一个偏移的亮点相位谱则保留波前倾斜信息。与实数输入相比复数输入的存储和计算量都翻倍FFT IP 核配置时需要显式选择复数模式。下表总结两种数据的差异。比较项实数输入复数输入频谱对称性共轭对称可只算一半无对称性必须算全平面内存占用双通道可压缩双通道全部保留典型来源图像、温度场相干成像、通信基带信号处理注意事项可视化常取 log 幅度幅相都要看归一化方式不同4. 参数怎么定补零、窗函数、谱域滤波的边界条件二维 FFT 的参数不止“变换长度”一个。补零长度、窗类型、滤波掩膜半径、数据位宽都会直接改变频谱图形态这一章给出每个参数的默认值和调整原则。4.1 补零与块大小分辨率上限和 DFT 插值补零是在矩阵尾部追加零行零列让参与 FFT 的尺寸变大。它带来的第一个效果是频率轴更密相邻 bin 间隔从 1/(N·dx) 降到 1/(N_padded·dx)。但这种细化只是插值不带来新的物理信息。真正限制分辨率的因素是原始数据的空间跨度跨度 L N·dx频率分辨率上限始终约等于 1/L。dataPad zeros(256, 256); dataPad(1:128, 1:128) data; % 数据放左上角其余补零 Fp fftshift(fft2(dataPad));注意把原始数据放在左上角而不是矩阵中心。FFT 默认起始点是 (0,0)数据位置移动相当于乘了一个线性相位因子放到中心会让相位谱多出一列线性项给相位解读带来额外负担。块大小选择上Xilinx FFT 和 STM32 的库通常要求 2 的幂但 MATLAB 没有这个限制用素数尺寸也能算。如果你的代码最终要移植到硬件建议从一开始就用 2 的幂。4.2 窗函数与频谱泄漏二维 Hann 窗怎么选任何有限空间范围的数据都相当于乘了一个矩形窗矩形窗在频域对应很宽的高频旁瓣这就是频谱泄漏。对二维数据一维窗函数要扩展成二维。MATLAB 里最直接的方法是外积win1D hann(N, periodic); % periodic 适合谱分析 win2D win1D * win1D; % 外积生成二维窗 dataWindowed data .* win2D; % 逐点相乘参数periodic表示周期型窗它首尾连线连续适合 FFT 谱估计另一个选项symmetric更适合 FIR 滤波器设计。二维外积窗是行方向和列方向窗的叠加所以频率选择性是各向同性的。对于多峰值信号我推荐下表选型。窗函数主瓣宽度第一旁瓣衰减适用场景矩形窗最窄-13 dB瞬态信号、硬件资源受限Hann 窗较宽-31 dB一般频谱分析默认选择Hamming 窗较宽-41 dB近距离多峰分离Blackman 窗更宽-58 dB动态范围要求极高窗函数会以 2 倍主瓣宽度为代价换取旁瓣衰减所以不要盲目追求高衰减。如果信号要做的不是谱显示而是严格反变换重建矩形窗往往是更安全的选择因为任何加窗都会破坏原信号的能量分布。4.3 谱域滤波从理想低通掩膜到边界反射二维谱域滤波的标准写法是先 fft2、再 fftshift、然后乘以和频谱同尺寸的 0/1 掩膜、ifftshift 后逆变换。最容易错的一步是掩膜是否需要ifftshift。如果你的频谱 S 已经过fftshift那么中心坐标 (M/21, N/21) 就是零频构建掩膜时用 meshgrid 直接从负 Nyquist 排到正 Nyquist乘法可以直接做但逆变换前必须ifftshift回去[M, N] size(data); fx (-N/2 : N/2-1) / (N * dx); fy (-M/2 : M/2-1) / (M * dx); [FX, FY] meshgrid(fx, fy); cutoff 0.2; % 截止频率单位与 fx 相同 maskLP double(sqrt(FX.^2 FY.^2) cutoff); % 圆形理想低通 F fftshift(fft2(data)); filtered ifft2(ifftshift(F .* maskLP)); % 先乘再 ifftshift这里maskLP是中心化布局F也是中心化布局二者可以直接逐点相乘。如果漏掉最后的ifftshift逆变换出来的图像会发生环形平移表现为目标整体搬移到图像边缘甚至四角。理想低通掩膜的边界在频域是一个硬切断会在空间域产生吉布斯振铃。想减轻振铃就把掩膜的非通带区域从 1 平滑过渡到 0常见做法是加一个余弦滚降带。4.4 在嵌入式与 FPGA 上做二维 FFT 的参数约束把二维 FFT 搬进单片机或 FPGA 时内存布局往往比算法本身更影响性能。STM32F4 的 DSP 库提供一维 fft 函数二维变换要分两遍运算数据先按行读入算完存到外部 RAM再以转置后的顺序读列算第二遍。转置操作的开销常常大于 FFT 本身这就是为什么很多嵌入式 FFT 实战文章强调存储顺序。// 伪代码用一维 FFT 核做二维 FFT 的经典两遍法 // 第一遍对每一行调用 1D FFT结果按列优先写入转置缓冲区 // 第二遍按行读取转置缓冲区再调用 1D FFT // 注意转置缓冲区的行大小和列大小互换Xilinx Vivado 里的 FFT IP 核配置页主要看三个参数变换长度、输入数据位宽和缩放调度scaling schedule。二维应用里 IP 核一次只处理一行你和系统时钟的换算关系是总时钟周期约等于行数乘列数再乘每个 FFT 的延迟。如果 IP 核配置成定点模式每级蝶形运算都会增加位宽必须开缩放否则数据会溢出。有一个经验值输入 16 位、变换长度 1024第一级缩放到 16 位、其余级不缩放信噪比通常能接受但最终要以实际频谱图与 MATLAB 双精度结果对比为准。5. 现场验证技巧频谱图和功率谱密度图的快速核对二维 FFT 项目里我最后一步不会直接看效果而是做一套可复现的数值验证把 MATLAB、嵌入式 DSP 库和 FPGA 仿真结果拉齐。下面这套流程能够回答“我算出来的频谱到底对不对”的问题也避免你把生成频谱图像误认为 FFT 正确。第一个验证是 Parseval 定理。二维 FFT 是正交变换能量在空间域和频率域保持守恒但有归一化常数Espace sum(abs(data(:)).^2); Efreq sum(abs(fft2(data)).^2) / (M*N); relError abs(Espace - Efreq) / Espace;relError在双精度下应该在 1e-15 量级。如果这个值达到 1e-3说明数据里混入了 NaN、Inf或者矩阵尺寸定义不一致先修数据再谈算法。第二个验证是峰值 bin 定位。构造一个已知周期 T 的二维正弦读出峰值所在矩阵坐标反算频率并与理论值比较。下表是我常用的三个测试案例和对应诊断。测试信号预期结果常见偏差原因中心冲激功率谱全平面平坦fftshift 使用次数不对水平正弦 8 周期亮点在 fx8/N 处频率轴公式用 N-1 代替 N无噪声常数平面只有直流点有值均值未提前减掉低频台阶高斯光斑输出仍是高斯宽度成反比补零后网格坐标没有同步更新第三个验证针对频域滤波把滤波后的结果再做一次正变换观察频谱是否真的把掩膜外的分量清掉。如果掩膜外还有亮点通常是因为maskLP里的坐标和频谱坐标没有对齐检查meshgrid的参数顺序是否写反。很多人在这一步能发现问题不在 vivado fft 核而在频率轴坐标定义。最后给你一个硬经验fftshift和ifftshift的配对使用永远成对出现。你在频域做任何乘法、开窗或剪裁之前先检查明显频率量是否已经被 shfit 到位做逆变换之前先对中心化数据ifftshift。当 N 为奇数时fftshift(fftshift(x))会把 DC 移回原来位置但相位关系已经被打乱所以不要试图用两次相同的函数相互抵消。记住这句话二维傅里叶变换路径上最大的坑你已经跨过去一半。本文还有配套的精品资源点击获取