二维傅里叶变换从原理到实践：fft2、频谱搬移与频域滤波解析

简介：这是一份基于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) = Σ_{x=0}^{M-1} Σ_{y=0}^{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 决定，与窗函数、补零长度都无关。这个关系在第四章还会用到。下面表格列出不同尺寸下频率轴的行列对应关系，帮助你快速核对。

FFT 的 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 核配置时需要显式选择复数模式。下表总结两种数据的差异。

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 滤波器设计。二维外积窗是行方向和列方向窗的叠加，所以频率选择性是各向同性的。对于多峰值信号，我推荐下表选型。

窗函数会以 2 倍主瓣宽度为代价换取旁瓣衰减，所以不要盲目追求高衰减。如果信号要做的不是谱显示而是严格反变换重建，矩形窗往往是更安全的选择，因为任何加窗都会破坏原信号的能量分布。

4.3 谱域滤波：从理想低通掩膜到边界反射

二维谱域滤波的标准写法是：先 fft2、再 fftshift、然后乘以和频谱同尺寸的 0/1 掩膜、ifftshift 后逆变换。最容易错的一步是掩膜是否需要ifftshift。如果你的频谱 S 已经过fftshift，那么中心坐标 (M/2+1, N/2+1) 就是零频，构建掩膜时用 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 的二维正弦，读出峰值所在矩阵坐标，反算频率并与理论值比较。下表是我常用的三个测试案例和对应诊断。

第三个验证针对频域滤波：把滤波后的结果再做一次正变换，观察频谱是否真的把掩膜外的分量清掉。如果掩膜外还有亮点，通常是因为maskLP里的坐标和频谱坐标没有对齐，检查meshgrid的参数顺序是否写反。很多人在这一步能发现，问题不在 vivado fft 核，而在频率轴坐标定义。

最后给你一个硬经验：fftshift和ifftshift的配对使用永远成对出现。你在频域做任何乘法、开窗或剪裁之前，先检查明显频率量是否已经被 shfit 到位；做逆变换之前，先对中心化数据ifftshift。当 N 为奇数时，fftshift(fftshift(x))会把 DC 移回原来位置，但相位关系已经被打乱，所以不要试图用两次相同的函数相互抵消。记住这句话，二维傅里叶变换路径上最大的坑你已经跨过去一半。

本文还有配套的精品资源，点击获取