字全息離軸干涉:從原理到代碼實(shí)現(xiàn))
簡(jiǎn)介本資源是一套面向光學(xué)工程、數(shù)字全息與計(jì)算成像方向的本科高年級(jí)或研究生實(shí)踐教學(xué)材料聚焦像面數(shù)字全息中離軸干涉圖樣的仿真建模與相位恢復(fù)全流程。資源完整實(shí)現(xiàn)從物光波前生成、參考光引入、干涉圖模擬、頻譜濾波、傅里葉逆變換到相位解包裹與強(qiáng)度再現(xiàn)的全部核心算法適用于課程設(shè)計(jì)、實(shí)驗(yàn)復(fù)現(xiàn)及算法原理驗(yàn)證。壓縮包共24個(gè)文件含6個(gè)Matlab源碼如fresnel.m、unwrap_LS.m、New_FFT2.m等關(guān)鍵函數(shù)、13張運(yùn)行結(jié)果圖jpg格式涵蓋原始干涉圖、頻譜圖、相位像、強(qiáng)度再現(xiàn)像等中間與最終可視化結(jié)果以及5幅bmp格式的原始/處理圖像如usaf.bmp、0相位像.bmp整體僅869KB輕量易用。已有476人學(xué)習(xí)下載配套代碼結(jié)構(gòu)清晰、變量命名規(guī)范、關(guān)鍵步驟附注釋可直接運(yùn)行并支持參數(shù)調(diào)整特別適合作為數(shù)字全息入門(mén)算法理解與Matlab光學(xué)仿真能力提升的實(shí)操范例。1. 項(xiàng)目概述從一束光到一幅圖如果你對(duì)光學(xué)、圖像處理或者計(jì)算成像感興趣手里又恰好有一份名為“模擬像面數(shù)字全息離軸干涉含Matlab源碼.zip”的文件那你算是來(lái)對(duì)地方了。這不僅僅是一個(gè)壓縮包它更像是一把鑰匙能幫你打開(kāi)一扇通往“計(jì)算光學(xué)”這個(gè)交叉領(lǐng)域的大門(mén)。簡(jiǎn)單來(lái)說(shuō)這個(gè)項(xiàng)目就是用計(jì)算機(jī)軟件Matlab來(lái)模擬和重現(xiàn)一個(gè)經(jīng)典的光學(xué)實(shí)驗(yàn)過(guò)程——離軸數(shù)字全息干涉并且是在“像面”這個(gè)特殊的位置上進(jìn)行。聽(tīng)起來(lái)有點(diǎn)繞別急我干了十幾年光學(xué)仿真和圖像處理今天就用最“人話”的方式帶你把這個(gè)項(xiàng)目從里到外、從原理到代碼徹底盤(pán)明白。數(shù)字全息你可以把它理解成給光波“拍一張包含所有信息的立體照片”。傳統(tǒng)相機(jī)只能記錄光的強(qiáng)度也就是明暗但全息技術(shù)能同時(shí)記錄下光的強(qiáng)度和相位你可以理解為光波的“形狀”或“起伏”。離軸干涉則是實(shí)現(xiàn)這種記錄的一種巧妙方法它通過(guò)引入一個(gè)傾斜的參考光把物體的相位信息“編碼”到一幅我們能直接拍攝到的干涉條紋圖中。而“模擬像面”指的是我們模擬的這個(gè)全息圖不是在物體所在的平面記錄的而是在經(jīng)過(guò)一個(gè)透鏡成像后的像平面上記錄的。這么做的好處是物體的像和全息圖在同一平面上重建過(guò)程更直觀也更貼近很多實(shí)際顯微成像的應(yīng)用場(chǎng)景。這個(gè)Matlab源碼包的價(jià)值在于它把一套復(fù)雜的光學(xué)理論和數(shù)學(xué)運(yùn)算封裝成了你可以一行行運(yùn)行、一步步修改的代碼。你不需要昂貴的激光器、精密的光學(xué)平臺(tái)和防震臺(tái)只需要一臺(tái)電腦就能親手“搭建”這個(gè)光學(xué)系統(tǒng)觀察干涉圖樣的生成并最終通過(guò)計(jì)算“再現(xiàn)”出物體的三維光場(chǎng)。無(wú)論是學(xué)生做課程設(shè)計(jì)、研究者驗(yàn)證算法還是工程師進(jìn)行系統(tǒng)前期仿真這都是一個(gè)極佳的起點(diǎn)。接下來(lái)我會(huì)帶你拆解它的核心思路、手把手過(guò)一遍代碼、分享我踩過(guò)的坑和調(diào)試技巧讓你不僅能運(yùn)行出結(jié)果更能理解每一步背后的“所以然”。2. 核心原理拆解光是如何被“計(jì)算”的要玩轉(zhuǎn)這個(gè)源碼我們不能只當(dāng)個(gè)“調(diào)包俠”或者“跑碼員”必須得弄清楚它到底在算些什么。這里面的核心物理是“干涉”核心數(shù)學(xué)是“傅里葉變換”和“角譜傳播”核心操作是“模擬”和“重建”。我們一層層來(lái)看。2.1 離軸全息干涉的編碼原理為什么離軸干涉這么重要在早期的同軸全息參考光和物光方向一致中重建的物體像會(huì)和共軛像、零級(jí)衍射光重疊在一起畫(huà)面非常混亂。離軸技術(shù)就像在合唱中讓領(lǐng)唱站到一邊他的聲音物光信息就能從和聲零級(jí)光和其他噪聲中被清晰地分辨出來(lái)。在模擬中我們通常用兩個(gè)復(fù)數(shù)矩陣來(lái)表示物光波 ( O(x, y) ) 和參考光波 ( R(x, y) )。物光波攜帶著物體信息它的振幅和相位由物體決定。參考光波通常被模擬成一個(gè)理想的平面波其數(shù)學(xué)表達(dá)式為 ( R(x, y) A_r \exp(i 2\pi (f_x x f_y y)) )。這里的 ( A_r ) 是常數(shù)振幅關(guān)鍵就在于 ( (f_x, f_y) )它們是參考光的空間頻率決定了光波的傾斜角度。正是這個(gè)傾斜使得物光和參考光干涉后物體的相位信息被轉(zhuǎn)換成了干涉條紋的疏密和走向。兩者在像面上疊加其光強(qiáng)分布也就是我們模擬能“拍”到的全息圖 ( I_H(x, y) )為 [ I_H |O R|^2 |O|^2 |R|^2 O^R OR^] 這個(gè)公式是全息的靈魂。( |O|^2 ) 是物光自身強(qiáng)度( |R|^2 ) 是參考光強(qiáng)度它們合起來(lái)構(gòu)成緩慢變化的背景零級(jí)項(xiàng)。而 ( O^R ) 和 ( OR^) 是交叉項(xiàng)其中就編碼了我們需要的物體相位信息。在離軸設(shè)置下這兩個(gè)交叉項(xiàng)在頻域傅里葉變換后的空間里會(huì)和零級(jí)項(xiàng)分離開(kāi)這是后續(xù)我們能干凈利落地提取出物體信息的前提。2.2 像面全息與角譜傳播為什么強(qiáng)調(diào)“像面”在菲涅耳全息或夫瑯禾費(fèi)全息中全息圖記錄面距離物體較遠(yuǎn)。而在像面全息中我們假設(shè)物體已經(jīng)通過(guò)一個(gè)理想透鏡完美成像在了記錄平面上。這意味著記錄平面上的物光復(fù)振幅 ( O(x, y) )直接就是物體透過(guò)率函數(shù)經(jīng)過(guò)透鏡成像后的結(jié)果無(wú)需考慮復(fù)雜的衍射傳播過(guò)程。這極大地簡(jiǎn)化了模擬的物理模型。但是當(dāng)我們從全息圖重建物體光場(chǎng)時(shí)我們通常需要“回溯”光場(chǎng)到物體平面或者至少到某個(gè)觀察平面。這時(shí)就需要用到“角譜傳播”理論。角譜傳播法是計(jì)算光波在自由空間中傳播的精確方法之一在標(biāo)量衍射理論范圍內(nèi)。其核心思想是任何一個(gè)光波場(chǎng)都可以分解成無(wú)數(shù)個(gè)不同方向傳播的平面波角譜。傳播一段距離 ( z ) 后每個(gè)平面波成分只會(huì)產(chǎn)生一個(gè)相位延遲。在頻域里這個(gè)操作變得異常簡(jiǎn)單。設(shè)輸入平面光場(chǎng)為 ( U(x, y; 0) )其二維傅里葉變換為 ( \mathcal{F}{U(x, y; 0)} G(f_x, f_y; 0) )其中 ( f_x, f_y ) 是空間頻率。那么傳播距離 ( z ) 后的光場(chǎng)頻譜為 [ G(f_x, f_y; z) G(f_x, f_y; 0) \cdot H(f_x, f_y) ] 其中傳遞函數(shù) ( H(f_x, f_y) \exp(i k z \sqrt{1 - (\lambda f_x)^2 - (\lambda f_y)^2}) )這里 ( k2\pi/\lambda ) 是波數(shù)( \lambda ) 是光波長(zhǎng)。這個(gè)公式是源碼中重建部分的核心。通過(guò)一次傅里葉變換、乘以一個(gè)相位因子、再做一次逆傅里葉變換我們就完成了光場(chǎng)的傳播模擬。這種方法特別適合在計(jì)算機(jī)上實(shí)現(xiàn)也是數(shù)字全息重建的主流算法。2.3 數(shù)字重建的四步流程源碼的算法骨架無(wú)論是模擬還是重建都遵循一個(gè)清晰的四步流程理解了它看代碼就事半功倍模擬生成全息圖構(gòu)建或加載物體圖像 - 賦予其相位分布模擬物體形貌- 將其作為物光波 - 模擬透鏡成像效應(yīng)可能包含衍射- 與傾斜的參考光干涉 - 計(jì)算光強(qiáng)得到模擬全息圖。這一步的輸出是一幅灰度圖像。頻域?yàn)V波與信息提取對(duì)模擬全息圖進(jìn)行二維傅里葉變換得到其頻譜。在頻譜圖中你會(huì)看到中心亮團(tuán)零級(jí)項(xiàng)和兩側(cè)的“邊帶”交叉項(xiàng)。通過(guò)設(shè)計(jì)一個(gè)合適的濾波器通常是圓形或矩形窗函數(shù)將其中一個(gè)邊帶比如 ( O^*R ) 項(xiàng)完整地截取出來(lái)并移動(dòng)到頻譜中心。這一步的目的是分離出攜帶物體信息的單項(xiàng)。角譜傳播重建對(duì)濾波并移中后的頻譜進(jìn)行逆傅里葉變換得到復(fù)數(shù)形式的像面光場(chǎng)分布。如果需要觀察不同深度的重建像就利用上面提到的角譜傳播公式將這個(gè)光場(chǎng)向不同距離 ( z ) 進(jìn)行“回溯”或“聚焦”計(jì)算。結(jié)果顯示與解包裹計(jì)算重建光場(chǎng)的振幅和相位。振幅圖就是物體的強(qiáng)度像比較直觀。相位圖則包含了物體的表面高度或折射率變化信息。但直接計(jì)算出的相位值被包裹在 ( [-\pi, \pi] ) 的主值區(qū)間內(nèi)呈鋸齒狀需要通過(guò)“相位解包裹”算法將其恢復(fù)成連續(xù)的相位分布這才是最終可用的三維形貌數(shù)據(jù)。注意在像面全息模擬中第一步的“透鏡成像效應(yīng)”有時(shí)會(huì)被簡(jiǎn)化甚至直接認(rèn)為物光波就是物體本身如果物體是二維圖像。而在重建時(shí)如果只關(guān)心像面本身傳播距離 ( z ) 可以設(shè)為0。源碼的靈活性往往就體現(xiàn)在這些地方。3. 源碼結(jié)構(gòu)深度解析與關(guān)鍵函數(shù)剖析拿到一個(gè)“含Matlab源碼”的壓縮包最忌諱的就是直接運(yùn)行main.m然后看個(gè)結(jié)果圖了事。我們必須像解剖一樣看清它的每一塊肌肉和骨骼。下面我以一個(gè)典型的、結(jié)構(gòu)良好的模擬像面離軸全息Matlab項(xiàng)目為例拆解其核心文件和作用。你的源碼包可能文件名不同但功能模塊大抵如此。3.1 主腳本 (main.m或demo_simulate_reconstruct.m)這是程序的入口通常負(fù)責(zé)設(shè)置全局參數(shù)、調(diào)用各個(gè)功能函數(shù)、控制流程和顯示結(jié)果。一個(gè)好的主腳本應(yīng)該像一份實(shí)驗(yàn)報(bào)告提綱。%% 1. 清空與初始化 clear; close all; clc; addpath(genpath(./utils)); % 添加工具函數(shù)路徑 %% 2. 系統(tǒng)參數(shù)設(shè)置 lambda 632.8e-9; % 波長(zhǎng)單位米 (He-Ne激光) k 2*pi/lambda; % 波數(shù) pixel_size 6.5e-6; % CCD像元尺寸單位米 Nx 1024; Ny 1024; % 全息圖尺寸像素 Lx Nx * pixel_size; % 全息圖物理寬度 Ly Ny * pixel_size; % 全息圖物理高度 % 參考光傾斜角度通常用離軸頻率表示 theta_x 1.0 * pi/180; % x方向傾斜角弧度 theta_y 0.5 * pi/180; % y方向傾斜角弧度 fx_r sin(theta_x) / lambda; % 參考光空間頻率x分量 fy_r sin(theta_y) / lambda; % 參考光空間頻率y分量 %% 3. 生成或加載物體 object_amp im2double(imread(test_object.png)); % 物體振幅 object_phase generate_phase_map(Nx, Ny); % 物體相位模擬形貌 object_field object_amp .* exp(1i * object_phase); % 物光復(fù)振幅 %% 4. 模擬像面全息圖記錄 [hologram, ref_field] simulate_imgholo(object_field, pixel_size, lambda, fx_r, fy_r); %% 5. 全息圖重建 recon_distance 0; % 像面重建距離為0 [recon_amp, recon_phase] reconstruct_imgholo(hologram, pixel_size, lambda, fx_r, fy_r, recon_distance); %% 6. 結(jié)果顯示 figure; subplot(2,3,1); imshow(object_amp, []); title(原始物體振幅); subplot(2,3,2); imshow(object_phase, []); title(原始物體相位); subplot(2,3,3); imshow(hologram, []); title(模擬像面全息圖); subplot(2,3,4); imshow(abs(fftshift(fft2(hologram))), []); title(全息圖頻譜對(duì)數(shù)); subplot(2,3,5); imshow(recon_amp, []); title(重建振幅); subplot(2,3,6); imshow(recon_phase, []); title(重建包裹相位);在主腳本中參數(shù)設(shè)置部分是重中之重。pixel_size和lambda的比例決定了系統(tǒng)的最大記錄角度。fx_r,fy_r的設(shè)置必須足夠大以確保頻譜中的各項(xiàng)充分分離但又不能太大否則干涉條紋過(guò)密超過(guò)CCD的奈奎斯特采樣頻率會(huì)導(dǎo)致混疊。一個(gè)經(jīng)驗(yàn)法則是參考光引入的載頻應(yīng)使邊帶中心與頻譜中心的距離大于物體帶寬的3倍左右。3.2 全息圖模擬函數(shù) (simulate_imgholo.m)這個(gè)函數(shù)是物理模型的核心。它精確地模擬了物光傳播到像面、與參考光干涉的過(guò)程。function [hologram, ref_field] simulate_imgholo(object_field, pixel_size, lambda, fx_r, fy_r) % 輸入物光場(chǎng)像素尺寸波長(zhǎng)參考光空間頻率 % 輸出全息圖強(qiáng)度參考光場(chǎng)可選 [Ny, Nx] size(object_field); Lx Nx * pixel_size; Ly Ny * pixel_size; % 生成坐標(biāo)網(wǎng)格 x linspace(-Lx/2, Lx/2, Nx); y linspace(-Ly/2, Ly/2, Ny); [X, Y] meshgrid(x, y); % 1. 模擬像面物光場(chǎng)這里簡(jiǎn)化處理認(rèn)為物體已完美成像在像面。 % 更復(fù)雜的模擬可以在此處加入透鏡傳遞函數(shù)或衍射傳播。 img_field object_field; % 2. 生成傾斜平面波作為參考光 % 關(guān)鍵參考光相位是線性相移對(duì)應(yīng)空間頻率 ref_field exp(1i * 2*pi * (fx_r * X fy_r * Y)); % 3. 干涉并記錄強(qiáng)度 interference_field img_field ref_field; hologram abs(interference_field).^2; % 4. (可選) 添加噪聲更接近真實(shí)實(shí)驗(yàn) % hologram imnoise(hologram, gaussian, 0, 0.01); end關(guān)鍵點(diǎn)剖析meshgrid生成的X,Y是物理坐標(biāo)網(wǎng)格單位是米。這是所有物理計(jì)算的基礎(chǔ)千萬(wàn)不能和像素索引i, j混淆。ref_field的生成是離軸干涉的靈魂代碼。fx_r * X fy_r * Y計(jì)算了每個(gè)坐標(biāo)點(diǎn)相對(duì)于原點(diǎn)的光程差乘以2π后就是相位延遲。這個(gè)線性相位項(xiàng)在傅里葉變換后會(huì)表現(xiàn)為頻譜的平移。干涉強(qiáng)度hologram是一個(gè)實(shí)數(shù)矩陣取值范圍從0到最大值。在顯示時(shí)通常需要?dú)w一化或調(diào)整對(duì)比度。3.3 全息圖重建函數(shù) (reconstruct_imgholo.m)這是數(shù)字全息的魔法所在將一幅干涉條紋圖還原為三維光場(chǎng)。function [recon_amp, recon_phase] reconstruct_imgholo(hologram, pixel_size, lambda, fx_r, fy_r, z) % 輸入全息圖像素尺寸波長(zhǎng)參考光頻重建距離z % 輸出重建平面的振幅和相位 [Ny, Nx] size(hologram); Lx Nx * pixel_size; Ly Ny * pixel_size; % 1. 全息圖頻譜分析 H fft2(hologram); H_shift fftshift(H); % 將零頻移到中心方便觀察 % 2. 設(shè)計(jì)頻域?yàn)V波器提取1級(jí)邊帶 % 生成頻域坐標(biāo) fx linspace(-1/(2*pixel_size), 1/(2*pixel_size), Nx); fy linspace(-1/(2*pixel_size), 1/(2*pixel_size), Ny); [Fx, Fy] meshgrid(fx, fy); % 預(yù)測(cè)1級(jí)邊帶中心位置對(duì)應(yīng)O*R項(xiàng) fx_center -fx_r; % 注意符號(hào)由干涉公式推導(dǎo)決定 fy_center -fy_r; % 創(chuàng)建圓形濾波器也可用矩形窗 radius min(1/(4*pixel_size), sqrt(fx_center^2 fy_center^2)/3); % 濾波器半徑經(jīng)驗(yàn)值 filter_mask ((Fx - fx_center).^2 (Fy - fy_center).^2) radius^2; % 3. 濾波并移中 H_filtered H_shift .* filter_mask; % 將濾波后的邊帶移回頻譜原點(diǎn)以消除參考光傾斜引入的線性相位 H_centered ifftshift(H_filtered); % 4. 逆傅里葉變換得到復(fù)振幅 recon_field_complex ifft2(H_centered); % 5. 角譜傳播如果z不為0 if z ~ 0 % 生成角譜傳遞函數(shù)的頻率坐標(biāo)注意與濾波時(shí)頻率坐標(biāo)的尺度可能不同 dfx 1 / Lx; dfy 1 / Ly; u ((0:Nx-1) - floor(Nx/2)) * dfx; % 注意此處的頻率坐標(biāo)排列需與fft2匹配 v ((0:Ny-1) - floor(Ny/2)) * dfy; [U, V] meshgrid(u, v); % 角譜傳遞函數(shù) H_as exp(1i * 2*pi/lambda * z .* sqrt(1 - (lambda*U).^2 - (lambda*V).^2)); % 傳播操作 recon_field_complex ifft2( fft2(recon_field_complex) .* ifftshift(H_as) ); end % 6. 提取振幅和相位 recon_amp abs(recon_field_complex); recon_phase angle(recon_field_complex); % 這是包裹相位 end這是最容易出錯(cuò)的模塊需要逐行理解頻譜坐標(biāo)生成linspace(-1/(2*pixel_size), 1/(2*pixel_size), Nx)生成的是物理頻率坐標(biāo)單位是m^{-1}。其范圍由采樣定理決定最大頻率奈奎斯特頻率是1/(2*pixel_size)。這是設(shè)計(jì)濾波器的基準(zhǔn)。邊帶中心預(yù)測(cè)fx_center -fx_r這個(gè)符號(hào)關(guān)系需要根據(jù)具體干涉公式推導(dǎo)來(lái)驗(yàn)證。一個(gè)快速驗(yàn)證方法是運(yùn)行模擬后查看全息圖的頻譜圖 (imshow(log(abs(H_shift)1), []))觀察亮斑位置然后調(diào)整符號(hào)使濾波器能對(duì)準(zhǔn)它。濾波器設(shè)計(jì)圓形濾波器的半徑radius是關(guān)鍵參數(shù)。太小會(huì)丟失物體高頻信息重建像模糊太大會(huì)包含零級(jí)或其他邊帶成分引入噪聲。通常需要根據(jù)頻譜圖手動(dòng)調(diào)整。高級(jí)方法可以用自適應(yīng)濾波。角譜傳播的頻率坐標(biāo)注意第5步中生成U, V的方式與第2步生成Fx, Fy的方式可能不同。這是因?yàn)閒ft2輸出的頻率順序是[0, 1, ..., N/2, -N/21, ..., -1]。使用ifftshift(H_as)是為了將傳遞函數(shù)的零頻調(diào)整到與fft2輸出匹配的位置。這一步是角譜傳播正確實(shí)現(xiàn)的核心細(xì)節(jié)很多開(kāi)源代碼在這里都有細(xì)微差別容易導(dǎo)致重建像散焦或出現(xiàn)偽影。3.4 工具函數(shù)集 (utils/文件夾)一個(gè)完整的項(xiàng)目通常還包含一些工具函數(shù)例如generate_phase_map.m: 生成模擬的物體相位分布如球形相位模擬微透鏡、斜坡相位模擬傾斜平面或隨機(jī)相位板。phase_unwrapping.m: 實(shí)現(xiàn)相位解包裹算法如 Goldstein枝切法、最小二乘法或基于質(zhì)量圖導(dǎo)向的路徑積分法。remove_tilt.m: 從重建相位中移除由于系統(tǒng)未對(duì)準(zhǔn)可能殘留的線性相位傾斜。calculate_metrics.m: 計(jì)算重建圖像與原始圖像的誤差指標(biāo)如均方誤差MSE、結(jié)構(gòu)相似性SSIM等。4. 從零開(kāi)始手把手實(shí)現(xiàn)與參數(shù)調(diào)試實(shí)戰(zhàn)理解了原理和代碼結(jié)構(gòu)我們現(xiàn)在可以自己動(dòng)手從一個(gè)簡(jiǎn)單的例子開(kāi)始逐步增加復(fù)雜度并學(xué)會(huì)如何調(diào)試。4.1 第一步創(chuàng)建一個(gè)簡(jiǎn)單的二維振幅物體我們從一個(gè)最簡(jiǎn)單的“雙縫”或“USAF分辨率板”圖像開(kāi)始。這能讓我們清晰地看到干涉和重建效果。% 生成一個(gè)簡(jiǎn)單的振幅物體 N 512; object_amp zeros(N, N); object_amp(200:220, 100:400) 1; % 一個(gè)水平亮條 object_amp(300:320, 100:400) 1; % 另一個(gè)水平亮條 % 或者加載一個(gè)標(biāo)準(zhǔn)分辨率板圖像 % object_amp im2double(imread(usaf1951.png)); % object_amp imresize(object_amp, [N, N]); % 初始相位設(shè)為0純振幅物體 object_phase zeros(N, N); object_field object_amp .* exp(1i * object_phase);運(yùn)行主腳本你應(yīng)該能看到清晰的干涉條紋全息圖其頻譜呈現(xiàn)對(duì)稱的三個(gè)亮斑中心零級(jí)和兩側(cè)邊帶。重建的振幅圖應(yīng)該能大致恢復(fù)出兩條亮帶。常見(jiàn)問(wèn)題1全息圖沒(méi)有條紋檢查參考光空間頻率fx_r,fy_r是否設(shè)置得太小嘗試增大角度例如從0.1度增加到1度。檢查物光場(chǎng)object_field是否為復(fù)數(shù)純實(shí)數(shù)只有振幅的物光與參考光干涉條紋對(duì)比度可能很低。確保有相位項(xiàng)exp(1i*phase)。4.2 第二步引入相位物體并觀察頻譜分離現(xiàn)在我們給物體加上相位信息模擬一個(gè)更真實(shí)的情況比如一個(gè)凸起的圓形相位物體。% 生成一個(gè)圓形相位物體 [X, Y] meshgrid(1:N, 1:N); center N/2; radius 50; phase_object zeros(N, N); ind sqrt((X-center).^2 (Y-center).^2) radius; phase_object(ind) 2 * pi * 2; % 給圓形區(qū)域一個(gè)2*2π的相位延遲一個(gè)“相位丘” object_field object_amp .* exp(1i * phase_object); % 結(jié)合振幅和相位再次運(yùn)行。觀察全息圖的頻譜 (fftshift(fft2(hologram)))。你會(huì)發(fā)現(xiàn)兩側(cè)的邊帶亮斑不再是干凈的一個(gè)點(diǎn)而是擴(kuò)散成了一個(gè)“云團(tuán)”。這個(gè)云團(tuán)的大小就代表了物體信息的帶寬。確保你的濾波器半徑大于這個(gè)云團(tuán)的半徑否則會(huì)丟失高頻細(xì)節(jié)重建的相位圖邊緣會(huì)變得平滑。4.3 第三步精細(xì)調(diào)整濾波器與重建距離這是獲得高質(zhì)量重建結(jié)果的關(guān)鍵步驟也是最需要耐心和經(jīng)驗(yàn)的地方??梢暬l譜并手動(dòng)選取濾波器H_spectrum log(abs(fftshift(fft2(hologram))) 1); figure; imshow(H_spectrum, []); title(全息圖對(duì)數(shù)頻譜); % 交互式選取濾波器中心和大小的代碼簡(jiǎn)化示意 % [fx_center, fy_center] ginput(1); % 用鼠標(biāo)點(diǎn)擊邊帶中心 % radius input(請(qǐng)輸入濾波器半徑像素: );在實(shí)際調(diào)試中我經(jīng)常先用ginput函數(shù)手動(dòng)在頻譜圖上點(diǎn)選邊帶中心獲取其像素坐標(biāo)再換算成物理頻率坐標(biāo)。這樣可以精確校準(zhǔn)fx_center,fy_center。優(yōu)化濾波器形狀圓形濾波器簡(jiǎn)單但可能會(huì)切掉方形成分??梢試L試高斯濾波器或巴特沃斯濾波器讓濾波邊緣更平滑減少吉布斯振鈴現(xiàn)象。% 高斯濾波器示例 sigma radius / 2; % 高斯標(biāo)準(zhǔn)差 gaussian_mask exp(-((Fx - fx_center).^2 (Fy - fy_center).^2) / (2*sigma^2)); filter_mask gaussian_mask;調(diào)整重建距離z對(duì)于像面全息z0通常就是最清晰的平面。但你可以嘗試微調(diào)z比如正負(fù)幾毫米觀察重建振幅圖像的清晰度變化。這可以幫你理解“數(shù)字聚焦”的過(guò)程。編寫(xiě)一個(gè)循環(huán)遍歷一系列z值計(jì)算每個(gè)z下重建圖像的銳度例如梯度平方和自動(dòng)找到最佳聚焦位置。4.4 第四步相位解包裹與結(jié)果分析直接angle()函數(shù)得到的相位是包裹的對(duì)于我們的圓形相位丘你會(huì)看到一個(gè)從-π到π反復(fù)跳變的同心圓環(huán)。% 使用簡(jiǎn)單的二維相位解包裹算法例如基于質(zhì)量圖導(dǎo)向的路徑積分法 % 假設(shè)有 phase_unwrapping 函數(shù) unwrapped_phase phase_unwrapping(recon_phase); % 可視化對(duì)比 figure; subplot(1,2,1); imshow(recon_phase, [-pi pi]); title(包裹相位); colorbar; subplot(1,2,2); imshow(unwrapped_phase, []); title(解包裹相位); colorbar; % 繪制中心線剖面 center_line unwrapped_phase(N/2, :); figure; plot(center_line); title(解包裹相位中心線剖面); xlabel(像素); ylabel(相位弧度);解包裹后的相位圖應(yīng)該呈現(xiàn)一個(gè)平滑的“丘”狀。你可以測(cè)量丘的高度相位值與我們最初設(shè)定的2 * 2π進(jìn)行對(duì)比驗(yàn)證重建的保真度。如果存在整體傾斜或彎曲可能是系統(tǒng)像差或參考光不完全準(zhǔn)直導(dǎo)致的可以使用remove_tilt函數(shù)進(jìn)行校正。5. 避坑指南與性能優(yōu)化來(lái)自一線的經(jīng)驗(yàn)在實(shí)際操作和代碼優(yōu)化中我積累了一些教科書(shū)和論文里很少提及但至關(guān)重要的經(jīng)驗(yàn)。5.1 頻譜混疊與采樣定理這是數(shù)字全息模擬中最常見(jiàn)的“坑”?,F(xiàn)象是重建圖像出現(xiàn)奇怪的波紋、扭曲或者高頻部分完全丟失。根源參考光傾斜角太大導(dǎo)致干涉條紋的空間頻率超過(guò)了奈奎斯特極限即每像素少于2個(gè)采樣點(diǎn)。診斷檢查全息圖本身。用plot(hologram(512, :))畫(huà)一行剖面線數(shù)一數(shù)條紋的周期。最小周期兩個(gè)亮紋或暗紋之間的距離對(duì)應(yīng)的像素?cái)?shù)應(yīng)大于2。解決降低參考光角度減小fx_r,fy_r。增大像素尺寸在模擬中這意味著你的CCD像元更大雖然不現(xiàn)實(shí)但可以作為驗(yàn)證。增加插值在模擬生成物光場(chǎng)和參考光場(chǎng)時(shí)使用更高的采樣率比如2倍生成全息圖后再降采樣回來(lái)。這相當(dāng)于過(guò)采樣。5.2 能量均衡與動(dòng)態(tài)范圍問(wèn)題全息圖頻譜中零級(jí)項(xiàng)中心亮斑的能量遠(yuǎn)高于邊帶導(dǎo)致濾波時(shí)邊帶信號(hào)微弱重建圖像信噪比低。解決調(diào)整物參光強(qiáng)比在模擬中可以給物光場(chǎng)乘以一個(gè)系數(shù)例如object_field 0.3 * object_field使物光與參考光強(qiáng)度匹配。經(jīng)驗(yàn)上物參光強(qiáng)比在1:1到1:10之間效果較好。頻譜濾波前進(jìn)行背景抑制從全息圖中減去其平均強(qiáng)度或進(jìn)行高通濾波可以一定程度上壓制零級(jí)。5.3 計(jì)算效率與內(nèi)存優(yōu)化當(dāng)圖像尺寸很大如2048x2048或需要批量處理時(shí)計(jì)算速度和內(nèi)存成為瓶頸。使用單精度Matlab默認(rèn)雙精度。對(duì)于大部分全息計(jì)算單精度浮點(diǎn)數(shù) (single) 足夠且能節(jié)省一半內(nèi)存提升計(jì)算速度。在數(shù)據(jù)初始化時(shí)使用single()。object_field single(object_amp) .* exp(1i * single(object_phase));使用GPU加速如果你的Matlab支持Parallel Computing Toolbox且有好顯卡將大數(shù)據(jù)轉(zhuǎn)移到GPU上計(jì)算是質(zhì)的飛躍。object_field_gpu gpuArray(object_field); hologram_gpu abs(object_field_gpu ref_field_gpu).^2; hologram gather(hologram_gpu);注意頻繁在CPU和GPU之間傳輸數(shù)據(jù)會(huì)成為新的瓶頸盡量讓整個(gè)計(jì)算流程都在GPU上完成。優(yōu)化循環(huán)與向量化盡量避免在大的for循環(huán)中進(jìn)行逐像素操作。Matlab的矩陣運(yùn)算向量化要快得多。我們之前的代碼已經(jīng)完全是向量化的。5.4 與實(shí)驗(yàn)數(shù)據(jù)的對(duì)接這個(gè)模擬項(xiàng)目的終極價(jià)值之一是指導(dǎo)真實(shí)實(shí)驗(yàn)。當(dāng)你用這套代碼處理實(shí)驗(yàn)采集的全息圖時(shí)要注意像素尺寸校準(zhǔn)實(shí)驗(yàn)中的pixel_size必須精確測(cè)量相機(jī)像元尺寸除以物鏡放大倍數(shù)。一個(gè)像素的誤差會(huì)導(dǎo)致重建尺度錯(cuò)誤。參考光頻率標(biāo)定實(shí)驗(yàn)中的fx_r,fy_r很難直接測(cè)量。最實(shí)用的方法是從實(shí)驗(yàn)全息圖的頻譜中直接讀取。用鼠標(biāo)點(diǎn)選頻譜圖中邊帶亮斑的中心像素坐標(biāo)(cx, cy)然后通過(guò)公式換算fx_r_exp -(cx - floor(Nx/2) - 1) / (Nx * pixel_size); % 注意坐標(biāo)偏移和符號(hào) fy_r_exp -(cy - floor(Ny/2) - 1) / (Ny * pixel_size);消除零級(jí)和共軛像實(shí)驗(yàn)全息圖的頻譜分離可能不如模擬的完美。除了濾波還可以采用相位偏移法Phase Shifting等記錄多幅全息圖來(lái)徹底消除零級(jí)和共軛像但這超出了本模擬項(xiàng)目的范圍。通過(guò)以上五個(gè)部分的拆解你應(yīng)該已經(jīng)從“知道這個(gè)源碼包”變成了“能駕馭、能修改、能debug、能應(yīng)用到實(shí)際問(wèn)題”的狀態(tài)。數(shù)字全息是一個(gè)深邃而有趣的領(lǐng)域這個(gè)像面離軸干涉的Matlab模擬是你探索計(jì)算成像世界一個(gè)非常扎實(shí)的起點(diǎn)。多動(dòng)手改參數(shù)多觀察中間結(jié)果的變化你會(huì)對(duì)光與計(jì)算的融合有更直覺(jué)的理解。本文還有配套的精品資源點(diǎn)擊獲取