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