同轴数字全息图及其重建的MATLAB程序

%-------------------LJCM11.m---------------------------------- % 功能: 四步相移法生成同轴数字全息图及1-IFFT重建模拟 %------------------------------------------------------------- % 1,读取一幅图像,加相位噪声 % 2,模拟该图像在距离z0处形成的同轴1-FFT的4步相移数字全息图及强度图像 % 3,利用4步相移公式获取到达CCD的物光复振幅 % 4,利用1-IFFT重建物光复振幅 % 主要变量： % h ——波长（mm）； Ih ——数字全息图； % L ——全息图宽度（mm）； z0——记录全息图的距离(mm)； %------------------------------------------------------------- clear;close all; chemin='D:\我的资料库\Pictures\'; [nom,chemin]=uigetfile([chemin,'*.*'],['调入模拟物体图像'],100,100); [XRGB,MAP]=imread([chemin,nom]); figure,imshow(XRGB); X0=rgb2gray(XRGB); figure,imshow(X0,[]); [M0,N0]=size(X0); N1=min(M0,N0); X1=zeros(N1,N1); X1(1:N1,1:N1)=X0(1:N1,1:N1); N=1024; %设定取样数 pix=0.00465 %设定CCD像素间隔 L=N*pix %CCD宽度 X1=imresize(X1,N/N1); z0=1000; % 衍射距离 mm(可以根据需要修改) h=0.532e-3; %光波长 mm(可以根据需要修改) k=2*pi/h; L0=h*z0*N/L; %S-FFT计算时初始场宽度mm Y=double(X1); bb=rand(N,N)*pi*2; X=Y.*exp(i.*bb);%叠加随机相位噪声模拟散射物体 Uh=[0:N-1]-N/2;Vh=[0:N-1]-N/2; [mh,nh]=meshgrid(Uh,Vh); T=L0/N; % 物光场像素间距 mm; figstr=strcat('物平面宽度=',num2str(L0),'mm'); figure(1), imagesc(Uh*T,Vh*T,abs(X));colormap(gray); xlabel(figstr);title('物平面图像'); U0=double(X); %---------------菲涅耳衍射的S-FFT计算起始 n=1:N; x=-L0/2+L0/N*(n-1); y=x; [yy,xx] = meshgrid(y,x); Fresnel=exp(i*k/2/z0*(xx.^2+yy.^2)); f2=U0.*Fresnel; Uf=fft2(f2,N,N);%对N*N点的离散函数f2作FFT计算 Uf=fftshift(Uf);%将FFT计算结果进行整序 x=-L/2+L/N*(n-1); y=x; [yy,xx] = meshgrid(y,x); phase=exp(i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));%菲涅耳衍射积分号前方的相位因子 Uf=Uf.*phase; Uf=Uf*T*T; %二维离散变换量值补偿 %---------------菲涅耳衍射的S-FFT计算结束 figure(2);imagesc(Uh*pix,Vh*pix,abs(Uf));colormap(gray);axis equal;axis tight;title('到达CCD的物光场振幅'); %---------------4步相移形成4幅全息图 Ar=mean(mean(abs(Uf)));%将衍射场振幅平均值设为参考光振幅 figstr1=strcat('全息图宽度=',num2str(L),'mm'); %---------------------------------0相移 Ur=Ar; Uh=Ur+Uf; Ih=Uh.*conj(Uh);%浮点型数字全息图 Imax=max(max(Ih)); I1=uint8(Ih./Imax*255);%变换为极大值255的整型数 figure(3);imshow(I1,[]);colormap(gray);xlabel(figstr1);axis equal;axis tight;title('全息图I1'); %---------------------------------pi/2相移 Ur=Ar*exp(i*pi/2); Uh=Ur+Uf; Ih=Uh.*conj(Uh);%浮点型数字全息图 Imax=max(max(Ih)); I2=uint8(Ih./Imax*255);%变换为极大值255的整型数 figure(4);imshow(I2,[]);colormap(gray);xlabel(figstr1);axis equal;axis tight;title('全息图I2'); %---------------------------------pi相移 Ur=Ar*exp(i*pi); Uh=Ur+Uf; Ih=Uh.*conj(Uh);%浮点型数字全息图 Imax=max(max(Ih)); I3=uint8(Ih./Imax*255);%变换为极大值255的整型数 figure(5);imshow(I3,[]);colormap(gray);xlabel(figstr1);axis equal;axis tight;title('全息图I3'); %---------------------------------3*pi/2相移 Ur=Ar*exp(i*3*pi/2); Uh=Ur+Uf; Ih=Uh.*conj(Uh);%浮点型数字全息图 Imax=max(max(Ih)); I4=uint8(Ih./Imax*255);%变换为极大值255的整型数 figure(6);imshow(I4,[]);colormap(gray);xlabel(figstr1);axis equal;axis tight;title('全息图I3'); %利用4幅全息图获取到达CCD的物光复振幅 I1=double(I1); I2=double(I2); I3=double(I3); I4=double(I4); U=((I4-I2)+i*(I1-I3))/4/Ar;%到达CCD的共轭物光复振幅 %1-IFFT重建 %---------------菲涅耳衍射的S-IFFT计算起始 n=1:N; x=-L/2+L/N*(n-1); y=x; [yy,xx] = meshgrid(y,x); Fresnel=exp(-i*k/2/z0*(xx.^2+yy.^2)); f2=U.*Fresnel; Uf=ifft2(f2,N,N);%对N*N点的离散函数f2作IFFT计算 x=-L0/2+L0/N*(n-1); y=x; [yy,xx] = meshgrid(y,x); phase=exp(-i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));%菲涅耳衍射积分号前的相位因子 Uf=Uf.*phase; T=L/N; %CCD取样间隔 Us=Uf*T*T; %二维离散变换量值补偿 %---------------菲涅耳衍射的S-IFFT计算结束 Uh=[0:N-1]-N/2;Vh=[0:N-1]-N/2; [mh,nh]=meshgrid(Uh,Vh); T=L0/N; % 物光场像素间距 mm; figstr=strcat('物平面宽度=',num2str(L0),'mm'); figure(6), imagesc(Uh*T,Vh*T,abs(Us));colormap(gray); xlabel(figstr);title('重建物平面图像');