我已经发过相关的帖子，是我以前做过的课程论文。

近来看见有很多人回帖说需要程序，故而总结一下方法和共享程序：通过MATLAB软件编程实现夫琅和费衍射的方法：(1)用衍射积分(2)傅立叶变换一、衍射积分

相关程序如下：

1.单缝衍射clear

lamba=500e-9;%波长

a=1e-3;D=1;

ym=3*lamba*D/a;%屏幕上y的范围

n=51;%屏幕上的点数

ys=linspace(-ym,ym,n);

n=51;%屏幕上的点数

yp=linspace(0,a,n);

fori=1:n

sinphi=ys(i)/D;

alpha=pi*yp*sinphi/lamba;

sumcos=sum(cos(alpha));

sumsin=sum(sin(alpha));

B(i,:)=(sumcos^2+sumsin^2)/n^2;

end

N=256;%确定灰度的等级

Br=(B/max(B))*N;

subplot(1,2,1)

image(ym,ys,Br);

colormap(gray(N));%色调处理

subplot(1,2,2)

plot(B,ys,'k');

2.多缝衍射

clear

lamda=500e-9; %波长

N=2; %缝数，可以随意更改变换

a=2e-4;D=5;d=5*a;

ym=2*lamda*D/a;xs=ym;

n=1001;

ys=linspace(-ym,ym,n);

for i=1:n

sinphi=ys(i)/D;

alpha=pi*a*sinphi/lamda;

beta=pi*d*sinphi/lamda;

B(i,:)=(sin(alpha)./alpha).^2.*(sin(N*beta)./sin(beta)).^2;

B1=B/max(B);

end

NC=256; %确定灰度的等级

Br=(B/max(B))*NC;

subplot(1,2,1)

image(xs,ys,Br);

colormap(gray(NC)); %色调处理

subplot(1,2,2)

plot(B1,ys,'k');

3.矩孔衍射

clear

lamda=500e-9;

a=1e-3;

b=1e-3;

f=1;

m=500;

ym=8000*lamda*f;

ys=linspace(-ym,ym,m);

xs=ys;

n=255;

for i=1:m

sinth1=xs(i)/sqrt(xs(i)^2+f^2);

sinth2=ys./sqrt(ys.^2+f^2);

angleA=pi*a*sinth1/lamda;

angleB=pi*b*sinth2./lamda;

B(:,i)=(sin(angleA).^2.*sin(angleB).^2.*5000./(angleA.^2.*angleB.^2));

end

subplot(1,2,1)

image(xs,ys,B)

colormap(gray(n))

subplot(1,2,2)

plot(B(m/2,:),ys,'k')

4.正弦光栅

clear all

xm=10*pi;ys=xm;

xs=linspace(-xm,xm,500);

B=cos(xs)+1;

N=255;

Br=B/2*N;

image(xs,ys,Br);

colormap(gray(N));

二、傅里叶变换

(1)基本思想：

在傅立叶变换光学中夫琅和费衍射场的强度分布就等于屏函数的功率谱。因此我们可以直接将光屏进行傅立叶变换，再处理得到衍射图样。

(2)具体操作：首先，物体图像的生成:(1)直接由矩阵运算生成(2)也可利用Windows下的画图工具，生成一幅黑白图像。其次，并调用命令函数imread()输入图像，输入的图像是一个巨大的二维矩阵，利用MATLAB函数库中的fft2()命令对该矩阵进行二位离散傅立叶变换，得到图像的频谱，该频谱是一个复数矩阵，然后用取模函数abs()对该复数矩阵取模，得到振幅谱矩阵，利用函数fftshift()对取模后的矩阵进行频谱位移，这是因为变换后的二维矩阵的直流分量位于图像的周边角，该函数交换矩阵的1、3象限和2、4象限，使直流分量移到频谱中心，从而使FFT频谱可视效果与实际图像相吻合。最后利用imshow()函数将图像显示出来。

***具体图像参考：https://www.ilovematlab.cn/thread-25982-1-1.html

基本的程序：

clear

a=imread('***');%你的图像的位置或直接编写孔生成

grid on

figure(1)

imshow(a,[])

afft=fft2(a);

aabs=abs(afft);

aabss=fftshift(aabs);

figure(2)

imshow(aabss,[])

colormap(gray);

colorbar

figure(3)

plot(aabss)

colormap(gray);

figure(4)

meshc(aabss)

maxx1=max(max(aabss));

如果还有什么不懂的直接发到我邮箱！！！