%模擬波浪圖像
%clear?all
%tic;
%L=2000;d=200;f0=1/16;w=2*pi*f0;
%x=1:512;y=x;
%[xy]=meshgrid（xy）;
%z=3*peaks（512）;
%figure（1）;mesh（z）;axis?on
%xlabel（‘x（pixel）‘）;ylabel（‘y（pixel）‘）;zlabel（‘h（mm）‘）

%基準光柵
%J0=128+127*cos（w*x）;
%J0=mat2gray（J0）;
%figure（2）;imshow（J0）;axis?on

%變形光柵
%J1=128+127*cos（w*（x+z*d./（L-z）））;
%J1=mat2gray（J1）;
%J1=imnoise（J1‘gaussian‘00.001）;
%figure（3）;imshow（J1）;axis?on

%參考李中偉博士論文實現相移法加三頻外差解相位仿真
%?程序開始

clc;

close?all;

clear;

?

%?圖片的初始化

width?=?1024;??

heigth?=?768;

?

%?三頻率

%?這個可以參見李中偉的博士論文

freq?=?[70?64?59];??%像素單位為個數，可以看做頻率；正弦函數為周期含義
f12?=abs（（freq（1）*freq（2））/（freq（1）-freq（2）））;
f23?=?abs（（freq（2）*freq（3））/（freq（2）-freq（3）））;
f13?=?abs（（freq（1）*freq（3））/（freq（1）-freq（3）））;
f123?=?abs（（f12*f23）/f12-f23）;

%?利用分塊矩陣C存儲3組共計12張圖

%?三種頻率，四組相位

C?=?cell（34）;??

for?i=1:3

????for?j=1:4

????????C{ij}?=?zeros（heigthwidth）;

????end

end

?

%?利用余弦函數計算12張圖的灰度值

%?圖像的生成

%?三種頻率，四組相位

for?i?=?1:3?%?對應三種不同的頻率

????for??j?=?0:3?%?對應四種相位

????????for?k?=?1:width

????????????C{ij+1}（:k）?=?128+127*sin（2*pi*k*freq（i）/width+j*pi/2）;

????????end

????end

end

?

%?對灰度值進行歸一化處理

for?i?=?1:3

????for?j?=?1:4

????????C{ij}?=?mat2gray（C{ij}）;

????end

end

?

%?顯示12張圖

?for?i?=?1:3

???for?j?=?1:4

???????n?=?4*（i-1）+j;

?????????h?=?figure（n）;

????????imshow（C{ij}）;

?????end

end

?

%?初始化三組處理后的圖片灰度矩陣

%?phi也是分塊矩陣

%?存儲相位主值圖像

phi?=?cell（31）;

for?i?=?1:3

????phi{i1}?=?zeros（heigthwidth）;

end

?

%?求取相位差

%?計算每種頻率對應的相位主值

%?輸出三種頻率的相位主值，用于相差計算

for?i?=?1:3?%?對于3組中的每一組圖片，每一組相同頻率的有四張圖片

?????I1?=?C{i1};

?????I2?=?C{i2};

?????I3?=?C{i3};

?????I4?=?C{i4};

??????for?g?=?1:heigth

??????????for?k?=?1:width??????????

????????????if?????I4（gk）==I2（gk）&&I1（gk）>I3（gk）?%四個特殊位置

????????????????????phi{i1}（gk）=0;

????????????elseif?I4（gk）==I2（gk）&&I1（gk）
????????????????????phi{i1}（gk）=pi;?

????????????elseif?I1（gk）==I3（gk）&&I4（gk）>I2（gk）?%四個特殊位置

????????????????????phi{i1}（gk）=pi/2;

????????????elseif?I1（gk）==I3（gk）&&I4（gk）
????????????????????phi{i1}（gk）=3*pi/2;

????????????elseif?I1（gk）
????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））+pi;

????????????elseif?I1（gk）>I3（gk）&&I4（gk）>I2（gk）?%第一象限

????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））;

????????????elseif?I1（gk）>I3（gk）&&I4（gk）
????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））+2*pi;??

????????????end

??????????end?????????????

?????end

end

??

%?計算相差

%?保存矩陣，用于多頻相差的計算

PH1?=?phi{11};???%頻率1

PH2?=?phi{21};???%頻率2

PH3?=?phi{31};???%頻率3

?

%?初始化相差變量

%?多頻相差

PH12?=?zeros（heigthwidth）;

PH23?=?zeros（heigthwidth）;

PH13?=?zeros（heigthwidth）;

PH123?=?zeros（

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????6449??2019-09-29?13:03??Multi_frequency_heterodyne_phase_solution.m

-----------?---------??----------?-----??----

?????????????????6449????????????????????1