PUNRAP.m demonstrates phase unwrapping using synthetic data for light

passing through a glass bead.  Neglect refraction, and simply compute
straight line distances.

by Chuck DiMarzio Northeastern University May 2009

!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.

Contents

Set parameters

x=(0:200)*1e-6; % location in meters
r=75e-6; % bead radius
dn=0.04;  % Index of refraction of bead minus that of background
lambda=0.6328e-6; % HeNe laser wavelength in microns
m=0.3;% prism slope

Compute thickness and phase

y=m*x;
field=exp(1i*2*pi*dn*y/lambda);  % Complex field
phi=angle(field);   % Phase

Add some random noise

rand('twister',5489);  % mlint doesn't like this line.
grvreal=sum((rand(10,length(x))-0.5))*0.65; % Gaussian random variable
grvimag=sum((rand(10,length(x))-0.5))*0.65; % Gaussian random variable
field1=field+grvreal+1i*grvimag;  % Add complex noise to field
phi1=angle(field1);   % Phase
phi2=unwrap(phi1);  % Matlab unwrap phase

fig1=figure;plot(x*1e6,phi,x*1e6,phi1,x*1e6,phi2);grid on;
xlabel('x, Location, \mu m');
ylabel('\phi, phase, rad');

Estimate y from unwrapped phase

yest=phi2*lambda/2/pi/dn;
fig2=figure;plot(x*1e6,yest*1e6,x*1e6,y*1e6);
 grid on;
xlabel('x, Location, \mu m');
ylabel('y, thickness, \mu m');