CX11_2.zip_全息D-FFT_全息角谱再现_角谱_角谱再现_角谱再现算法

clear;clc;close all II=double(imread(1.bmp'));%读入数字全息图 II=II(1:500,1:700); figure(1),imshow(II,[]),title('原图'); FF=fftshift(fft2(II)); %计算全息图的频谱 figure(2),imshow(FF),title('频谱图')%显示全息图的频谱 lamda=5328*10^(-10);k=2*pi/lamda; %赋值波长、波数 [r,c]=size(II); %计算全息图的大小(像素) Lox=c*4.65*10^(-6) %赋值全息图的尺寸(纵向),单位:米 Loy=r*4.65*10^(-6) %赋值全息图的尺寸(横向),单位:米 xo=linspace(-Lox/2,Lox/2,c); %生成全息图的x坐标 yo=linspace(-Loy/2,Loy/2,r); %生成全息图的y坐标 [xo,yo]=meshgrid(xo,yo); %生成全息图的坐标网格 %% s1=zeros(1024,1024); s2=80; %设定高通滤波器的窗口大小(像素) H1=zeros(r,c); %预设高通滤波器 H1(314-s2:314+s2,333-s2:333+s2)=1;%理想高通滤波器 HFF=FF(314-s2:314+s2,333-s2:333+s2); %滤波 figure(2),imshow(HFF,[0,max(max(HFF))/20]) s1(512-s2:512+s2,512-s2:512+s2)=HFF; figure,imshow(IIFF,[0,max(max(FF))/200]) %figure(5),imshow(abs(IIFF),[0,20000]),title('滤波前频谱图');%显示全息图的频谱 %figure(6),imshow(abs(s1),[0,max(max(s1))/200]),title('滤波后频谱图'); %显示理想dai通滤波器 %break IIyp=ifft2(fftshift(s1)); %计算滤波后的全息图 %figure,imshow(abs(IIyp),[]) %下面用S-FFT算法重构再现像 zi=0.336; %全息记录面到像面的距离,单位:米 Lix=c*lamda*zi/Lox %给出像面的尺寸(x方向),单位:米 Liy=r*lamda*zi/Loy %给出像面的尺寸(y方向),单位:米 x=linspace(-Lix/2,Lix/2,c);y=linspace(-Liy/2,Liy/2,r); [x,y]=meshgrid(x,y); F0=exp(j*k*zi)/(j*lamda*zi)*exp(j*k/2/zi*(x.^2+y.^2)); F=exp(j*k/2/zi*(xo.^2+yo.^2)); % 取再现照明光垂直入射C=1 holo=Lox/c*Loy/r*fftshift(fft2(IIyp.*F.*1)); holo=holo.*F0; Ii=holo.*conj(holo); figure,imshow(Ii,[0,max(max(Ii))./2]),colormap(gray),title('S-FFT再现像') %break %% d=0.255; %赋值观察屏到衍射面的距离,单位:米 fx=linspace(-1./2./Lox,1./2./Lox,c).*c; %给出频域坐标 fy=linspace(-1./2./Loy,1./2./Loy,r).*r; [fx,fy]=meshgrid(fx,fy); H=exp(j*k*d.*sqrt(1-(lamda*fx).^2-(lamda*fy).^2));%传递函数 ft=fft2(IIyp); ft=ifft2(fftshift(ft.*H)); I=ft.*conj(ft);%subplot(1,2,1), figure,imshow(I,[0,max(max(I))./2]),,title('角谱算法再现像') kethi=linspace(-1./2./Lox,1./2./Lox,c).*c; %给出频域坐标 nenta=linspace(-1./2./Loy,1./2./Loy,r).*r; [kethi,nenta]=meshgrid(kethi,nenta); H=exp(j*k*zi.*(1-lamda.*lamda.*(kethi.*kethi+nenta.*nenta)./2)); %传递函数H %[U3,I3]=DFFT(IIyp,R,H); % alpha=pi*0.43; %参考光与x轴间的夹角 % beita=pi*0.43; % R=exp(j*k*(xo*cos(alpha)+yo*cos(beita))); %参考光 fa2=fftshift(fft2(IIyp.*1)); %衍射面上光场的傅里叶变换 Fuf2=fa2.*H; %光场的频谱与传递函数相乘 U3=ifft2(Fuf2); %在观察屏上的光场分布 I3=U3.*conj(U3); %在观察屏上的光强分布 phyp2=angle(U3); figure,imshow(I3,[0,max(max(I3))./2.5]),title('角谱再现像'); %break figure(11),imshow(phyp2,[]),title('角谱再现像相位图');