4F相关器系统Matlab仿真：从透镜傅里叶变换到相关峰定位

简介：基于透镜傅立叶变换特性的4F相关器系统的Matlab仿真项目，聚焦光学信息处理与傅立叶光学中的频谱分析，适用于光电、通信、自动化等专业学生完成毕业设计、课程设计或大作业，也可作为相关科研人员的入门参考。压缩包共12个文件，包含1个Matlab源码文件（.m）、10张仿真结果图（.jpg）和1个README说明文档，整体大小约953KB。仿真结果覆盖无滤波、针孔、菲涅耳波带板、水平/垂直单缝与双缝、网格等多种输入场景，每张图片对应不同滤波器或物体的输出光场分布，便于对照验证。全部代码均经过测试并成功运行，答辩评审平均分达96分，可信度高。用户可通过源码与结果图像快速理解4F相关器的仿真实现，并在此基础上修改参数、更换滤波器或扩展功能，适合从基础学习到项目演进的多种需求，目前已有124人浏览学习。

1. 4F相关器系统的Matlab仿真，落点不在FFT而在物理坐标

4F相关器是光学信息处理和图像模式识别里的经典结构：两个焦距相同的透镜，中间放一个频谱面，把输入图像和参考图像的互相关运算在光路中一次算完。毕设或课设要求用Matlab仿真这个系统时，多数人第一反应是把fft2连起来跑一遍，结果频谱图出来了，相关峰也亮了，仔细一问却解释不清峰的位置为什么出现在那里、频谱面的坐标轴对应多少周期每毫米。原因是把仿真当纯数字运算，丢了透镜傅里叶变换特性里的关键映射：空间频率与后焦面物理坐标之间的比例关系。

这篇文章从透镜傅里叶变换原理出发，给出一套可直接运行的4F相关器系统Matlab仿真方案。内容包括坐标系统设计、匹配滤波器构造、相关峰定位，以及仿真参数调整和排错方法。读完之后你不仅能跑通代码，还能在答辩时把“透镜做了什么、仿真步对应哪个光学面”讲清楚。本文适合正在做光电方向毕设的学生，也适合需要用光学相关思路快速验证算法的工程师。

2. 透镜傅里叶变换和4F系统：先把光学结构翻译成仿真语言

2.1 薄透镜的相位变换与后焦面场分布

薄透镜对入射光场的作用，可以写成一个相位变换函数。设透镜焦距为 f，波长为 λ，忽略透镜有限孔径和厚度，则出射复振幅与入射复振幅满足：

U_out(x,y) = U_in(x,y) · exp(-jπ(x²+y²)/(λf))

这个二次相位项把平面波前变成球面波前，使得平行光会聚到焦点。把输入面放在透镜前焦面，在后焦面上观察光场，借助菲涅尔衍射公式可以得到：

U_f(u,v) = (e^(jkf))/(jλf) · exp(jπ(u²+v²)/(λf)) · ∬ U_in(x,y) exp(-j2π(xu+yv)/(λf)) dx dy

积分项正是二维傅里叶变换的标准形式。后焦面上的复振幅分布，除了一个不影响强度的球面相位因子，就是入射光场的傅里叶变换。这里最重要的不是公式本身，而是坐标代换关系：

fx = u/(λf)， fv = v/(λf)

后焦面上的物理坐标 u 除以 λf 才是空间频率。λf 这个乘积决定频谱面尺度：波长越长、焦距越长，同样的空间频率在频谱面上展得越开。仿真代码里的坐标轴必须按这个式子标注，否则画出来的频谱图的 u 轴是像素编号，不是物理尺寸。

2.2 4F系统的基本结构与仿真映射

4F系统由两个焦距相同的透镜 L1、L2 组成，四个关键平面之间的间距都是 f。输入面 P1 在 L1 前焦面，频谱面 P2 在 L1 后焦面同时也是 L2 前焦面，输出面 P3 在 L2 后焦面。光从 P1 到 P3 一共走四个焦距的距离，所以叫 4F。

这个表格可以直接映射为仿真代码的运算顺序。注意 P3 面上的数学操作在光学中是第二次傅里叶变换，等价于逆变换加坐标翻转。仿真中常用 ifft2 替代，省去输出面坐标翻转的麻烦，后面代码部分会具体说明。

2.3 频域共轭相乘与空域互相关的等价性

现在把参考图像 g 的频谱取共轭，构成匹配滤波器 H(fx,fy) = conj(G(fx,fy))。将 H 放在频谱面上，入射光场 f 经过 L1 后的频谱 F(fx,fy) 与 H 相乘，再经过 L2 的第二次变换，输出面得到的就是 f 与 g 的互相关。这一步的数学依据就是互相关定理：空域互相关的傅里叶变换等于一方频谱的共轭与另一方频谱的乘积。

在Matlab里用一行代码验证这个等价性：

corr_check = ifft2(fft2(input_img) .* conj(fft2(ref_img)), 'symmetric');

symmetric参数利用了实数输入经共轭相乘后结果仍近似实数的特点，把计算中残留的微小虚部直接清零，避免输出被复数误差干扰。这段代码没有引入任何光学参数，只验证频域乘法与空域相关性的一致关系。真正的4F仿真要在这一步的基础上补上坐标映射、滤波器幅度控制和透镜孔径约束。

3. 用Matlab实现4F相关器系统仿真：坐标、FFT、匹配滤波三步走

3.1 仿真坐标参数的设计

Matlab仿真4F相关器，第一步不是写傅里叶变换，而是把坐标系统定下来。需要设定四个量：采样边长 L、采样点数 N、波长 λ 和焦距 f。常见的一组取值如下：

N = 256; % 采样点数，取2的幂便于FFT L = 5e-3; % 输入面物理边长，单位m，对应5mm lambda = 632.8e-9; % 波长，氦氖激光经典值632.8nm f = 0.2; % 透镜焦距，单位m，对应200mm dx = L / N; % 空间采样间隔 fx = (-N/2 : N/2-1) / L; % 空间频率轴，单位cycles/m u = fx * lambda * f; % 频谱面物理坐标，单位m

这里频率轴 fx 的分母是 L 而不是 dx，因为离散傅里叶变换的频率分辨率由总观察范围决定，Δfx = 1/L。最大可表示频率由采样间隔决定，fmax = 1/(2dx)。这两个公式是所有傅里叶光学仿真的基础，很多人把频率轴写错就是在这里混淆了分辨率与最大频率。

频谱面物理坐标 u 的计算对应 2.1 节的公式。如果改成其他波长或焦距，只需要替换 lambda 和 f，u 轴的刻度会自动缩放。代码里所有显示坐标都用物理单位，避免出现“像素频率”这种没法跟实验对照的量。

3.2 第一步：输入图像的透镜傅里叶变换

生成一个简单的二值图像作为输入，计算它经过 L1 后在频谱面的分布：

input_img = zeros(N); input_img(80:176, 64:192) = 1; % 矩形目标，模拟透光物体 F1 = fftshift(fft2(input_img)); % 傅里叶变换并零频居中

fft2计算二维离散傅里叶变换，结果零频在 (1,1) 位置。fftshift把零频移动到数组中心，这一步对应光学频谱面上零频位于中心的物理事实。显示频谱幅度时，直接画 abs(F1) 会被中心极大值压缩掉全部细节，通常取对数后再显示：

figure; imagesc(u*1e3, u*1e3, log10(1 + abs(F1))); axis image; colormap bone; xlabel('u (mm)'); ylabel('v (mm)'); title('频谱面强度（对数刻度）');

log10(1+abs(F1))压缩动态范围让高频旁瓣可见，u*1e3把米转成毫米。矩形目标的频谱是 sinc 函数形状，中心主瓣加十字旁瓣。如果旁瓣看不到，检查 N 是否太小或矩形边缘恰好落在网格边界上。

3.3 第二步：构造匹配滤波器并做频谱面调制

4F相关器的核心操作在频谱面。对参考图像做同样的傅里叶变换，取共轭就得到匹配滤波器：

ref_img = zeros(N); ref_img(96:160, 96:160) = 1; % 参考目标，正方形 F2 = fftshift(fft2(ref_img)); H = conj(F2) / max(abs(F2(:))); % 共轭加幅度归一化

匹配滤波器的数学形式是参考频谱的共轭。max(abs(F2(:)))是零频分量的幅度，也是数组中的最大值，用它归一化可以把 H 的幅度限制在 1 以内，避免零频分量在频谱面相乘后产生巨大的直流背景淹没相关峰。如果目标是做严格的匹配滤波理论验证，可以不做归一化，但要意识到输出面会出现一个强度远高于其他区域的中心亮斑。

纯相位型滤波器在光学实现中更常见，做法是只保留共轭频谱的相位：

H_phase = exp(-1j * angle(F2));

纯相位滤波器对光照变化更鲁棒，但相关峰宽度和旁瓣水平与匹配滤波不同。建议先跑通幅度版本，再换相位版本对比。

频谱面乘法就是一次逐点相乘：

F_matched = F1 .* H;

这一步对应光学中透过 P2 面滤波器的光场。如果要在频谱面同时模拟透镜孔径截断，可以在这一步乘以一个圆形或矩形孔径函数：

[U, V] = meshgrid(u, u); pupil = double(sqrt(U.^2 + V.^2) <= D/2); % D为透镜口径 F_matched = (F1 .* H) .* pupil;

孔径截断会滤掉高于截止频率的频谱分量，实际效果是输出面的相关峰变宽、旁瓣增多。仿真时不加这个约束也能跑，但和真实光路的差距会变大。

3.4 第三步：L2透镜变换得到相关输出

频谱面调制完成后，经过 L2 在输出面得到相关结果：

output = ifft2(ifftshift(F_matched)); % 第二次变换，得到相关面复振幅 I_out = abs(output).^2; % 光强分布

这里用ifft2而不是fft2做第二次变换。原因在 2.2 节提过：光学透镜做的是正傅里叶变换，输出面相关函数带有坐标翻转。ifft2等价于逆变换，避免输出面的峰位置上下颠倒，读坐标时不用再手动翻转数组。如果要严格模拟光学系统的坐标翻转，可以用fft2(ifftshift(F_matched))再rot90(output, 2)，两种做法幅度分布一致，只是坐标方向不同。

abs(output).^2是把复振幅转成光强。光学探测器记录的是光强而不是复振幅，仿真输出也应该取强度。如果只看实部或幅度，相关峰的背景形状会不同。

3.5 输出面解读：相关峰的位置就是目标位置

把输入图像和参考图像设计成有明确空间关系的情况，可以验证输出是否正确：

target = zeros(N); target(96:160, 96:160) = 1; % 参考目标 input_img = circshift(target, [30, -20]); % 把目标移到新位置 F_input = fftshift(fft2(input_img)); H_filter = conj(fftshift(fft2(target))); F_matched = F_input .* (H_filter / max(abs(H_filter(:)))); I_out = abs(ifft2(ifftshift(F_matched))).^2; [maxv, idx] = max(I_out(:)); [r, c] = ind2sub(size(I_out), idx); peak = [r, c] - (N/2 + 1); fprintf('检测到目标偏移: 行 %d, 列 %d\n', peak(1), peak(2));

circshift把目标沿行方向移 30 像素、列方向移 -20 像素。相关峰在输出面上的位置应该对应这个偏移量，即峰出现在中心点偏移 (30, -20) 处。ind2sub把线性索引转回行列坐标，减去中心索引得到相对零频的偏移。如果输出峰位置与期望值符号相反，检查第二次变换用的是 fft2 还是 ifft2；如果峰值固定在中心不动，检查 fftshift/ifftshift 是否配对出错。

4. 4F仿真参数怎么调、相关峰为什么不锐：常见坑与排错

4.1 参数选择对照表

仿真参数直接影响相关峰的质量，调整时的参考依据如下：

一组适合毕设起步的推荐值是 N=256、L=5mm、λ=632.8nm、f=200mm。先跑通再逐项调整，每次只改一个参数观察输出面变化。

4.2 相关峰不锐利的三个原因

第一个原因是输入图像没有填满整个仿真面。输入图在观察域边缘有硬边界，FFT默认图像周期延拓，硬边界在频谱中产生高频拖尾。表现为相关峰周围出现十字亮线和周期性重复峰。处理方法是在输入图像外围补零或者用窗函数平滑边缘：

window = hann(N) * hann(N)'; % 二维汉宁窗 input_win = input_img .* window;

加窗会略微展宽相关峰，但能大幅压低旁瓣。实际光学实验里，有限尺寸的物体天然带边缘效应，所以加窗是仿真向实验靠拢的合理手段。

第二个原因是滤波器直流分量没有归一化。参考图像如果是全白的大方块，它的零频分量数值会是图像面积的量级，H 中其他频率分量相对很小。相乘后频谱面的直流分量过强，逆变换输出面会出现一个高耸的中心峰，真正的互相关峰被淹没在背景里。解决方法是 3.3 节里的幅度归一化，或者在频谱面把零频分量单独置零再参与运算。

第三个原因是频谱面采样不够。当目标在图像中只有几个像素大时，它的频谱主瓣很宽，而仿真频率轴 Δfx=1/L 由整幅观察域决定。如果 L 太大，Δfx 过小，频谱主瓣的采样点数会稀疏，相关峰形状不再平滑。检查方法是对输出面的相关峰做插值显示，看峰是否存在“振铃”。如果出现，减小 L 或增大 N。

4.3 用已知平移量验证仿真链路的正确性

最可靠的验证方法是构造一个完全可控的输入：参考图是确定图案，输入图是参考图的精确平移。跑完整个4F仿真后，检查输出峰的位置是否等于预设平移量。

N = 256; ref = zeros(N); ref(96:160, 96:160) = 1; delta = [25, -18]; input = circshift(ref, delta); F_in = fftshift(fft2(input)); F_re = fftshift(fft2(ref)); H = conj(F_re) / max(abs(F_re(:))); out = abs(ifft2(ifftshift(F_in .* H))).^2; [~, idx] = max(out(:)); [r, c] = ind2sub(size(out), idx); measured = [r, c] - (N/2 + 1); if isequal(measured, delta) disp('仿真链路正确'); else fprintf('期望(%d, %d) 实测(%d, %d)\n', delta(1), delta(2), measured(1), measured(2)); end

如果实测偏移与预设不一致，按下述顺序排查：先去掉所有 fftshift 相关调用，用最朴素的 ifft2(fft2(f).*conj(fft2(g))) 验证数学链路；确认无误后再加回坐标映射和归一化。顺带检查 N 为奇数时 fftshift 与 ifftshift 行为不同，N 为偶数时二者等价。建议全程保持 N 为 2 的幂，省掉这一类边界问题。

5. 扩展方向：从4F相关器到联合变换相关器与纯相位滤波

5.1 联合变换相关器的仿真实现

联合变换相关器（JTC，Joint Transform Correlator）与4F的区别在于不做匹配滤波器，而是把参考图像和输入图像并排放在输入面，先记录联合功率谱，再做一次傅里叶变换得到相关峰。仿真代码比4F更短：

N = 512; input_img = zeros(N); ref_img = zeros(N); input_img(50:150, 50:150) = 1; % 场景图 ref_img(300:400, 300:400) = 1; % 参考图放另一侧 joint = input_img + ref_img; % 联合输入面 J = fft2(joint); JPS = abs(J).^2; % 联合功率谱 out_JTC = abs(ifft2(JPS)).^2; % 逆变换得到相关输出

输出面上中心是零频亮斑，两侧对称分布两个相关峰。峰值位置分别对应两幅图的相对位移。JTC 的优势是不需要预先计算滤波器复数模板，实时更新参考图只需重新放置输入面即可，因此实战中用得比4F多。代价是峰背景比较高，且需要额外处理中心亮斑。仿真时可以把 JPS 减去其均值再逆变换，中心亮斑能明显减弱。

5.2 峰值背景比的量化验证

答辩或报告里如果只有相关峰图，说服力不够。建议加一个量化指标：峰值背景比（Peak-to-Background Ratio，PBR），定义为相关峰强度与输出面背景均值的比值。

peak_val = max(I_out(:)); mask_bg = true(size(I_out)); mask_bg(peak_r-2:peak_r+2, peak_c-2:peak_c+2) = false; % 挖掉峰周围 bg_mean = mean(I_out(mask_bg)); PBR = peak_val / bg_mean;

这个数值可以写进结论：PBR 越高代表识别越可靠。实际应用中 PBR 低于 10 时目标基本不可辨识，高于 100 则是很清晰的相关峰。对不同噪声水平的输入重复计算 PBR，画一条 PBR 随噪声变化的曲线，整份毕设的仿真数据就完整了。

5.3 在频谱面做纯相位滤波的实验对比

把4F仿真的滤波器替换为纯相位版，观察相关峰形态的变化：

H_phase = exp(-1j * angle(F2)); % 纯相位滤波器 out_phase = abs(ifft2(ifftshift(F1 .* H_phase))).^2;

纯相位滤波丢失了幅度信息，因此对光照不均匀和对比度变化更稳健；但旁瓣水平会比匹配滤波高一点。仿真中对比两组输出面的 PBR 和峰宽，是很好的分析实验。如果时间充裕，还可以扩展梅林变换做尺度不变识别，或者换成三波长照明做彩色图像的4F相关。每一个扩展都是毕设里能独立成一节的素材，且都建立在本文这套坐标系统和仿真链路上。

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