簡(jiǎn)介本資源是一套面向光學(xué)工程、信息光學(xué)及計(jì)算成像方向初學(xué)者與高校實(shí)驗(yàn)教學(xué)的數(shù)字全息仿真實(shí)驗(yàn)MATLAB實(shí)現(xiàn)方案聚焦數(shù)字全息圖生成、零級(jí)像抑制、波前再現(xiàn)等核心原理的編程驗(yàn)證。壓縮包共2個(gè)文件1個(gè)BMP格式原始全息圖數(shù)據(jù)、1個(gè)holographic.m主程序腳本總大小3.42MB結(jié)構(gòu)精煉便于快速運(yùn)行與代碼剖析其中MATLAB腳本完整涵蓋圖像讀取、傅里葉變換、空域高斯濾波去零級(jí)、逆變換及衍射再現(xiàn)全流程可直接用于課堂演示或課后復(fù)現(xiàn)。已有2423人學(xué)習(xí)下載適用于《信息光學(xué)》《計(jì)算全息》課程實(shí)驗(yàn)環(huán)節(jié)幫助學(xué)習(xí)者打通“光學(xué)原理—數(shù)值建模—MATLAB實(shí)現(xiàn)—圖像分析”的完整鏈路切實(shí)掌握從干涉記錄到三維物場(chǎng)重建的關(guān)鍵技術(shù)細(xì)節(jié)與調(diào)試邏輯。1. 項(xiàng)目概述從理論到屏幕的數(shù)字全息之旅數(shù)字全息仿真實(shí)驗(yàn)聽起來像是光學(xué)實(shí)驗(yàn)室里高深莫測(cè)的玩意兒但實(shí)際上它正是一扇連接經(jīng)典光學(xué)理論與現(xiàn)代計(jì)算成像的絕佳窗口。簡(jiǎn)單來說這項(xiàng)目就是用MATLAB這把“數(shù)字瑞士軍刀”在計(jì)算機(jī)里完整地模擬一套全息記錄與再現(xiàn)系統(tǒng)。你不用真的去搭建昂貴且嬌貴的光學(xué)平臺(tái)不用擔(dān)心激光器的穩(wěn)定性更不用在暗房里小心翼翼地處理全息干板。所有過程從生成模擬的物光波到模擬參考光干涉形成全息圖再到最終的數(shù)字再現(xiàn)全部在代碼和矩陣運(yùn)算中完成。這解決了什么問題對(duì)于學(xué)生和研究者它降低了學(xué)習(xí)全息原理的門檻讓你可以直觀地、可重復(fù)地觀察每一個(gè)參數(shù)變化對(duì)最終成像的影響。對(duì)于工程師它成為了一個(gè)強(qiáng)大的設(shè)計(jì)驗(yàn)證工具可以在實(shí)際搭建光路前預(yù)先仿真不同系統(tǒng)配置比如改變光波長(zhǎng)、記錄距離、探測(cè)器像素尺寸的成像效果節(jié)省大量時(shí)間和成本。無論你是光學(xué)工程專業(yè)的學(xué)生想深入理解《信息光學(xué)》課本里的公式還是從事計(jì)算成像、顯微成像或三維顯示研發(fā)的工程師需要快速驗(yàn)證一個(gè)新想法這個(gè)基于MATLAB的數(shù)字全息仿真實(shí)驗(yàn)都能提供一個(gè)清晰、可控且功能強(qiáng)大的沙盤。2. 仿真實(shí)驗(yàn)的核心思路與框架設(shè)計(jì)數(shù)字全息仿真的核心思路是對(duì)物理全息過程的嚴(yán)格數(shù)學(xué)建模和離散化計(jì)算。整個(gè)過程可以清晰地拆解為三個(gè)核心階段對(duì)應(yīng)著三個(gè)主要的MATLAB函數(shù)模塊。2.1 第一階段模擬物光波前生成全息記錄的是物光波的振幅和相位信息。在仿真中我們首先要“創(chuàng)造”一個(gè)虛擬的物體及其發(fā)出的光波。這里的關(guān)鍵在于如何用數(shù)學(xué)描述一個(gè)復(fù)雜的光場(chǎng)。最直接的方法是采用角譜傳播理論。我們假設(shè)物體是一個(gè)二維的透射或反射率分布圖例如一個(gè)簡(jiǎn)單的字母“A”的圖片。這個(gè)分布圖可以看作是一個(gè)平面上的復(fù)振幅分布其中振幅代表物體的透射或反射強(qiáng)度初始相位可以設(shè)為零或一個(gè)隨機(jī)相位板用于模擬粗糙表面。然后我們需要計(jì)算這個(gè)初始平面光場(chǎng)傳播一定距離即記錄距離d0后到達(dá)全息記錄平面如CCD靶面的復(fù)振幅分布。這個(gè)過程通過角譜傳播函數(shù)實(shí)現(xiàn)其本質(zhì)是求解標(biāo)量衍射的積分方程在頻域里它表現(xiàn)為一個(gè)傳遞函數(shù)的乘積運(yùn)算。注意為什么不直接用菲涅爾衍射或卷積法角譜理論在數(shù)學(xué)上是最嚴(yán)格的標(biāo)量衍射近似只要采樣滿足奈奎斯特頻率它對(duì)任何距離的傳播計(jì)算都是準(zhǔn)確的避免了菲涅爾近似在極近場(chǎng)時(shí)的誤差。這對(duì)于構(gòu)建一個(gè)基礎(chǔ)扎實(shí)的仿真框架至關(guān)重要。在MATLAB中這意味著我們需要對(duì)物體的二維矩陣進(jìn)行二維快速傅里葉變換2D-FFT乘以一個(gè)對(duì)應(yīng)于傳播距離的相位傳遞函數(shù)exp函數(shù)構(gòu)成再進(jìn)行逆傅里葉變換。這個(gè)傳遞函數(shù)H是仿真精度的心臟其表達(dá)式為H exp(1i*2*pi*d0/lambda * sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2))其中fx,fy是空間頻率坐標(biāo)lambda是光波長(zhǎng)。這里就涉及到第一個(gè)關(guān)鍵參數(shù)選擇如何根據(jù)模擬的物理尺寸和像素?cái)?shù)正確構(gòu)建這個(gè)頻率坐標(biāo)網(wǎng)格。2.2 第二階段全息圖記錄干涉模擬得到物光波O(x,y)后我們需要模擬它與參考光R(x,y)的干涉。參考光通常模擬為平面波或球面波。平面波最簡(jiǎn)單其復(fù)振幅可表示為R Ar * exp(1i * 2*pi/lambda * (sin(theta_x)*x sin(theta_y)*y))其中Ar是振幅常設(shè)為1theta_x和theta_y是參考光的傾斜角。這個(gè)傾斜角引入了載頻對(duì)于后續(xù)的離軸全息分離衍射級(jí)至關(guān)重要。兩者干涉后記錄平面的光強(qiáng)分布即為全息圖I_hologram abs(O R).^2。這里得到的I_hologram是一個(gè)實(shí)數(shù)值矩陣模擬了CCD相機(jī)記錄到的強(qiáng)度信息。它丟失了光波的相位但編碼了物光波的振幅和相位信息于干涉條紋中。實(shí)操心得參考光與物光的光強(qiáng)比IR/IO是一個(gè)需要仔細(xì)調(diào)節(jié)的參數(shù)。比值太大參考光過強(qiáng)全息圖條紋對(duì)比度低再現(xiàn)像信噪比差比值太小物光過強(qiáng)可能導(dǎo)致干涉條紋超過探測(cè)器的動(dòng)態(tài)范圍產(chǎn)生非線性畸變。通常將這個(gè)比值設(shè)置在3:1到10:1之間進(jìn)行仿真嘗試觀察再現(xiàn)效果。2.3 第三階段數(shù)字全息再現(xiàn)這是從全息圖中“解壓”出物體信息的過程。數(shù)字再現(xiàn)的核心是模擬參考光照射全息圖后的衍射過程。最常用的方法是菲涅爾衍射法盡管生成用角譜但再現(xiàn)常用菲涅爾近似因?yàn)橛?jì)算更直觀。再現(xiàn)過程在數(shù)學(xué)上表示為U_recon IFFT2( FFT2(I_hologram .* R_conj) .* H_prop )。這里R_conj是模擬再現(xiàn)照明光通常與參考光共軛即conj(R)用于抵消記錄時(shí)的傾斜相位使像回到中心。H_prop是菲涅爾衍射的傳遞函數(shù)形式為exp(1i*pi/(lambda*d_recon)*(fx.^2fy.^2))其中d_recon是再現(xiàn)距離通常等于記錄距離d0。計(jì)算得到的U_recon是一個(gè)復(fù)矩陣其振幅abs(U_recon)即為再現(xiàn)物體的強(qiáng)度像其相位angle(U_recon)包含了物體的三維形貌信息。對(duì)于離軸全息在頻譜上會(huì)存在三個(gè)分離的衍射級(jí)零級(jí)、正負(fù)一級(jí)我們需要通過頻域?yàn)V波提取出包含物體信息的那個(gè)一級(jí)衍射項(xiàng)再進(jìn)行逆?zhèn)鞑ヒ垣@得清晰的再現(xiàn)像。3. 關(guān)鍵參數(shù)解析與MATLAB實(shí)現(xiàn)細(xì)節(jié)一個(gè)仿真能否成功、結(jié)果是否物理可信完全取決于一系列關(guān)鍵參數(shù)的正確設(shè)置和匹配。這些參數(shù)構(gòu)成了連接數(shù)字世界與物理世界的橋梁。3.1 空間采樣與模擬尺度這是最容易出錯(cuò)的地方。在MATLAB中一切都是以像素為單位的離散數(shù)組。我們必須為這些像素賦予物理尺寸。像素尺寸delta這模擬的是CCD相機(jī)像元的物理大小例如6.45e-6 m6.45微米。它決定了仿真系統(tǒng)的空間截止頻率。網(wǎng)格大小Nx, Ny這是圖像矩陣的行列數(shù)如1024 x 1024?？偰M的物理尺寸為L(zhǎng)x Nx * delta。波長(zhǎng)lambda模擬激光的波長(zhǎng)如氦氖激光的632.8e-9 m。記錄距離d0物體平面到記錄平面的距離。這個(gè)距離不能隨便設(shè)必須滿足菲涅爾近似或角譜傳播的采樣條件以避免混疊。它們之間的約束關(guān)系由采樣定理決定。對(duì)于角譜傳播需要滿足d0 delta * Lx / lambda以避免頻域混疊。在編程時(shí)我們首先根據(jù)lambda、delta和期望的視場(chǎng)Lx來估算最大允許的d0或者先確定d0再反推所需的delta。MATLAB實(shí)現(xiàn)時(shí)構(gòu)建坐標(biāo)網(wǎng)格的代碼至關(guān)重要lambda 632.8e-9; % 波長(zhǎng) delta 6.45e-6; % 像素尺寸 N 1024; % 像素?cái)?shù) L N * delta; % 總物理尺寸 % 空間坐標(biāo) x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 頻率坐標(biāo) fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy);注意使用meshgrid生成網(wǎng)格并且頻率坐標(biāo)的構(gòu)建方式這是后續(xù)所有FFT運(yùn)算的基礎(chǔ)。3.2 參考光設(shè)計(jì)與載頻控制對(duì)于離軸全息參考光傾斜角的選擇直接決定了全息圖頻譜中各級(jí)次的分離程度。參考光波矢在x方向的投影為k_x 2*pi/lambda * sin(theta_x)。在頻譜上這會(huì)使得物光信息即1級(jí)的中心從零頻點(diǎn)移動(dòng)到(f_x0, f_y0) (sin(theta_x)/lambda, sin(theta_y)/lambda)。為了在再現(xiàn)時(shí)能完美分離出1級(jí)需要滿足分離條件f_x0必須大于物光頻譜的帶寬B約等于物體尺寸除以lambda*d0的1.5倍即f_x0 1.5 * B。否則各級(jí)頻譜會(huì)重疊產(chǎn)生串?dāng)_。采樣條件f_x0 B必須小于奈奎斯特頻率1/(2*delta)否則會(huì)發(fā)生混疊。在MATLAB中我們通過調(diào)整theta_x來滿足這些條件。通常先估算物光帶寬B然后設(shè)置f_x0 2 * B左右再反推theta_x asin(lambda * f_x0)。3.3 相位解包裹與像質(zhì)評(píng)價(jià)數(shù)字全息再現(xiàn)得到的是包裹相位值域在[-π, π]對(duì)于測(cè)量物體三維形貌需要進(jìn)行相位解包裹。MATLAB中有unwrap函數(shù)但對(duì)于噪聲大或不連續(xù)的相位圖需要更穩(wěn)健的算法如最小二乘法、質(zhì)量圖導(dǎo)引法。仿真中因?yàn)閿?shù)據(jù)干凈一維或二維的unwrap通常就夠用。評(píng)價(jià)再現(xiàn)像質(zhì)量除了主觀觀察常用客觀指標(biāo)均方誤差MSE比較再現(xiàn)像振幅與原始物體圖像的差異。峰值信噪比PSNR基于MSE計(jì)算值越高越好。結(jié)構(gòu)相似性SSIM從亮度、對(duì)比度、結(jié)構(gòu)三方面評(píng)價(jià)圖像相似性更符合人眼感知。在仿真中我們可以通過計(jì)算這些指標(biāo)定量分析不同噪聲水平、不同參數(shù)誤差對(duì)成像質(zhì)量的影響。4. 完整MATLAB仿真流程與代碼實(shí)現(xiàn)下面我們將上述思路整合成一個(gè)可運(yùn)行的、模塊化的MATLAB仿真示例。我們將模擬一個(gè)簡(jiǎn)單的方形孔徑作為物體進(jìn)行離軸菲涅爾全息記錄與再現(xiàn)。4.1 步驟一初始化參數(shù)與創(chuàng)建物體%% 1. 參數(shù)初始化 clear; close all; clc; % 物理參數(shù) lambda 632.8e-9; % 波長(zhǎng)單位米 (He-Ne激光) k 2 * pi / lambda; % 波數(shù) delta 6.45e-6; % CCD像素尺寸單位米 N 1024; % 像素?cái)?shù) (假設(shè)為正方形) L N * delta; % 總模擬尺寸單位米 d0 0.5; % 記錄距離單位米 (需滿足采樣條件) % 參考光參數(shù) (離軸角) theta_x 0.5 * pi / 180; % x方向傾斜角單位弧度 (0.5度) theta_y 0; % y方向無傾斜 Ar 1.0; % 參考光振幅 % 坐標(biāo)網(wǎng)格 x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 頻率坐標(biāo) (用于角譜傳播) fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy); %% 2. 創(chuàng)建模擬物體 % 生成一個(gè)方形孔徑 obj_size 2e-3; % 物體尺寸2mm obj double(abs(X) obj_size/2 abs(Y) obj_size/2); % 可以添加相位信息模擬一個(gè)傾斜的相位物體 phase_obj 0.5 * pi * X / (obj_size/2); % 線性相位傾斜 U_obj obj .* exp(1i * phase_obj); % 物體平面復(fù)振幅 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, abs(U_obj)); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物體振幅分布‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, angle(U_obj)); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物體相位分布包裹‘);這段代碼定義了所有核心物理參數(shù)并創(chuàng)建了一個(gè)帶有線性相位變化的方形物體。坐標(biāo)網(wǎng)格的構(gòu)建是后續(xù)所有運(yùn)算的基石。注意我們將單位從米轉(zhuǎn)換到毫米進(jìn)行顯示更符合視覺習(xí)慣。4.2 步驟二角譜傳播與全息圖記錄%% 3. 角譜傳播計(jì)算物體到記錄平面的光場(chǎng) % 角譜傳遞函數(shù) H_as exp(1i * 2*pi*d0/lambda * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 對(duì)物體場(chǎng)進(jìn)行FFT乘以傳遞函數(shù)再IFFT U_obj_fft fft2(fftshift(U_obj)); % 注意fftshift將零頻移到中心與我們的頻率坐標(biāo)匹配 U_rec_fft U_obj_fft .* H_as; U_rec ifftshift(ifft2(U_rec_fft)); % ifftshift將結(jié)果移回標(biāo)準(zhǔn)順序 % 物光波在記錄平面的振幅 Ao abs(U_rec); % 為了模擬實(shí)際情況可以給物光添加一個(gè)衰減使其強(qiáng)度弱于參考光 Ao Ao / max(Ao(:)) * 0.3; % 歸一化后調(diào)整相對(duì)強(qiáng)度 %% 4. 生成參考光波并干涉記錄全息圖 % 生成平面參考光波帶有離軸角 R Ar * exp(1i * k * (sin(theta_x)*X sin(theta_y)*Y)); % 記錄平面總光場(chǎng) U_total U_rec .* (Ao ./ abs(U_rec)) R; % 保持U_rec的相位但使用調(diào)整后的振幅Ao % 全息圖光強(qiáng)分布 I_hologram abs(U_total).^2; % 顯示全息圖 subplot(1,3,3); imagesc(x*1e3, y*1e3, I_hologram); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘記錄的全息圖‘);這里有幾個(gè)關(guān)鍵點(diǎn)fftshift與ifftshift由于我們構(gòu)建的頻率坐標(biāo)FX,FY是以零頻為中心的所以在對(duì)空間域信號(hào)做FFT前需要用fftshift將信號(hào)零頻也移到中心與傳遞函數(shù)對(duì)齊。運(yùn)算完成后再用ifftshift移回來。物光強(qiáng)度調(diào)整通過Ao Ao / max(Ao(:)) * 0.3將物光峰值振幅設(shè)為參考光振幅的0.3倍大致符合IR/IO ≈ 10:1的強(qiáng)度比以獲得高對(duì)比度干涉條紋。全息圖I_hologram是模擬CCD實(shí)際采集到的數(shù)據(jù)它是一個(gè)實(shí)值矩陣丟失了相位信息但包含了重建所需的一切。4.3 步驟三數(shù)字再現(xiàn)與像分離%% 5. 數(shù)字全息再現(xiàn) % 5.1 頻域?yàn)V波分離衍射級(jí) I_hologram_fft fft2(I_hologram); I_hologram_fft_shifted fftshift(I_hologram_fft); % 將零頻移到中心以便觀察 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(log(1 abs(I_hologram_fft_shifted))); axis image; colormap(‘jet‘); title(‘全息圖頻譜對(duì)數(shù)顯示‘); xlabel(‘空間頻率 f_x‘); ylabel(‘空間頻率 f_y‘); % 可以觀察到三個(gè)亮斑中心是零級(jí)兩側(cè)是正負(fù)一級(jí)。 % 創(chuàng)建濾波器提取1級(jí) [fxx, fyy] meshgrid(1:N, 1:N); % 估算1級(jí)中心位置對(duì)應(yīng)參考光載頻 f0_x round(N/2 sin(theta_x) * d0 / (lambda * delta)); % 近似計(jì)算 f0_y round(N/2); filter_radius 50; % 濾波器半徑需小于載頻與零級(jí)的距離 % 生成圓形帶通濾波器 filter_mask double((fxx - f0_x).^2 (fyy - f0_y).^2 filter_radius^2); % 應(yīng)用濾波器 I_filtered_fft I_hologram_fft .* filter_mask; subplot(1,3,2); imagesc(filter_mask); axis image; title(‘頻域?yàn)V波器‘); subplot(1,3,3); imagesc(log(1 abs(fftshift(I_filtered_fft)))); axis image; colormap(‘jet‘); title(‘濾波后的頻譜1級(jí)‘); % 5.2 菲涅爾衍射法再現(xiàn) % 構(gòu)建菲涅爾衍射傳遞函數(shù)再現(xiàn)距離為-d0即共軛再現(xiàn) d_recon -d0; % 負(fù)號(hào)表示反向傳播 H_fresnel exp(1i * pi/(lambda * d_recon) * (FX.^2 FY.^2) * (delta^2 * N^2)); % 注意這里FX,FY是歸一化頻率需要轉(zhuǎn)換為實(shí)際頻率并考慮離散采樣效應(yīng) % 更標(biāo)準(zhǔn)的寫法是使用空間坐標(biāo)構(gòu)建傳遞函數(shù) % H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 使用空間坐標(biāo)構(gòu)建傳遞函數(shù)更直觀 H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 再現(xiàn)過程濾波后的全息圖乘以共軛參考光再進(jìn)行菲涅爾衍射 R_conj conj(R); % 共軛參考光用于消除載頻 U_temp ifft2(I_filtered_fft) .* R_conj; % 回到空域并消除傾斜相位 % 菲涅爾衍射積分通過卷積計(jì)算先FFT乘傳遞函數(shù)再IFFT U_recon_fft fft2(U_temp) .* fftshift(H_fresnel); % 注意傳遞函數(shù)需要fftshift對(duì)齊 U_recon ifft2(U_recon_fft); % 提取再現(xiàn)像的振幅和相位 amp_recon abs(U_recon); phase_recon angle(U_recon); figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, amp_recon); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再現(xiàn)像振幅‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, phase_recon); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再現(xiàn)像相位包裹‘);這一步是整個(gè)仿真的核心。頻域?yàn)V波是關(guān)鍵操作濾波器的大小和位置直接影響再現(xiàn)像的質(zhì)量和分辨率。濾波器半徑filter_radius需要足夠大以包含全部物體頻譜信息但又不能太大以至于包含零級(jí)或其他級(jí)的成分這需要根據(jù)全息圖頻譜圖手動(dòng)調(diào)整或通過算法自動(dòng)估計(jì)。注意事項(xiàng)菲涅爾衍射傳遞函數(shù)H_fresnel的構(gòu)建有兩種常見方式一種在頻率域使用FX, FY一種在空間域使用X, Y。兩者在數(shù)學(xué)上等價(jià)但離散化計(jì)算時(shí)要注意坐標(biāo)縮放因子。使用空間域形式exp(1i*k/(2*d)*(X.^2Y.^2))通常更直觀且不易出錯(cuò)但計(jì)算量稍大。在仿真中我們更關(guān)注正確性因此推薦空間域形式。4.4 步驟四相位解包裹與結(jié)果分析%% 6. 相位解包裹與結(jié)果分析 % 相位解包裹 (使用MATLAB內(nèi)置的unwrap對(duì)于仿真簡(jiǎn)單相位通常有效) phase_unwrapped unwrap(phase_recon, [], 1); % 先按行解包裹 phase_unwrapped unwrap(phase_unwrapped, [], 2); % 再按列解包裹 % 去除傾斜相位因?yàn)槲覀兡M的物體相位本身就是傾斜的這里減去一個(gè)平面擬合值作為演示 % 實(shí)際上這一步在定量相位測(cè)量中用于消除系統(tǒng)誤差。 subplot(1,3,3); imagesc(x*1e3, y*1e3, phase_unwrapped); axis image; colormap(‘jet‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再現(xiàn)像相位解包裹后‘); colorbar; %% 7. 像質(zhì)評(píng)價(jià) % 裁剪出中心區(qū)域與原始物體進(jìn)行比較 crop_ratio 0.3; % 裁剪比例 crop_N round(N * crop_ratio); center_idx N/2 (-crop_N/2 : crop_N/2-1); center_idx round(center_idx); amp_original_crop abs(U_obj(center_idx, center_idx)); amp_recon_crop amp_recon(center_idx, center_idx); % 歸一化 amp_original_crop amp_original_crop / max(amp_original_crop(:)); amp_recon_crop amp_recon_crop / max(amp_recon_crop(:)); % 計(jì)算均方誤差(MSE)和峰值信噪比(PSNR) mse mean((amp_original_crop(:) - amp_recon_crop(:)).^2); max_val 1; % 歸一化后最大值為1 psnr 10 * log10(max_val^2 / mse); fprintf(‘圖像質(zhì)量評(píng)價(jià)\n‘); fprintf(‘ 均方誤差 (MSE): %.4e\n‘ mse); fprintf(‘ 峰值信噪比 (PSNR): %.2f dB\n‘ psnr); % 顯示對(duì)比 figure(‘Position‘ [100 100 800 400]); subplot(1,2,1); imagesc(amp_original_crop); axis image; colormap(‘gray‘); title(‘原始物體裁剪后‘); subplot(1,2,2); imagesc(amp_recon_crop); axis image; colormap(‘gray‘); title(‘再現(xiàn)像裁剪后‘); sgtitle(sprintf(‘PSNR %.2f dB‘ psnr));相位解包裹是獲取連續(xù)相位分布的必要步驟。MATLAB的unwrap函數(shù)對(duì)仿真生成的、噪聲低的相位圖效果很好。但在實(shí)際實(shí)驗(yàn)數(shù)據(jù)中由于噪聲、陰影和相位跳變可能需要更復(fù)雜的算法如phase_unwrap工具箱中的算法。像質(zhì)評(píng)價(jià)環(huán)節(jié)讓我們能定量評(píng)估仿真系統(tǒng)的性能。PSNR值越高說明再現(xiàn)像與原始物體越接近。在理想仿真中無噪聲參數(shù)完美匹配PSNR可以非常高60 dB。通過引入噪聲或參數(shù)誤差可以觀察PSNR如何下降從而理解系統(tǒng)對(duì)各因素的敏感度。5. 仿真中的典型問題、調(diào)試技巧與進(jìn)階應(yīng)用即使按照上述流程初學(xué)者在仿真中仍會(huì)遇到各種問題。下面是一些常見“坑”及其排查思路。5.1 問題一再現(xiàn)像一片模糊或出現(xiàn)鬼影可能原因1頻譜濾波不徹底零級(jí)或共軛像干擾。排查仔細(xì)檢查全息圖的頻譜圖log(1abs(fftshift(fft2(I_hologram)))。你是否能看到三個(gè)明顯分離的亮斑如果零級(jí)和1級(jí)靠得太近說明參考光載頻theta_x太小。解決增大參考光傾斜角theta_x重新計(jì)算。確保f_x0 1.5 * B。排查檢查你應(yīng)用的頻域?yàn)V波器。用imagesc(filter_mask)顯示濾波器看其位置是否準(zhǔn)確覆蓋了1級(jí)頻譜且沒有包含零級(jí)中心。解決調(diào)整濾波器的中心坐標(biāo)(f0_x, f0_y)和半徑filter_radius?？梢試L試先手動(dòng)在頻譜圖上選取區(qū)域。可能原因2再現(xiàn)距離d_recon設(shè)置錯(cuò)誤。排查再現(xiàn)距離理論上應(yīng)等于記錄距離d0共軛再現(xiàn)。如果使用菲涅爾衍射法嘗試微調(diào)d_recon的值觀察再現(xiàn)像是否變得清晰。可以寫一個(gè)循環(huán)讓d_recon在d0附近微小變化尋找圖像最清晰的點(diǎn)聚焦。解決使用自動(dòng)聚焦算法。常用方法是定義一個(gè)清晰度評(píng)價(jià)函數(shù)如圖像梯度平方和遍歷一系列d_recon取函數(shù)值最大的距離作為最佳再現(xiàn)距離。5.2 問題二再現(xiàn)像邊緣出現(xiàn)周期性條紋或振鈴效應(yīng)可能原因頻譜泄露與吉布斯現(xiàn)象。分析當(dāng)物體是理想的矩形陡峭邊緣時(shí)其頻譜是無限的sinc函數(shù)。我們用有限大小的頻域?yàn)V波器去截?cái)嗨喈?dāng)于在空域與一個(gè)sinc函數(shù)卷積導(dǎo)致邊緣出現(xiàn)振蕩。解決對(duì)原始物體加窗在生成物體U_obj時(shí)對(duì)其振幅分布乘以一個(gè)緩變的窗函數(shù)如高斯窗、漢寧窗使邊緣平滑過渡。增大濾波器尺寸適當(dāng)增加filter_radius包含更多高頻分量但要注意不要引入其他級(jí)的干擾。使用更優(yōu)的濾波器將圓形二值濾波器改為高斯衰減濾波器即filter_mask exp(-((fxx-f0_x).^2(fyy-f0_y).^2)/(2*sigma^2))可以平滑截?cái)鄿p少振鈴。5.3 問題三計(jì)算速度慢特別是對(duì)大尺寸圖像分析角譜傳播和菲涅爾衍射涉及大量FFT運(yùn)算N1024時(shí)很快但當(dāng)N4096或更大時(shí)計(jì)算和內(nèi)存消耗顯著增加。優(yōu)化技巧使用單精度如果精度要求可接受將數(shù)據(jù)轉(zhuǎn)換為單精度single。U_obj single(U_obj);FFT在單精度下更快內(nèi)存減半。預(yù)計(jì)算傳遞函數(shù)如果參數(shù)不變可以將H_as或H_fresnel計(jì)算一次并保存避免在循環(huán)中重復(fù)計(jì)算。利用GPU如果MATLAB安裝了Parallel Computing Toolbox且擁有支持CUDA的NVIDIA GPU可以使用gpuArray將數(shù)據(jù)轉(zhuǎn)移到GPU上計(jì)算。FFT在GPU上對(duì)大規(guī)模數(shù)據(jù)有巨大加速。U_obj_gpu gpuArray(U_obj); H_as_gpu gpuArray(H_as); U_rec_gpu ifft2(fft2(U_obj_gpu) .* H_as_gpu); U_rec gather(U_rec_gpu); % 將結(jié)果取回CPU減少不必要的可視化在調(diào)試完成后關(guān)閉中間的圖形顯示 (close all;)或使用set(0,‘DefaultFigureVisible‘,‘off‘)禁止圖形彈出可以節(jié)省大量時(shí)間。5.4 進(jìn)階應(yīng)用引入噪聲與像差仿真一個(gè)更貼近現(xiàn)實(shí)的仿真需要引入噪聲和像差。添加噪聲模擬CCD讀出噪聲、散粒噪聲等。SNR_dB 20; % 信噪比 I_hologram_noiseless I_hologram; signal_power mean(I_hologram(:).^2); noise_power signal_power / (10^(SNR_dB/10)); noise sqrt(noise_power) * randn(size(I_hologram)); % 高斯白噪聲 I_hologram I_hologram_noiseless noise; I_hologram(I_hologram 0) 0; % 確保強(qiáng)度非負(fù)通過改變SNR_dB可以研究噪聲對(duì)再現(xiàn)像質(zhì)量PSNR的影響。引入像差模擬光學(xué)系統(tǒng)的不完美如球差、彗差、像散等。這可以在角譜傳遞函數(shù)H_as或參考光波前R上乘以一個(gè)像差相位板W。% 例如引入初級(jí)球差 r2 (X.^2 Y.^2) / (L/2)^2; % 歸一化半徑 W_spherical 2 * pi / lambda * 0.1e-6 * r2.^2; % 0.1微米的球差 H_as_aberrated H_as .* exp(1i * W_spherical);觀察像差如何導(dǎo)致再現(xiàn)像模糊、變形從而理解像差校正如數(shù)字相位補(bǔ)償?shù)闹匾浴Ｍㄟ^這個(gè)完整的MATLAB數(shù)字全息仿真框架你不僅能夠復(fù)現(xiàn)教科書中的經(jīng)典現(xiàn)象更能將其作為一個(gè)靈活的工具箱用于探索更復(fù)雜的全息成像問題如相移全息、彩色全息、顯微全息等為真正的光學(xué)實(shí)驗(yàn)或工程應(yīng)用打下堅(jiān)實(shí)的理論和實(shí)踐基礎(chǔ)。本文還有配套的精品資源點(diǎn)擊獲取