移相干涉技术1-多种干涉条纹仿真模拟生成（原理转载+代码实现 包括模拟生成干涉条纹图）

过去的干涉测量技术是通过人的肉眼或者相机拍摄，来直观判断干涉图中条纹特征进而完成测量，该方法的不稳定因素（比如人的主观意志）很多，其精度误差在/10左右38]；现代干涉测量技术通过将电子技术、计算机技术、光电图像处理技术、算法模型和高精度光学元件结合，大大推动了干涉测量领域的发展，使其进入高精度时代。移相干涉技术(P$)作为现代干涉测量技术的代表，已经成为干涉测量领域学者的研究重点。移相干涉技术测量的原理：移相式激光干涉仪通过控制压电陶瓷驱动器(PZT)移动参考镜位置或者改变激光器波长的方式，令参考光和测试光的光程差改变，从干涉图像中可看出干涉条纹进行了相应的位移。在移相过程中，通过光电探测器(CCD或者CMOS)采集不同移相量下的干涉图，图2-1所示为一组移相间隔为90°的移相干涉图组。计算机根据特定的数学算法和模型对干涉图组的像素数据进行系列运算，进而可以求得被测元件的面形信息，并根据计算好的相关参数进行质量评价[39（本文移相干涉系统通过驱动PZT模块完成移相，因此不对基于改变激光器波长以完成移相原理的波长调谐移相技术作介绍)

基于PZT实现移相的过程为：当PZT控制模块发出指定命令时，PZT伸长量进行相应改变，即参考镜和被测镜的相对距离改变，进而完成相位调制以达到移相目的。在上述移相干涉测量过程中，干涉图像中的条纹分布是由参考光束和测试光束叠加形成的干涉场导致，干涉场的光强信息与干涉图中像素点的灰度信息相关，其分布函数可表示为：

其中：A(X,y)是干涉图的背景光强，B(x,y)是调制度，p(X,y)是待测物的波面初始相位信息，6：是移相量，(x,y)是干涉图中像素点的坐标信息。参考镜和被测镜间的面形差P(X,y)与位相分布的关系可表示为： 

 λ为激光器波长（本文系统为He-Ne激光器一一632.8nm),在移相干涉测量技术中，测量单位为入。当通过解相算法处理干涉图组的像素数据I:(X,y),可计算出p(x,y),并进一步得到面形信息P(x,y)以完成检测工作。下面简单介绍三种本文所设计系统应用到的解相算法（两种传统解相算法+一种新颖的平均算法)：

下面是模拟生成干涉条纹图的代码

%zernike
clear;clc;close all;

[x,y]=meshgrid(-1:0.01:1);

k=36;%多项式项数

Z=cell(1,36);

Z{1,1}=x;

Z{1,2}=y;

Z{1,3}=-1+2*(x.^2+y.^2);

Z{1,4}=x.^2-y.^2;

Z{1,5}=2*x.*y;

Z{1,6}=-2*y+3*x.*(x.^2+y.^2);

Z{1,7}=-2*y+3*y.*(x.^2+y.^2);

Z{1,8}=1-6*(x.^2+y.^2)+6*(x.^2+y.^2).^2;

Z{1,9}=x.^3-3*x.*y.^2;

Z{1,10}=3*x.^2.*y-y.^3;

Z{1,11}=-3*x.^2+3*y.^2+4*x.^4-4*y.^4;

Z{1,12}=-6*x.*y+8*x.*y*(x.^2+y.^2);

Z{1,13}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*x;

Z{1,14}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*y;

Z{1,15}=20*(x.^2+y.^2).^3-30*(x.^2+y.^2).^2+12*(x.^2+y.^2)-1;

Z{1,16}=x.^4-6*x.^2.*y.^2+y.^4;

Z{1,17}=4*x.^3.*y-4*x.*y.^3;

Z{1,18}=15*(x.^2+y.^2).^2.*x+20*(x.^2+y.^2).*x.^3+12*(x.^2+y.^2).*x-16*x.^3;

Z{1,19}=(5*(x.^2+y.^2)-4).*(3*x.^2.*y-y.^3);

Z{1,20}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(x.^2-y.^2);

Z{1,21}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(2*x.*y);

Z{1,22}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*x;

Z{1,23}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*y;

Z{1,24}=70*(x.^2+y.^2).^4-140*(x.^2+y.^2).^3+90*(x.^2+y.^2).^2-20*(x.^2+y.^2)+1;

Z{1,25}=5*x*(x.^2+y.^2).^2-20*x.^3.*(x.^2+y.^2)+16*x.^5;

Z{1,26}=-6*(x.^2+y.^2).*x.^2.*y+14*x.^4.*y-3*x.^2.*y.*(x.^2+y.^2)-2*x.^2.*y.^3+(x.^2+y.^2).*y.^3;

Z{1,27}=(6*(x.^2+y.^2)-5).*((x.^2+y.^2).^2-8*(x.^2+y.^2).*x.^2+8*x.^4);

Z{1,28}=(6*(x.^2+y.^2)-5).*(4*x.^3.*y-4*x.*y.^3);

Z{1,29}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(-3*x*(x.^2+y.^2)+4*x.^3);

Z{1,30}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(3*x.^2.*y-y.^3);

Z{1,31}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(x.^2-y.^2);

Z{1,32}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(2*x.*y);

Z{1,33}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*x;

Z{1,34}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*y;

Z{1,35}=252*(x.^2+y.^2).^5-630*(x.^2+y.^2).^4+560*(x.^2+y.^2).^3-210*(x.^2+y.^2).^2+30*(x.^2+y.^2)-1;

Z{1,36}=924*(x.^2+y.^2).^6-2772*(x.^2+y.^2).^5+3150*(x.^2+y.^2).^4-1680*(x.^2+y.^2).^3+420*(x.^2+y.^2).^2-42*(x.^2+y.^2)+1;

%单位圆内

mask=double((x.^2+y.^2)<=1);

for i=1:k

    Z{1,i}=Z{1,i}.*mask;

end

%显示泽尔尼克系数

for i=1:k

    subplot(6,6,i),imshow(Z{1,i},[]);

end

%a+(b-a)rand(m,n) ,[a,b]区间m，n维随机数生成

%计算相位P(x,y)，ai根据论文取随机数

% 输出

num=100;%输出数量

for n=1:num

    for i=1:k %ai是泽尔尼克基函数的系数，根据论文给出的范围来定

      if i==1 %a1=0,P0=0 

        P=0;

      end

      if i==2||i==3 %ai 0,5

        P=5*rand(1)*Z{1,i}+P;

      end

      if i==4 %ai -5,5

        P=-5+10*rand(1)*Z{1,i}+P;

      end

      if i>=5&&i<=9 %ai -0.3,0.3

        P=-0.2+0.6*rand(1)*Z{1,i}+P;

      end

      if i>=10&&i<=36 %ai -0.03,0.03

        P=-0.03+0.006*rand(1)*Z{1,i}+P;

      end

    end

  P=pi*P.*mask;

  w=0.2+0.8*rand(1);%w是0.2,1上的随机数

  A0min=0.5/(1+w+2*sqrt(w));%A0范围下限

  A0max=0.95/(1+w+2*sqrt(w));%A0范围上限

  A0=A0min+(A0max-A0min)*rand(1);%生成振幅强度A0

  x0=-0.2+0.4*rand(201,201);

  y0=-0.2+0.4*rand(201,201);%x0,y0表示光点在x和y方向上的偏移，这两个方向都是在[0.2,0.2]范围内随机生成的

  derta=0.9+0.2*rand(1);%光束宽度0.9,1.1

  %激光源参考光强度

  Ir=A0*exp(-(((x-x0).^2+(y-y0).^2)/derta^2));%越靠中间(x-x0).^2+(y-y0).^2)越小，越亮

  % c=cos(P);%cos相位

  A=(w+1)*Ir;%背景光强

  B=2*sqrt(w)*Ir;%调制光强

  I=A+B.*cos(P);%条纹图

  I=I.*mask;

  name=num2str(n);

  I=imresize(I,[256,256]);

  imwrite(I,['E:\my_self_progect\matlab\fringes\',name,'.jpg']);

end

figure(5),imshow(A,[]);

 figure(6),imshow(B,[]);

figure(7),imshow(A+B.*cos(P),[]);

这段Matlab代码是一个生成干涉图像的脚本，其中使用了泽尼克多项式（Zernike polynomials）来模拟光波的相位分布。泽尼克多项式是一组定义在单位圆盘上的正交多项式，它们在光学领域被用于表示波前的相位误差。这段代码包含以下几个主要部分：

1. 泽尼克多项式的定义：

   • 脚本首先创建了一个坐标网格(x, y)，这个网格覆盖了[-1, 1]的范围，步长为0.01。
   • 接下来，定义了36个泽尼克多项式，每个多项式都是x和y的函数，表示不同的波前误差模式。

2. 单位圆内的应用：

   • 使用逻辑索引创建了一个掩膜mask，它定义了单位圆内的区域。这意味着多项式只在这个圆内部被考虑。
   • 所有的泽尼克多项式都与这个掩膜相乘，从而确保它们只在单位圆内非零。

3. 相位计算：

   • 代码计算了一个相位图P，这是通过为每个泽尼克多项式随机分配系数并将它们相加来实现的。
   • 这些系数ai从不同的范围内随机选取，这些范围在注释中有说明，例如a1=0，ai对于第2和第3项是0到5之间的随机数，等等。

4. 干涉图像生成：

   • 生成了一个强度分布Ir，模拟了一个高斯分布的激光光束。
   • 振幅A模拟了背景光强度，而B模拟了调制光强度。
   • 最后，生成的干涉图像I由背景光强度和调制光强度的余弦调制组成，余弦项包含了上面计算的相位图P。

5. 图像保存：

   • 脚本将生成的干涉图像保存到指定的文件夹中。保存的图像数量由num变量确定。

6. 显示图像：

   • 最后，脚本显示了背景光强度A，调制光强度B，以及最终的干涉图像I。

% 清除环境并关闭所有图像
clear;
clc;
close all;

% 定义x和y的网格
[x, y] = meshgrid(linspace(-1, 1, 256), linspace(-1, 1, 256));
mask = double(x.^2 + y.^2 <= 1); % 单位圆的掩膜

% 定义泽尼克多项式的系数
% 这里为了简化，我们只使用前9个泽尼克多项式
coeffs = rand(1, 9) - 0.5; % 随机生成系数，范围在[-0.5, 0.5]

% 计算泽尼克多项式的线性组合
P = zeros(size(x));
for i = 1:length(coeffs)
    Z = zernike_polynomial(i, x, y); % 计算第i个泽尼克多项式
    P = P + coeffs(i) * Z; % 将泽尼克多项式加权后累加
end
P = P .* mask; % 应用掩膜

% 设置保存图像的文件夹
outputFolder = 'E:\gansheyi\matlab\output\folder'; % 更改为您的输出文件夹路径
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder); % 如果输出文件夹不存在，则创建它
end

% 生成和保存9张相移干涉图像
for n = 1:9
    phase_shift = (n - 1) * pi / 2; % 相位移动
    I = 0.5 + 0.5 * cos(P + phase_shift); % 生成干涉图像
    I = I .* mask; % 再次应用掩膜

    % 保存图像
    filename = fullfile(outputFolder, sprintf('interference_%d.bmp', n));
    imwrite(I, filename);

    % 显示图像
    subplot(3, 3, n);
    imshow(I, []);
    title(sprintf('Interference Pattern %d', n));
end

function Z = zernike_polynomial(n, x, y)
    % 这个函数根据索引n返回对应的泽尼克多项式
    r = sqrt(x.^2 + y.^2);
    theta = atan2(y, x);
    % 添加更多的泽尼克多项式，如果需要的话
    switch n
        case 1
            Z = 1;
        case 2
            Z = 2 * r .* cos(theta);
        case 3
            Z = 2 * r .* sin(theta);
        % 添加更多的泽尼克多项式，如果需要的话
        % ...
        otherwise
            Z = 0;
    end
end

User

 

 增加条纹数

% 初始化环境
clear;
clc;
close all;

% 定义x和y的网格
[x, y] = meshgrid(-1:0.01:1);

% 定义泽尼克多项式

k=36;%多项式项数

Z=cell(1,36);

Z{1,1}=x;

Z{1,2}=y;

Z{1,3}=-1+2*(x.^2+y.^2);

Z{1,4}=x.^2-y.^2;

Z{1,5}=2*x.*y;

Z{1,6}=-2*y+3*x.*(x.^2+y.^2);

Z{1,7}=-2*y+3*y.*(x.^2+y.^2);

Z{1,8}=1-6*(x.^2+y.^2)+6*(x.^2+y.^2).^2;

Z{1,9}=x.^3-3*x.*y.^2;

Z{1,10}=3*x.^2.*y-y.^3;

Z{1,11}=-3*x.^2+3*y.^2+4*x.^4-4*y.^4;

Z{1,12}=-6*x.*y+8*x.*y*(x.^2+y.^2);

Z{1,13}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*x;

Z{1,14}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*y;

Z{1,15}=20*(x.^2+y.^2).^3-30*(x.^2+y.^2).^2+12*(x.^2+y.^2)-1;

Z{1,16}=x.^4-6*x.^2.*y.^2+y.^4;

Z{1,17}=4*x.^3.*y-4*x.*y.^3;

Z{1,18}=15*(x.^2+y.^2).^2.*x+20*(x.^2+y.^2).*x.^3+12*(x.^2+y.^2).*x-16*x.^3;

Z{1,19}=(5*(x.^2+y.^2)-4).*(3*x.^2.*y-y.^3);

Z{1,20}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(x.^2-y.^2);

Z{1,21}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(2*x.*y);

Z{1,22}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*x;

Z{1,23}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*y;

Z{1,24}=70*(x.^2+y.^2).^4-140*(x.^2+y.^2).^3+90*(x.^2+y.^2).^2-20*(x.^2+y.^2)+1;

Z{1,25}=5*x*(x.^2+y.^2).^2-20*x.^3.*(x.^2+y.^2)+16*x.^5;

Z{1,26}=-6*(x.^2+y.^2).*x.^2.*y+14*x.^4.*y-3*x.^2.*y.*(x.^2+y.^2)-2*x.^2.*y.^3+(x.^2+y.^2).*y.^3;

Z{1,27}=(6*(x.^2+y.^2)-5).*((x.^2+y.^2).^2-8*(x.^2+y.^2).*x.^2+8*x.^4);

Z{1,28}=(6*(x.^2+y.^2)-5).*(4*x.^3.*y-4*x.*y.^3);

Z{1,29}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(-3*x*(x.^2+y.^2)+4*x.^3);

Z{1,30}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(3*x.^2.*y-y.^3);

Z{1,31}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(x.^2-y.^2);

Z{1,32}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(2*x.*y);

Z{1,33}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*x;

Z{1,34}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*y;

Z{1,35}=252*(x.^2+y.^2).^5-630*(x.^2+y.^2).^4+560*(x.^2+y.^2).^3-210*(x.^2+y.^2).^2+30*(x.^2+y.^2)-1;

Z{1,36}=924*(x.^2+y.^2).^6-2772*(x.^2+y.^2).^5+3150*(x.^2+y.^2).^4-1680*(x.^2+y.^2).^3+420*(x.^2+y.^2).^2-42*(x.^2+y.^2)+1;
% 应用单位圆掩膜
mask = double((x.^2 + y.^2) <= 1);

% 选择使用一个径向对称的泽尼克多项式生成干涉图
% 这里我们使用Z{1,4} = x^2 - y^2，并增加其系数以增加条纹数
coeff_Z4 = 50; % 调整这个系数可以增加条纹数
P = coeff_Z4 * Z{1,2}; % 使用增大系数的泽尼克多项式
P = P .* mask; % 应用掩膜

% 设置保存图像的文件夹
outputFolder = 'E:\gansheyi_daima\matlab\path_to_your_output_folder2'; % 更改为您的输出文件夹路径
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder); % 如果输出文件夹不存在，则创建它
end

% 生成和保存9张相移干涉图像
A = 0.5;  % 设置振幅，增加对比度
for n = 1:9
    % 相位移动
    phase_shift = (n - 1) * pi / 2;
    
    % 生成干涉图像，增加对比度
    I = 0.5 + A * cos(P + phase_shift) / 2;
    I = I .* mask; % 再次应用掩膜
    
    % 保存图像
    filename = fullfile(outputFolder, sprintf('interference_%d.bmp', n));
    imwrite(I, filename);
end

 

 生成斜的条纹

直接生成均匀斜向条纹，这通常不需要复杂的泽尼克多项式，因为它们用于模拟更复杂的波前畸变。然而，如果想要使用泽尼克多项式来调整或添加特定的波前畸变，您可以选择合适的泽尼克多项式并将其加入到相位P中。

% 初始化环境
clear;
clc;
close all;

% 定义x和y的网格
[x, y] = meshgrid(-1:0.01:1);

% 应用单位圆掩膜
mask = double((x.^2 + y.^2) <= 1);

% 创建斜着的条纹相位变化
P = x + y;
P = P .* mask; % 应用掩膜

% 设置保存图像的文件夹
outputFolder = 'E:\gansheyi_daima\matlab\path_to_your_output_folder2'; % 更改为您的输出文件夹路径
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder); % 如果输出文件夹不存在，则创建它
end

% 生成和保存9张相移干涉图像
A = 0.5;  % 设置振幅，增加对比度
B = 30;   % 系数，用于调整条纹的间距和数量
for n = 1:9
    % 相位移动
    phase_shift = (n - 1) * pi / 2;
    
    % 生成干涉图像，增加对比度
    I = 0.5 + A * cos(B * P + phase_shift) / 2;
    I = I .* mask; % 再次应用掩膜
    
    % 保存图像
    filename = fullfile(outputFolder, sprintf('interference_%d.bmp', n));
    imwrite(I, filename);
end

 

上面舍去了泽妮克多项式的方法，我们下面加上泽妮克系数多项式来

% 初始化环境
clear;
clc;
close all;

% 定义x和y的网格
[x, y] = meshgrid(-1:0.01:1);

% 定义泽尼克多项式

k=36;%多项式项数

Z=cell(1,36);

Z{1,1}=x;

Z{1,2}=y;

Z{1,3}=-1+2*(x.^2+y.^2);

Z{1,4}=x.^2-y.^2;

Z{1,5}=2*x.*y;

Z{1,6}=-2*y+3*x.*(x.^2+y.^2);

Z{1,7}=-2*y+3*y.*(x.^2+y.^2);

Z{1,8}=1-6*(x.^2+y.^2)+6*(x.^2+y.^2).^2;

Z{1,9}=x.^3-3*x.*y.^2;

Z{1,10}=3*x.^2.*y-y.^3;

Z{1,11}=-3*x.^2+3*y.^2+4*x.^4-4*y.^4;

Z{1,12}=-6*x.*y+8*x.*y*(x.^2+y.^2);

Z{1,13}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*x;

Z{1,14}=(10*(x.^2+y.^2).^2-12*(x.^2+y.^2)+3).*y;

Z{1,15}=20*(x.^2+y.^2).^3-30*(x.^2+y.^2).^2+12*(x.^2+y.^2)-1;

Z{1,16}=x.^4-6*x.^2.*y.^2+y.^4;

Z{1,17}=4*x.^3.*y-4*x.*y.^3;

Z{1,18}=15*(x.^2+y.^2).^2.*x+20*(x.^2+y.^2).*x.^3+12*(x.^2+y.^2).*x-16*x.^3;

Z{1,19}=(5*(x.^2+y.^2)-4).*(3*x.^2.*y-y.^3);

Z{1,20}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(x.^2-y.^2);

Z{1,21}=(15*(x.^2+y.^2).^2-20*(x.^2+y.^2)+6).*(2*x.*y);

Z{1,22}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*x;

Z{1,23}=(35*(x.^2+y.^2).^3-60*(x.^2+y.^2).^2+30*(x.^2+y.^2)-4).*y;

Z{1,24}=70*(x.^2+y.^2).^4-140*(x.^2+y.^2).^3+90*(x.^2+y.^2).^2-20*(x.^2+y.^2)+1;

Z{1,25}=5*x*(x.^2+y.^2).^2-20*x.^3.*(x.^2+y.^2)+16*x.^5;

Z{1,26}=-6*(x.^2+y.^2).*x.^2.*y+14*x.^4.*y-3*x.^2.*y.*(x.^2+y.^2)-2*x.^2.*y.^3+(x.^2+y.^2).*y.^3;

Z{1,27}=(6*(x.^2+y.^2)-5).*((x.^2+y.^2).^2-8*(x.^2+y.^2).*x.^2+8*x.^4);

Z{1,28}=(6*(x.^2+y.^2)-5).*(4*x.^3.*y-4*x.*y.^3);

Z{1,29}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(-3*x*(x.^2+y.^2)+4*x.^3);

Z{1,30}=(21*(x.^2+y.^2).^2-30.*(x.^2+y.^2)+10).*(3*x.^2.*y-y.^3);

Z{1,31}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(x.^2-y.^2);

Z{1,32}=(56*(x.^2+y.^2).^3-105*(x.^2+y.^2).^2+60*(x.^2+y.^2)-10).*(2*x.*y);

Z{1,33}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*x;

Z{1,34}=(126*(x.^2+y.^2).^4-280*(x.^2+y.^2).^3+210*(x.^2+y.^2).^2-60*(x.^2+y.^2)+5).*y;

Z{1,35}=252*(x.^2+y.^2).^5-630*(x.^2+y.^2).^4+560*(x.^2+y.^2).^3-210*(x.^2+y.^2).^2+30*(x.^2+y.^2)-1;

Z{1,36}=924*(x.^2+y.^2).^6-2772*(x.^2+y.^2).^5+3150*(x.^2+y.^2).^4-1680*(x.^2+y.^2).^3+420*(x.^2+y.^2).^2-42*(x.^2+y.^2)+1;
% 应用单位圆掩膜
mask = double((x.^2 + y.^2) <= 1);

% 选择使用一些泽尼克多项式来调整条纹
Zernike_indices = [2, 3]; % 选择要使用的泽尼克多项式的索引
Zernike_coefficients = [0.5, 0.1]; % 设置选定泽尼克多项式的系数

% 创建斜着的条纹相位变化并添加泽尼克畸变
B = 50; % 系数，用于调整条纹的间距和数量
P = B * (x + y); % 基础斜向条纹
for i = 1:length(Zernike_indices)
    P = P + Zernike_coefficients(i) * Z{Zernike_indices(i)};
end
P = P .* mask; % 应用掩膜

% 设置保存图像的文件夹
outputFolder = 'E:\gansheyi_daima\matlab\path_to_your_output_folder2'; % 确保这是正确的路径
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder); % 如果输出文件夹不存在，则创建它
end

% 生成和保存9张相移干涉图像
A = 0.5; % 设置振幅，增加对比度
for n = 1:9
    % 相位移动
    phase_shift = (n - 1) * pi / 2;
    
    % 生成干涉图像，增加对比度
    I = 0.5 + A * cos(P + phase_shift);
    I = I .* mask; % 再次应用掩膜
    
    % 保存图像
    filename = fullfile(outputFolder, sprintf('interference_%d.bmp', n));
    imwrite(I, filename);
end

 

整体来说，这段代码通过模拟光波的相位分布并将其应用于高斯光束，来生成具有波前误差的干涉图像。这种方法可以用于模拟实际的光学系统中的波前，并为光学测试和分析提供理论模型。

下面 是读取9帧移相干涉条纹图

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import matplotlib.patches as mpatches

# 设置中文兼容字体
plt.rcParams['font.sans-serif'] = ['SimHei']

# 创建一个3x3网格的子图
fig, axs = plt.subplots(3, 3, figsize=(10, 8))

# 图像文件路径
image_files = [f'C:/Users/Administrator/Desktop/Untitled Folder/matlab/simulation/{i}.bmp' for i in range(9)]

# 为每张图像添加标题和图例
for ax, image_file, i in zip(axs.ravel(), image_files, range(1, 10)):
    img = mpimg.imread(image_file)
    ax.imshow(img, cmap='gray')
    ax.axis('off')

    # 设置标题
    ax.set_title(f'图{i}')

    # 创建代理艺术家用于图例
    proxy = mpatches.Rectangle((0, 0), 1, 1, fc="w", transform=ax.transAxes, alpha=0)
    ax.legend([proxy], [f'图例{i}'], loc='upper right')

plt.tight_layout()
plt.show()

matlab 的代码：

% 设置图像文件夹和文件名
folder = 'C:/Users/Administrator/Desktop/Untitled Folder/matlab/simulation/'; % 更改为您的图像文件夹路径
filePattern = fullfile(folder, '*.bmp');
bmpFiles = dir(filePattern);

% 创建一个新图形窗口
figure;

% 遍历所有图像文件
for k = 1:9
    baseFileName = bmpFiles(k).name;
    fullFileName = fullfile(folder, baseFileName);
    fprintf(1, 'Now reading %s\n', fullFileName);

    % 读取图像
    imageArray = imread(fullFileName);
    imageSize = size(imageArray);
    
    % 创建掩膜（圆形）
    [columnsInImage, rowsInImage] = meshgrid(1:imageSize(2), 1:imageSize(1));
    % 中心点
    centerX = imageSize(2) / 2;
    centerY = imageSize(1) / 2;
    % 半径
    radius = min([centerX centerY]) - 5; % 减5像素以保持圆形在图像内部
    % 创建二维逻辑索引
    circlePixels = (rowsInImage - centerY).^2 ...
        + (columnsInImage - centerX).^2 <= radius.^2;
    
    % 初始化一个全黑的图像
    maskedImage = uint8(zeros(size(imageArray)));
    
    % 使用圆形掩膜来定义ROI区域内的像素
    maskedImage(circlePixels) = imageArray(circlePixels);
    
    % 显示图像
    subplot(3,3,k);
    imshow(maskedImage); % 显示应用掩膜后的图像
    title(sprintf('Image #%d', k));
end

四步移相法：

作为经典的解相算法，四步移相法4]可以有效计算出理想移相图组的相位分布。先获取无移相的干涉图像，再进行三次移相间隔为π/2的移相，每次移相完成采集1帧图像，共得到4帧相位相差π/2的干涉图像，其图像光强信息分别为：

 (2.3)~(2.6)进行公式变换消去A和B可得到p(x,y)和四帧图像像素信息的关系式：

 该算法求解速度快，适合对理想干涉图组进行求解。当有外部干扰时，相位变化呈无规则态，此时四步移相法的解相效果非常糟糕，无任何抗振能力。与四步移相法类似的单周期定步长算法还有三步算法等。

其中原理部分转载与下方硕士论文：

参考文献：翟天保 硕士论文 基于干涉图像质量分析的激光干涉仪抗振技术研究