数字全息matlab模拟程序

clear all; close all;clc lamda=532e-6; delta_x=8e-3; k=2*pi/lamda; delta_y=8e-3; H=imread('1.bmp'); H=H(:,640-511:640+512); I1=double(imread('1.bmp')); I1=I1(:,640-511:640+512); % I1=I1(:,:,2); % I1=I1/255; % figure(2);imshow(I1); I2=double(imread('2.bmp')); I2=I2(:,640-511:640+512); % I2=I2(:,:,2); % I2=I2/255; I3=double(imread('3.bmp')); I3=I3(:,640-511:640+512); % I3=I3(:,:,2); % I3=I3/255; % I4=double(imread('4.bmp')); % I4=I4(:,:,2); % I4=I4/255; % H=(I1-I3)+1i*(I2-I4); % H=H(:,:,2); a1=0; a2=0; sita1=0;sita2=2/3*pi+a1;sita3=4/3*pi+a2; % H=I1*(exp(1i*sita3)-exp(1i*sita2))+I2*(exp(1i*sita1)-exp(1i*sita3))+I3*(exp(1i*sita2)-exp(1i*sita1)); H=double(H); % H=double(H(:,640-511:640+512))/4; [N,M]=size(H); % figure,imshow(H,[]),title('截取的全息图'); % break strum = New_FFT2(H,delta_x,delta_x,0,0 ); figure(),imshow(abs(strum)./20),title('频谱'); % % figure(1); % % plot(i,r); % i=i+1; % end % end % imwrite(abs(strum)./20,'频谱图4.bmp'); % break strum1=strum(612-50:612+50,594-50:594+50); % strum1=strum(553:701,522:656); figure,imshow(abs(strum1)./0.01,[]),title('截取的频谱'); [F,G]=size(strum1); strum2=zeros(N,N); strum2(N/2-F/2+1:N/2+F/2,N/2-G/2+1:N/2+G/2)=strum1; figure,imshow(abs(strum2)./0.10,[]),title('补0后的频谱'); % return delta_fx=1/N/delta_x; H2=New_IFFT2(strum2,delta_fx,delta_fx,0,0); figure,imshow(abs(H2),[]),title('物光波在记录面上的复振幅'); d=390;zr=d;dx=delta_x; Lox=1024*delta_x; Loy=1024*delta_y; xo=linspace(-Lox/2,Lox/2,1024);yo=linspace(-Loy/2,Loy/2,1024); [xo,yo]=meshgrid(xo,yo); alpha=pi/2.000096; %参考光与x轴间的夹角 beita=pi/2.0001; %参考光与y轴间的夹角 R=exp(j*k*(xo*cos(alpha)+yo*cos(beita))); %参考光 [u,dx2,dy2,x2,y2] = fresnel(H2.*conj(R),M,N,dx,dx,zr,lamda); % u=u(707:841,262:397); phase1=(angle(u)); figure,imshow(phase1,[]),title('相位像'); % imwrite(phase1,'相位再现像p1.bmp'); u=abs(u).^2; u=u/max(max(u))./0.6; % figure();imshow(u); % imwrite(abs(u)./1,'强度再现像z1.bmp'); break I=Single_FFT(H2,d,lamda,delta_x,delta_y,0,0 ); I=I(530:806,380:700); phase2=(angle(I)); figure,mesh(phase2),title('相位像'); figure(d),imshow(abs(I)./0.00009),title('强度'); % imwrite(abs(I)./0.00009,'强度再现像2.bmp'); break; i=1; for sita2=1.1:0.01:1.4 for sita3=3.6:0.01:3.9 Hx=I1*(exp(1i*sita3)-exp(1i*sita2))+I2*(exp(1i*sita1)-exp(1i*sita3))+I3*(exp(1i*sita2)-exp(1i*sita1)); [f2,dx2,dy2,x2,y2] = fresnel(Hx,M,M,dx,dx,zr,lamda); f2=abs(f2); f2=f2/max(max(f2)); w=f2(186:220,698:738); b=f2(398:442,714:442); % w=f2(726:768,867:897); % b=f2(608:658,738:778); wmax=max(max(w)); wmin=min(min(w)); bmax=max(max(b)); bmin=min(min(b)); rw=(wmax-wmin)/(wmax+wmin); rb=(bmax-bmin)/(bmax+bmin); r=(rw+rb)/2; a(i,1)=sita2; a(i,2)=sita3; a(i,3)=r; % % r=rw; % % figure(1); % % plot(i,r); i=i+1; % % hold on; end end m=min(min(a)); [k1,k2]=find(a==m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % sita2=a(k1,1);sita3=a(k1,2); Hx=I1*(exp(1i*sita3)-exp(1i*sita2))+I2*(exp(1i*sita1)-exp(1i*sita3))+I3*(exp(1i*sita2)-exp(1i*sita1)); [f2,dx2,dy2,x2,y2] = fresnel(Hx,M,M,dx,dx,zr,lamda); f2=abs(f2); f2=f2/max(max(f2)); % w=f2(408:424,297:306); % b=f2(395:409,238:251); % wmax=max(max(w)); % wmin=min(min(w)); % bmax=max(max(b)); % bmin=min(min(b)); % rw=(wmax-wmin)/(wmax+wmin); % rb=(bmax-bmin)/(bmax+bmin); % r=(rw+rb)/2; % imwrite(f2,'USAF_after.jpg'); figure;imshow((f2)); break strumnew = New_FFT2(Hx,delta_x,delta_x,0,0 ); figure,imshow(abs(strumnew)./20),title('频谱') % imwrite(abs(strum)./20,'频谱图4.bmp'); % break % strum1=strum(512:1024,512:1024); % strum11=strumnew(553:701,522:656); figure,imshow(abs(strum11)./0.01,[]),title('截取的频谱'); [F,G]=size(strum11); strum2=zeros(N,N); strum2(N/2-F/2+1:N/2+F/2,N/2-G/2+1:N/2+G/2)=strum11; figure,imshow(abs(strum2)./0.10,[]),title('补0后的频谱'); % return delta_fx=1/N/delta_x; H3=New_IFFT2(strum2,delta_fx,delta_fx,0,0); figure,imshow(abs(H3),[]),title('物光波在记录面上的复振幅'); d=390;zr=d;dx=delta_x; [u,dx2,dy2,x2,y2] = fresnel(H3,M,N,dx,dx,zr,lamda); % u=u(707:841,262:397); phase3=(angle(u)); figure,mesh(phase3),title('相位像') u=abs(u).^2; u=u/max(max(u))./0.6; figure;imshow(u);title('校正再现像'); return for ii=250:1:290 d=1*ii; % d=254; I=Precise_Angular_Spectrum(H2,lamda,d,delta_x,delta_x,0,0 ); figure(d),imshow(abs(I),[]),title('强度') end break phase1=angle(H2); figure(6),imshow(phase1,[]),title('相位像') %I2 = Single_FFT( H2,d,lamda,delta_x,delta_x,0,0 ); %figure(7),imshow(abs(I2),[]),title('强度'); %break % clear C ph_unwrap1=unwrap_LS(phase1); figure(7),imshow(ph_unwrap1,[]) figure(8),mesh(ph_unwrap1) [t,t]=size(ph_unwrap1); x=1:t;y=1:t; [x,y]=meshgrid(x,y); B=[ones(1,t);1:t;(1:t).^2]'; G=[ones(1,t);1:t;(1:t).^2]'; C=inv(B'*B)*B'*ph_unwrap1*G*inv(G'*G); % W=C(1,1)+C(2,1)*y+C(3,1)*y.^2+C()*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2; % W=inv(B*B')*B*(B'*B)*C(G'*G)*G'*(G*G'); % W=C(1,1)+C(1,2)*y+C(1,3)*y.^2+C(2,1)*x+C(2,2)*x.*y+C(2,3)*x.*y.^2+C(3,1)*x.^2+C(3,2)*x.^2.*y+C(3,3)*x.^2.*y.^2; % W=C(1)+C(2)*y+C(3)*y.^2+C(4)*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2; % W=C(1)+C(2)*y+C(3)*y.^2+C(4)*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2; W=C(1)+C(2)*y+C(3)*y.^2+C(4)*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2;%+C(4,1)*x.^3+C(4,2).*x.^3.*y+C(4,3).*x.^3.*y.^2+C(4,4).*x^3.*y.^3+C(1,4).*y.^3+C(2,4).*x.*y.^3+C(3,4).*x.^2.*y.^3; ph_unwrap=ph_unwrap1-W; ph_unwrap=ph_unwrap./10; %figure,imshow(W,[]); %figure,imshow(ph_unwrap,[]); figure(9),imshow(ph_unwrap,[]) % figure,mesh(ph_unwrap) return ph_unwrap2=ph_unwrap(106:501,486:881); figure,imshow(ph_unwrap2,[]) yy=ph_unwrap2(230,104:240); % t=size(yy); phase_c2=zeros(396,396); B=[ones(size(yy)) ;yy;(yy).^2]; AA=B/yy; xx=ph_unwrap2(230-68:230+68,172); B=[ones(size(xx));xx;(xx).^2]; BB=B/xx; m=230-68:230+68; n=104:240; for ii=230-68:230+68; for jj=104:240; phase_c2(ii,jj)=AA(1).*ii+BB(1).*(jj)+AA(2)^2*(ii)^2+BB(2)^2*(jj)^2; end end % phase_c2=AA(1).*m+BB(1).*n+AA(2)^2.*m.^2+BB(2)^2.*n.^2; figure imshow(phase_c2,[]) % [t,t]=size(ph_unwrap2); % x=1:t;y=1:t; % [x,y]=meshgrid(x,y); % B=[ones(1,t);1:t;(1:t).^2]'; % G=[ones(1,t);1:t;(1:t).^2]'; % C=inv(B'*B)*B'*ph_unwrap2*G*inv(G'*G); % W=C(1)+C(2)*y+C(3)*y.^2+C(4)*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2; %划线法去畸变 ph_unwrap3=ph_unwrap2-phase_c2; figure,imshow(ph_unwrap3,[]) % ph_unwrap=ph_unwrap./10; return % H2=fftshift(ifft2(strum2)); % delta_fx=1/N/delta_x; % H2=New_IFFT2(strum2,delta_fx,delta_fx,0,0); % % H2=mat2gray(abs(H2)); figure(4),imshow(abs(H2),[]),title('new hologram'); phase1=angle(H2); figure(5),imshow(phase1,[]),title('相位像') ph_unwrap1=unwrap_LS(phase1); figure,imshow(ph_unwrap1,[]) [t,t]=size(ph_unwrap1); x=1:t;y=1:t; [x,y]=meshgrid(x,y); B=[ones(1,t);1:t;(1:t).^2]'; G=[ones(1,t);1:t;(1:t).^2]'; C=inv(B'*B)*B'*ph_unwrap1*G*inv(G'*G); W=C(1)+C(2)*y+C(3)*y.^2+C(4)*x+C(5)*x.*y+C(6)*x.*y.^2+C(7)*x.^2+C(8)*x.^2.*y+C(9)*x.^2.*y.^2; ph_unwrap=ph_unwrap1-W; ph_unwrap=ph_unwrap./10; %figure,imshow(W,[]); %figure,imshow(ph_unwrap,[]); figure,imshow(-ph_unwrap,[]) figure,mesh(-ph_unwrap) return return %两步相减法 W=imread('E:\实验数据\5.9\20x_800.bmp'); W=double(W(1:1024,129:1152)); figure(12),imshow(W,[]),title('截取的全息图1') strum = New_FFT2(H,delta_x,delta_x,0,0 ); figure(13),imshow(abs(strum))