企业官网建设流程全解析

简介：以MATLAB为工具的傅里叶变换与4F系统仿真资源，面向图像处理入门者和光学信息处理方向学生，解决频域分析、滤波操作及光学系统模拟中的动手实践问题。压缩包共3个文件，包含MATLAB脚本、示例图像和FIG图窗文件，整体82.54MB；脚本覆盖从读取图像、fft2正变换、频域掩模处理到ifft2逆变换的完整流程，示例图片供测试效果，FIG文件可直接查看滤波结果。该资源已有751人学习/下载。运行后可掌握二维傅里叶变换的核心函数调用、滤波器设计与频谱可视化方法，并借助4F系统模拟理解光学图像处理的频域机制，参照代码修改掩模即可实现高通、低通等不同滤波实验，适合课程实验、期末设计或自学进阶。

1. 4F 系统仿真为什么绕不开傅里叶变换

4F 系统是光学信息处理里最经典的级联结构：两个焦距相等的透镜相距 2f，物面放在第一个透镜的前焦面，像面落在第二个透镜的后焦面，中间那个公共焦面就是频谱面。物面到像面的总距离正好是 4 个焦距，所以叫 4F。它之所以重要，是因为两个透镜在这里分别承担“傅里叶变换”和“逆傅里叶变换”的角色——你在频谱面上加一个挡片或相位板，就等价于对图像的二维频谱做了一次滤波。

真实光路上调 4F 系统很麻烦：透镜要对准、焦面位置要用剪切板反复确认、频谱面上的针孔孔径稍微偏一点，像面就完全变形。而在 MATLAB 里，这件事可以用 fft2 和 ifft2 两句核心调用复现。你只需要把物光场铺成离散网格，做一次正变换进频谱域，再在频谱面乘一个滤波 mask，最后逆变换回像面，就能看到实验里基本一致的滤波结果。这篇文章就沿着这条路走一遍：先解决网格与采样，再讲透镜怎么“变成”傅里叶变换器，然后给出完整的 4F 滤波仿真代码，最后用角谱法做交叉验证。

2. 用 MATLAB 搭 4F 仿真前的光场网格：采样间隔与 fftshift

FFT 本身不知道你算的是光学还是信号，它只处理一个复数矩阵。光学仿真里每个矩阵元素代表一个微小面元上的复振幅：模长是幅度，幅角是相位，单位是 V/m 或直接归一化都行。要让 fft2 算出来的结果能对应到 4F 系统里真实的频谱面，第一步是你得把坐标网格造对。

2.1 光场网格与物理参数

我一般先用一组固定参数搭骨架，后面所有代码都复用这套变量。波长为 632.8 nm（He-Ne 激光），焦距 f=0.5 m，网格取 1024×1024，采样间隔 8 μm。网格总宽 8.192 mm。

lambda = 632.8e-9; % 波长 (m) f = 0.5; % 透镜焦距 (m) N = 1024; % 一维采样点数 pitch = 8e-6; % 物面采样间隔 (m) L = N * pitch; % 物面总宽度 (m) x = (-N/2 : N/2-1) * pitch; [X, Y] = meshgrid(x); k = 2*pi / lambda; % 波数

坐标从-L/2排到L/2之前，原点在矩阵中心，这符合光学里“光轴沿 z、横截面以光轴为原点”的约定。meshgrid生成二维坐标网格，后面生成物函数、透镜相位都要在X,Y上求值。

注意meshgrid的输出形状：第一个输出X的每一行相同，第二个输出Y的每一列相同。如果后面把X当横坐标、Y当纵坐标就不会搞反。采样间隔 8 μm 对可见光来说是偏粗的，但仿真物是微米到毫米尺度的特征，这个尺度足够，而且减小 N 能明显加快 FFT。

2.2 fftshift 的必要性：直流分量不在中心

MATLAB 的fft2输出遵循 DFT 的约定，零频（直流）分量在矩阵的四个角上，而不是在中心。光学里频谱面关注的是零频在中心、高频向外扩散的分布，所以每一轮fft2之后都要做一次fftshift：

Uf = fft2(U0); % 直接变换，直流在四角 Uf_shifted = fftshift(Uf); % 移中心，直流在正中

对应地，逆变换之前要先ifftshift再ifft2。一个容易出错的地方是fftshift和ifftshift在 N 为奇数时行为不同：前者把原索引 0 移到中心，后者把中心元素移回索引 0，二者互为逆操作但实现细节不同。用 1024 这类偶数点数时两者等价，所以我建议网格尺寸固定取 2 的幂，能少踩一个坑。

2.3 空间频率坐标和物理频谱面坐标

fft2输出的每个点对应一个空间频率(fx, fy)，单位是 cycles/m（周/米）。频率坐标必须和网格匹配，否则滤波 mask 不知道自己的带宽在哪里。

fx = (-N/2 : N/2-1) / L; % 频率分辨率 1/L [FX, FY] = meshgrid(fx);

fftshift之后的Uf_shifted矩阵中，中心点对应(0,0)，向右一格频率增加1/L。最高频率是1/(2*pitch)，即 62500 cycles/m，这就是该网格能表示的频谱范围边界。把这些参数的关系整理成一个表，调参时对着查就行。

频谱面物理坐标（单位 m）由空间频率换算得到：

xi = FX * lambda * f; eta = FY * lambda * f;

这里FX是空间频率，乘以lambda*f得到的是 4F 系统里真实频谱面上的位置。后面做滤波 mask 时，mask 半径若按物理尺寸写，就必须用这套坐标。

2.4 物函数从哪来：合成图案优先

4F 仿真最适合的输入是二元几何图案和灰度图像，两者我都常用。合成图案的好处是不依赖外部文件，任何机器上都能跑：

% 两个圆孔，中心分别在 x = ±1.5 mm r = 0.5e-3; U0 = double(sqrt((X - 1.5e-3).^2 + Y.^2) <= r) ... + double(sqrt((X + 1.5e-3).^2 + Y.^2) <= r);

如果物函数来自实验测量数据，比如 CCD 采集的振幅分布存成了 CSV，用readmatrix读进来后确保它是一个二维矩阵，再插值或裁剪到N×N，转成 double 就能进入同一条流水线。灰度图也同理，imread后取单通道并转 double 即可。这类输入的重点是确保矩阵尺寸与N一致，否则频谱坐标会错位。

3. 4F 系统的两次傅里叶变换：透镜相位、频谱面与像面

FFT 只是数值工具，真正让 4F 系统成立的是透镜的相位变换作用。理解了透镜相位，你才知道为什么频谱面能放滤波器，也知道纯fft2仿真里丢掉了哪个全局相位。

3.1 透镜是把光场变换到频域的器件

薄透镜近似的复振幅透过率是：

phiLens = exp(-1j * k / (2*f) * (X.^2 + Y.^2));

它的作用是对入射光场加一个与半径平方成正比的相位延迟。入射场U_in经过透镜后变成：

U_out = U_in .* phiLens;

这个二次相位因子不是随便写的。入射光场带上球面相位后继续向前传播，在透镜后焦面上，菲涅尔衍射积分里的二次相位项恰好被透镜相位抵消，剩下的积分式正好是入射光场的二维傅里叶变换：

U_f(fx, fy) ∝ ∫∫ U_in(x, y) * exp(-j*2*pi*(fx*x + fy*y)) dx dy

其中fx = xi/(lambda*f)，fy = eta/(lambda*f)。这就是“透镜做傅里叶变换”的来历。焦距相同的两个透镜级联，第一个把物场变到频域，第二个把频域场变回空域，于是构成 4F 系统。

3.2 频谱面上的球面相位与 4F 的相位补偿

严格推导时，透镜后焦面上的场不是纯傅里叶变换，而是要乘一个球面相位exp(j*k/(2f)*(xi^2 + eta^2))。这个相位的存在意味着频谱面上的复振幅带有弯曲波前，直接再放一个透镜做逆变换时，前面累积的相位和第二级透镜的相位会再次相消，最后像面得到的是近似几何像（倒像）。4F 的精妙之处就在于两段对称的结构自动完成了二次相位的补偿。

离散仿真里这个相位怎么处理？我见过不少写法是直接fft2拿频谱、ifft2还原，完全不加透镜相位，这在“只看强度”的场景下完全够用，因为频谱面上的球面相位只影响相位分布，不影响abs(Uf).^2的功率谱。但如果你要在频谱面放相位型滤波器（涡旋相位板、相位二分板等），或者关心复振幅的干涉结果，就应该把这一项补上：

Uf_phys = fftshift(fft2(U0)) .* exp(1j * k/(2*f) * (xi.^2 + eta.^2));

xi,eta就是 2.3 节里换算出的物理坐标。这里乘上的相位模拟的是“真实频谱面上的波前弯曲”，最终ifft2还原时会自然消掉，因此像面强度结果基本不变，但中间过程的相位结构更接近实验。写滤波器代码时我通常把这两条路径都保留：显示时用Uf_phys，滤波时直接用fftshift(fft2(U0))，避免多一次复数乘法。

3.3 频谱面的可视化

频谱动态范围往往跨越好几个数量级，线性显示只能看到中心亮点，所以可视化用对数强度：

S = log10(1 + abs(Uf_phys).^2); imagesc(xi, eta, S); axis image; colormap(gray);

xi和eta作为imagesc的坐标轴传入，x 轴、y 轴的单位就是米。这一步建议每次滤波前后都看一眼：它能确认 mask 中心是否与零频对齐，也能直观看到截断半径对应的频谱范围。

4. MATLAB 里的 4F 滤波链：低通去噪、高通边缘与像面重建

这一章给出完整可运行的 4F 滤波仿真。物函数用双圆孔，目的是观察衍射条纹和滤波后的细节变化。实际替换成灰度图同样成立。

4.1 完整仿真代码

close all; clear; % ---------- 物理参数 ---------- lambda = 632.8e-9; f = 0.5; N = 1024; pitch = 8e-6; L = N * pitch; x = (-N/2 : N/2-1) * pitch; [X, Y] = meshgrid(x); k = 2*pi / lambda; % ---------- 物函数：两个圆孔 ---------- r = 0.5e-3; U0 = double(sqrt((X - 1.5e-3).^2 + Y.^2) <= r) ... + double(sqrt((X + 1.5e-3).^2 + Y.^2) <= r); % ---------- 频谱面 ---------- Uf = fftshift(fft2(U0)); fx = (-N/2 : N/2-1) / L; [FX, FY] = meshgrid(fx); % ---------- 低通滤波 ---------- rho_c = 3162; % 截止频率，单位 cycles/m Hlp = double(sqrt(FX.^2 + FY.^2) <= rho_c); U_lp = ifft2(ifftshift(Uf .* Hlp)); % ---------- 高通滤波 ---------- Hhp = 1 - Hlp; U_hp = ifft2(ifftshift(Uf .* Hhp)); % ---------- 显示 ---------- figure('Color','w'); subplot(1,3,1); imagesc(x*1e3, x*1e3, abs(U0)); axis image; colormap(gray); title('物面'); subplot(1,3,2); imagesc(fx, fx, log10(1 + abs(Uf).^2)); axis image; title('频谱（对数强度）'); subplot(1,3,3); imagesc(x*1e3, x*1e3, abs(U_lp).^2); axis image; title('低通像面');

代码的核心逻辑是：fftshift(fft2(U0))把频谱中心移到坐标原点，Hlp按空间频率半径做二值 mask，乘到频谱上相当于把高于截止频率的分量清零，ifftshift把频谱还原成 FFT 的原始排列，最后ifft2回到空域。注意abs(U_lp).^2是强度，实验里探测器只能看到这个量。

4.2 截止频率怎么换算

rho_c的单位是 cycles/m，它和真实频谱面上针孔半径的换算关系是：

rho_physical = lambda * f * rho_c

比如rho_c = 3162cycles/m 时，对应物理半径632.8e-9 * 0.5 * 3162 ≈ 1e-3m，也就是 1 mm 的针孔。这个换算在做实验预演时特别重要：你在代码里设的是 1 mm，光路里就应该准备直径 2 mm 的针孔。反过来，你知道实验上只想放一个 0.5 mm 的小孔，代码里的截止频率就是：

rho_c = 0.5e-3 / (lambda * f);

算出来约 1581 cycles/m。常见截止半径对应的效果参考下表。

注意频谱面上能表示的最大物理半宽是0.5 * N * lambda*f / L，本例约为 19.3 mm，所以 2 mm 仍远小于频谱面尺寸，属于窄带滤波。rho_c最大不能超过最大频率1/(2*pitch)，即 62500，否则 mask 失去滤波意义。

4.3 高通滤波与边缘提取

高通滤波把零频附近的低频分量去掉，只保留高频边缘信息，效果等价于图像处理里的边缘提取。代码如下：

U_hp = ifft2(ifftshift(Uf .* (1 - Hlp))); imagesc(x*1e3, x*1e3, abs(U_hp).^2); axis image;

高通输出里圆孔内部变暗，而圆孔的边界会变成一条亮环。这是因为边界对应频谱中的高频分量，被原样保留甚至增强。实验上高通滤波常用一个小黑屏挡住频谱中心，而不是用大孔径光阑，正是对应这里的1 - Hlpmask。

一个值得说的是：高通滤波去除直流后，像面强度的平均值为零，所以abs(U_hp).^2里会出现大量暗区，这是正常现象，不是数值误差。如果你想让边界环更亮，可以把 mask 改成带增益的形式Hhp = 1 + a * Hlp，通过调节a控制边缘增强强度。

4.4 把物函数换成灰度图

换输入只需替换U0那一句：

cam = imread('cameraman.tif'); % 灰度图，256x256 cam = imresize(cam, [N, N]); % 重采样到网格 U0 = double(cam) / 255; % 归一化到 [0,1]

imresize要凑到N×N，否则频率坐标全部错位。摄像师图本身的频谱能量集中在中低频，rho_c=3162时像面会明显变平滑，细节如支架、人脸轮廓被抹掉；rho_c=6324时基本恢复原图。这一步能直观看到低通在图像上的效果。

由于仿真里没有加入噪声，低通的去噪效果看不出来。想模拟实验场景，可以在物面加高斯噪声：

U0 = double(imread(...)) / 255 + 0.05 * randn(N, N);

噪声的能量分布在整个频谱面上，低通 mask 切掉高频后，像面的颗粒感明显减少，这就更进一步接近真实实验了。

5. 用角谱传播交叉验证 4F 仿真，三个必查项

纯fft2的两步变换算得再顺，也只能说明数值自洽，不能证明它等价于真实 4F 光路。最稳的验证方式是用角谱传播把 4F 拆成五段：传播到透镜、乘透镜相位、再传播到频谱面、乘第二个透镜相位、最后传播到像面。两种路径的结果如果误差在可忽略范围，就说明简化模型没有用错。

5.1 角谱交叉验证代码

% 空间频率网格，沿用 4.1 节 fx = (-N/2 : N/2-1) / L; [FX, FY] = meshgrid(fx); % 自由空间传播 z 米的角谱算子（菲涅尔近似） prop = @(U, z) ifft2(ifftshift(fftshift(fft2(U)) .* ... exp(-1j * pi * lambda * z * (FX.^2 + FY.^2)))); phiLens = exp(-1j * k / (2*f) * (X.^2 + Y.^2)); Uv = prop(U0, f); % 物面到透镜1 Uv = Uv .* phiLens; % 过透镜1 Uv = prop(Uv, f); % 透镜1到频谱面 Uv = Uv .* phiLens; % 过透镜2 Uv = prop(Uv, f); % 透镜2到像面 err = norm(abs(Uv).^2 - abs(U_lp).^2) / norm(abs(U_lp).^2);

prop里先fftshift把频谱中心化，乘角谱转移函数，再ifftshift还原后逆变换。这里值要强调的是，角谱路径里每段都是真实的物理传播，而简化路径只有两次 FFT，所以误差主要来自简化路径省略的球面相位项。对强度而言，err通常在 1e-6 量级，基本可以忽略。

5.2 三个必查项

第一个必查项是像面方向。真实 4F 系统成倒像，而代码里fft2之后再ifft2得到的是正立像，因为两次变换的方向约定抵消了。所以判断仿真是否正确不能看方向，要看结构是否与原物一致，方向问题属于坐标约定，留到和实验对照时再处理。

第二个必查项是全通滤波恢复。把 mask 设成全 1，像面强度应该与原物强度完全一致，误差同样是 1e-6 量级。如果你发现像面有明显条纹或光晕，先检查fftshift/ifftshift是否成对使用，这是最常见的错误源。

第三个必查项是频谱 center 对齐。滤波 mask 的圆心必须对准频谱中心(N/2+1, N/2+1)，有偏移时像面会出现明显的整体倾斜调制。精确做法是把 mask 构造在频域网格上，而不是用imcrop之类的图像工具手动画圆。

最后给一个实用技巧：如果你要在频谱面模拟带方向性的滤波器，比如让图案只在竖直方向模糊，可以使用椭圆 mask：

H_ell = double((FX.^2 / a^2 + FY.^2 / b^2) <= 1);

其中a,b分别是两个主轴方向的截止频率。这样一套代码就能覆盖 4F 系统的大部分频域滤波实验，从各向同性低通到方向滤波，参数只需要改一个 mask 表达式。

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